// PROGRAM fieldline
// draw electric field lines in two dimensions

#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 field lines
    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;                             // "radius" of visual image of charges
    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_];

    // input charge values and location
    (*N) = 0;                               // # of point charges
    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); // location of charge
            q[(*N)] = charge;
            (*qtotal) = (*qtotal) + charge;
        }
    } while(!(charge == 0));
    // redraw charges
    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)
// number of lines per unit charge leaving positive charges
{
    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;                                // lines per charge
    Emin = 0.01;                            // if E < Emin stop drawing fieldline
    sign = sgn(qtotal);
    if(sign == 0) sign = 1;
    ds_small = 0.01*sign;
    ds_big = sign;
    for(i = 1; i <= N; ++i)                 // loop over all charges
    {
        if(q[i]*sign > 0)                   // start fields at positive charges
        {
            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)
                    {
                        // x-distance from point to charge j
                        dx = xline - x[j];
                        // y-distance from point to charge j
                        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
                        {
                            // field line reached another charge
                            strcpy(stop_plot, "yes");
                        }
                    }
                    E = sqrt(Ex*Ex + Ey*Ey);
                    if(E > Emin || r < 20)
                    {
                        // new position on fieldline
                        xline = xline + ds_small*Ex/E;
                        yline = yline + ds_small*Ey/E;
                    }
                    else
                    {
                        // new position on fieldline
                        xline = xline + ds_big*Ex/E;
                        yline = yline + ds_big*Ey/E;
                    }
                    GWline2((float)(xline), (float)(yline));
                    if(TBkeyinput())
                    {
                        // user can stop drawing field line
                        key = TBgetkey();
                        strcpy(stop_plot, "yes");
                    }
                } while(!(strcmp(stop_plot, "yes") == 0));
                // turn beam off
            }
        }
    }
}