/*
 * dynamics.c
 */

/**********************************************************************
 **********************************************************************
 *
 *   2 次元の定数係数線形常微分方程式
 *       x'(t) = a x + b y
 *       y'(t) = c x + d y
 *  に初期条件
 *       x(0)=x0
 *       y(0)=y0
 *  を課した常微分方程式の初期値問題を解いて、相図を描く。
 *
 *  このプログラムは次の4つの部分から出来ている。
 *     main()
 *       主プログラム
 *       行列の係数の入力、ウィンドウのオープン等の初期化をした後に、
 *       ユーザーにメニュー形式でコマンドを入力してもらう。
 *       実際の計算のほとんどは他の副プログラムに任せている。
 *     draworbit(x0,y0,h,tlimit)
 *       (x0,y0) を初期値とする解の軌道を描く。
 *       刻み幅を h、追跡時間を tlimit とする。
 *       近似解の計算には Runge-Kutta 法を用いる。
 *     fx(x,y)
 *       微分方程式の右辺の x 成分
 *     fy(x,y)
 *       微分方程式の右辺の y 成分
 **********************************************************************/

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

/* 係数行列 A の成分（common 文により function fx,fy と共有する） */
double a, b, c, d;
/* ウィンドウに表示する範囲（common 文により draworbit と共有する） */
double xleft, xright, ybottom, ytop;

main()
{
    /* 初期値 */
    double x0, y0;
    /* 時間の刻み幅、追跡時間 */
    double h, tlimit, newh, newT;
    /* メニューに対して入力されるコマンドの番号 */
    int cmd;
    /* マウスのボタンの状態 */
    int lbut, mbut, rbut;
    /* ウィンドウに表示する範囲の設定 */
    xleft = -1.0;
    xright = 1.0;
    ybottom = -1.0;
    ytop = 1.0;
    /* 時間刻み幅、追跡時間（とりあえず設定） */
    h = 0.01;
    tlimit = 10.0;
    /* 行列の成分を入力 */
    printf(" a,b,c,d=");
    scanf("%lf %lf %lf %lf", &a, &b, &c, &d);
    /* ウィンドウを開く */
    openpl();
    fspace2(xleft, ybottom, xright, ytop);
    /* x 軸、y 軸を描く */
    fline(xleft, 0.0, xright, 0.0);
    fline(0.0, ybottom, 0.0, ytop);
    xsync();
    /* メイン・ループの入口 */
    while (1) {
        /* メニューを表示して、何をするか、番号で選択してもらう */
        printf("したいことを番号で選んで下さい。\n");
        printf(" -1:メニュー終了, 0:初期値のキー入力, 1:初期値のマウス入力,");
        printf("  2:刻み幅 h,追跡時間 T 変更(現在 h=%7.4f, T=%7.4f)\n",
            h, tlimit);
        scanf("%d", &cmd);
        /* 番号 cmd に応じて、指示された仕事をする */
        if (cmd == 0) {
            /* 初期値の入力 */
            printf(" 初期値 x0,y0=");
            scanf("%lf %lf", &x0, &y0);
            draworbit(x0, y0, h, tlimit);
        } else if (cmd == 1) {
            while (1) {
                printf("マウスの左ボタンで初期値を指定して下さい (右ボタンで中止)\n");
                fmouse(&lbut, &mbut, &rbut, &x0, &y0);
                if (lbut == 1) {
                    printf("(x0,y0)=(%g,%g)\n", x0, y0);
                    draworbit(x0, y0, h, tlimit);
                } else if (mbut == 1) {
                    printf("(x0,y0)=(%g,%g)\n", x0, y0);
                    draworbit(x0, y0, -h, tlimit);
                } else {
                    printf("マウスによる初期値の入力を打ち切ります。\n");
                    break;
                }
            }
        } else if (cmd == 2) {
            /* 時間刻み幅、追跡時間の変更 */
            printf(" 時間刻み幅 h, 追跡時間 T= ");
            scanf("%lf %lf", &newh, &newT);
            if (newh != 0 && newT != 0) {
                h = newh;
                tlimit = newT;
                printf("新しい時間刻み幅 h = %g, 新しい追跡時間 T = %g\n",
                    h, tlimit);
            } else {
                printf(" h=%g, T=%g は不適当です。\n", newh, newT);
            }
        } else if (cmd == -1) {
            /* 終了 -- メイン・ループを抜ける */
            break;
        }
    }
    mkplot("dynamics.plot");
    printf("fplotウィンドウを左ボタンでクリックして下さい\n");
    closepl();
    return 0;
}

/***********************************************************************/
/***********************************************************************/
/*      指示された初期値に対する解軌道を描く */
draworbit(x0, y0, h, tlimit)
double x0, y0, h, tlimit;
{
    double x, y, fx(), fy();
    double k1x, k1y, k2x, k2y, k3x, k3y, k4x, k4y, t;
    /* 時刻を 0 にセットする */
    t = 0.0;
    /* 初期値のセット */
    x = x0;
    y = y0;
    /* 初期点を描く */
    fpoint(x, y);
    /* ループの入口 */
    do {
        /* Runge-Kutta 法による計算 */
        /* k1 の計算 */
        k1x = h * fx(x, y);
        k1y = h * fy(x, y);
        /* k2 の計算 */
        k2x = h * fx(x + k1x / 2.0, y + k1y / 2.0);
        k2y = h * fy(x + k1x / 2.0, y + k1y / 2.0);
        /*  k3 の計算 */
        k3x = h * fx(x + k2x / 2.0, y + k2y / 2.0);
        k3y = h * fy(x + k2x / 2.0, y + k2y / 2.0);
        /*  k4 の計算 */
        k4x = h * fx(x + k3x, y + k3y);
        k4y = h * fy(x + k3x, y + k3y);
        /*  (Xn+1,Yn+1) の計算 */
        x += (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0;
        y += (k1y + 2.0 * k2y + 2.0 * k3y + k4y) / 6.0;
        /* 解軌道を延ばす */
        fcont(x, y);
        /* 時刻を 1 ステップ分進める */
        t += h;
        /* まだ範囲内かどうかチェック */
    }    while (abs(t) <= fabs(tlimit));
    /* サブルーチンを抜ける */
    /* 確実に描画が終るのを待つ */
    xsync();
}

/**********************************************************************/
/**********************************************************************/
/* 微分方程式の右辺のベクトル値関数 f の x 成分 */
double fx(x, y)
double x, y;
{
    return a * x + b * y;
}

/* 微分方程式の右辺のベクトル値関数 f の y 成分 */
double fy(x, y)
double x, y;
{
    return c * x + d * y;
}