/*
 * wave1n-naive.c --- 1次元波動方程式 (Neumann 境界条件)
 *
 *   ccx wave1n-naive.c -lmatrix
 * ノートパソコンでは
 *   ccx wave1n-naive.c -I/usr/local/include -L/usr/local/lib -lmatrix
 *
 *   入力の目安:
 *     n は適度に大きめに (500 とか 1000)
 *     λは 1 (安定性のために≦1, 精度のために=1 が望ましいので)
 *     Tmax は適当に (10 とか 100 とか 1000 とか)
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fplot.h>

typedef double *vector;

void copy_vector(int, vector, vector);
int msleep(int);

double pi;

int main()
{
  /* N は区間の分割数 */
  int N;
  int i, n, nmax, nfunc;
  /* lambda はλ, tau はτ */
  double lambda, Tmax, h, tau, lambda2;
  /* u1[i]=u_{i,j-1}
   * u2[i]=u_{i,j}
   * u3[i]=u_{i,j+1]
   */
  vector u1, u2, u3;
  double phi(double, int), psi(double, int);

  pi = 4.0 * atan(1.0);

  /* n は 100 とか 1000 とか */
  printf("N=");
  scanf("%d", &N);
  u1 = malloc(sizeof(double) * (N + 1));
  u2 = malloc(sizeof(double) * (N + 1));
  u3 = malloc(sizeof(double) * (N + 1));

  /* λについては本に詳しいが、λ=1 が良い結果を出す */
  printf("λ=");
  scanf("%lf", &lambda);
  /* Tmax は最終時刻だから、長く計算したければ大きな値を入れる */
  printf("Tmax=");
  scanf("%lf", &Tmax);

  lambda2 = lambda * lambda;
  h = 1.0 / N;
  tau = lambda * h;

  /* 初期値を選ぶ */
  printf("1: 定常波, 2: 割れる山, 3: 滑かな山, 4: ホイヘンスの原理反例\n");
  printf("nfunc(初期値の種類 1..5)=");
  scanf("%d", &nfunc);
  /* 初期値 */
  /* u_{i,0}=φi */
  for (i = 0; i <= N; i++)
    u1[i] = phi(i * h, nfunc);
  /* u_{i,1}=(1-λ2)... */
  for (i = 1; i < N; i++)
    u2[i] = (1 - lambda2) * u1[i]
      + 0.5 * lambda2 * (u1[i - 1] + u1[i + 1]) + tau * psi(i * h, nfunc);
  u2[0] = u2[1]; u2[N] = u2[N-1];

  /* グラフを描くための準備 */
  openpl();
  fspace2(-0.2, -1.2, 1.2, 1.2);

  /* t=0 でのグラフ */
  fmove(0.0, u1[0]);
  for (i = 1; i <= N; i++)
    fcont(i * h, u1[i]);

  /* t=t1=τ でのグラフ */
  erase();
  fmove(0.0, u2[0]);
  for (i = 1; i <= N; i++)
    fcont(i * h, u2[i]);

  nmax = rint(Tmax / tau);
  /* 時間に関するループ */
  for (n = 1; n <= nmax; n++) {
    for (i = 1; i < N; i++)
      u3[i] = 2 * (1.0 - lambda2) * u2[i]
	+ lambda2 * (u2[i - 1] + u2[i + 1]) - u1[i];
    u3[0] = u3[1]; u3[N] = u3[N-1];

#ifdef USESLEEP
    /* 5 倍の秒数待つ */
    msleep(5 * (int) (tau * 1000));
#endif
    /* t=t_{n+1}=(n+1)τ でのグラフ */
    erase();
    fmove(0.0, u3[0]);
    for (i = 1; i <= N; i++)
      fcont(i * h, u3[i]);
    /* グラフ描画が終わるまで待つ */
    xsync();

    /* u1 <- u2, u2 <- u3 */
    copy_vector(N + 1, u1, u2);
    copy_vector(N + 1, u2, u3);
  }
  closepl();

  return 0.0;
}

/* 初期値 (nfunc=3,4) を作るための道具 */
double sqr(double x)
{
  return x * x;
}

double mountain(double x)
{
  double a = 0.4, b = 0.6;
  if (x <= a || x >= b)
    return 0;
  else {
    double c = sqr(sqr(2 / (b - a)));
    return c * sqr((x - a) * (x - b));
  }
}

/* 初期値 (nfunc=5) を作るための道具 */

/* これは初期値ではない！ */
double f(double x)
{
  if (x <= 0)
    return 0;
  else
    return exp(- 1.0 / x);
}

/* これは初期値ではない！ */
double g(double x)
{
  return f(x) / (f(x) + f(1 - x));
}

/*
 *             1   (|x|<a)
 *  h_ab(x) ={ 
 *             0   (|x|>b)
 *
 */
double h_ab(double x, double a, double b)
{
  return g((x + a) / (a - b)) * g((- x + a) / (a - b));
}

/* 初期値 φ=u(・,0) */
double phi(double x, int nfunc)
{
  if (nfunc == 1)
    return sin(pi * x);
  else if (nfunc == 2) {
    if (x > 0.375 && x <= 0.5)
      return 8.0 * x - 3.0;
    else if (x > 0.5 && x < 0.625)
      return -8.0 * x + 5.0;
    else
      return 0.0;
  }
  else if (nfunc == 3)
    return mountain(x);
  else if (nfunc == 4)
    return mountain(x);
  else if (nfunc == 5) {
    double a = 0.3, b = 0.2;
    return h_ab(x - 0.5, a, b);
  }
  else
    return 0.0;
}

/* 初期値 ψ=ut(・,0) */
double psi(double x, int nfunc)
{
  if (nfunc == 1)
    return 0.0;
  else if (nfunc == 2)
    return 0.0;
  else if (nfunc == 3)
    return 0.0;
  else if (nfunc == 4)
    return 5 * mountain(x);
  else if (nfunc == 5)
    return 0.0;
  else
    return 0.0;
}

/*
 * ベクトル b をベクトル a にコピー
 */
void copy_vector(int N, vector a, vector b)
{
  int i;
  for (i = 0; i < N; i++)
    a[i] = b[i];
  return;
}