// gcc pl.c p.c -Wall -std=c99 -lm -O3 // Derived from http://cap-lore.com/MathPhys/Orbits/c.c #include "h.h" #define C const static C R pi = 3.14159265358979323846L; typedef struct{R x, y, t;} sl; typedef struct{R x[3], xd[3];} co; typedef struct{R m[3];} v; static C R cpi = 3.14159265358979323846L/180; static R fm(R a, R m){return a - m*floorl(a/m);} static v cp(v x, v y){return (v){{ x.m[1]*y.m[2] - x.m[2]*y.m[1], x.m[2]*y.m[0] - x.m[0]*y.m[2], x.m[0]*y.m[1] - x.m[1]*y.m[0]}};} static co Kepler(C R e, C R L){sl w = {1, 0, 0}; static R Cos[64], Sin[64]; {static char first = 1; if(first){R b = 4; first = 0; for(int j=0; j<64; ++j) {Cos[j] = cosl(b); Sin[j] = sinl(b); b /=2;}}} C R a = 1/(1 - e*e), ys = sqrtl(a), M = fm(L, 2*pi); {R dt = 4; for(int j=0; j<64; ++j){C sl q = {w.x*Cos[j] - w.y*Sin[j], w.x*Sin[j] + w.y*Cos[j], w.t + dt}; dt *= .5; if(q.t - e*q.y < M) w = q;}} assert(fabsl((w.t - e*sinl(w.t)) - M) < 1.e-18); {C R sp = 1/(a*ys*(1 - e*w.x)); return (co){{(w.x - e)*a, w.y*ys, 0}, {-w.y*a*sp, w.x*ys*sp, 0}};}} void pco(char * m, co x){printf("co %s:\n", m); for(int j=0; j<3; ++j) printf("%26.19Lf %26.19Lf\n", x.x[j], x.xd[j]);} //typedef struct{R a, e, i, Omega, omega, L;} elm; typedef R elm[6]; static v bx(R C * C x){v a; for(int j=0; j<3; ++j) a.m[j] = x[j]; return a;} static R ip(v x, v y){R s = 0; for(int j=0; j<3; ++j) s += x.m[j]*y.m[j]; return s;} static void el(C co p, elm E){C v x = bx(p.x), xd = bx(p.xd), L = cp(x, xd); C R A = ip(L, L), ri = sqrtl(1/ip(x, x)), T = 0.5*ip(xd, xd) - ri, dis = 1+2*T*A, vdvd = ip(x, xd), e = A/ph - 1; C R ph = (-1 + sqrtl(dis))/(2*T), a = -1/(2*T); if(dis < 1.e-8) {E[1] = sqrtl(s(vdvd)+s(A*ri-1)); E[5] = fm(atan2l(vdvd, A*ri-1), 2*pi)/cpi;} else {E[1] = e; E[5] = acosl((ri/A-1)/e)/cpi;} printf("dis = %Le, T = %Le\n", dis, T); printf("v×v̇ = %24.19Lf ri = %24.19Lf,\n v∙v̇ = %24.19Lf atan = %24.19Lf, %Le\n", A, ri, vdvd, atan2l(vdvd, A*ri-1)/cpi, sqrtl(s(vdvd)+s(A*ri-1))); E[0] = a; E[2] = asinl(sqrtl((s(L.m[0])+s(L.m[1]))/A))/cpi; E[3] = atan2l(L.m[0], -L.m[1])/cpi; E[4] = atan2l(L.m[1], L.m[0])/cpi - E[3];} co trco(R m[3][3], co a){ v mm(R y[3]){v z; for(int k=0; k<3; ++k){ R s=0; for(int j=0; j<3; ++j) s += y[j]*m[k][j]; z.m[k] = s;} return z;} {v z = mm(a.x), zd = mm(a.xd); co c; for(int j=0; j<3; ++j) {c.x[j] = z.m[j]; c.xd[j] = zd.m[j];} return c;}} static co initp(elm p){R e = p[1]; co n = Kepler(e, (p[5] - p[4])*cpi); C R in = p[2]*cpi, Om = p[3]*cpi, om = p[4]*cpi, da = om-Om; R uo[3][3] = {{cosl(da), -sinl(da), 0},{sinl(da), cosl(da), 0},{0,0,1}}; R ui[3][3] = {{cosl(in), 0, -sinl(in)}, {0,1,0}, {sinl(in), 0, cosl(in)}}; R uO[3][3] = {{cosl(Om), -sinl(Om), 0},{sinl(Om), cosl(Om), 0},{0,0,1}}; {R A = p[0]*(1-e*e), P = 1/sqrtl(A); for(int j=0; j<3; ++j) {n.x[j] *= A; n.xd[j] *= P;}} return trco(uO, trco(ui, trco(uo, n)));} #define np 9 // For 2000 Jan 1, JD = 2451544.500000 (noon, midnight?, GMT!) static elm planets[np] = { { 0.38709893, 0.20563069, 7.00487, 48.33167, 77.45645, 252.25084}, //{1, 0, 0, 0, 0, 0}, // => {1 0, 0 1, 0 0} //{4.L/3, 0.5, 0, 0, 0, 0}, // => {2/3 0, 0 3/2, 0 0} //{ 2, 0, 0, 0, 0, 0}, // => {2 0, 0 (√2)/2, 0 0} { 0.72333199, 0.00677323, 3.39471, 76.68069, 131.53298, 181.97973}, { 1.00000011, 0.01671022, 0.00005, -11.26064, 102.94719, 100.46435}, { 1.52366231, 0.09341233, 1.85061, 49.57854, 336.04084, 355.45332}, { 5.20336301, 0.04839266, 1.30530, 100.55615, 14.75385, 34.40438}, { 9.53707032, 0.05415060, 2.48446, 113.71504, 92.43194, 49.94432}, {19.19126393, 0.04716771, 0.76986, 74.22988, 170.96424, 313.23218}, {30.06896348, 0.00858587, 1.76917, 131.72169, 44.97135, 304.88003}, {39.48168677, 0.24880766, 17.14175, 110.30347, 224.06676, 238.92881}}; static elm planetsD[np] = { { 0.00000066, 0.00002527, -23.51, -446.30, 573.57, 538101628.29}, { 0.00000092, -0.00004938, -2.86, -996.89, -108.80, 210664136.06}, {-0.00000005, -0.00003804, -46.94, -18228.25, 1198.28, 129597740.63}, {-0.00007221, 0.00011902, -25.47, -1020.19, 1560.78, 68905103.78}, { 0.00060737, -0.00012880, -4.15, 1217.17, 839.93, 10925078.35}, {-0.00301530, -0.00036762, 6.11, -1591.05,-1948.89, 4401052.95}, { 0.00152025, -0.00019150, -2.09, -1681.40, 1312.56, 1542547.79}, {-0.00125196, 0.0000251, -3.64, -151.25, -844.43, 786449.21}, {-0.00076912, 0.00006465, 11.07, -37.33, -132.25, 522747.90}}; int main(){ { // Demonstrate 1/r^2 acceleration towards origin. void ry(C R e){ C R dL = 1e-6, a = 1/(1 - e*e), xf2 = -2*e*a, ys = sqrtl(a), ts = 1/(a*ys); void rx(C R t){ C co a = Kepler(e, t-dL*ts), b = Kepler(e, t), c = Kepler(e, t+dL*ts); printf("t=%5.2Lf, ", t); for(int j=0; j<2; ++j) assert(fabs((c.x[j]-a.x[j])/dL-2*b.xd[j]) < 1e-6); assert(fabs(b.x[0]*b.xd[1] - b.x[1]*b.xd[0] - 1) < 1e-12); assert(fabs(sqrtl(s(b.x[0]) + s(b.x[1])) + sqrtl(s(b.x[0] - xf2) + s(b.x[1])) - 2/(1-e*e)) < 1.e-13); // On ellipse of orbit for(int j=0; j<2; ++j) printf("x=%12.9Lf xd=%12.9Lf ", b.x[j], b.xd[j]); printf("\n"); {C R rm3 = ({C R r2 = s(b.x[0]) + s(b.x[1]); 1/(sqrtl(r2)*r2);}); for(int j=0; j<2; ++j) printf("a=%12.9Lf p=%12.9Lf ", (a.x[j]+c.x[j] - 2*b.x[j])/(dL*dL), rm3*b.x[j]); printf("ac\n");}} printf("xf2 = %Le, e = %Le\n", xf2, e); rx(0); rx(1); rx(2.31); rx(3.14); rx(4.7);} ry(0.42); ry(0); ry(0.5); ry(0.9); if(0) ry(1.1);} { // Demonstrate that Kepler produces correct orbit for special points. C R e = 0.5; // a = 1/(1-e*e), ts = a*sqrtl(a), af = -2*e*a; R g(C R x){return acosl(x) - sqrtl(1-x*x)*e;} void tt(R t){C co c = Kepler(e, t); printf("%22.19Lf %22.19Lf %Le\n%22.19Lf %22.19Lf .\n", c.x[0], c.x[1], t, c.xd[0], c.xd[1]);} tt(g(1)); tt(g(e)); tt(g(0)); tt(g(-e)); tt(g(-1)); tt(2*pi-g(-e)); /* 0.6666666666666666667 0.0000000000000000000 0.000000e+00 -0.0000000000000000000 1.5000000000000000000 . 0.0000000000000000004 0.9999999999999999998 6.141848e-01 -0.9999999999999999999 0.5000000000000000004 . -0.6666666666666666666 1.1547005383792515292 1.070796e+00 -0.8660254037844386471 0.0000000000000000001 . -1.3333333333333333330 1.0000000000000000001 1.661382e+00 -0.6000000000000000002 -0.2999999999999999999 . -2.0000000000000000002 0.0000000000000000002 3.141593e+00 -0.0000000000000000001 -0.5000000000000000000 . -1.3333333333333333334 -1.0000000000000000000 4.621803e+00 0.6000000000000000001 -0.3000000000000000000 . */ } if(0) {C R d = 0.0001; // Drive el for(int i=0; i<4; ++i) for(int j=0; j<4; ++j) { co c = {{1, 0, 0}, {0, 1, 0}}; if(i) c.x[i-1] += d; if(j) c.xd[j-1] += d; {elm e; el(c, e); printf("%d %d\n", i, j); for(int j=0; j<6; ++j) printf("%23.19Lf\n", e[j]);}}} // Drive el {void tin(elm e, int ex){ // Illustrate initp and el: co c = initp(e); elm f; el(c, f); for(int j=0; j<6; ++j) {int t = fabsl(e[j] - f[j]) < 1.e-12 || (1&(ex>>(5-j))); if (0) assert(t); else if(1||!t) {pco("tz", c); printf("er %d, ex %d\n", j, ex); for(int j=0; j<6; ++j) printf("%Le %26.18Le,\n", e[j], f[j]); break;}}} for(int j = 0; j<360; j += 10) tin((elm){1, 1.e-8, 0, 0, 0, j}, 6); tin((elm){2, 1.e-8, 0, 0, 0, 0}, 6); tin((elm){2, 4.e-9, 0, 0, 0, 0}, 6); tin((elm){2, 0, 0, 0, 0, 0}, 6); tin((elm){1, 0.5, 0, 0, 0, 0}, 4); tin((elm){3, 0.5, 0, 0, 0, 0}, 4); tin((elm){1, 0.7, 0, 0, 0, 75}, 6); tin((elm){1, 0, 30, 0, 0, 0}, 2); tin((elm){1, 0, 0, 45, 0, 0}, 6); tin((elm){1, 0, 30, 45, 0, 0}, 2); tin((elm){1, 0, 0, 0, 0, 60}, 6); tin((elm){4, 0.6, 15, 60, 20, 140}, 0);} if(0) for(R t=0; t<5; t+=0.02) {C co x = Kepler(0.4, t); // Print an orbit. printf("%3.2Lf %8.5Lf %8.5Lf %8.5Lf %8.5Lf ", t, x.x[0], x.x[1], x.xd[0], x.xd[1]); printf("%13.10Lf\n", 1./sqrtl(s(x.x[0]) + s(x.x[1])) - 0.5*(s(x.xd[0]) + s(x.xd[1])));} void loc (C R years, pt z[N]) {co state[np]; elm orbs[np]; // computes Cartesian coordinates for time. // in centuries since January 2000; Stores into z. // Next block adjusts for time. for(int j = 0; j