// PROGRAM radiation
// 加速度運動する電荷のつくる電場の計算と電気力線の描画

#include "TrueBASIC.h"
#include "csgraphics.h"

int main();
void initial(double *L, int *nsnap, double *a, int *lpc, double *dt, 
             double *dtheta);
void draw_field(double L, int isnap, double a, double dt, double dtheta);
void motion(double R[], double t_ret, double r_ret[], double v_ret[], 
            double a_ret[]);
void field(double R[], double t, double E[]);
void retarded_time(double R[], double t, double *t_ret);
void fanddf(double *f, double *dfdt, double R[], double t, double t_ret);
void zero_of_f(double *f, double *dfdt, double R[], double t, double *t_ret, 
               double *ta, double *tb);
void animate(int nsnap);
void crossproduct(double a[], double b[], double c[]);
double dotproduct(double a[], double b[]);

int main()
{
    int isnap, lpc, nsnap;
    double L, dt, dtheta, radius;

    GWopen(0);

    initial(&L, &nsnap, &radius, &lpc, &dt, &dtheta);
    for(isnap = 1; isnap <= nsnap; ++isnap)
        draw_field(L, isnap, radius, dt, dtheta);
    animate(nsnap);

    GWquit();
    return 0;
}

void initial(double *L, int *nsnap, double *a, int *lpc, double *dt, 
             double *dtheta)
{
    float xwin, ywin;

    (*L) = 20;
    compute_aspect_ratio((float)(2*(*L)), &xwin, &ywin);
    GWindow(-xwin, -ywin, xwin, ywin);
    (*a) = 0.3;                             // 画面上の電荷の"半径"
    (*lpc) = 10;                            // 電気力線の数
    (*dt) = 0.4;                            // 表示と表示の間の時間
    (*dtheta) = 2*pi/(*lpc);
    (*nsnap) = 15;                          // スナップ・ショットの数
}

void draw_field(double L, int isnap, double a, double dt, double dtheta)
// プログラム fieldline より
{
    double E[4], Emag, R[4], a_ret[4], ds, r_ret[4], t, theta, v_ret[4], x, 
           xline, y, yline;
    int Itmp1_;
    for(Itmp1_ = 0; Itmp1_ < 4; ++Itmp1_)
        E[Itmp1_] = R[Itmp1_] = a_ret[Itmp1_] = r_ret[Itmp1_] = 
        v_ret[Itmp1_] = 0;

    GWsetogn(0, -isnap);
    GWsetpen(-1, -1, -1, -1);
    ds = 1;
    t = isnap*dt;                           // 観測時刻
    R[1] = 0;
    R[2] = 0;
    motion(R, t, r_ret, v_ret, a_ret);
    // 時刻 t の電荷の位置を求める
    x = R[1] - r_ret[1];
    y = R[2] - r_ret[2];
    GWellipse((float)(x-a), (float)(y-a), (float)(x+a), (float)(y+a));
    for(theta = 0; theta <= 2*pi; theta += dtheta)
    {
        xline = x + a*cos(theta);
        yline = y + a*sin(theta);
        GWline2((float)(xline), (float)(yline));
        do
        {
            R[1] = xline;
            R[2] = yline;
            field(R, t, E);                 // 電場の計算
            Emag = sqrt(dotproduct(E, E));
            // 電気力線の新しい位置
            xline += ds*E[1]/Emag;
            yline += ds*E[2]/Emag;
            GWline2((float)(xline), (float)(yline));
        } while(!(fabs(xline) > L || fabs(yline) > L));
        // 描画の停止
    }
    GWrect((float)(-L), (float)(-L), (float)(L), (float)(L));
}

void motion(double R[], double t_ret, double r_ret[], double v_ret[], 
            double a_ret[])
// 電荷の運動を計算
{
    r_ret[1] = R[1] - 0.2*cos(t_ret);       // 調和振動
    r_ret[2] = R[2];
    r_ret[3] = R[3];
    v_ret[1] = -0.2*sin(t_ret);             // 粒子の速度
    v_ret[2] = 0;
    v_ret[3] = 0;
    a_ret[1] = -0.2*cos(t_ret);             // 粒子の加速度
    a_ret[2] = 0;
    a_ret[3] = 0;
}

