/*
 * wave1d-glsc.c --- 1次元波動方程式 (Dirichlet 境界条件)
 *
 *   ccmg wave1d-glsc.c
 *   ccmg wave1d-glsc.c msleep.c -DUSESLEEP
 *
 *  http://nalab.mind.meiji.ac.jp/~mk/program/fdm から入手可能
 *
 *   入力の目安:
 *     n は適度に大きめに (500 とか 1000)
 *     λは 1 (安定性のために≦1, 精度のために=1 が望ましいので)
 *     Tmax は適当に (10 とか 100 とか 1000 とか)
 */

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

#ifndef G_DOUBLE
#define G_DOUBLE
#endif
#include <glsc.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);
  double win_width, win_height, w_margin, h_margin;
  char message[100];

  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[N] = 0.0;

  /* ***************** グラフィックスの準備 ***************** */
  /* メタファイル名は "WAVE",
   * ウィンドウのサイズは、
      横 win_width  + 2 * w_margin, 縦 win_height + 2 * h_margin */
  win_width = 150.0; win_height = 150.0; w_margin = 10.0; h_margin = 10.0;
  g_init("WAVE", win_width  + 2 * w_margin, win_height + 2 * h_margin);
  /* 画面とメタファイルの両方に記録する */
  g_device(G_BOTH);
  /* 座標系の定義: [-0.1,1.1]×[-1.1,1.1] という閉領域を表示する */
  g_def_scale(0,
              -0.1, 1.1, -1.1, 1.1,
              w_margin, h_margin, win_width, win_height);
  /* 線を二種類用意する */
  g_def_line(0, G_BLACK, 0, G_LINE_SOLID);
  g_def_line(1, G_BLACK, 0, G_LINE_DOTS);
  /* 表示するための文字列の属性を定義する */
  g_def_text(0, G_BLACK, 3);
  /* 定義したものを選択する */
  g_sel_scale(0); g_sel_line(0); g_sel_text(0);

  /* タイトルと入力パラメーターを表示する */
  g_text(30.0, 30.0,
         "wave equation, homogeneous Dirichlet boundary condition");
  snprintf(message, sizeof(message),
          "N=%d, lambda=%g, Tmax=%g", N, lambda, Tmax);
  g_text(30.0, 60.0, message);

  /* 座標軸を表示する */
  g_sel_line(1);
  g_move(-0.1, 0.0); g_plot(1.1, 0.0);
  g_move(0.0, -0.1); g_plot(0.0, 1.1);
  g_sel_line(0);

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

  /* t=t1=τ でのグラフ */
  g_cls();
  g_move(0.0, u2[0]);
  for (i = 1; i <= N; i++)
    g_plot(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[N] = 0.0;

#ifdef USESLEEP
    /* 5 倍の秒数待つ */
    msleep(5 * (int) (tau * 1000));
#endif
    /* t=t_{n+1}=(n+1)τ でのグラフ */
    g_cls();
    g_move(0.0, u3[0]);
    for (i = 1; i <= N; i++)
      g_plot(i * h, u3[i]);

    /* u1 <- u2, u2 <- u3 */
    copy_vector(N + 1, u1, u2);
    copy_vector(N + 1, u2, u3);
  }
  printf("終りました。X の場合はウィンドウをクリックして下さい。\n");
  g_sleep(-1.0);
  /* ウィンドウを閉じる */
  g_term();
  return 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;
}