/*
 * heat1d-i.c -- 陰的スキームで熱方程式を解く
 *
 *   菊地・山本, 微分方程式と計算機演習, 山海堂 (1991) 第11章
 *   p237 のプログラムを修正・拡張したもの
 *
 C     ** IMPLICIT FINITE DIFFERENCE METHOD **
 C     **         FOR HEAT EQUATION         **
 C     nfunc : 関数の番号（5種類の関数を用意している）
 C     theta : スキームのパラメーター
 C     N : 分割の数
 C     t : 時刻
 C     tau : 時間の刻み幅
 C     Tmax : 時刻の上限
 C     h : 空間の刻み幅
 C     AL,AD,AU : 係数行列
 C     u : DIMENSION FOR UNKNOWN FUNCTION
 *
 *  このプログラムは
 *    http://nalab.mind.meiji.ac.jp/%7Emk/program/
 *  から入手可能
 */

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

/* 円周率 --- main() で pi = 4.0 * atan(1.0); と代入 */

double pi;

void trilu(int, double *, double *, double *);
void trisol(int, double *, double *, double *, double *);
double f(double, int);
void print100(int, double, double *);
void axis(double, double, double, double);
void fsymbol(double, double, char *);
void print_parameters(double, double, double, double,
                      int, double, int,
                      double, double, double, double);
void save_picture(void);

int main()
{
    int N, nfunc, i, n, nMax;
    double a, b, theta, h, Tmax, tau, c1, c2, c3, c4, lambda, t;
    double *AL, *AD, *AU, *ff, *u;
    double xleft, ybottom, xright, ytop;
    /* グラフを書き換える時間間隔 */
    double dt;
    /* 何ステップおきにグラフを書き換えるか */
    int skip;
    /* 以前描いたグラフを消すか？ (0=No, 1=Yes) */
    int erase_always = 0;
    /* 数値を表示するか? (0=No, 1=Yes) */
    int print_numerical_data = 0;

    /* 円周率 */
    pi = 4.0 * atan(1.0);
    /* 区間の左端、右端 */
    a = 0.0;
    b = 1.0;

    /* 入力 */
    printf("入力して下さい : nfunc(1..5)=");
    scanf("%d", &nfunc);
    printf("入力して下さい : θ=");
    scanf("%lf", &theta);
    printf("入力して下さい : N=");
    scanf("%d", &N);
    printf("入力して下さい : λ=");
    scanf("%lf", &lambda);

    /* ベクトルの準備 */
    if ((AL = malloc(N * sizeof(double))) == NULL) {
        fprintf(stderr, "cannot allocate AL[]\n");
        exit(1);
    }
    if ((AD = malloc(N * sizeof(double))) == NULL) {
        fprintf(stderr, "cannot allocate AD[]\n");
        exit(1);
    }
    if ((AU = malloc(N * sizeof(double))) == NULL) {
        fprintf(stderr, "cannot allocate AU[]\n");
        exit(1);
    }
    if ((ff = malloc(N * sizeof(double))) == NULL) {
        fprintf(stderr, "cannot allocate ff[]\n");
        exit(1);
    }
    if ((u = malloc((N + 1) * sizeof(double))) == NULL) {
        fprintf(stderr, "cannot allocate u[]\n");
        exit(1);
    }
    /* */
    h = (b - a) / N;
    tau = lambda * h * h;
    printf(" 時間の刻み幅τ= %g になりました。\n", tau);

    /* 入力 */
    printf("入力して下さい : 最終時刻 Tmax=");
    scanf("%lf", &Tmax);
    printf("入力して下さい : グラフ書き換え時間間隔 Δt=");
    scanf("%lf", &dt);
    skip = (int) rint(dt / tau);
    /* skip が 0 になっていないかチェック */
    if (skip == 0) {
	fprintf(stderr, "Δtが小さすぎるので、Δt=τ=%g にします。\n", tau);
	skip = 1;
	dt = tau;
    }

    /* 係数行列の準備 */
    c1 = - theta * lambda;
    c2 = 1.0 + 2.0 * theta * lambda;
    c3 = 1.0 - 2.0 * (1.0 - theta) * lambda;
    c4 = (1.0 - theta) * lambda;
    for (i = 1; i < N; i++) {
        AL[i] = c1;
        AD[i] = c2;
        AU[i] = c1;
    }

    /* 初期値 */
    for (i = 0; i <= N; i++)
        u[i] = f(a + i * h, nfunc);

    /* 出力（t=0） */
    /* 出力（ここでは画面に数値を表示すること）は FORTRAN と C で
       かなり異なるので、直接の翻訳は出来ない。ごたごたするので、
       print100 という関数にまとめた */
    t = 0.0;
    print100(N, t, u);

    /* ***************** グラフィックスの準備 ***************** */
    /* グラフィックス画面の準備 -- まず画面の範囲を決めてから */
    xleft = a - (b - a) / 10;
    xright = b + (b - a) / 10;
    ybottom = - 0.1;
    ytop = 1.1;

    openpl();
    fspace2(xleft, ybottom, xright, ytop);

    /* パラメーター、入力データ等の表示 */
    print_parameters(xleft, ybottom, xright, ytop,
                     nfunc, theta, N, lambda, tau, Tmax, t);

    /* 解 u の t=0 におけるグラフを描く */
    fmove(a, u[0]);
    for (i = 1; i <= N; i++)
        fcont(a + i * h, u[i]);
    xsync();

    /* 係数行列の LU 分解 */
    trilu(N - 1, AL + 1, AD + 1, AU + 1);

    /* 繰り返し計算 */
    nMax = (int) rint(Tmax / tau);
    for (n = 1; n <= nMax; n++) {
        /* 連立1次方程式を作って解く */
        for (i = 1; i < N; i++)
            ff[i] = c3 * u[i] + c4 * (u[i - 1] + u[i + 1]);
        trisol(N - 1, AL + 1, AD + 1, AU + 1, ff + 1);
        for (i = 1; i < N; i++)
            u[i] = ff[i];
        /* 境界条件 */
        u[0] = 0.0;
        u[N] = 0.0;

        /* 時刻 */
        t = n * tau;

        if (n % skip == 0) {
            /* 数値データの表示 */
            if (print_numerical_data)
                print100(N, t, u);

            /* 解のグラフを描く */
            if (erase_always) {
                /* これまで描いたものを消し、パラメーターを再表示する */
                erase();
                print_parameters(xleft, ybottom, xright, ytop,
                                 nfunc, theta, N, lambda, tau, Tmax, t);
            }
            fmove(a, u[0]);
            for (i = 1; i <= N; i++)
                fcont(a + i * h, u[i]);
            xsync();
        }
    }

    save_picture();
    printf("マウスでウィンドウをクリックして下さい。\n");
    closepl();

    return 0;
}

