[152] | 1 | package solarsim.sochacki; |
---|
| 2 | |
---|
| 3 | import solarsim.renderer.Renderer2D; |
---|
| 4 | import solarsim.common.State; |
---|
| 5 | import solarsim.common.Particle; |
---|
| 6 | import solarsim.common.Vector; |
---|
| 7 | import java.util.LinkedList; |
---|
| 8 | import java.util.List; |
---|
| 9 | |
---|
| 10 | public class Sochacki { |
---|
| 11 | int nt; |
---|
| 12 | int Np; |
---|
| 13 | static double AU = 149597870691F; // Astronomical Unit |
---|
| 14 | static double D = 86400; // Astronomical Day |
---|
| 15 | static double K = 0.01720209895; // Gauss' gravitational constant |
---|
| 16 | static double G = K*K; |
---|
| 17 | static double S = 1.9891E30; // Astronomical Unit for mass - the mass of the Sun |
---|
| 18 | //static double Gm = 6.6743E-11; |
---|
| 19 | |
---|
| 20 | double [] mass; |
---|
| 21 | double [][][] xx; |
---|
| 22 | double [][][] vv; |
---|
| 23 | |
---|
| 24 | double [][][] aa; |
---|
| 25 | double [][][] u1; |
---|
| 26 | double [][][] u2; |
---|
| 27 | double [][][] u3; |
---|
| 28 | |
---|
| 29 | public Sochacki(int nt, int Np) { |
---|
| 30 | this.nt = nt; |
---|
| 31 | this.Np = Np; |
---|
| 32 | mass = new double[2]; |
---|
| 33 | xx = new double[3][Np][nt]; |
---|
| 34 | vv = new double[3][Np][nt]; |
---|
| 35 | |
---|
| 36 | aa = new double[Np][Np][nt]; |
---|
| 37 | u1 = new double[Np][Np][nt]; |
---|
| 38 | u2 = new double[Np][Np][nt]; |
---|
| 39 | u3 = new double[Np][Np][nt]; |
---|
| 40 | } |
---|
| 41 | |
---|
| 42 | void initState() { |
---|
| 43 | // Sun |
---|
| 44 | mass[0] = G; |
---|
| 45 | xx[0][0][0] = 0; |
---|
| 46 | xx[1][0][0] = 0; |
---|
| 47 | xx[2][0][0] = 0; |
---|
| 48 | vv[0][0][0] = 0; |
---|
| 49 | vv[1][0][0] = 0; |
---|
| 50 | vv[2][0][0] = 0; |
---|
| 51 | |
---|
| 52 | // Earth |
---|
| 53 | mass[1] = G / 332840;//S |
---|
| 54 | xx[0][1][0] = 0.98328989;//AU |
---|
| 55 | xx[1][1][0] = 0; |
---|
| 56 | xx[2][1][0] = 0; |
---|
| 57 | vv[0][1][0] = 0; |
---|
| 58 | vv[1][1][0] = 30286.7 * D / AU; |
---|
| 59 | vv[2][1][0] = 0; |
---|
| 60 | |
---|
| 61 | for( int j = 0; j < Np; j++ ) { |
---|
| 62 | for( int k = 0; k < Np; k++ ) { |
---|
| 63 | if( j == k ) |
---|
| 64 | continue; |
---|
| 65 | |
---|
| 66 | double s = 0; |
---|
| 67 | for( int i = 0; i < 3; i++ ) { |
---|
| 68 | double dx = (xx[i][j][0] - xx[i][k][0]); |
---|
| 69 | s += dx * dx; |
---|
| 70 | } |
---|
| 71 | u1[j][k][0] = 1 / Math.sqrt(s); |
---|
| 72 | u2[j][k][0] = u1[j][k][0] * u1[j][k][0]; |
---|
| 73 | u3[j][k][0] = u2[j][k][0] * u1[j][k][0]; |
---|
| 74 | |
---|
| 75 | s = 0; |
---|
| 76 | for( int i = 0; i < 3; i++ ) { |
---|
| 77 | double dx = (xx[i][k][0] - xx[i][j][0]); |
---|
| 78 | double dv = (vv[i][k][0] - vv[i][j][0]); |
---|
| 79 | s += dx * dv; |
---|
| 80 | } |
---|
| 81 | aa[j][k][0] = s; |
---|
| 82 | } |
---|
| 83 | } |
---|
| 84 | } |
---|
| 85 | |
---|
| 86 | void genPoly() { |
---|
| 87 | double a, b; |
---|
| 88 | int i, j, k, L, m; |
---|
| 89 | |
---|
| 90 | for( m = 1; m < nt; m++ ) { |
---|
| 91 | for( j = 0; j < Np; j++ ) { |
---|
| 92 | for( i = 0; i < 3; i++ ) { |
---|
| 93 | xx[i][j][m] = vv[i][j][m-1] / m; |
---|
| 94 | a = 0; |
---|
| 95 | for( k = 0; k < Np; k++ ) { |
---|
| 96 | b = 0; |
---|
| 97 | for( L = 0; L <= m-1; L++ ) { |
---|
| 98 | b = b + (xx[i][k][L] - xx[i][j][L])*u3[j][k][m-1-L]; |
---|
| 99 | } |
---|
| 100 | a = a + mass[k] * b; |
---|
| 101 | } |
---|
| 102 | vv[i][j][m] = a / m; |
---|
| 103 | } |
---|
| 104 | for( k = 0; k <= j-1; k++ ) { |
---|
| 105 | a = 0; |
---|
| 106 | for( L = 0; L <= m-1; L++ ) { |
---|
| 107 | a = a - u3[j][k][L] * aa[j][k][m-1-L]; |
---|
| 108 | } |
---|
| 109 | u1[j][k][m] = a / m; |
---|
| 110 | u1[k][j][m] = a / m; |
---|
| 111 | a = 0; |
---|
| 112 | for( L = 0 ; L <= m; L++ ) { |
---|
| 113 | a = a + u1[j][k][L] * u1[j][k][m-L]; |
---|
| 114 | } |
---|
| 115 | u2[j][k][m] = a; |
---|
| 116 | u2[k][j][m] = a; |
---|
| 117 | a = 0; |
---|
| 118 | b = 0; |
---|
| 119 | for( L = 0; L <= m; L++ ) { |
---|
| 120 | b = b + u2[j][k][L] * u1[j][k][m-L]; |
---|
| 121 | for( i = 0; i < 3; i++ ) { |
---|
| 122 | a = a + (xx[i][j][L]-xx[i][k][L])*(vv[i][j][m-L]-vv[i][k][m-L]); |
---|
| 123 | } |
---|
| 124 | } |
---|
| 125 | aa[j][k][m] = a; |
---|
| 126 | aa[k][j][m] = a; |
---|
| 127 | u3[j][k][m] = b; |
---|
| 128 | u3[k][j][m] = b; |
---|
| 129 | } |
---|
| 130 | aa[j][j][m] = 0; |
---|
| 131 | u1[j][j][m] = 0; |
---|
| 132 | u2[j][j][m] = 0; |
---|
| 133 | u3[j][j][m] = 0; |
---|
| 134 | } |
---|
| 135 | } |
---|
| 136 | } |
---|
| 137 | |
---|
| 138 | double evalPoly( double t, double [] c ) { |
---|
| 139 | double p = 0; |
---|
| 140 | |
---|
| 141 | for( int i = c.length-1; i >= 0; i-- ) { |
---|
| 142 | p = p * t + c[i]; |
---|
| 143 | } |
---|
| 144 | |
---|
| 145 | return p * AU; |
---|
| 146 | } |
---|
| 147 | |
---|
| 148 | State getState( double time ) { |
---|
| 149 | List<Particle> particles = new LinkedList<Particle>(); |
---|
| 150 | for( int i = 0; i < Np; i++) { |
---|
| 151 | double x = evalPoly(time, xx[0][i]); |
---|
| 152 | double y = evalPoly(time, xx[1][i]); |
---|
| 153 | double z = evalPoly(time, xx[2][i]); |
---|
| 154 | Vector r = new Vector(x, y, z); |
---|
| 155 | |
---|
| 156 | double vx = evalPoly(time, vv[0][i]); |
---|
| 157 | double vy = evalPoly(time, vv[1][i]); |
---|
| 158 | double vz = evalPoly(time, vv[2][i]); |
---|
| 159 | Vector v = new Vector(vx, vy, vz); |
---|
| 160 | |
---|
| 161 | Particle p = new Particle(mass[i], r, v); |
---|
| 162 | particles.add(p); |
---|
| 163 | } |
---|
| 164 | |
---|
| 165 | return new State(particles); |
---|
| 166 | } |
---|
| 167 | |
---|
| 168 | public static void main(String [] args) { |
---|
| 169 | Renderer2D renderer = new Renderer2D(); |
---|
| 170 | new Thread(renderer).start(); |
---|
| 171 | |
---|
| 172 | Sochacki sch = new Sochacki(100,2); |
---|
| 173 | sch.initState(); |
---|
| 174 | sch.genPoly(); |
---|
| 175 | for( double t = -228; t < 229; t += 1 ) { |
---|
| 176 | State s = sch.getState(t); |
---|
| 177 | //System.out.println(s); |
---|
| 178 | renderer.drawState(s); |
---|
| 179 | } |
---|
| 180 | try { |
---|
| 181 | Thread.currentThread().sleep(60 * 1000); |
---|
| 182 | } catch( InterruptedException e ) { |
---|
| 183 | e.printStackTrace(); |
---|
| 184 | } |
---|
| 185 | renderer.stop(); |
---|
| 186 | } |
---|
| 187 | } |
---|