/*
 * heat1d-e-glsc3d.c -- 1次元熱伝導方程式の初期値境界値問題
 *   コンパイル:  heat1d-e-glsc3d.c
 *
 *  u_t(x,t)=u_{xx}(x,t)        x∈(0,1), t∈(0,∞)
 *  u(0,t)=u(1,t)=0             t∈(0,∞)
 *  u(x,0)=f(x)                 x∈[0,1]
 *
 * を陽解法で解く。
 */

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

int main(void)
{
  int i, n, nMax, N;
  double tau, h, lambda, Tmax;
  double *u, *newu;
  double f(double);
  double win_width, win_height, w_margin, h_margin;
  double xleft, xright, ybottom, ytop;

  /* N, λ を入力する */
  printf("区間の分割数 N = "); scanf("%d", &N);
  printf("λ (=τ/h^2) = "); scanf("%lf", &lambda);

  /* 時間の刻み幅 h, τ を計算する */
  h = 1.0 / N;
  tau = lambda * h * h;
  printf("τ=%g\n", tau);

  /* 最終時刻を入力する */
  printf("最終時刻 Tmax = "); scanf("%lf", &Tmax);

  /* ベクトル u, newu を用意する */
  u = malloc(sizeof(double) * (N+1));
  newu = malloc(sizeof(double) * (N+1));
  if (u == NULL || newu == NULL) {
    fprintf(stderr, "メモリーが準備できませんでした。\n");
    exit(1);
  }
  /* 初期値の代入 */
  for (i = 0; i <= N; i++)
    u[i] = f(i * h);

  /* ***************** グラフィックスの準備 ***************** */
  /* 表示する範囲 [xleft,xright]×[ybottom,ytop] を決定 */
  xleft = - 0.1; xright = 1.1; ybottom = - 0.1; ytop = 1.1;
  /* メタファイル名は "HEAT1DE",
   * ウィンドウのサイズは、
      横 win_width  + 2 * w_margin, 縦 win_height + 2 * h_margin
   */
  win_width = 600.0; win_height = 600.0; w_margin = 30.0; h_margin = 30.0;
  g_init("HEAT1DE", win_width  + 2 * w_margin, win_height + 2 * h_margin);
  /* 座標系の定義: [xleft,xright]×[ybottom,ytop] という閉領域を表示する */
  g_def_scale_2D(0,
		 xleft, xright, ybottom, ytop,
		 w_margin, h_margin, win_width, win_height);
  g_cls(); // ないとマズイ？

  /* 線を二種類用意する */
  // 色は黒 0,0,0,1 (R,G,B,透明度)
  // 線のtypeは用意されていない 0と1
#define G_LINE_SOLID (0)
#define G_LINE_DOTS (1)
  g_def_line(0, 0,0,0,1, 1, G_LINE_SOLID);
  g_def_line(1, 0,0,0,1, 2, G_LINE_DOTS); // 太くしないと点線のように見える
  /* 表示するための文字列の属性を定義する */
  g_def_text(0, 0,0,0,1, 18); // 文字のサイズの単位が全然違う
  /* 定義したものを選択する */
  g_sel_scale(0); g_sel_line(0); g_sel_text(0);

  /* タイトルと入力パラメーターを表示する */
  g_text_standard(90.0, 90.0,
		  "heat equation, homogeneous Dirichlet boundary condition");
  g_text_standard(90.0, 180.0, "N=%d, λ=%g, Tmax=%g", N, lambda, Tmax);
  /* 座標軸を描く */
  g_sel_line(1);
  g_move_2D(xleft, 0.0); g_plot_2D(xright, 0.0);
  g_move_2D(0.0, ybottom); g_plot_2D(0.0, ytop);
  g_sel_line(0);

  /* t=0 の状態を表示する */
  g_move_2D(0.0, u[0]);
  for (i = 1; i <= N; i++)
    g_plot_2D(i * h, u[i]);

  /* ループを何回まわるか計算する (四捨五入) */
  nMax = rint(Tmax / tau);

  /* 時間に関するステップを進めるループ */
  for (n = 0; n < nMax; n++) {
    /* 差分方程式 (n -> n+1) */
    for (i = 1; i < N; i++)
      newu[i] = (1.0 - 2 * lambda) * u[i] + lambda * (u[i+1] + u[i-1]);
    /* 計算値を更新 */
    for (i = 1; i < N; i++)
      u[i] = newu[i];
    /* Dirichlet 境界条件 */
    u[0] = u[N] = 0.0;
    /* この時刻 (t=(n+1)τ) の状態を表示する */
    g_move_2D(0.0, u[0]);
    for (i = 1; i <= N; i++)
      g_plot_2D(i * h, u[i]);
    // 描画する (これをループの外に出すと速くなるが…)
    g_finish(); g_sleep(0.005);
  }

  printf("終りました。X の場合はウィンドウをクリックして下さい。\n");
  g_sleep(-1.0);
  /* ウィンドウを閉じる */
  return 0;
}

/* 初期条件 --- ここでは min(x,1-x) という山形のグラフを持つ関数 */
double f(double x)
{
  if (x <= 0.5)
    return x;
  else
    return 1.0 - x;
}