// PROGRAM planet2
// d = 2 solar system with major and minor planet

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

int main();
void initial(double x[], double y[], double vx[], double vy[], double *t, 
             double *GM, double ratio[], double *dt, int *nshow);
void Euler(double x[], double y[], double vx[], double vy[], double *t, 
           double GM, double ratio[], double dt);
void output(double x[], double y[], double t);

int main()
{
    double GM, dt, ratio[3], t, vx[3], vy[3], x[3], y[3];
    int counter, nshow;

    GWopen(0);

    initial(x, y, vx, vy, &t, &GM, ratio, &dt, &nshow);
    output(x, y, t);
    counter = 0;
    do
    {
        Euler(x, y, vx, vy, &t, GM, ratio, dt);
        if((++counter % nshow) == 0)
            output(x, y, t);
    } while(!TBkeyinput());

    GWquit();
    return 0;
}

void initial(double x[], double y[], double vx[], double vy[], double *t, 
             double *GM, double ratio[], double *dt, int *nshow)
{
    float xwin, ywin;

    (*t) = 0;
    (*GM) = 4.0*pi*pi;                      // gravitational constant
    (*dt) = 0.001;                          // time step (yrs)
    (*nshow) = 20;                          // # of iterations between output
    // planet one initial conditions
    x[1] = 2.52;
    y[1] = 0;
    vx[1] = 0;
    vy[1] = sqrt((*GM)/x[1]);
    // mass of planet 2 is 0.04 mass of sun 
    ratio[1] = 0.04*(*GM);
    // planet two initial conditions
    x[2] = 5.24;
    y[2] = 0;
    vx[2] = 0;
    vy[2] = sqrt((*GM)/x[2]);
    // mass of planet 1 is 0.001 mass of sun 
    ratio[2] = -0.001*(*GM);
    // set up window   
    compute_aspect_ratio((float)(2*x[2]), &xwin, &ywin);
    GWindow(-xwin, -ywin, xwin, ywin);
}

void Euler(double x[], double y[], double vx[], double vy[], double *t, 
           double GM, double ratio[], double dt)
{
    int i;
    double ax, ay, dr, dr2, dr3, dx, dy, r, r2, r3;

    // compute distance between planets 1 and 2
    dx = x[2] - x[1];
    dy = y[2] - y[1];
    dr2 = dx*dx + dy*dy;
    dr = sqrt(dr2);
    dr3 = dr2*dr;
    for(i = 1; i <= 2; ++i)                 // sum over planets
    {
        r2 = x[i]*x[i] + y[i]*y[i];
        r = sqrt(r2);
        r3 = r2*r;
        ax = -GM*x[i]/r3 + ratio[i]*dx/dr3;
        ay = -GM*y[i]/r3 + ratio[i]*dy/dr3;
        vx[i] += ax*dt;
        vy[i] += ay*dt;
        x[i] += vx[i]*dt;
        y[i] += vy[i]*dt;
    }
    (*t) += dt;
}

void output(double x[], double y[], double t)
// plot orbits
{
    if(t == 0)
    {
        GWsetpen(RED, -1, -1, -1);
        GWsetmrk(6, 0.2f, RED, -1, 4);
        GWputmrk(0.0f, 0.0f);   // sun at origin
    }
    GWsetpxl((float)(x[1]), (float)(y[1]), BLUE);   // planet one
    GWsetpxl((float)(x[2]), (float)(y[2]), GREEN);  // planet two
}