//** calculate a orbits of the planets copyright Adam Treat 2/29/00 **/
public class Orbits2 implements Differential
{
public static void main(String argv[])
{
double x[] = new double[44]; // creating arrays
x[0] = 0; // suns x position
x[1] = 0; // suns y position
x[2] = 0; // suns vx velocity
x[3] = 0; // suns vy velocity
x[4] = 0; // eart x position
x[5] = 1; // eart y position
x[6] = 1.72125e - 2; // eart vx velocity in AU*days^-1
x[7] = 0; // eart vy velocity
x[8] = 0; // moon x position
x[9] = 1 + 2.55e - 3; // moon y position
x[10] = 1.72125e - 2 + 5.90680e - 4; // moon vx velocity in AU*days^-1
x[11] = 0; // moon vy velocity
x[12] = 0; // merc x position
x[13] = -0.3870975; // merc y position
x[14] = 2.77e - 2; // merc vx velocity in AU*days^-1
x[15] = 0; // merc vy velocity
x[16] = 0; // venus x position
x[17] = -0.7233235; // venus y position
x[18] = 2.02e - 2; // venus vx velocity in AU*days^-1
x[19] = 0; // venus vy velocity
x[20] = 0; // mars x position
x[21] = -1.5237131; // mars y position
x[22] = 1.39e - 2; // mars vx velocity in AU*days^-1
x[23] = 0; // mars vy velocity
x[24] = 0; // jup x position
x[25] = 5.202427; // jup y position
x[26] = 7.55e - 3; // jup vx velocity in AU*days^-1
x[27] = 0; // jup vy velocity
x[28] = 0; // sat x position
x[29] = -9.561943; // sat y position
x[30] = 5.57e - 3; // sat vx velocity in AU*days^-1
x[31] = 0; // sat vy velocity
x[32] = 0; // ura x position
x[33] = 19.30425; // ura y position
x[34] = 3.92e - 3; // ura vx velocity in AU*days^-1
x[35] = 0; // ura vy velocity
x[36] = 0; // nep x position
x[37] = 30.27750; // nep y position
x[38] = 3.13e - 3; // nep vx velocity in AU*days^-1
x[39] = 0; // nep vy velocity
x[40] = 0; // plu x position
x[41] = -39.8067; // plu position
x[42] = 2.73e - 3; // plu vx velocity in AU*days^-1
x[43] = 0; // plu vy velocity
Orbits2 traj = new Orbits2();
double timestop = Integer.parseInt(argv[0]);// time in days
int steps = (int)(timestop * (10000 / 365)); // makes the plot accurate to about a meter
double dt = timestop / steps; // time step
double ts = 0;
RungaKutta rk = new RungaKutta(traj);
for (int time = 1; time < timestop; time++)
{
rk.solve(x, ts, time, dt);
ts = time;
//System.out.println(time+" "+x[10]+" "+x[11]+" ");
System.out.println(time + " "+x[0] + " "+x[1] + " "+x[4] + " "+x[5] + " "+
x[8] + " "+x[9] + " "+x[12] + " "+x[13] + " "+x[16] + " "+x[17] + " "+
x[20] + " "+x[21] + " "+x[24] + " "+x[25] + " "+x[28] + " "+x[29] + " "+
x[32] + " "+x[33] + " "+x[36] + " "+x[37] + " "+x[40] + " "+x[41] + " ");
}
}
public void deriv(double t, double x[], double dx[])
{
int n = 11; // number of planets
double G = 8.89704e - 10; // gravitational constant in AU^3*em^-1*days^-2
double m[] = new double[n];
m[0] = 3.33e5; // mass sun in em
m[1] = 1; // mass earth in em
m[2] = 1.23e - 2; // mass moon in em
m[3] = 5.52e - 2; // mass merc in em
m[4] = 8.14e - 1; // mass venu in em
m[5] = 1.07e - 1; // mass mars in em
m[6] = 318; // mass jup in em
m[7] = 95.2; // mass sat in em
m[8] = 14.5; // mass ura in em
m[9] = 17.1; // mass nep in em
m[10] = 2.16e - 3; // mass plu in em
//accel( int u, int n, int pn, double G, double m[], double x[])
Acceleration ac = new Acceleration();
dx[0] = x[2];
dx[1] = x[3];
dx[2] = ac.accel(0, n, 0, G, m, x);
dx[3] = ac.accel(1, n, 0, G, m, x);
dx[4] = x[6];
dx[5] = x[7];
dx[6] = ac.accel(0, n, 1, G, m, x);
dx[7] = ac.accel(1, n, 1, G, m, x);
dx[8] = x[10];
dx[9] = x[11];
dx[10] = ac.accel(0, n, 2, G, m, x);
dx[11] = ac.accel(1, n, 2, G, m, x);
dx[12] = x[14];
dx[13] = x[15];
dx[14] = ac.accel(0, n, 3, G, m, x);
dx[15] = ac.accel(1, n, 3, G, m, x);
dx[16] = x[18];
dx[17] = x[19];
dx[18] = ac.accel(0, n, 4, G, m, x);
dx[19] = ac.accel(1, n, 4, G, m, x);
dx[20] = x[22];
dx[21] = x[23];
dx[22] = ac.accel(0, n, 5, G, m, x);
dx[23] = ac.accel(1, n, 5, G, m, x);
dx[24] = x[26];
dx[25] = x[27];
dx[26] = ac.accel(0, n, 6, G, m, x);
dx[27] = ac.accel(1, n, 6, G, m, x);
dx[28] = x[30];
dx[29] = x[31];
dx[30] = ac.accel(0, n, 7, G, m, x);
dx[31] = ac.accel(1, n, 7, G, m, x);
dx[32] = x[34];
dx[33] = x[35];
dx[34] = ac.accel(0, n, 8, G, m, x);
dx[35] = ac.accel(1, n, 8, G, m, x);
dx[36] = x[38];
dx[37] = x[39];
dx[38] = ac.accel(0, n, 9, G, m, x);
dx[39] = ac.accel(1, n, 9, G, m, x);
dx[40] = x[42];
dx[41] = x[43];
dx[42] = ac.accel(0, n, 10, G, m, x);
dx[43] = ac.accel(1, n, 10, G, m, x);
//System.out.println(dx[2]+" "+dx[3]+" ");
}
}
//** Runga Kutta copyright Adam Treat 2/29/00 **/
public class RungaKutta
{
Differential eqns;
public RungaKutta(Differential eqns)
{
this.eqns = eqns;
}
public void solve(double x[], double tstart, double tfinish, double dt)
{
int n = x.length;
double k1[] = new double [n];
double k2[] = new double [n];
double k3[] = new double [n];
double k4[] = new double [n];
double tx[] = new double [n];
double t = tstart;
double tt;
while ((t + dt / 2) < tfinish)
{
eqns.deriv(t, x, k1);
tt = t + dt / 2;
for (int i = 0; i < x.length; i++)
{
tx[i] = x[i] + k1[i] * dt / 2;
}
eqns.deriv(tt, tx, k2);
tt = t + dt / 2;
for (int i = 0; i < x.length; i++)
{
tx[i] = x[i] + k2[i] * dt / 2;
}
eqns.deriv(tt, tx, k3);
t = t + dt;
for (int i = 0; i < x.length; i++)
{
tx[i] = x[i] + k3[i] * dt;
}
eqns.deriv(t, tx, k4);
for (int i = 0; i < x.length; i++)
{
x[i] = x[i] + (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6.*dt;
}
}
}
}
//** Differential copyright Adam Treat 2/29/00 **/
interface Differential
{
void deriv(double t, double x[], double dx[]);
}
//** Acceleration due to Gravity in three dimensions... copyright Adam Treat 2/29/00 **/
public class Acceleration
{
public double accel(int u, int n, int pn, double G, double m[], double x[])
{
double kxx = 0, kyy = 0, kx = 0, ky = 0;
for (int i = 0; i < n; i++)
{
if (i != pn)
{
double xdist = x[pn * 4] - x[i * 4];
double ydist = x[(pn * 4) + 1] - x[(i * 4) + 1];
double r = Math.sqrt((xdist * xdist) + (ydist * ydist));
kx = kxx + (G * m[i] * xdist / (r * r * r));
kxx = kx;
ky = kyy + (G * m[i] * ydist / (r * r * r));
kyy = ky;
}
else
{}
}
double accelx = -1 * kx;
double accely = -1 * ky;
if (u == 0)
{
return accelx;
}
else
{
return accely;
}
}
}
|