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