/*
 * euler-graph-eggx.c -- Euler 法の例の解曲線を描く (EGGXバージョン)
 *
 * 微分方程式の初期値問題
 *    x'(t)=f(t,x), x(a)=xi
 * を Euler 法で解いて図示する。
 */

#include <stdio.h>
#include <math.h>   /* exp() を使うため */
#include <eggxlib.h>  /* EGGX */

/* 微分方程式の右辺に現われる関数の宣言 */
double f(double, double);

int main()
{
  int win;
  int j, N;
  double a, b, xi, h, t, x, new_t, new_x;
  /* 描画範囲 */
  double xmin, xmax, ymin, ymax;

  /* 問題の設定 */
  a = 0; b = 1; xi = 1;
  /* 刻み幅の決定 */
  printf("区間を何等分するか: "); scanf("%d", &N);
  h = (b - a) / N;

  /* EGGXの初期化、画面のクリア、座標系の設定 */
  win = eggx_gopen(500, 500);
  eggx_gsetbgcolor(win, "black"); eggx_gclr(win);
  xmin = a - (b - a) / 10; xmax = b + (b - a) / 10;
  ymin = - 0.2; ymax = 3.2;
  eggx_window(win, xmin, ymin, xmax, ymax);
  /* 座標軸を描く */
  eggx_newcolor(win, "cyan"); /* CYAN = 水色 */
  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); 

  /* 比較用に厳密解 */
  eggx_newcolor(win, "green");
  t = a; x = exp(a);
  eggx_line(win, t, x, PENUP);
  for (j = 1; j <= N; j++) {
    t = a + j * h;
    x = exp(t);
    printf("%f %f\n", t, x);
    eggx_line(win, t, x, PENDOWN);
  }
  eggx_drawstr(win, 0.2, 2.5, 16, 0, "exact solution");

  /* 初期値の設定 */
  t = a; x = xi;
  printf("%f %20.15f\n", t, x);

  /* 線の色を赤にする */
  eggx_newcolor(win, "red");
  /* (t_0,x_0) に移動 */
  eggx_line(win, t, x, PENUP);

  /* 繰り返し */
  for (j = 0; j < N; j++) {
    new_x = x + h * f(t, x);
    new_t = a + (j + 1) * h;
    /* 値の出力 */
    printf("%f %20.15f\n", new_t, new_x);
    /* 値の更新 */
    x = new_x; t = new_t;
    /* 現在点から (t_{j+1}, x_{j+1}) まで線を引くj*/
    eggx_line(win, t, x, PENDOWN);
  }
  eggx_drawstr(win, 0.2, 2.8, 16, 0, "Euler method");

  /* PNG形式で保存。"display euler.png" で再表示可能。 */
  eggx_saveimg(win, 0, xmin, ymin, xmax, ymax,
               "pnmtopng", 256, "euler.png");
  /* PostScript形式で保存。"ggv euler.eps" で再表示可能。 */
  eggx_saveimg(win, 0, xmin, ymin, xmax, ymax,
               "convert", 256, "euler.eps");
  /* ウィンドウ内でのキー入力を待つ */
  eggx_ggetch(win);
  /* EGGX の終了 */
  eggx_gclose(win);
  return 0;
}

double f(double t, double x)
{
  return x;
}