/***************************************************************************
 *       etpro.c                                                           *
 *       自励系常微分方程式の解のグラフ                                    *
 *       I ported etpro .c (for Turbo C) to Unix X's environment.          *
 *       (C) Kenrou Adachi 2004                                            *
 *                                                                         *
 *                                                                         *
 *                                                                         *
 *     dx/dt = fx(x,y), dy/dt = gy(x,y)                                    *
 *                                                                         *
 *                                                                         *
 * cc -O2 -o etpro etpro.c func.c -I/usr/X11R6/include \                   *
 * -I/usr/local/include -L/usr/X11R6/lib -L/usr/local/lib -lX11 -leggx -lm *
 *                                                                         *
 *  You need to compile etpro.c with eggx graphics library                 * 
 ***************************************************************************/

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

void data_in(double *a, double *b, double *eps, int *m, int *l);
void drawing(double a, double b, double eps, double h, int l);
void axes(void);
void repeat(double *x, double *y,
	    double a, double b, double eps, double h, double e);
extern double fx(double x, double y, double r);
extern double gy(double x, double y, double r);

int win;

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


int main(void)
{
  double a, b, eps, h;
  int    m, l;
  int    n;

  /* データの設定 */

  data_in(&a, &b, &eps, &m, &l);

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

  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, eps, h, l); /* 線素の描画 */
  }
  fprintf(stderr, "終了しました。\n");
  fprintf(stderr, "Enter key を押してください。\n\n");
  getchar();
  getchar();
  gclose(win);
  return EXIT_SUCCESS;
}


void drawing(double a, double b, double eps, double h, int l)
{
  int    j, count;
  double rnd, e1, e2, e, x[40], y[40], xx[40], yy[40];

  /* 描画範囲に４０個の初期点をランダムに生成する。 */
  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);  /* 過去 */
    }

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


void repeat(double *x, double *y,
	    double a, double b, double eps, double h, double e)
{
  register int j;
  double x0, y0, xr0, x1, y1, xr1, wx, wy, wx1, wy1;
  double u, v, u1, v1, q, q_h;

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

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

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

      /* 近似計算 */

      if (xr0 >= h)
      {
	u = fx(x0, y0, xr0);
        v = gy(x0, y0, xr0);
	q = sqrt(u * u + v * v);
	if (q >= eps)
	{
	  q_h = h / q;
	  x1 = x0 + 0.5 * e * u * q_h;
	  y1 = y0 + 0.5 * e * v * q_h;
	  xr1 = sqrt(x1 * x1 + y1 * y1);
	  u1 = fx(x1, y1, xr1);
	  v1 = gy(x1, y1, xr1);

	  /* j 番目の点を少し動かす */

	  *(x + j) = x0 + e * u1 * q_h;
	  *(y + j) = y0 + e * v1 * q_h;
	  
	  /* Window に描画範囲を割当るように座標変換する */

	  wx = 320.0 * x0 / a;
	  wy = 200.0 * y0 / b;
	  wx1 = 320.0 * (*(x + j)) / a;
	  wy1 = 200.0 * (*(y + j)) / 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 data_in(double *a, double *b, double *eps, int *n, int *l)
{
   fprintf(stderr, "描画範囲を設定して下さい。\n");
   fprintf(stderr, " x 座標 = [-a, a]: a = ? ");
   scanf("%lf", a);
   fprintf(stderr, " y 座標 = [-b, b]: b = ? ");
   scanf("%lf", b);
   fprintf(stderr, "計算限界精度(0.01~0.1): eps = ? "); 
   scanf("%lf", eps);
   fprintf(stderr, "描画回数 n は、n = ? ");
   scanf("%d", n);
   fprintf(stderr, "計算処理回数は、l = ? ");
   scanf("%d", l);
}


/* 場の関数 */

/* double fx(double x, double y, double r)
 * {
 *   return (-y + 5.0 * x * exp(-1.0 / sqrt(r)) * sin(1.0 / r));
 * }
 *
 *
 * double gy(double x, double y, double r)
 * {
 *   return (x + 5.0 * y * exp(-1.0 / sqrt(r)) * sin(1.0 / r));
 * }
 */