source: proiecte/SolarSim/Java/SolarSim/src/solarsim/sochacki/Sochacki.java @ 152

package solarsim.sochacki;

import solarsim.renderer.Renderer2D;
import solarsim.common.State;
import solarsim.common.Particle;
import solarsim.common.Vector;
import java.util.LinkedList;
import java.util.List;

public class Sochacki {
    int nt;
    int Np;
    static double AU = 149597870691F;   // Astronomical Unit
    static double D = 86400;            // Astronomical Day
    static double K = 0.01720209895;    // Gauss' gravitational constant
    static double G = K*K;
    static double S = 1.9891E30;        // Astronomical Unit for mass - the mass of the Sun
    //static double Gm = 6.6743E-11;

    double [] mass;
    double [][][] xx;
    double [][][] vv;

    double [][][] aa;
    double [][][] u1;
    double [][][] u2;
    double [][][] u3;

    public Sochacki(int nt, int Np) {
        this.nt = nt;
        this.Np = Np;
        mass = new double[2];
        xx = new double[3][Np][nt];
        vv = new double[3][Np][nt];

        aa = new double[Np][Np][nt];
        u1 = new double[Np][Np][nt];
        u2 = new double[Np][Np][nt];
        u3 = new double[Np][Np][nt];
    }

    void initState() {
        // Sun
        mass[0] = G;
        xx[0][0][0] = 0;
        xx[1][0][0] = 0;
        xx[2][0][0] = 0;
        vv[0][0][0] = 0;
        vv[1][0][0] = 0;
        vv[2][0][0] = 0;

        // Earth
        mass[1] = G / 332840;//S
        xx[0][1][0] = 0.98328989;//AU
        xx[1][1][0] = 0;
        xx[2][1][0] = 0;
        vv[0][1][0] = 0;
        vv[1][1][0] = 30286.7 * D / AU;
        vv[2][1][0] = 0;

        for( int j = 0; j < Np; j++ ) {
            for( int k = 0; k < Np; k++ ) {
                if( j == k )
                    continue;
                
                double s = 0;
                for( int i = 0; i < 3; i++ ) {
                    double dx = (xx[i][j][0] - xx[i][k][0]);
                    s += dx * dx;
                }
                u1[j][k][0] = 1 / Math.sqrt(s);
                u2[j][k][0] = u1[j][k][0] * u1[j][k][0];
                u3[j][k][0] = u2[j][k][0] * u1[j][k][0];

                s = 0;
                for( int i = 0; i < 3; i++ ) {
                    double dx = (xx[i][k][0] - xx[i][j][0]);
                    double dv = (vv[i][k][0] - vv[i][j][0]);
                    s += dx * dv;
                }
                aa[j][k][0] = s;
            }
        }
    }

    void genPoly() {
        double a, b;
        int i, j, k, L, m;

        for( m = 1; m < nt; m++ ) {
            for( j = 0; j < Np; j++ ) {
                for( i = 0; i < 3; i++ ) {
                    xx[i][j][m] = vv[i][j][m-1] / m;
                    a = 0;
                    for( k = 0; k < Np; k++ ) {
                        b = 0;
                        for( L = 0; L <= m-1; L++ ) {
                            b = b + (xx[i][k][L] - xx[i][j][L])*u3[j][k][m-1-L];
                        }
                        a = a + mass[k] * b;
                    }
                    vv[i][j][m] = a / m;
                }
                for( k = 0; k <= j-1; k++ ) {
                    a = 0;
                    for( L = 0; L <= m-1; L++ ) {
                        a = a - u3[j][k][L] * aa[j][k][m-1-L];
                    }
                    u1[j][k][m] = a / m;
                    u1[k][j][m] = a / m;
                    a = 0;
                    for( L = 0 ; L <= m; L++ ) {
                        a = a + u1[j][k][L] * u1[j][k][m-L];
                    }
                    u2[j][k][m] = a;
                    u2[k][j][m] = a;
                    a = 0;
                    b = 0;
                    for( L = 0; L <= m; L++ ) {
                        b = b + u2[j][k][L] * u1[j][k][m-L];
                        for( i = 0; i < 3; i++ ) {
                            a = a + (xx[i][j][L]-xx[i][k][L])*(vv[i][j][m-L]-vv[i][k][m-L]);
                        }
                    }
                    aa[j][k][m] = a;
                    aa[k][j][m] = a;
                    u3[j][k][m] = b;
                    u3[k][j][m] = b;
                }
                aa[j][j][m] = 0;
                u1[j][j][m] = 0;
                u2[j][j][m] = 0;
                u3[j][j][m] = 0;
            }
        }
    }

    double evalPoly( double t, double [] c ) {
        double p = 0;

        for( int i = c.length-1; i >= 0; i-- ) {
            p = p * t + c[i];
        }

        return p * AU;
    }

    State getState( double time ) {
        List<Particle> particles = new LinkedList<Particle>();
        for( int i = 0; i < Np; i++) {
            double x = evalPoly(time, xx[0][i]);
            double y = evalPoly(time, xx[1][i]);
            double z = evalPoly(time, xx[2][i]);
            Vector r = new Vector(x, y, z);

            double vx = evalPoly(time, vv[0][i]);
            double vy = evalPoly(time, vv[1][i]);
            double vz = evalPoly(time, vv[2][i]);
            Vector v = new Vector(vx, vy, vz);

            Particle p = new Particle(mass[i], r, v);
            particles.add(p);
        }

        return new State(particles);
    }

    public static void main(String [] args) {
        Renderer2D renderer = new Renderer2D();
        new Thread(renderer).start();

        Sochacki sch = new Sochacki(100,2);
        sch.initState();
        sch.genPoly();
        for( double t = -228; t < 229; t += 1 ) {
            State s = sch.getState(t);
            //System.out.println(s);
            renderer.drawState(s);
        }
        try {
            Thread.currentThread().sleep(60 * 1000);
        } catch( InterruptedException e ) {
            e.printStackTrace();
        }
        renderer.stop();
    }
}