sun.c (2267B)
1 #include "astro.h" 2 3 void 4 sun(void) 5 { 6 double mven, merth, mmars, mjup, msat; 7 double dmoon, mmoon, gmoon; 8 double pturbb, pturbl, pturbr, lograd; 9 10 ecc = .01675104 - 4.180e-5 * capt - 1.26e-7*capt2; 11 incl = 0; 12 node = 0; 13 argp = 281.220833 + .0000470684*eday + .000453*capt2 14 + .000003*capt3; 15 mrad = 1; 16 anom = 358.475845 + .9856002670*eday - .000150*capt2 17 - .000003*capt3; 18 motion = .9856473354; 19 20 dmoon = 350.737681+12.1907491914*eday-.001436*capt2; 21 gmoon = 11.250889 + 13.2293504490*eday - .003212*capt2; 22 mmoon = 296.104608 + 13.0649924465*eday + 9.192e-3*capt2; 23 mven = 212.448 + 1.602121635*eday; 24 merth = 358.476 + 0.985600267*eday; 25 mmars = 319.590 + .524024095*eday; 26 mjup = 225.269 + .083082362*eday; 27 msat = 175.593 + .033450794*eday; 28 29 dmoon = fmod(dmoon, 360.)*radian; 30 gmoon = fmod(gmoon, 360.)*radian; 31 mmoon = fmod(mmoon, 360.)*radian; 32 mven *= radian; 33 merth *= radian; 34 mmars *= radian; 35 mjup *= radian; 36 msat *= radian; 37 38 icosadd(sunfp, suncp); 39 anom += cosadd(4, mmars, merth, mven, mjup)/3600.; 40 anom += sinadd(5, mmars, merth, mven, mjup, .07884*capt)/3600.; 41 42 incl *= radian; 43 node *= radian; 44 argp *= radian; 45 anom = fmod(anom, 360.)*radian; 46 47 /* 48 * computation of elliptic orbit 49 */ 50 51 lambda = anom + argp; 52 53 pturbl = (6910.057 - 17.240*capt - 0.052*capt2)*sin(anom) 54 + (72.338 - 0.361*capt) * sin(2.*anom) 55 + (1.054 - 0.001*capt) * sin(3.*anom) 56 + 0.018 * sin(4.*anom); 57 58 lambda += pturbl*radsec; 59 60 beta = 0.; 61 62 lograd = (30.57e-6 - 0.15e-6*capt) 63 - (7274.12e-6 - 18.14e-6*capt - 0.05e-6*capt2)*cos(anom) 64 - (91.38e-6 - 0.46e-6*capt) * cos(2.*anom) 65 - (1.45e-6 - 0.01e-6*capt) * cos(3.*anom) 66 - 0.02e-6 * cos(4.*anom); 67 68 pturbl = cosadd(5, mmars, merth, mven, mjup, msat); 69 pturbl += sinadd(3, dmoon, mmoon, merth) + .9; 70 pturbl *= radsec; 71 72 pturbb = cosadd(3, merth, mven, mjup); 73 pturbb += sinadd(3, gmoon, mmoon, dmoon); 74 pturbb *= radsec; 75 76 pturbr = cosadd(5, mmars, merth, mven, mjup, msat); 77 pturbr += cosadd(3, dmoon, mmoon, merth); 78 79 lambda += pturbl; 80 if(lambda > pipi) 81 lambda -= pipi; 82 83 beta += pturbb; 84 85 lograd = (lograd+pturbr) * 2.30258509; 86 rad = 1 + lograd * (1 + lograd * (.5 + lograd/6)); 87 88 motion *= radian*mrad*mrad/(rad*rad); 89 90 semi = 961.182; 91 if(flags['o']) 92 semi = 959.63; 93 mag = -26.5; 94 }