void field(double R[], double t, double E[])
// 電場のベクトルを計算
{
    int i;
    double E0, a_ret[4], dist_ret, r_ret[4], ru, t_ret, u_ret[4], uxa[4], v2, 
           v_ret[4], w_ret[4];

    retarded_time(R, t, &t_ret);            // 遅延時刻を求める
    // 運動する粒子の力学変数
    motion(R, t_ret, r_ret, v_ret, a_ret);
    v2 = dotproduct(v_ret, v_ret);
    dist_ret = sqrt(dotproduct(r_ret, r_ret));
    for(i = 1; i <= 3; ++i)
        u_ret[i] = r_ret[i]/dist_ret - v_ret[i];
    crossproduct(u_ret, a_ret, uxa);
    crossproduct(r_ret, uxa, w_ret);
    ru = dotproduct(r_ret, u_ret);
    E0 = dist_ret/(ru*ru*ru);
    for(i = 1; i <= 3; ++i)
        E[i] = E0*(u_ret[i]*(1 - v2) + w_ret[i]);
}

void retarded_time(double R[], double t, double *t_ret)
{
    double dfdt, f, fa, fb, ta, tb;

    tb = t;                                 // 遅延時刻の予想値の上限
    fanddf(&fb, &dfdt, R, t, tb);
    ta = -1/1.6;                            // 遅延時刻の予想値の下限
    // f(ta) > 0 と f(tb) < 0 を保証する
    do
    {
        ta *= 1.6;
        fanddf(&fa, &dfdt, R, t, ta);
    } while(!(fa > 0));
    (*t_ret) = 0.5*(ta + tb);
    fanddf(&f, &dfdt, R, t, *t_ret);
    zero_of_f(&f, &dfdt, R, t, t_ret, &ta, &tb);
}

void fanddf(double *f, double *dfdt, double R[], double t, double t_ret)
// f と df/dt_r を計算する
{
    double a_ret[4], dist_ret, r_ret[4], v_ret[4];

    motion(R, t_ret, r_ret, v_ret, a_ret);
    dist_ret = sqrt(dotproduct(r_ret, r_ret));
    (*f) = t - t_ret - dist_ret;
    // 遅延時刻における微分の計算
    (*dfdt) = -1 + dotproduct(r_ret, v_ret)/dist_ret;
}

void zero_of_f(double *f, double *dfdt, double R[], double t, double *t_ret, 
               double *ta, double *tb)
// f = t - t_ret - dist_ret = 0 となる時刻 t_ret を求めるために
// 100回以上の反復は行わない
{
    double dt, dt_old, dt_r, eps, j, sign;

    eps = 1e-6;
    dt_r = (*tb) - (*ta);
    for(j = 1; j <= 100; ++j)
    {
        // ニュートン法における前の t_ret の計算値との差
        dt = (*f)/(*dfdt);
        sign = (((*t_ret) - (*ta)) - dt)*(((*t_ret) - (*tb)) - dt);
        dt_old = dt_r;
        if(sign >= 0 || fabs(2*dt) > fabs(dt_old))
        {
            // ニュートン法の次の反復計算の結果が ta と tb の間に存在しないか
            // dt が dt_r の前の値の半分より小さくないときには二分法を使う
            dt_r = 0.5*((*ta) - (*tb));
            (*t_ret) = (*tb) + dt_r;
        }
        else                      // ニュートン法を使う
        {
            dt_r = dt;
            (*t_ret) = (*t_ret) - dt_r;
        }
        if(fabs(dt_r) < eps)      // 収束性のテスト
            return;
        fanddf(f, dfdt, R, t, *t_ret);
        if((*f) < 0)
            (*tb) = (*t_ret);
        else
            (*ta) = (*t_ret);
    }
    printf("too many iterations\n");
    GWquit();
    exit(0);
}

void animate(int nsnap)
{
    int isnap = 1, lsnap, loop = 10, i;

    GWflush(isnap);
    for(i = 2; i <= nsnap*loop ; ++i) // 動画として表示
    {
        Sleep(100);
        lsnap = isnap;
        if(++isnap > nsnap) isnap = 1;
        GWsetogn(lsnap, -lsnap);
        GWflush(isnap);
    }
}

void crossproduct(double a[], double b[], double c[])
{
    c[1] = a[2]*b[3] - a[3]*b[2];
    c[2] = a[3]*b[1] - a[1]*b[3];
    c[3] = a[1]*b[2] - a[2]*b[1];
}

double dotproduct(double a[], double b[])
{
    double dotproduct_;

    dotproduct_ = a[1]*b[1] + a[2]*b[2] + a[3]*b[3];
    return dotproduct_;
}