/*
 *   wave2d.c --- 2次元波動方程式 (Dirichlet 境界条件)
 *   コンパイル
 *   ccmg wave2d.c
 */

#include <stdio.h>
#include <math.h>
#include <matrix.h>
#ifndef G_DOUBLE
#define G_DOUBLE
#endif
#include <glsc.h>

double sqr(double x) { return x * x; }

double pi;

int main(void)
{

  int     Nx, Ny, i, j, n, skip, nMax, nfunc;
  double  hx,hy,lambdax,lambday,lambdax2,lambday2,tau,Tmax,t,c,dt;
  matrix  u1,u2,u3;
  double  phi( double, double,int), psi( double, double,int);
  double  exact(double x, double y, double t);
  char    message[100];

  pi = 4.0 * atan(1.0);  
  
  printf("Nx = , Ny = : "); scanf("%d %d", &Nx, &Ny);

  hx = 1.0 / Nx; hy = 1.0/ Ny;
  printf("τ (限界は%g) : ", 1.0 / sqrt(1/(hx*hx) + 1/(hy*hy)));
  scanf("%lf", &tau);
  
  lambdax = tau / hx; lambday = tau / hy;
  printf("λx = %g になりました。\n", lambdax); 
  printf("λy = %g になりました。\n", lambday); 
  lambdax2 = lambdax * lambdax; lambday2 = lambday * lambday;
  printf("λx * λx + λy * λy = %g になりました。\n", lambdax2 + lambday2);  
  
  printf("Tmax = : "); scanf("%lf", &Tmax);
  printf("Δt = : ");  scanf("%lf", &dt);
  printf("nfunc = "); scanf("%d", &nfunc);
  
  /* メモリの確保 */
  if ((u1 = new_matrix(Nx + 1, Ny + 1)) == NULL) {
    fprintf(stderr, "配列u1を確保できませんでした,\n");
    exit(1);
  }
  if ((u2 = new_matrix(Nx + 1, Ny + 1)) == NULL) {
    fprintf(stderr, "配列u2を確保できませんでした,\n");
    exit(1); 
  }
  if ((u3 = new_matrix(Nx + 1, Ny + 1)) == NULL) {
    fprintf(stderr, "配列u3を確保できませんでした,\n");
    exit(1);
  }
  
  c = 1.0 - lambdax2 - lambday2;
  
  skip = dt / tau;

  /* GLSC の開始 */
  g_init("wave2d", 250.0, 180.0);
  g_device(G_BOTH);
  g_def_scale(0,
              0.0, 1.0, -4.0, 4.0,
              150.0, 35.0, 80.0, 140.0);
  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, 2);
  g_sel_scale(0);
  
  /* t=0 の時グラフを描く*/
  for (i=0; i <= Nx; i++)
    for (j=0; j <= Ny; j++)
      u1[i][j] = phi(i * hx, j * hy, nfunc);
  g_cls();
  g_text(50.0, 10.0,
         "The numerical simulation of a 2-dimensional wave equation");
  snprintf(message, sizeof(message),
           "Nx=%d,Ny=%d,tau=%g,lambdax=%g,lambday=%g,Tmax=%g, dt=%g",
           Nx, Ny, tau, lambdax, lambday, Tmax, dt);
  g_text(50.0, 20.0, message);
  snprintf(message, sizeof(message), "t = %g", 0 * tau);
  g_text(50.0, 30.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, 2.0);
  /* 断面図を描く */
  g_sel_line(0);
  g_move(0.0, u1[0][Ny/2]);
  for (i=1; i<= Nx; i++)
        g_plot(i * hx, u1[i][Ny/2]);
  /* 鳥瞰図を描く */
  g_hidden2(5.0, 5.0, 1.0, -1.0, 1.0,
            15.0, 25.0, 30.0,
            20.0, 60.0, 100.0, 80.0,
            u1, Nx + 1, Ny + 1, 1,
            G_SIDE_NONE, 2, 1);
  printf("u1 まで描けた\n");

  /* t=τ の時のグラフを描く*/
  for (i = 1; i < Nx; i++)
    for (j = 1; j < Ny; j++)
      u2[i][j] = c * u1[i][j] 
                 + 0.5 * lambdax2 * (u1[i+1][j] + u1[i-1][j])
                 + 0.5 * lambday2 * (u1[i][j+1] + u1[i][j-1])
                 + tau * psi(i * hx, j * hy, nfunc);

  /* Dirichlet 境界条件 */
  for (i = 0; i <= Nx; i++)
    u2[i][0] = u2[i][Ny] = 0.0;
  for (j = 0; j <= Ny; j++)
    u2[0][j] = u2[Nx][j] = 0.0;  

  g_cls();
  snprintf(message, sizeof(message), "t = %g", tau);
  g_text(50.0, 30.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, 2.0);
  /* 断面図を描く */
  g_sel_line(0);
  g_move(0.0, u2[0][Ny/2]);
  for (i = 1; i <= Nx; i++)
    g_plot(i * hx, u2[i][Ny/2]);
  /* 鳥瞰図を描く */
  g_hidden2(5.0, 5.0, 1.0, -1.0, 1.0,
            15.0, 25.0, 30.0,
            20.0, 60.0, 100.0, 80.0,
            u2, Nx + 1, Ny + 1, 1,
            G_SIDE_NONE, 2, 1);
  printf("u2 まで描けた\n");

  nMax = rint(Tmax /tau);

  /*　t=nτ (n≧2)　の時のグラフを描く*/
  for (n = 2; n <= nMax; n++) {
    for (i = 1; i < Nx; i++)
      for (j = 1; j < Ny; j++)
        u3[i][j] = 2 * c * u2[i][j] 
                   + lambdax2 * (u2[i+1][j] + u2[i-1][j])
                   + lambday2 * (u2[i][j+1] + u2[i][j-1])
                   - u1[i][j];
    /* Dirichlet 境界条件 */
    for (i = 0; i <= Nx; i++) 
      u3[i][0] = u3[i][Ny] = 0.0;   
    for (j = 0; j <= Ny; j++)
      u3[0][j] = u3[Nx][j] = 0.0;

    /* 配列の交換 */
    for (i = 0; i <= Nx; i++)
      for (j = 0; j <= Ny; j++) {
        u1[i][j] = u2[i][j];  
        u2[i][j] = u3[i][j];
      }
    if (n % skip == 0) {
       g_cls();
       /* 座標軸を描く */
       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, 2.0);
       /* 断面図を描く */
       g_sel_line(0);
       g_move(0.0, u2[0][Ny/2]);
       for (i = 1; i <= Nx; i++)
         g_plot(i * hx, u2[i][Ny/2]);
       t = n * tau;
       snprintf(message, sizeof(message), "t = %g", t);
       g_text(50.0, 30.0, message);
       /* 鳥瞰図を描く */
       g_hidden2(5.0, 5.0, 1.0, -1.0, 1.0,
                 15.0, 25.0, 30.0,
                 20.0, 60.0, 100.0, 80.0,
                 u3, Nx + 1, Ny + 1, 1,
                 G_SIDE_NONE, 2, 1);
       if (nfunc == 1)
	 printf("t=%g, u[Nx/2][Ny/2]=%g, 厳密解=%g\n",
		t, u2[Nx/2][Ny/2], exact(0.5,0.5,t));
    }
  }
  printf("終了したらウィンドウをクリックして終了してください。\n");
  g_sleep(-1.0);
  g_term();
  
  return 0;
}

