/*
 * dynamics2-eggx.c --- 2次元の力学系 例題2
 *
 *   x'(t)=y                x(0)=x0
 *   y'(t)=-x+μ(1-x^2)y    y(0)=y0
 *
 * (0,0) という一つの平衡点を持つ
 */

#include <stdio.h>
#include <eggxlib.h>

/* 微分方程式の右辺 */
double mu = 1.0;

double fx(double t, double x, double y)
{
  return y;
}

double fy(double t, double x, double y)
{
  return - x + mu * (1- x * x) * y;
}

int main()
{
  /* 初期値 */
  double x0, y0;
  /* 時間刻み幅, 最終時刻 (0≦t≦Tmax の範囲で解くことにする) */
  double h, Tmax;
  /* 変数 */
  double t, x, y, new_t, new_x, new_y;
  double k1x, k1y, k2x, k2y, k3x, k3y, k4x, k4y;
  /* 図を描く範囲 [xmin,xmax]×[ymin,ymax] */
  double xmin, xmax, ymin, ymax;
  int N, j;
  /* EGGX のウィンドウ */
  int win;
  /* */
  char buf[BUFSIZ];
  /* */
  int verbose = 0;

  /* EGGX の用意 */
  win = eggx_gopen(500, 500);
  eggx_gclr(win);
  printf("xmin, xmax, ymin, ymax: ");
  while (1) {
    fgets(buf, sizeof(buf), stdin);
    if (buf[0] != '#') break;
  }
  sscanf(buf, "%lf %lf %lf %lf", &xmin, &xmax, &ymin, &ymax);
  printf("xmin=%g, xmax=%g, ymin=%g, ymax=%g\n", xmin, xmax, ymin, ymax);
  eggx_window(win, xmin, ymin, xmax, ymax);
  /* 座標軸を緑色で描く */
  eggx_newcolor(win, "green");
  eggx_line(win, xmin, 0.0, PENUP); eggx_line(win, xmax, 0.0, PENDOWN);
  eggx_line(win, 0.0, ymin, PENUP); eggx_line(win, 0.0, ymax, PENDOWN);

  /* 0≦t≦100 の範囲を 1000 等分して計算 */
  N = 1000;

  printf("初期値の入力を終了するならば行頭に \"$\" を入力してください。\n");
  while (1) {
    printf("x0, y0, Tmax: ");
    /* 一行読みファイルの終端や行頭が $ であったら while ループを抜ける */
    if (fgets(buf, sizeof(buf), stdin) == NULL) break;
    if (buf[0] == '$') break;
    /* 行頭が # だったら注釈としてパスする */
    if (buf[0] == '#') {
      printf("注釈行\n");
      continue;
    }
    sscanf(buf, "%lf %lf %lf", &x0, &y0, &Tmax);
    /* Runge-Kutta 法で計算 */
    t = 0; x = x0; y = y0; h = Tmax / N;
    printf("%f %f %f\n", x, y, Tmax);
    eggx_newcolor(win, "red");
    eggx_line(win, x, y, PENUP);
    for (j = 0; j < N; j++) {
      k1x = h * fx(t, x, y);
      k1y = h * fy(t, x, y);
      k2x = h * fx(t + h / 2, x + k1x / 2, y + k1y / 2);
      k2y = h * fy(t + h / 2, x + k1x / 2, y + k1y / 2);
      k3x = h * fx(t + h / 2, x + k2x / 2, y + k2y / 2);
      k3y = h * fy(t + h / 2, x + k2x / 2, y + k2y / 2);
      k4x = h * fx(t + h, x + k3x, y + k3y);
      k4y = h * fy(t + h, x + k3x, y + k3y);
      /* t=t_{j+1} での値の算出 */
      new_t = (j + 1) * h;
      new_x = x + (k1x + 2 * k2x + 2 * k3x + k4x) / 6;
      new_y = y + (k1y + 2 * k2y + 2 * k3y + k4y) / 6;
      /* t_{j+1} での値の表示 */
      if (verbose)
	printf("t=%f (x,y)=(%20.15f,%20.15f)\n", new_t, new_x, new_y);
      /* 値の更新 */
      x = new_x; y = new_y; t = new_t;
      eggx_line(win, x, y, PENDOWN);
    }
  }
  /* PNG形式で保存。"display dynamics2.png" で再表示可能。 */
  eggx_saveimg(win, 0, xmin, ymin, xmax, ymax,
               "pnmtopng", 256, "dynamics2.png");
  /* PostScript形式で保存。"ggv dynamics2.eps" で再表示可能。 */
  eggx_saveimg(win, 0, xmin, ymin, xmax, ymax,
               "convert", 256, "dynamics2.eps");
  printf("プログラムを終了するにはウィンドウ内でキーを入力して下さい。\n");
  /* キー入力を待つ */
  eggx_ggetch(win);
  /* EGGX の終了 */
  eggx_gclose(win);
  return 0;
}