source: Sophya/trunk/SophyaExt/XephemAstroLib/mooncolong.c@ 1719

Last change on this file since 1719 was 1719, checked in by cmv, 24 years ago

Adapted to version 3.5 xephem cmv 22/10/2001
File size: 6.3 KB
RevLine 
[1457]1/* code to compute lunar sunrise position and local sun angle.
2 */
3
4#include <stdio.h>
5#include <math.h>
6#if defined(__STDC__)
7#include <stdlib.h>
8#endif
9
10#include "P_.h"
11#include "astro.h"
12
13static void Librations P_((double RAD, double LAMH, double BH, double OM,
14 double F, double L, double L1, double *L0, double *B0));
15static void Moon P_((double RAD, double T, double T2, double LAM0, double R,
16 double M, double *F, double *L1, double *OM, double *LAM, double *B,
17 double *DR, double *LAMH, double *BH));
18static void Sun P_((double RAD, double T, double T2, double *L, double *M,
19 double *R, double *LAM0));
20
21/* given a Julian date and a lunar location, find selenographic colongitude of
22 * rising sun, lunar latitude of subsolar point, illuminated fraction, and alt
23 * of sun at the given location. Any pointer may be 0 if not interested.
24 * From Bruning and Talcott, October 1995 _Astronomy_, page 76.
25 * N.B. lunar coordinates use +E, but selenograhic colongs are +W.
26 */
27void
28moon_colong (jd, lt, lg, cp, kp, ap, sp)
29double jd; /* jd */
30double lt, lg; /* lat/long of location on moon, rads +N +E */
31double *cp; /* selenographic colongitude (-lng of rising sun), rads */
32double *kp; /* illuminated fraction of surface from Earth */
33double *ap; /* sun altitude at location, rads */
34double *sp; /* lunar latitude of subsolar point, rads */
35{
36 double RAD = .0174533;
37 double T;
38 double T2;
39 double L, M, R, LAM0;
40 double F, L1, OM, LAM, B, DR, LAMH, BH;
41 double L0, B0;
42 double TEMP;
43 double C0;
44 double PSI;
45 double NUM, DEN;
46 double I, K;
47 double THETA, ETA;
48 double H;
49
50 T = (jd - 2451545)/36525.0;
51 T2 = T * T;
52
53 Sun(RAD, T, T2, &L, &M, &R, &LAM0);
54 Moon(RAD, T, T2, LAM0, R, M, &F, &L1, &OM, &LAM, &B, &DR, &LAMH, &BH);
55 Librations(RAD, LAMH, BH, OM, F, L, L1, &L0, &B0);
56 if (sp)
57 *sp = B0;
58
59 TEMP = L0 / 360;
60 L0 = ((TEMP) - (int)(TEMP)) * 360;
61 if (L0 < 0) L0 = L0 + 360;
62 if (L0 <= 90) C0 = 90 - L0; else C0 = 450 - L0;
63 if (cp) {
64 *cp = degrad(C0);
65 range (cp, 2*PI); /* prefer 0..360 +W */
66 }
67
68 if (kp) {
69 TEMP = cos(B * RAD) * cos(LAM - LAM0 * RAD);
70 PSI = acos(TEMP);
71 NUM = R * sin(PSI);
72 DEN = DR - R * TEMP;
73 I = atan(NUM / DEN);
74 if (NUM * DEN < 0) I = I + 3.14159;
75 if (NUM < 0) I = I + 3.14159;
76 K = (1 + cos(I)) / 2;
77 *kp = K;
78 }
79
80 if (ap) {
81 THETA = lt;
82 ETA = lg;
83 C0 = C0 * RAD;
84 TEMP = sin(B0) * sin(THETA) + cos(B0) * cos(THETA) * sin(C0+ETA);
85 H = asin(TEMP);
86 *ap = H;
87 }
88}
89
90static void
91Librations (RAD, LAMH, BH, OM, F, L, L1, L0, B0)
92double RAD;
93double LAMH;
94double BH;
95double OM;
96double F;
97double L;
98double L1;
99double *L0;
100double *B0;
101{
102 double I, PSI, W, NUM, DEN, A, TEMP;
103
104 /* inclination of lunar equator */
105 I = 1.54242 * RAD;
106
107 /* nutation in longitude, in arcseconds */
108 PSI = -17.2 * sin(OM) - 1.32 * sin(2 * L) - .23 * sin(2 * L1) +
109 .21 * sin(2 * OM);
110 PSI = PSI * RAD / 3600;
111
112 /* optical librations */
113 W = (LAMH - PSI) - OM;
114 NUM = sin(W) * cos(BH) * cos(I) - sin(BH) * sin(I);
115 DEN = cos(W) * cos(BH);
116 A = atan(NUM / DEN);
117 if (NUM * DEN < 0) A = A + 3.14159;
118 if (NUM < 0) A = A + 3.14159;
119 *L0 = (A - F) / RAD;
120 TEMP = -sin(W) * cos(BH) * sin(I) - sin(BH) * cos(I);
121 *B0 = asin(TEMP);
122}
123
124static void
125Moon (RAD, T, T2, LAM0, R, M, F, L1, OM, LAM, B, DR, LAMH, BH)
126double RAD;
127double T;
128double T2;
129double LAM0;
130double R;
131double M;
132double *F;
133double *L1;
134double *OM;
135double *LAM;
136double *B;
137double *DR;
138double *LAMH;
139double *BH;
140{
141 double T3, M1, D2, SUMR, SUML, DIST;
142
143 T3 = T * T2;
144
145 /* argument of the latitude of the Moon */
146 *F = (93.2721 + 483202 * T - .003403 * T2 - T3 / 3526000) * RAD;
147
148 /* mean longitude of the Moon */
149 *L1 = (218.316 + 481268. * T) * RAD;
150
151 /* longitude of the ascending node of Moon's mean orbit */
152 *OM = (125.045 - 1934.14 * T + .002071 * T2 + T3 / 450000) * RAD;
153
154 /* Moon's mean anomaly */
155 M1 = (134.963 + 477199 * T + .008997 * T2 + T3 / 69700) * RAD;
156
157 /* mean elongation of the Moon */
158 D2 = (297.85 + 445267 * T - .00163 * T2 + T3 / 545900) * 2 * RAD;
159
160 /* Lunar distance */
161 SUMR = -20954 * cos(M1) - 3699 * cos(D2 - M1) - 2956 * cos(D2);
162 *DR = 385000 + SUMR;
163
164 /* geocentric latitude */
165 *B = 5.128 * sin(*F) + .2806 * sin(M1 + *F) + .2777 * sin(M1 - *F) +
166 .1732 * sin(D2 - *F);
167 SUML = 6.289 * sin(M1) + 1.274 * sin(D2 - M1) + .6583 * sin(D2) +
168 .2136 * sin(2 * M1) - .1851 * sin(M) - .1143 * sin(2 * *F);
169 *LAM = *L1 + SUML * RAD;
170 DIST = *DR / R;
171 *LAMH = (LAM0 + 180 + DIST * cos(*B) * sin(LAM0 * RAD - *LAM) / RAD)
172 * RAD;
173 *BH = DIST * *B * RAD;
174}
175
176static void
177Sun (RAD, T, T2, L, M, R, LAM0)
178double RAD;
179double T;
180double T2;
181double *L;
182double *M;
183double *R;
184double *LAM0;
185{
186 double T3, C, V, E, THETA, OM;
187
188 T3 = T2 * T;
189
190 /* mean longitude of the Sun */
191 *L = 280.466 + 36000.8 * T;
192
193 /* mean anomaly of the Sun */
194 *M = 357.529 + 35999 * T - .0001536 * T2 + T3 / 24490000;
195 *M = *M * RAD;
196
197 /* correction for Sun's elliptical orbit */
198 C = (1.915 - .004817 * T - .000014 * T2) * sin(*M) +
199 (.01999 - .000101 * T) * sin(2 * *M) + .00029 * sin(3 * *M);
200
201 /* true anomaly of the Sun */
202 V = *M + C * RAD;
203
204 /* eccentricity of Earth's orbit */
205 E = .01671 - .00004204 * T - .0000001236 * T2;
206
207 /* Sun-Earth distance */
208 *R = .99972 / (1 + E * cos(V)) * 145980000;
209
210 /* true geometric longitude of the Sun */
211 THETA = *L + C;
212
213 /* apparent longitude of the Sun */
214 OM = 125.04 - 1934.1 * T;
215 *LAM0 = THETA - .00569 - .00478 * sin(OM * RAD);
216}
217
218#ifdef TESTCOLONG
219
220/* insure 0 <= *v < r.
221 */
222void
223range (v, r)
224double *v, r;
225{
226 *v -= r*floor(*v/r);
227}
228
229/* To be sure the program is functioning properly, try the test case
230 * 2449992.5 (1 Oct 1995): the colongitude should be 3.69 degrees.
231 */
232int
233main (int ac, char *av[])
234{
235 double jd, lt, lg;
236 double c, k, a;
237
238 if (ac != 2) {
239 fprintf (stderr, "%s: JD\n", av[0]);
240 exit (1);
241 }
242
243 jd = atof(av[1]);
244
245 printf ("Latitude of lunar feature: ");
246 fscanf (stdin, "%lf", &lt);
247 lt = degrad(lt);
248 printf ("Longitude: ");
249 fscanf (stdin, "%lf", &lg);
250 lg = degrad(lg);
251
252 moon_colong (jd, lt, lg, &c, &k, &a);
253
254 printf ("Selenographic colongitude is %g\n", raddeg(c));
255 printf ("The illuminated fraction of the Moon is %g\n", k);
256 printf ("Altitude of Sun above feature is %g\n", raddeg(a));
257
258 return (0);
259}
260
261#endif
262
263/* For RCS Only -- Do Not Edit */
[1719]264static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: mooncolong.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.