#include <stdio.h>
#include <math.h>
#include <conio.h>
#include <string.h>
#include <iostream.h>
#include <fstream.h>

double zd, A[101], B[101], C[101], D[101], 
  dz = 2,
  t = 0,
  Tmax = 30,
  Tmin = 10,
  omega = 2 * 3.14 / 24.0;

void get_analytical_solution(double T[], int n)
{
  int i;
  double z;
  for (i=0; i<=n; i++)
  {
    z = (i - 0.5) * dz;
    T[i] = 0.5 * (Tmax + Tmin)
         + 0.5 * (Tmax - Tmin) * exp(- z / zd) * sin(omega * t - z / zd);
  }
}

void solve_tridiagobal_matrix(double x[], int st, int n)
{
  int i;
  double d;
  for (i = st; i < n; i ++) {
    d = A[i+1] / B[i];
    B[i+1] -= d * C[i];
    D[i+1] -= d * D[i];
  }
  x[n] = D[n] / B[n];
  for (i = n - 1; i >= st; i --)
    x[i] = (D[i] - C[i] * x[i+1]) / B[i];
}

main()
{
  clrscr();
  int i, n;
  FILE *fp;
  double dt, T[101], Tan[101], ratio,
  Ch = 1.7,
  kh = 12,
  Dh = kh / Ch,
  z_bot = 30,
  eta = 0.7,
  t_end = 72;
  //cout << "Please input ratio of the time increment to the maximum allowed step:"; cin >> ratio;
  ratio = 2;
  n = z_bot / dz;
  dt = ratio * (dz * dz) / (2 * Dh);
  zd = sqrt(2 * Dh / omega);
  get_analytical_solution(T, n + 1);
  //------------------------//
  ofstream fout;
  do {
    do {
      fclose(fp);
      fp = fopen("Ready.flg", "r");
    }
    while(fp == NULL); // "Ready.flg"が見つかるまで読み込み試行
    fclose(fp);
    t += dt;
    get_analytical_solution(Tan, n + 1);
    T[0] = Tan[0];
    T[n+1] = Tan[n+1];
    for (i = 1; i <= n; i ++) {
      A[i] = eta * kh / (dz * dz);
      B[i] = - 2 * eta * kh  / (dz * dz) - Ch / dt;
      C[i] = A[i];
      D[i] = -(1 - eta) * kh  / (dz * dz) * T[i-1]
             + (2 * (1 - eta) * kh  / (dz * dz) - Ch / dt) * T[i]
             -(1 - eta) * kh  / (dz * dz) * T[i+1];
    }
    D[1] -= A[1] * T[0];
    D[n] -= C[n] * T[n+1];
    solve_tridiagobal_matrix(T, 1, n);
    clrscr();
    cout << "t = " << t << endl;
  //  cout << "depth(cm)  Temperature(C)" << endl;
    fout.open("Soil Temperature.csv");
    for (i = 1; i < n; i ++) {
      fout << (i - 0.5) * dz << ", " << T[i] << ", " << Tan[i] << endl; 
    //  cout << (i - 0.5) * dz << ", " << T[i] << ", " << Tan[i] << endl; 
    }
    fout.close();
    fout.open("XYGrapher.msg");
      fout << t;
    fout.close();
    unlink("Ready.flg"); // "Ready.flg"を削除
  } while (t < t_end);
  
  cout << "Press any key to exit.\n";
  getch();
  return 0;
}
