/*
 * wave1-glsc.c --- 1次元波動方程式
 *
 *   ccmg wave1-glsc.c wave1-initial.c
 *
 *   入力の目安:
 *     n は適度に大きめに (500 とか 1000)
 *     λは 1 (安定性のために≦1, 精度のために=1 が望ましいので)
 *     Tmax は適当に (10 とか 100 とか 1000 とか)
 */

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

#ifndef G_DOUBLE
#define G_DOUBLE
#endif
#include <glsc.h>

typedef double *vector;

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

#include <sys/types.h>
#include <sys/time.h>
#include <unistd.h>

/* ミリ秒単位のタイマー */
int msleep(int ms)
{
	struct timeval timeout;
	timeout.tv_sec = ms / 1000;
	timeout.tv_usec = (ms % 1000) * 1000;
	if (select(0, (fd_set *) 0, (fd_set *) 0, (fd_set *) 0, &timeout) < 0) {
		perror("usleep");
		return -1;
	}
	return 0;
}

void usage(char *prog)
{
  fprintf(stderr, "usage: %s [-n] [-input]\n", prog);
  fprintf(stderr, "  -n: Neumann, -input: input parameters manually\n");
  exit(0);
}

double pi;

int main(int argc, char **argv)
{
  /* 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];
  int neumann;
  double dt = 0.01;
  int skip;

  pi = 4.0 * atan(1.0);

  N = 360;
  lambda = 1;
  Tmax = 1000;
  neumann = 0;

  for (i = 1; i < argc; i++) {
    if (strcmp(argv[i], "-input") == 0) {
      /* N は 100 とか 1000 とか */
      printf("N=");
      scanf("%d", &N);
      /* λについては本に詳しいが、λ=1 が良い結果を出す */
      printf("λ=");
      scanf("%lf", &lambda);
      /* Tmax は最終時刻だから、長く計算したければ大きな値を入れる */
      printf("Tmax=");
      scanf("%lf", &Tmax);
    }
    else if (strcmp(argv[1], "-n") == 0) {
      neumann = 1;
    }
    else
      usage(argv[0]);
  }
  prompt();
  scanf("%d", &nfunc);

  u1 = malloc(sizeof(double) * (N + 1));
  u2 = malloc(sizeof(double) * (N + 1));
  u3 = malloc(sizeof(double) * (N + 1));

  lambda2 = lambda * lambda;
  h = 1.0 / N;
  tau = lambda * h;
  skip = rint(dt / tau);
  if (skip <= 0)
    skip = 1;
  dt = skip * tau;

  /* 初期値 */
  /* 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);
  if (neumann) {
    u2[0] = (1 - lambda2) * u1[0] + lambda2 * u1[1] + tau * psi(0.0, nfunc);
    u2[N] = (1 - lambda2) * u1[N] + lambda2 * u1[N - 1] + tau * psi(1.0, nfunc);
  }
  else {
    u2[0] = u2[N] = 0.0;
  }

  /* ***************** グラフィックスの準備 ***************** */
  /* メタファイル名は "WAVE",
   * ウィンドウのサイズは、
      横 win_width  + 2 * w_margin, 縦 win_height + 2 * h_margin */
  win_width = 200.0; win_height = 200.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];
    if (neumann) {
      u3[0] = 2 * (1.0 - lambda2) * u2[0] + 2 * lambda2 * u2[1] - u1[0];
      u3[N] = 2 * (1.0 - lambda2) * u2[N] + 2 * lambda2 * u2[N - 1] - u1[N];
    }
    else {
      u3[0] = u3[N] = 0.0;
    }

    if (n % skip == 0) {
      msleep(2 * (int) (dt * 1000));
      /* 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;
}