// PROGRAM magnetism

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

int main();
void wire(int *N, double *a, double *delta, double x[], double y[], double z[], 
          double dLx[], double dLy[], double dLz[]);
void screen(double a);
void draw(int N, double delta, double x[], double y[], double z[], double dLx[], 
          double dLy[], double dLz[], float xcursor, float ycursor);

int main()
{
    int N;
    float xcursor, ycursor;
    double a, dLx[51], dLy[51], dLz[51], delta, x[51], y[51], z[51];

    GWopen(0);

    wire(&N, &a, &delta, x, y, z, dLx, dLy, dLz);
    screen(a);
    do
    {
        GWcappnt(&xcursor, &ycursor, NULL);
        draw(N, delta, x, y, z, dLx, dLy, dLz, xcursor, ycursor);
    } while(!((fabs(xcursor) < .1) && (fabs(ycursor) < .1)));

    GWquit();
    return 0;
}

void wire(int *N, double *a, double *delta, double x[], double y[], double z[], 
          double dLx[], double dLy[], double dLz[])
// 電流セグメントの分割数の読み込み
// 単位電流を仮定する
{
    int i;
    double angle, dL, delta_angle;
    char Stmp1_[_LBUFF_];

    printf("radius of loop = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", a);
    printf("number of segments = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", N);
    // 磁力線の分割の大きさ
    printf("delta = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", delta);
    delta_angle = 2*pi/(*N);
    angle = 0;
    dL = 2*pi*(*a)/(*N);
    for(i = 1; i <= (*N); ++i)
    {
        x[i] = (*a)*cos(angle);
        y[i] = 0;
        z[i] = (*a)*sin(angle);
        // (x,y,z) における電流セグメントの向き
        dLx[i] = -dL*sin(angle);
        dLy[i] = 0;
        dLz[i] = dL*cos(angle);
        angle += delta_angle;
    }
}

void screen(double a)
// スクリーン幅の 1/10 を距離の単位にとる
{
    float xwin, ywin;

    compute_aspect_ratio((float)(3*a), &xwin, &ywin);
    GWindow(-xwin, -ywin, xwin, ywin);
    GWellipse((float)(-a), (float)(-0.4), (float)(a), (float)(0.4));
    GWsetmrk(1, 0.2f, BLUE, -1, -1);
    GWputmrk(0.0f, 0.0f);
    printf("Click with mouse to start line\n");
    printf("Hit any key to stop line\n");
    printf("Click on filled box and hit any key to end\n");
}

void draw(int N, double delta, double x[], double y[], double z[], double dLx[], 
          double dLy[], double dLz[], float xcursor, float ycursor)
{
    int i, kk;
    double B, B0, Bx, By, Bz, dx, dy, dz, r, rx, ry, rz = 0, zcursor = 0;

    rx = xcursor;
    ry = ycursor;
    GWmove2((float)(rx), (float)(ry));
    do
    {
        Bx = 0.0;
        By = 0.0;
        Bz = 0.0;
        for(i = 1; i <= N; ++i)
        {
            // 電流セグメントまでの変位の計算
            dx = rx - x[i];
            dy = ry - y[i];
            dz = zcursor - z[i];
            r = sqrt(dx*dx + dy*dy + dz*dz);
            B0 = 1/(r*r*r);                 // 単位電流を仮定する
            // dL x r に比例する B
            Bx += B0*(dLy[i]*dz - dLz[i]*dy);
            By += B0*(dLz[i]*dx - dLx[i]*dz);
            Bz += B0*(dLx[i]*dy - dLy[i]*dx);
        }
        B = sqrt(Bx*Bx + By*By + Bz*Bz);
        rx += delta*Bx/B;                   // 磁力線上の新しい位置
        ry += delta*By/B;
        rz += delta*Bz/B;
        GWline2((float)(rx), (float)(ry));
    } while(!TBkeyinput());
    kk = TBgetkey();                        // キー入力の初期化
}