// PROGRAM fieldline
// 2 次元の電気力線を描く

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

int main();
void screen(double *a);
void charges(int *N, float x[], float y[], double q[], double *qtotal);
void draw_charges(int N, float x[], float y[], double q[]);
void draw_lines(int N, float x[], float y[], double q[], double qtotal, double a);

int main()
{
    int N;
    float x[11], y[11];
    double a, q[11], qtotal = 0;

    GWopen(0);

    screen(&a);
    charges(&N, x, y, q, &qtotal);
    // 電気力線の描画
    draw_lines(N, x, y, q, qtotal, a);

    GWquit();
    return 0;
}

void screen(double *a)
{
    float xwin, ywin;
    int L;

    L = 10;
    compute_aspect_ratio((float)L, &xwin, &ywin);
    GWindow(-xwin, -ywin, xwin, ywin);
    (*a) = 0.2;                             // 画面上の電荷の"半径"
    GWsetmrk(6, (float)((*a)*2), -1, -1, -1);
}

void charges(int *N, float x[], float y[], double q[], double *qtotal)
{
    double charge;
    char Stmp1_[_LBUFF_];

    // 電荷の大きさと位置を入力
    (*N) = 0;                               // 点電荷の数
    do
    {
        draw_charges(*N, x, y, q);
        GWinput("charge (0 to exit input mode) = ", Stmp1_, _LBUFF_);
        charge = atof(Stmp1_);
        if(sscanf(Stmp1_, "%lg", &charge) == 0) charge = 0;
        if(charge != 0)
        {
            GWsetmsg("place mouse at charge location and click.");
            (*N) = (*N) + 1;
            GWcappnt(x+(*N), y+(*N), NULL); // 電荷の位置
            q[(*N)] = charge;
            (*qtotal) = (*qtotal) + charge;
        }
    } while(!(charge == 0));
    // 再度電荷を描画
    draw_charges(*N, x, y, q);
}

void draw_charges(int N, float x[], float y[], double q[])
{
    int i;

    for(i = 1; i <= N; ++i)
    {
        if(q[i] > 0)
            GWsetmrk(-1, -1.0f, BLUE, -1, -1);
        else
            GWsetmrk(-1, -1.0f, RED, -1, -1);
        GWputmrk((float)(x[i]), (float)(y[i]));
    }
}

void draw_lines(int N, float x[], float y[], double q[], double qtotal, double a)
// 正電荷から出る単位電荷当たりの電気力線の本数
{
    int i, j, key, lpc, sign;
    double E, E0, Emin, Ex, Ey, ds_big, ds_small, dtheta, dx, dy, r = 0, 
           theta, xline, yline;
    char stop_plot[_LBUFF_];

    lpc = 8;                                // 単位電荷当たりの電気力線の本数
    Emin = 0.01;                            // E < Emin なら描画を停止する
    sign = sgn(qtotal);
    if(sign == 0) sign = 1;
    ds_small = 0.01*sign;
    ds_big = sign;
    for(i = 1; i <= N; ++i)                 // すべての電荷について行う
    {
        if(q[i]*sign > 0)                   // 電気力線が正の電荷から出る
        {
            dtheta = 2*pi/(lpc*fabs(q[i]));
            for(theta = 0; theta < 2*pi; theta += dtheta)
            {
                GWmove2(x[i], y[i]);
                strcpy(stop_plot, "no");
                xline = x[i] + a*cos(theta);
                yline = y[i] + a*sin(theta);
                do
                {
                    Ex = 0;
                    Ey = 0;
                    for(j = 1; j <= N; ++j)
                    {
                        // 点から電荷jまでのx方向の距離
                        dx = xline - x[j];
                        // 点から電荷jまでのy方向の距離
                        dy = yline - y[j];
                        r = sqrt(dx*dx + dy*dy);
                        if(r > 0.9*a)
                        {
                            E0 = q[j]/(r*r*r);
                            Ex = Ex + E0*dx;
                            Ey = Ey + E0*dy;
                        }
                        else
                        {
                            // 電気力線が他の電荷の位置に到達
                            strcpy(stop_plot, "yes");
                        }
                    }
                    E = sqrt(Ex*Ex + Ey*Ey);
                    if(E > Emin || r < 20)
                    {
                        // 電気力線上の新しい位置
                        xline = xline + ds_small*Ex/E;
                        yline = yline + ds_small*Ey/E;
                    }
                    else
                    {
                        // 電気力線上の新しい位置
                        xline = xline + ds_big*Ex/E;
                        yline = yline + ds_big*Ey/E;
                    }
                    GWline2((float)(xline), (float)(yline));
                    if(TBkeyinput())
                    {
                        // 電気力線の描画を停止することができる
                        key = TBgetkey();
                        strcpy(stop_plot, "yes");
                    }
                } while(!(strcmp(stop_plot, "yes") == 0));
                // 描画の停止
            }
        }
    }
}