/**************************************************************************
 *       selflinear.c                                                     *
 *       定数係数線形自励系常微分方程式の解のグラフ                       *
 *       (C) Kenrou Adachi 2004/11/22                                     *
 *                                                                        *
 *     dv/dt = Av                                                         *
 *     v = (x, y), A = [[s, t], [u, v]]                                   *
 *     solution  v(t) = v(0) + exp(tA)v(0)                                *
 *                                                                        *
 *  You needs to compile etude_x.c with eggx graphics library             * 
 **************************************************************************/

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

#define RED   2
#define GREEN 3
#define BLUE  4

void message(void);
void drawing(double a, double b, double h, int l,
	     double s, double t, double u, double v);
void axes(void);
void repeat(double *x, double *y,
	    double a, double b, double h,
	    double A, double B, double C, double D);
void inputcoeff(double *s, double *t, double *u, double *v);
void detectcoeff(double s, double t, double u, double v,
                 double *A, double *B, double *C, double *D, double h);

int win;

int main(void)
{
  double a, b, h;
  double s, t, u, v;
  int    m, l;
  int    n;

  /* 挨拶 */

  message();

  /* 係数行列の入力 */
  /*   d  [x]    [ s  t ][x]
      --- [ ]  = [      ][ ]
      dt  [y]    [ u  v ][y]   */

  inputcoeff(&s, &t, &u, &v);

  /* データの設定 */

  a = 6.4;
  b = 4.0;
  m = 3;
  l = 50;

  /* グラフィックの初期化 */

  win = gopen(640, 400);
  gsetbgcolor(win, "black");
  gclr(win);

  /* きざみ幅の設定 */

  h = 2.0 * sqrt((a / 320.0) *(a / 320.0) +(b / 200.0) * (b / 200.0));

  /* 座標軸の描画 */

  axes();

  /* 繰り返し処理 */

  for (n = 0; n < m; n++)
  {
    drawing(a, b, h, l, s, t, u, v); /* 線素の描画 */
  }
  fprintf(stderr, "終了しました。\n");
  fprintf(stderr, "Enter key を押してください。\n\n");
  getchar();
  getchar();
  gclose(win);
  return EXIT_SUCCESS;
}


void drawing(double a, double b, double h, int l,
	     double s, double t, double u, double v)
{
  int    j, count;
  double rnd, e1, e2, e, x[40], y[40], xx[40], yy[40];
  double eh, A, B, C, D;

  /* 描画範囲に４０個の初期点をランダムに生成する。 */

  j = 0;
  for (e1 = -1.0; e1 <= 0.9; e1 += 0.25)
  {
    for (e2 = -1.0; e2 <= 0.9; e2 += 0.4)
    {
      rnd = (double)(rand()) / (double)RAND_MAX;
      *(xx + j) = a * (e1 + 0.25 * rnd);
      rnd = (double)(rand()) / (double)RAND_MAX;
      *(yy + j) = b * (e2 + 0.4 * rnd);
      j++;
    }
  }

  /* 軌道計算 routine を呼出す。*/

  e = 1.0; /* 先ず未来を計算 */
  do
  {
    /* 初期点情報を変数にコピーする。 */
    for (j = 0; j <= 39; j++)
    {
      *(x + j) = *(xx + j);
      *(y + j) = *(yy + j);
    }
    /* 描画色の変更 */
    if (e > 0)
    {
      newpen(win, GREEN);  /* 未来 */
    } 
    else 
    {
      newpen(win, RED);  /* 過去 */
    }

    /* Ah + (Ah)^2/2 + (Ah)^3/6 + (Ah)^4/24 + (Ah)^5/120 */

    eh = e * h;
    detectcoeff(s, t, u, v, &A, &B, &C, &D, eh);

    /* 軌道計算 */
    for (count = 0; count < l; count++)
    {
      repeat(x, y, a, b, h, A, B, C, D);
    }
    e = -e; /* 時間を逆転して過去を計算 */
  } while (e < 0);
}


