/******************************************************************************* ** Compute the position of the major planets and return the epli value. ** *******************************************************************************/ #include #include #define RAD 0.01745329251994329651 /******************************************************************************* ** Crude magnitudes of planets (x100) for output filtering by brightness. ** ***************************************************************************** */ #define MAGSOL -9.9 #define MAGMER -1.5 #define MAGVEN -4.0 #define MAGMAR -1.0 #define MAGJUP -2.0 #define MAGSAT 1.0 #define MAGURA 5.9 #define MAGNEP 8.0 double planet_pos(jd,T0) double jd,T0; { double M0,M1,M2,M4,M5,M6,M7,M8,D,G,H,M,N,P,Q,S,V,W,RA,DEC,ECC; double sinQ,cosQ,sin2Q,cos2Q,sin3Q,cos3Q,sin4Q,cos4Q; double sinze,cosze,sin2ze,cos2ze,sin3ze,cos3ze,sin4ze,cos4ze; double sinps,cosps,sin2ps,cos2ps,sin3ps,cos3ps; double sinH,cosH,sin2H,cos2H,sinG,cosG; double sinV,cosV,sin2V,cos2V,sinS,cosS,sin2S,cos2S; double w1,w2,a,b,e,i,l,ll,r,u,nu,nu2,psi,eta,th,ze,epli,thapp,omeg; double l1pert,w1pert; double esun,Cen,Stheta,Sr; double kepler(),truean(),longi(),lati(),poly(),range(); void helio_trans(),speak(); /******************************************************************************* ** Calculate orbital elements for Mercury. ** *******************************************************************************/ l = range(poly(178.179078,149474.07078,0.0003011,0.,T0)); a = 0.3870986; e = poly(0.20561421,0.00002046,-0.000000030,0.,T0); i = poly(7.00288100,0.00186080,-0.000018300,0.,T0); w1 = range(poly(28.753753,0.3702806,0.0001208,0.,T0)); w2 = range(poly(47.145944,1.1852083,0.0001739,0.,T0)); M1 = range(poly(102.27938,149472.51529,0.000007,0.,T0)); /******************************************************************************* ** Solving Kepler find the eccentric anomaly ECC. ** *******************************************************************************/ ECC = kepler(e,M1); nu = truean(e,ECC); r = a * (1. - e * cos(ECC*RAD)); u = l + nu - M1 - w2; ll = longi(w2,i,u); b = lati(u,i); /******************************************************************************* ** Now make corrections due to perturbations. ** *******************************************************************************/ M2 = range(poly(212.60322,58517.80387, 0.001286,0.,T0)); M4 = range(poly(319.51913,19139.85475, 0.000181,0.,T0)); M5 = range(poly(225.32833, 3034.69202,-0.000722,0.,T0)); M6 = range(poly(175.46622, 1221.55147,-0.000502,0.,T0)); ll += 0.00204 * cos((-2. * M1 + 5. * M2 + 12.220) * RAD) + 0.00103 * cos(( -M1 + 2. * M2 - 160.692) * RAD) + 0.00091 * cos(( -M1 + 2. * M5 - 37.003) * RAD) + 0.00078 * cos((-3. * M1 + 5. * M2 + 10.137) * RAD); r += 0.000007525 * cos(( -M1 + 2. * M5 + 53.013) * RAD) + 0.000006802 * cos((-3. * M1 + 5. * M2 - 259.918) * RAD) + 0.000005457 * cos((-2. * M1 + 2. * M2 - 71.188) * RAD) + 0.000003569 * cos(( -M1 + 5. * M2 - 77.750) * RAD); /******************************************************************************* ** Now find the Longi and radius vector of the sun. ** *******************************************************************************/ M = range(poly(358.47583,35999.0497500,-0.0001500,-0.0000033,T0)); esun = poly(0.01675104,-0.0000418,-0.000000126,0.,T0); Cen = poly(1.919460,-0.004789,-0.000014,0.,T0) * sin( M * RAD) + poly(0.020094,-0.000100, 0. ,0.,T0) * sin(2. * M * RAD) + poly(0.000293, 0. , 0. ,0.,T0) * sin(3. * M * RAD); Stheta = range(Cen + range(poly(279.69668,36000.76892,0.0003025,0.,T0))); Sr = 1.0000002 * (1. - esun * esun) / (1. + esun * cos((M + Cen) * RAD)); omeg = 259.18 - 1934.142 * T0; thapp = Stheta - 0.00569 - 0.00479 * sin(omeg * RAD); epli = poly(23.45229400,-0.013012500,-0.000001640,0.000000503,T0); N = cos(epli * RAD) * sin(thapp * RAD); D = cos(thapp * RAD); RA = atan2(N,D) / RAD; DEC = asin(sin(epli * RAD) * sin(thapp * RAD)) / RAD; #ifdef EUTSCH speak(RA,DEC,Sr,MAGSOL,"PS","Sonne"); #else speak(RA,DEC,Sr,MAGSOL,"PS","Sol"); #endif /******************************************************************************* ** Tansformation of coordinates on Mercury and output. ** *******************************************************************************/ #ifdef EUTSCH helio_trans(r,b,ll,Stheta,Sr,epli,MAGMER,"PM","Merkur"); #else helio_trans(r,b,ll,Stheta,Sr,epli,MAGMER,"PM","Mercury"); #endif /******************************************************************************* ** Now start on Venus. ** *******************************************************************************/ l = range(poly(342.767053,58519.21191,0.0003097,0.,T0)); a = 0.7233316; e = poly(0.00682069,-0.00004774, 0.000000091,0.,T0); i = poly(3.39363100, 0.00100580,-0.000001000,0.,T0); w1 = range(poly(54.384186,0.5081861,-0.0013864,0.,T0)); w2 = range(poly(75.779647,0.8998500, 0.0004100,0.,T0)); /******************************************************************************* ** Venus has a long period pert that needs to be added before Kelper. ** *******************************************************************************/ M0 = M2 + 0.00077 * sin((237.24 + 150.27 * T0) * RAD); l += M0 - M2; ECC = kepler(e,M0); nu = truean(e,ECC); r = a * (1. - e * cos(ECC * RAD)); u = l + nu - M0 - w2; ll = longi(w2,i,u); b = lati(u,i); /******************************************************************************* ** Now Venus perturbations. ** *******************************************************************************/ ll += 0.00313 * cos((2. * M - 2. * M2 - 148.225) * RAD) + 0.00198 * cos((3. * M - 3. * M2 + 2.565) * RAD) + 0.00136 * cos(( M - M2 - 119.107) * RAD) + 0.00096 * cos((3. * M - 2. * M2 - 135.912) * RAD) + 0.00082 * cos(( M5 - M2 - 208.087) * RAD); r += 0.000022501 * cos((2. * M - 2. * M2 - 58.208) * RAD) + 0.000019045 * cos((3. * M - 3. * M2 + 92.577) * RAD) + 0.000006887 * cos(( M5 - M2 - 118.090) * RAD) + 0.000005172 * cos(( M - M2 - 29.110) * RAD) + 0.000003620 * cos((5. * M - 4. * M2 - 104.208) * RAD) + 0.000003283 * cos((4. * M - 4. * M2 + 63.513) * RAD) + 0.000003074 * cos((2. * M5 - 2. * M2 - 55.167) * RAD); /******************************************************************************* ** Tansformation of coordinates on Venus and output. ** *******************************************************************************/ helio_trans(r,b,ll,Stheta,Sr,epli,MAGVEN,"PV","Venus"); /******************************************************************************* ** Now start the planet Mars. ** *******************************************************************************/ l = range(poly(293.737334,19141.69551,0.0003107,0.,T0)); a = 1.5236883; e = poly(0.09331290, 0.000092064,-0.000000077,0.,T0); i = poly(1.85033300,-0.000675000, 0.000012600,0.,T0); w1 = range(poly(285.431761,1.0697667, 0.0001313, 0.00000414,T0)); w2 = range(poly( 48.786442,0.7709917,-0.0000014,-0.00000533,T0)); /******************************************************************************* ** Mars has a long period perturbation. ** *******************************************************************************/ M0 = M4 + -0.01133 * sin((3. * M5 - 8. * M4 + 4. * M) * RAD) -0.00933 * cos((3. * M5 - 8. * M4 + 4. * M) * RAD); l += M0 - M4; ECC = kepler(e,M0); nu = truean(e,ECC); r = a * (1. - e * cos(ECC * RAD)); u = l + nu - M0 - w2; ll = longi(w2,i,u); b = lati(u,i); /******************************************************************************* ** Now Mars Perturbations. ** *******************************************************************************/ ll += 0.00705 * cos(( M5 - M4 - 48.958) * RAD) + 0.00607 * cos((2. * M5 - M4 - 188.350) * RAD) + 0.00445 * cos((2. * M5 - 2. * M4 - 191.897) * RAD) + 0.00388 * cos(( M - 2. * M4 + 20.495) * RAD) + 0.00238 * cos(( M - M4 + 35.097) * RAD) + 0.00204 * cos((2. * M - 3. * M4 + 158.638) * RAD) + 0.00177 * cos(( -M2 + 3. * M4 - 57.602) * RAD) + 0.00136 * cos((2. * M - 4. * M4 + 154.093) * RAD) + 0.00104 * cos(( M5 + 17.618) * RAD); r += 0.000053227 * cos(( M5 - M4 + 41.1306) * RAD) + 0.000050989 * cos((2. * M5 - 2. * M4 - 101.9847) * RAD) + 0.000038278 * cos((2. * M5 - M4 - 98.3292) * RAD) + 0.000015996 * cos(( M - M4 - 55.5550) * RAD) + 0.000014764 * cos((2. * M - 3. * M4 + 68.6220) * RAD) + 0.000008966 * cos(( M5 - 2. * M4 + 43.6150) * RAD) + 0.000007914 * cos((3. * M5 - 2. * M4 - 139.7370) * RAD) + 0.000007004 * cos((2. * M5 - 3. * M4 - 102.8880) * RAD) + 0.000006620 * cos(( M - 2. * M4 + 113.2020) * RAD) + 0.000004930 * cos((3. * M5 - 3. * M4 - 76.2430) * RAD) + 0.000004693 * cos((3. * M - 5. * M4 + 190.6030) * RAD) + 0.000004571 * cos((2. * M - 4. * M4 + 244.7020) * RAD) + 0.000004409 * cos((3. * M5 - M4 - 115.8280) * RAD); /******************************************************************************* ** Tansformation of coordinates on Mars and output. ** *******************************************************************************/ helio_trans(r,b,ll,Stheta,Sr,epli,MAGMAR,"Pm","Mars"); /******************************************************************************* ** Now start Jupiter. ** *******************************************************************************/ l = range(poly(238.049257,3036.301986,0.0003347,-0.00000165,T0)); a = 5.202561; e = poly(0.04833475, 0.000164180,-0.0000004676,-0.0000000017,T0); i = poly(1.30873600,-0.005696100, 0.0000039000, 0.0000000000,T0); w1 = range(poly(273.277558,0.5994317,0.00070405, 0.00000508,T0)); w2 = range(poly( 99.443414,1.0105300,0.00035222,-0.00000851,T0)); /******************************************************************************* ** Now start perturbation calclations. ** *******************************************************************************/ nu2 = T0 / 5. + 0.1; P = 237.47555 + 3034.9061 * T0; Q = 265.91650 + 1222.1139 * T0; S = 243.51721 + 428.4677 * T0; V = range(5. * Q - 2. * P) * RAD; W = range(2. * P - 6. * Q + 3. * S) * RAD; ze = range(Q - P) * RAD; psi = range(S - Q) * RAD; P = range(P) * RAD; Q = range(Q) * RAD; S = range(S) * RAD; /******************************************************************************* ** An extreme acceleration can be achieved by using the addtion theorems. ** *******************************************************************************/ sinQ = sin(Q); cosQ = cos(Q); sin2Q = 2. * sinQ * cosQ; cos2Q = cosQ*cosQ - sinQ*sinQ; sin3Q = ( 3. - 4. * sinQ * sinQ) * sinQ; cos3Q = (- 3. + 4. * cosQ * cosQ) * cosQ; sin4Q = (8. * cosQ * cosQ - 4.) * cosQ * sinQ; cos4Q = 8. * cosQ * cosQ * ( cosQ * cosQ - 1. ) + 1.; /*******************************************************************************/ sinze = sin(ze); cosze = cos(ze); sin2ze = 2. * sinze * cosze; cos2ze = cosze*cosze - sinze*sinze; sin3ze = ( 3. - 4. * sinze * sinze) * sinze; cos3ze = (- 3. + 4. * cosze * cosze) * cosze; sin4ze = (8. * cosze * cosze - 4.) * cosze * sinze; cos4ze = 8. * cosze * cosze * (cosze * cosze - 1. ) + 1.; /*******************************************************************************/ sinV = sin(V); cosV = cos(V); sin2V = 2. * sinV * cosV; cos2V = cosV * cosV - sinV * sinV; /*******************************************************************************/ sinS = sin(S); cosS = cos(S); sin2S = 2. * sinS * cosS; cos2S = cosS * cosS - sinS * sinS; /*******************************************************************************/ sinps = sin(psi); cosps = cos(psi); sin2ps = 2. * sinps * cosps; cos2ps = cosps*cosps - sinps*sinps; sin3ps = ( 3. - 4. * sinps * sinps) * sinps; cos3ps = (- 3. + 4. * cosps * cosps) * cosps; /******************************************************************************* ** Calculating the perturbation terms. ** *******************************************************************************/ l1pert = poly(0.331364,-0.010281,-0.004692,0.,nu2) * sinV + poly(0.003228,-0.064436, 0.002075,0.,nu2) * cosV - poly(0.003083, 0.000275,-0.000489,0.,nu2) * sin2V + 0.002472 * sin(W) + 0.013619 * sinze + 0.018472 * sin2ze + 0.006717 * sin3ze + 0.002775 * sin4ze + (0.007275 - 0.001253 * nu2) * sinze * sinQ + 0.006417 * sin2ze * sinQ + 0.002439 * sin3ze * sinQ - (0.033839 + 0.001125 * nu2) * cosze * sinQ - 0.003767 * cos2ze * sinQ - (0.035681 + 0.001208 * nu2) * sinze * cosQ - 0.004261 * sin2ze * cosQ + 0.002178 * cosQ - (0.006333 - 0.001161 * nu2) * cosze * cosQ - 0.006675 * cos2ze * cosQ - 0.002664 * cos3ze * cosQ - 0.002572 * sinze * sin2Q - 0.003567 * sin2ze * sin2Q + 0.002094 * cosze * cos2Q + 0.003342 * cos2ze * cos2Q; w1pert = (0.007192 - 0.003147 * nu2) * sinV - poly(0.020428,0.000675,-0.000197,0.,nu2) * cosV + (0.007269 + 0.000672 * nu2) * sinze * sinQ - 0.004344 * sinQ + 0.034036 * cosze * sinQ + 0.005614 * cos2ze * sinQ + 0.002964 * cos3ze * sinQ + 0.037761 * sinze * cosQ + 0.006158 * sin2ze * cosQ - 0.006603 * cosze * cosQ - 0.005356 * sinze * sin2Q + 0.002722 * sin2ze * sin2Q + 0.004483 * cosze * sin2Q - 0.002642 * cos2ze * sin2Q + 0.004403 * sinze * cos2Q - 0.002536 * sin2ze * cos2Q + 0.005547 * cosze * cos2Q - 0.002689 * cos2ze * cos2Q; l += l1pert; w1 += w1pert; M0 = M5 + l1pert - (w1pert / e); e += poly(0.0003606, 0.0000130,-0.0000043,0.,nu2) * sinV + poly(0.0001289,-0.0000580,0. ,0.,nu2) * cosV - 0.0006764 * sinze * sinQ - 0.0001110 * sin2ze * sinQ - 0.0000224 * sin3ze * sinQ - 0.0000204 * sinQ + (0.0001284 + 0.0000116 * nu2) * cosze * sinQ + 0.0000188 * cos2ze * sinQ + (0.0001460 + 0.0000130 * nu2) * sinze * cosQ + 0.0000224 * sin2ze * cosQ - 0.0000817 * cosQ + 0.0006074 * cosze * cosQ + 0.0000992 * cos2ze * cosQ + 0.0000508 * cos3ze * cosQ + 0.0000230 * cos4ze * cosQ + 0.0000108 * cos(5. * ze) * cosQ - (0.0000956 + 0.0000073 * nu2) * sinze * sin2Q + 0.0000448 * sin2ze * sin2Q + 0.0000137 * sin3ze * sin2Q - (0.0000997 - 0.0000108 * nu2) * cosze * sin2Q + 0.0000480 * cos2ze * sin2Q + 0.0000148 * cos3ze * sin2Q - (0.0000956 - 0.0000099 * nu2) * sinze * cos2Q + 0.0000490 * sin2ze * cos2Q + 0.0000158 * sin3ze * cos2Q + 0.0000179 * cos2Q + (0.0001024 + 0.0000075 * nu2) * cosze * cos2Q - 0.0000437 * cos2ze * cos2Q - 0.0000132 * cos3ze * cos2Q; a += -0.000263 * cosV + 0.000205 * cosze + 0.000693 * cos2ze + 0.000312 * cos3ze + 0.000147 * cos4ze + 0.000299 * sinze * sinQ + 0.000181 * cos2ze * sinQ + 0.000204 * sin2ze * cosQ + 0.000111 * sin3ze * cosQ - 0.000337 * cosze * cosQ - 0.000111 * cos2ze * cosQ; ECC = kepler(e,M0); nu = truean(e,ECC); r = a * (1. - e * cos(ECC * RAD)); u = l + nu - M0 - w2; ll = longi(w2,i,u); b = lati(u,i); /******************************************************************************* ** Tansformation of coordinates on Jupiter and output. ** *******************************************************************************/ helio_trans(r,b,ll,Stheta,Sr,epli,MAGJUP,"PJ","Jupiter"); /******************************************************************************* ** Now start Saturn. ** *******************************************************************************/ l = range(poly(266.564377,1223.509884,0.0003245,-0.0000058,T0)); a = 9.554747; e = poly(0.05589232,-0.00034550,-0.000000728,0.00000000074,T0); i = poly(2.49251900,-0.00391890,-0.000015490,0.00000004000,T0); w1 = range(poly(338.307800,1.0852207, 0.00097854, 0.00000992,T0)); w2 = range(poly(112.790414,0.8731951,-0.00015218,-0.00000531,T0)); /******************************************************************************* ** Now Saturn's perturbations. ** *******************************************************************************/ l1pert = poly(-0.814181,0.018150, 0.016714,0.,nu2) * sinV + poly(-0.010497,0.160906,-0.004100,0.,nu2) * cosV + 0.007581 * sin2V - 0.007986 * sin(W) - 0.148811 * sinze - 0.040786 * sin2ze - 0.015208 * sin3ze - 0.006339 * sin4ze - 0.006244 * sinQ + ( 0.008931 + 0.002728 * nu2) * sinze * sinQ - 0.016500 * sin2ze * sinQ - 0.005775 * sin3ze * sinQ + ( 0.081344 + 0.003206 * nu2) * cosze * sinQ + 0.015019 * cos2ze * sinQ + ( 0.085581 + 0.002494 * nu2) * sinze * cosQ + ( 0.025328 - 0.003117 * nu2) * cosze * cosQ + 0.014394 * cos2ze * cosQ + 0.006319 * cos3ze * cosQ + 0.006369 * sinze * sin2Q + 0.009156 * sin2ze * sin2Q + 0.007525 * sin3ps * sin2Q - 0.005236 * cosze * cos2Q - 0.007736 * cos2ze * cos2Q - 0.007528 * cos3ps * cos2Q; w1pert = poly(0.077108, 0.007186,-0.001533,0.,nu2) * sinV + poly(0.045803,-0.014766,-0.000536,0.,nu2) * cosV - 0.007075 * sinze - 0.075825 * sinze * sinQ - 0.024839 * sin2ze * sinQ - 0.008631 * sin3ze * sinQ - 0.072586 * cosQ - 0.150383 * cosze * cosQ + 0.026897 * cos2ze * cosQ + 0.010053 * cos3ze * cosQ - ( 0.013597 + 0.001719 * nu2) * sinze * sin2Q + (-0.007742 + 0.001517 * nu2) * cosze * sin2Q + ( 0.013586 - 0.001375 * nu2) * cos2ze * sin2Q + (-0.013667 + 0.001239 * nu2) * sinze * cos2Q + 0.011981 * sin2ze * cos2Q + ( 0.014861 + 0.001136 * nu2) * cosze * cos2Q - ( 0.013064 + 0.001628 * nu2) * cos2ze * cos2Q; l += l1pert; w1 += w1pert; M0 = M6 + l1pert - (w1pert / e); e += poly(-0.0007927,0.0002548, 0.0000091,0.,nu2) * sinV + poly( 0.0013381,0.0001226,-0.0000253,0.,nu2) * cosV + ( 0.0000248 - 0.0000121 * nu2) * sin2V - ( 0.0000305 + 0.0000091 * nu2) * cos2V + 0.0000412 * sin2ze + 0.0012415 * sinQ + ( 0.0000390 - 0.0000617 * nu2) * sinze * sinQ + ( 0.0000165 - 0.0000204 * nu2) * sin2ze * sinQ + 0.0026599 * cosze * sinQ - 0.0004687 * cos2ze * sinQ - 0.0001870 * cos3ze * sinQ - 0.0000821 * cos4ze * sinQ - 0.0000377 * cos(5. * ze) * sinQ + 0.0000497 * cos2ps * sinQ + ( 0.0000163 - 0.0000611 * nu2) * cosQ - 0.0012696 * sinze * cosQ - 0.0004200 * sin2ze * cosQ - 0.0001503 * sin3ze * cosQ - 0.0000619 * sin4ze * cosQ - 0.0000268 * sin(5. * ze) * cosQ - ( 0.0000282 + 0.0001306 * nu2) * cosze * cosQ + (-0.0000086 + 0.0000230 * nu2) * cos2ze * cosQ + 0.0000461 * sin2ps * cosQ - 0.0000350 * sin2Q + ( 0.0002211 - 0.0000286 * nu2) * sinze * sin2Q - 0.0002208 * sin2ze * sin2Q - 0.0000568 * sin3ze * sin2Q - 0.0000346 * sin4ze * sin2Q - ( 0.0002780 + 0.0000222 * nu2) * cosze * sin2Q + ( 0.0002022 + 0.0000263 * nu2) * cos2ze * sin2Q + 0.0000248 * cos3ze * sin2Q + 0.0000242 * sin3ps * sin2Q + 0.0000467 * cos3ps * sin2Q - 0.0000490 * cos2Q - ( 0.0002842 + 0.0000279 * nu2) * sinze * cos2Q + ( 0.0000128 + 0.0000226 * nu2) * sin2ze * cos2Q + 0.0000224 * sin3ze * cos2Q + (-0.0001594 + 0.0000282 * nu2) * cosze * cos2Q + ( 0.0002162 - 0.0000207 * nu2) * cos2ze * cos2Q + 0.0000561 * cos3ze * cos2Q + 0.0000343 * cos4ze * cos2Q + 0.0000469 * sin3ps * cos2Q - 0.0000242 * cos3ps * cos2Q - 0.0000205 * sinze * sin3Q + 0.0000262 * sin3ze * sin3Q + 0.0000208 * cosze * cos3Q - 0.0000271 * cos3ze * cos3Q - 0.0000382 * cos3ze * sin4Q - 0.0000376 * sin3ze * cos4Q; a += 0.000572 * sinV - 0.001590 * sin2ze * cosQ + 0.002933 * cosV - 0.000647 * sin3ze * cosQ + 0.033629 * cosze - 0.000344 * sin4ze * cosQ - 0.003081 * cos2ze + 0.002885 * cosze * cosQ - 0.001423 * cos3ze + (0.002172 + 0.000102 * nu2) * cos2ze * cosQ - 0.000671 * cos4ze + 0.000296 * cos3ze * cosQ - 0.000320 * cos(5. * ze) - 0.000267 * sin2ze * sin2Q + 0.001098 * sinQ - 0.000778 * cosze * sin2Q - 0.002812 * sinze * sinQ + 0.000495 * cos2ze * sin2Q + 0.000688 * sin2ze * sinQ + 0.000250 * cos3ze * sin2Q - 0.000393 * sin3ze * sinQ - 0.000228 * sin4ze * sinQ + 0.002138 * cosze * sinQ - 0.000999 * cos2ze * sinQ - 0.000642 * cos3ze * sinQ - 0.000325 * cos4ze * sinQ - 0.000890 * cosQ + 0.002206 * sinze * cosQ - 0.000856 * sinze * cos2Q + 0.000441 * sin2ze * cos2Q + 0.000296 * cos2ze * cos2Q + 0.000211 * cos3ze * cos2Q - 0.000427 * sinze * sin3Q + 0.000398 * sin3ze * sin3Q + 0.000344 * cosze * cos3Q - 0.000427 * cos3ze * cos3Q; ECC = kepler(e,M0); nu = truean(e,ECC); r = a * (1. - e * cos(ECC * RAD)); u = l + nu - M0 - w2; ll = longi(w2,i,u); b = lati(u,i) + 0.000747 * cosze * sinQ + 0.001069 * cosze * cosQ + 0.002108 * sin2ze * sin2Q + 0.001261 * cos2ze * sin2Q + 0.001236 * sin2ze * cos2Q - 0.002075 * cos2ze * cos2Q; /******************************************************************************* ** Tansformation of coordinates on Saturn and output. ** *******************************************************************************/ helio_trans(r,b,ll,Stheta,Sr,epli,MAGSAT,"Ps","Saturn"); /******************************************************************************* ** Now Start on Uranus. ** *******************************************************************************/ l = range(poly(244.197470,429.863546,0.0003160,-0.00000060,T0)); a = 19.21814; e = poly(0.0463444,-0.00002658,0.000000077,0.,T0); i = poly(0.7724640, 0.00062530,0.000039500,0.,T0); w1 = range(poly(98.071581,0.9857650,-0.0010745,-0.00000061,T0)); w2 = range(poly(73.477111,0.4986678, 0.0013117, 0.00000000,T0)); M7 = range(poly(72.64878,428.37911,0.000079,0.,T0)); /******************************************************************************* ** Now perturbation corrections for Uranus. ** *******************************************************************************/ G = (83.76922 + 218.4901 * T0) * RAD; H = range((2. * G - S) / RAD) * RAD; ze = range(( S - P) / RAD) * RAD; eta = range(( S - Q) / RAD) * RAD; th = range(( G - S) / RAD) * RAD; G = range(G / RAD) * RAD; /******************************************************************************* ** Addition theorems. ** *******************************************************************************/ sinze = sin(ze); cosze = cos(ze); sinH = sin(H); cosH = cos(H); sin2H = 2. * sinH * cosH; cos2H = cosH * cosH - sinH * sinH; sinG = sin(G); cosG = cos(G); /******************************************************************************/ l1pert = (0.864319 - 0.001583 * nu2) * sinH + (0.082222 - 0.006833 * nu2) * cosH + 0.036017 * sin2H - 0.003019 * cos2H + 0.008122 * sin(W); w1pert = 0.120303 * sinH + (0.019472 - 0.000947 * nu2) * cosH + 0.006197 * sin2H; l += l1pert; w1 += w1pert; M0 = M7 + l1pert - (w1pert / e); e += (-0.0003349 + 0.0000163 * nu2) * sinH + 0.0020981 * cosH + 0.0001311 * cosH; a -= 0.003825 * cosH; ECC = kepler(e,M0); nu = truean(e,ECC); r = a * (1. - e * cos(ECC * RAD)); u = l + nu - M0 - w2; ll = longi(w2,i,u); b = lati(u,i); ll += (0.010122 - 0.000988 * nu2) * sin( eta + S) + poly(-0.038581, 0.002031,-0.001910,0.,nu2) * cos( eta + S) + poly( 0.034964,-0.001038, 0.000868,0.,nu2) * cos( eta + 2. * S) + 0.005594 * sin(3. * th + S) - 0.014808 * sinze - 0.005794 * sin( eta) + 0.002347 * cos( eta) + 0.009872 * sin( th) + 0.008803 * sin(2. * th) - 0.004308 * sin(3. * th); b += (0.000458 * sin( eta) - 0.000642 * cos( eta) - 0.000517 * cos(4. * th)) * sinS - (0.000347 * sin( eta) + 0.000853 * cos( eta) + 0.000517 * sin(4. * eta)) * cosS + 0.000403 * (cos(2. * th) * sin2S + sin(2. * th) * cos2S); r += -0.025948 + 0.004985 * cosze - 0.001230 * cosS + 0.003354 * cos( eta) + (0.005795 * cosS - 0.001165 * sinS + 0.001388 * sin( eta) * cos2S) + (0.001351 * cosS + 0.005702 * sinS + 0.001388 * cos( eta) * sin2S) + 0.000904 * cos(2. * th) + 0.000894 * (cos( th) - cos(3. * th)); /******************************************************************************* ** Tansformation of coordinates on Uranus and output. ** *******************************************************************************/ helio_trans(r,b,ll,Stheta,Sr,epli,MAGURA,"PU","Uranus"); /******************************************************************************* ** Now start Neptune. ** *******************************************************************************/ l = range(poly(84.457994,219.885914,0.0003205,-0.00000060,T0)); a = 30.10957; e = poly(0.00899704, 0.000006330,-0.000000002,0.,T0); i = poly(1.77924200,-0.009543600,-0.000009100,0.,T0); w1 = range(poly(276.045975,0.3256394,0.00014095, 0.000004113,T0)); w2 = range(poly(130.681389,1.0989350,0.00024987,-0.000004718,T0)); M8 = poly(37.73063,218.46134,-0.000070,0.,T0); /******************************************************************************* ** Now perturbation corrections for neptune. ** *******************************************************************************/ ze = range((G - P) / RAD) * RAD; eta = range((G - Q) / RAD) * RAD; th = range((G - S) / RAD) * RAD; sinze = sin(ze); cosze = cos(ze); l1pert = (-0.589833 + 0.001089 * nu2) * sinH + (-0.056094 + 0.004658 * nu2) * cosH - 0.024286 * sin2H; w1pert = 0.024039 * sinH - 0.025303 * cosH + 0.006206 * sin2H - 0.005992 * cos2H; l += l1pert; w1 += w1pert; M0 = M8 + l1pert - (w1pert / e); e += 0.0004389 * sinH + 0.0004262 * cosH + 0.0001129 * sin2H + 0.0001089 * cos2H; a += -0.000817 * sinH + 0.008189 * cosH + 0.000781 * cos2H; ECC = kepler(e,M0); nu = truean(e,ECC); r = a * (1. - e * cos(ECC * RAD)); u = l + nu - M0 - w2; ll = longi(w2,i,u); b = lati(u,i); ll += -0.009556 * sinze - 0.005178 * sin( eta) + 0.002572 * sin(2. * th) - 0.002972 * cos(2. * th) * sinG - 0.002833 * sin(2. * th) * cosG; b += 0.000336 * cos(2. * th) * sinG + 0.000364 * sin(2. * th) * cosG; r += -0.040596 + 0.004992 * cosze + 0.002744 * cos( eta) + 0.002044 * cos( th) + 0.001051 * cos(2. * th); /******************************************************************************* ** Tansformation of coordinates on Neptune and output. ** *******************************************************************************/ #ifdef EUTSCH helio_trans(r,b,ll,Stheta,Sr,epli,MAGNEP,"PN","Neptun"); #else helio_trans(r,b,ll,Stheta,Sr,epli,MAGNEP,"PN","Neptune"); #endif return(epli); }