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