| 2 |
- |
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 |
}
|