/*
 * heat1d-e.c -- 1次元熱伝導方程式の初期値境界値問題を陽解法で解く。
 *   コンパイル:  ccx heat1d-e.c
 *
 *  http://nalab.mind.meiji.ac.jp/~mk/program/fdm/heat1d-e.c
 */

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

int main()
{
  int i, j, Jmax, N;
  double tau, h, lambda, Tmax;
  double *u, *newu;
  double f(double);

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

  /* h, τ を計算する */
  h = 1.0 / N;
  tau = lambda * h * h;
  printf("τ=%f\n", tau);

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

  /* ベクトル u, newu を用意する */
  u = malloc(sizeof(double) * (N+1));
  newu = malloc(sizeof(double) * (N+1));
  
  /* 初期値の代入 */
  for (i = 0; i <= N; i++)
    u[i] = f(i * h);

  /* ウィンドウを出す */
  openpl();
  fspace2(-0.1, -0.1, 1.1, 1.1);

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

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

  for (j = 0; j < Jmax; j++) {

    /* 差分方程式 (j -> j+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=(j+1)τ) の状態を表示する */
    fmove(0.0, u[0]);
    for (i = 1; i <= N; i++)
      fcont(i * h, u[i]);
  }
  mkplot("heat1d-e.plot");
  printf("終りました。ウィンドウをクリックして下さい。\n");
  /* ウィンドウを閉じる */
  closepl();
  return 0;
}

/* 初期条件 --- ここでは min(x,1-x)=1/2-|x-1/2| という中央が高い山形の関数 */
double f(double x)
{
  if (x <= 0.5)
    return x;
  else
    return 1.0 - x;
}