void repeat(double *x, double *y,
	    double a, double b, double h,
            double A, double B, double C, double D)
{
  register int j;
  double x0, y0, r0, wx, wy, wx1, wy1;

  for (j = 0; j <= 39; j++)
  {
    /* j 番目の初期点 */

    x0 = *(x + j);
    y0 = *(y + j);

    if ((fabs(x0) <= 2.0 * a) && (fabs(y0) <= 2.0 * b))
    {
      r0 = sqrt(x0 * x0 + y0 * y0);

      /* x(t) = x(0) + exp(tA)x(0) */

      if (r0 >= h)
      {
	wx1 = *(x + j) = x0 + (A * x0 + B * y0);
	wy1 = *(y + j) = y0 + (C * x0 + D * y0);

	/* Window に描画範囲を割当るように座標変換する */

	wx = 320.0 * x0 / a;
	wy = 200.0 * y0 / b;
	wx1 *= 320.0 / a;
	wy1 *= 200.0 / b;

	/* 線素の描画 */

        line(win, (float)wx + 320.0, (float)(-wy) + 200.0, PENUP); 
	line(win, (float)wx1 + 320.0, (float)(-wy1) + 200.0, PENDOWN);
      }
    }
  }
}


/* 座標軸の描画 */

void axes(void)
{
  newpen(win, BLUE);
  line(win, 0.0, 200.0, PENUP);
  line(win, 639.0, 200.0, PENDOWN);
  line(win, 320.0, 0.0, PENUP);
  line(win, 320.0, 399.0, PENDOWN);
}


/* 係数行列の入力 */

void inputcoeff(double *s, double *t, double *u, double *v)
{
  fprintf(stderr, " 係数行列 A = [[s, t], [u, v]] を入力して下さい。\n");
  fprintf(stderr, "\n s = ? ");
  scanf("%lf", s);
  fprintf(stderr, "\n t = ? ");
  scanf("%lf", t);
  fprintf(stderr, "\n u = ? ");
  scanf("%lf", u);
  fprintf(stderr, "\n v = ? ");
  scanf("%lf", v);
}

/* 挨拶 */

void message(void)
{
  fprintf(stderr, " このプログラムは2次元定数係数線形自励系微分方程式:\n");
  fprintf(stderr, "\t dx/dt = A * x \n の解軌道群を描画するものです。\n");
  fprintf(stderr, " 計算には解公式 x(t) = x(0) + exp(t * A) * x(0) の近似式を使います。\n");
}


void detectcoeff(double s, double t, double u, double v,
                 double *A, double *B, double *C, double *D, double h)
{
  double s2, t2, u2, v2, s3, t3, u3, v3, s4, t4, u4, v4, s5, t5, u5, v5;

  /* dx/dt 
     = x + hAx + ((hA)^2/2)x + ((hA)^3/6)x + ((hA^4)/24)x + ((hA)^5/120)x */

  /* A^2 */

  s2 = s * s + t * u;
  t2 = s * t + t * v;
  u2 = u * s + v * u;
  v2 = u * t + v * v;

  /* A^3 */

  s3 = s2 * s + t2 * u;
  t3 = s2 * t + t2 * v;
  u3 = u2 * s + v2 * u;
  v3 = u2 * t + v2 * v; 

  /* A^4 */

  s4 = s2 * s2 + t2 * u2;
  t4 = s2 * t2 + t2 * v2;
  u4 = u2 * s2 + v2 * u2;
  v4 = u2 * t2 + v2 * v2;

  /* A^5 */

  s5 = s2 * s3 + t2 * u3;
  t5 = s2 * t3 + t2 * v3;
  u5 = u2 * s3 + v2 * u3;
  v5 = u2 * t3 + v2 * v3; 

  /* exp(hA) - I の Maclaurin 近似 */
  
  *A = h * (s + h * (s2 / 2.0 + h * (s3 / 6.0 + h * (s4 / 24.0 + h * s5 / 120.0))));
  *B = h * (t + h * (t2 / 2.0 + h * (t3 / 6.0 + h * (t4 / 24.0 + h * t5 / 120.0))));
  *C = h * (u + h * (u2 / 2.0 + h * (u3 / 6.0 + h * (u4 / 24.0 + h * u5 / 120.0))));
  *D = h * (v + h * (v2 / 2.0 + h * (v3 / 6.0 + h * (v4 / 24.0 + h * v5 / 120.0))));

}