/**********************************************************/

/* 三重対角行列の LU 分解 (pivoting なし) */
void trilu(int n, double *al, double *ad, double *au)
{
  int i, nm1 = n - 1;
  /* 前進消去 (forward elimination) */
  for (i = 0; i < nm1; i++) {
    al[i + 1] /= ad[i];
    ad[i + 1] -= au[i] * al[i + 1];
  }
}

/* LU 分解済みの三重対角行列を係数に持つ三項方程式を解く */
void trisol(int n, double *al, double *ad, double *au, double *b)
{
  int i, nm1 = n - 1;
  /* 前進消去 (forward elimination) */
  for (i = 0; i < nm1; i++) b[i + 1] -= b[i] * al[i + 1];
  /* 後退代入 (backward substitution) */
  b[nm1] /= ad[nm1];
  for (i = n - 2; i >= 0; i--) b[i] = (b[i] - au[i] * b[i + 1]) / ad[i];
}

/**********************************************************/

/* 初期条件を与える関数。複数登録して番号 nfunc で選択する。 */

double f(double x, int nfunc)
{
    /* f(x)=min(x,1-x) */
    if (nfunc == 1) {
        if (x <= 0.5)
            return x;
        else
            return 1.0 - x;
    }
    /* f(x)=1 */
    else if (nfunc == 2)
        return 1.0;
    else if (nfunc == 3)
        return sin(pi * x);
    /* f(x)= 変な形をしたもの(by M.Sakaue) */
    else if (nfunc == 4) {
        if (x <= 0.1)
            return 5 * x;
        else if (x <= 0.3)
            return -2 * x + 0.7;
        else if (x <= 0.5)
            return 4.5 * x - 1.25;
        else if (x <= 0.7)
            return -x + 1.5;
        else if (x <= 0.9)
            return x + 0.1;
        else
            return -10 * x + 10.0;
    }
    /* やはり非対称な形をしたもの(by M.Sakaue) */
    else if (nfunc == 5) {
        if (x <= 0.2)
            return -x + 0.8;
        else if (x <= 0.3)
            return 4 * x - 0.2;
        else if (x <= 0.4)
            return -4 * x + 2.2;
        else if (x <= 0.7)
            return -x + 1.0;
        else if (x <= 0.8)
            return 4 * x - 2.6;
        else
            return -x + 1.5;
    } else
        return 0;
}

/**********************************************************/

/* 配列 u[] の内容(u[0],u[1],...,u[n])を一行 5 個ずつ表示する。 */

void print100(int N, double t, double *u)
{
    int i;
    printf("T= %12.4e\n", t);
    for (i = 0; i < 5; i++)
        printf("  I     u(i)    ");
    printf("\n");
    for (i = 0; i <= N; i++) {
        printf("%3d%12.4e ", i, u[i]);
        if (i % 5 == 4)
            printf("\n");
    }
    printf("\n");
}

/**********************************************************/

/* ウィンドウに座標軸を表示する */
void axis(double xleft, double ybottom, double xright, double ytop)
{
    linemod("dotted");
    fline(xleft, 0.0, xright, 0.0);
    fline(0.0, ybottom, 0.0, ytop);
    linemod("solid");
}

/* ウィンドウ内の位置 (x,y) から文字列 s を表示する */
void fsymbol(double x, double y, char *s)
{
    fmove(x, y);
    label(s);
}

#define X(x)    (xleft+(x)*(xright-xleft))
#define Y(y)    (ybottom+(y)*(ytop-ybottom))

/* パラメーターの値等をウィンドウに表示する */
void print_parameters(double xleft, double ybottom, double xright, double ytop,
                      int nfunc, double theta, int N,
                      double lambda, double tau, double Tmax, double t)
{
    char message[80];

    axis(xleft, ybottom, xright, ytop);
    snprintf(message, sizeof(message), "heat equation, u(0)=u(1)=0");
    fsymbol(X(0.2), Y(0.95), message);
    snprintf(message, sizeof(message),
            "nfunc=%d, theta=%g, N=%d, lambda=%g, tau=%g, Tmax=%g",
            nfunc, theta, N, lambda, tau, Tmax);
    fsymbol(X(0.2), Y(0.9), message);
    snprintf(message, sizeof(message), "t = %g", t);
    if (t != 0.0)
        fsymbol(X(0.2), Y(0.85), message);
}

/* 現在ウィンドウに表示されている図をファイルに保存する */
void save_picture(void)
{
    char filename[200];
    printf("図を保存するファイル名: ");
    scanf("%s", filename);
    mkplot(filename);
}