tnept.c - plan9port - [fork] Plan 9 from user space
HTML git clone git://src.adamsgaard.dk/plan9port
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tnept.c (3333B)
---
1 #include "astro.h"
2
3 static double elem[] =
4 {
5 36525.0, /* [0] eday of epoc */
6
7 30.06896348, /* [1] semi major axis (au) */
8 0.00858587, /* [2] eccentricity */
9 1.76917, /* [3] inclination (deg) */
10 131.72169, /* [4] longitude of the ascending node (deg) */
11 44.97135, /* [5] longitude of perihelion (deg) */
12 304.88003, /* [6] mean longitude (deg) */
13
14 -0.00125196, /* [1+6] (au/julian century) */
15 0.0000251, /* [2+6] (e/julian century) */
16 -3.64, /* [3+6] (arcsec/julian century) */
17 -151.25, /* [4+6] (arcsec/julian century) */
18 -844.43, /* [5+6] (arcsec/julian century) */
19 786449.21, /* [6+6] (arcsec/julian century) */
20 };
21
22 void
23 nept(void)
24 {
25 double pturbl, pturbb, pturbr;
26 double lograd;
27 double dele, enom, vnom, nd, sl;
28
29 double capj, capn, eye, comg, omg;
30 double sb, su, cu, u, b, up;
31 double sd, ca, sa;
32
33 double cy;
34
35 cy = (eday - elem[0]) / 36525.; /* per julian century */
36
37 mrad = elem[1] + elem[1+6]*cy;
38 ecc = elem[2] + elem[2+6]*cy;
39
40 cy = cy / 3600; /* arcsec/deg per julian century */
41 incl = elem[3] + elem[3+6]*cy;
42 node = elem[4] + elem[4+6]*cy;
43 argp = elem[5] + elem[5+6]*cy;
44
45 anom = elem[6] + elem[6+6]*cy - argp;
46 motion = elem[6+6] / 36525. / 3600;
47
48 incl *= radian;
49 node *= radian;
50 argp *= radian;
51 anom = fmod(anom,360.)*radian;
52
53 enom = anom + ecc*sin(anom);
54 do {
55 dele = (anom - enom + ecc * sin(enom)) /
56 (1. - ecc*cos(enom));
57 enom += dele;
58 } while(fabs(dele) > converge);
59 vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
60 cos(enom/2.));
61 rad = mrad*(1. - ecc*cos(enom));
62
63 lambda = vnom + argp;
64 pturbl = 0.;
65 lambda += pturbl*radsec;
66 pturbb = 0.;
67 pturbr = 0.;
68
69 /*
70 * reduce to the ecliptic
71 */
72
73 nd = lambda - node;
74 lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
75
76 sl = sin(incl)*sin(nd) + pturbb*radsec;
77 beta = atan2(sl, pyth(sl));
78
79 lograd = pturbr*2.30258509;
80 rad *= 1. + lograd;
81
82
83 lambda -= 1185.*radsec;
84 beta -= 51.*radsec;
85
86 motion *= radian*mrad*mrad/(rad*rad);
87 semi = 83.33;
88
89 /*
90 * here begins the computation of magnitude
91 * first find the geocentric equatorial coordinates of Saturn
92 */
93
94 sd = rad*(cos(beta)*sin(lambda)*sin(obliq) +
95 sin(beta)*cos(obliq));
96 sa = rad*(cos(beta)*sin(lambda)*cos(obliq) -
97 sin(beta)*sin(obliq));
98 ca = rad*cos(beta)*cos(lambda);
99 sd += zms;
100 sa += yms;
101 ca += xms;
102 alpha = atan2(sa,ca);
103 delta = atan2(sd,sqrt(sa*sa+ca*ca));
104
105 /*
106 * here are the necessary elements of Saturn's rings
107 * cf. Exp. Supp. p. 363ff.
108 */
109
110 capj = 6.9056 - 0.4322*capt;
111 capn = 126.3615 + 3.9894*capt + 0.2403*capt2;
112 eye = 28.0743 - 0.0128*capt;
113 comg = 168.1179 + 1.3936*capt;
114 omg = 42.9236 - 2.7390*capt - 0.2344*capt2;
115
116 capj *= radian;
117 capn *= radian;
118 eye *= radian;
119 comg *= radian;
120 omg *= radian;
121
122 /*
123 * now find saturnicentric ring-plane coords of the earth
124 */
125
126 sb = sin(capj)*cos(delta)*sin(alpha-capn) -
127 cos(capj)*sin(delta);
128 su = cos(capj)*cos(delta)*sin(alpha-capn) +
129 sin(capj)*sin(delta);
130 cu = cos(delta)*cos(alpha-capn);
131 u = atan2(su,cu);
132 b = atan2(sb,sqrt(su*su+cu*cu));
133
134 /*
135 * and then the saturnicentric ring-plane coords of the sun
136 */
137
138 su = sin(eye)*sin(beta) +
139 cos(eye)*cos(beta)*sin(lambda-comg);
140 cu = cos(beta)*cos(lambda-comg);
141 up = atan2(su,cu);
142
143 /*
144 * at last, the magnitude
145 */
146
147
148 sb = sin(b);
149 mag = -8.68 +2.52*fabs(up+omg-u)-
150 2.60*fabs(sb) + 1.25*(sb*sb);
151
152 helio();
153 geo();
154 }