double h(double r)
{
  double a = 0.1, b = 0.2;
  if (r >= a && r <= b)
    return 4 * sqr(r - a) * sqr(r - b) / sqr(sqr(b - a));
  else
    return 0.0;
}

double dh(double r)
{
  double a = 0.1, b = 0.2;
  if (r >= a && r <= b)
    return 12 * (r - a) * (r - b) * (2 * r - a - b) / sqr(sqr(b - a));
  else
    return 0.0;
}

/* 初期条件：φ、ψ */
double phi(double x, double y, int nfunc)
{
  if (nfunc == 1)
    return sin(pi * x) * sin(pi * y);
  else if (nfunc ==2 || nfunc == 3) {
    double r;
    r = sqrt(sqr(x - 0.5) + sqr(y - 0.5));
    if (r > 1e-4)
      return h(r) / r;
    else
      return 0.0;
  }
  else
    return 0.0;
}

double psi( double x, double y, int nfunc)
{
  if (nfunc == 1) {
    return 0.0;
  }
  else if (nfunc == 2) {
    double r;
    r = sqrt(sqr(x - 0.5) + sqr(y - 0.5));
    if (r > 1e-4)
      return - 1.0 * dh(r) / r;
    else
      return 0.0;
  }
  else if (nfunc == 3) {
    return 0.0;
  }
  else
    return 0.0;
}

double exact(double x, double y, double t)
{
  return cos(sqrt(2.0) * pi * t) * sin(pi * x) * sin(pi * y);
}