/*
 * constlinear2d-eggx.c --- 2次元の定数係数線形常微分方程式
 *
 *   x'(t)=a x+b y       x(0)=x0
 *   y'(t)=c x+d y       y(0)=y0
 *
 * ここでは単振動の方程式 x''(t)=-ω^2 x, x(0)=x0, x'(0)=y0 を一階化した
 * 問題を解く。
 */

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

/* 係数行列の成分 */
double a, b, c, d;

/* 微分方程式の右辺 */
double fx(double t, double x, double y)
{
  return a * x + b * y;
}

double fy(double t, double x, double y)
{
  return c * x + d * y;
}

int main()
{
  int win;
  double omega;
  /* 初期値 */
  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 の用意 */
  win = eggx_gopen(500, 500);
  eggx_gsetbgcolor(win, "black"); eggx_gclr(win);
  xmin = - 1; xmax = 1;
  ymin = - 1; ymax = 1;
  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); 

  /* 問題の決定 (単振動のパラメーター ω=1/2, 初期値 (0,1/2) */
  omega = 0.5;
  a =  0; b = 1; c = - omega * omega; d = 0;
  x0 = 0.0; y0 = 0.5;
  /* 0≦t≦100 の範囲を 1000 等分して計算 */
  Tmax = 100;
  N = 1000;
  h = Tmax / N;

  /* Runge-Kutta 法で計算 */
  t = 0; x = x0; y = y0;
  printf("%f %20.15f %20.15f\n", t, x, y);
  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 = a + (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} での値の表示 */
    printf("%f %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 constlinear2d.png" で再表示可能。 */
  eggx_saveimg(win, 0, xmin, ymin, xmax, ymax,
               "pnmtopng", 256, "constlinear2d.png");
  /* PostScript形式で保存。"ggv constlinear2d.eps" で再表示可能。 */
  eggx_saveimg(win, 0, xmin, ymin, xmax, ymax,
               "convert", 256, "constlinear2d.eps");
  /* キー入力を待つ */
  eggx_ggetch(win);
  /* EGGX の終了 */
  eggx_gclose(win);
  return 0;
}