1 #include <math.h>
2 #include "MoonRise.h"
3 #include "Moon.h"
4
5 #define DegPerRad 57.29577951308232087680
6 #define RadPerDeg 0.01745329251994329576
7
8 extern double Glon, SinGlat, CosGlat, TimeZone;
9
10
MoonRise(int year,int month,int day,double LocalHour,double * UTRise,double * UTSet)11 void MoonRise(int year, int month, int day, double LocalHour, double *UTRise, double *UTSet){
12
13 double UT, ym, y0, yp, SinH0;
14 double xe, ye, z1, z2, SinH(), hour24();
15 int Rise, Set, nz;
16
17 SinH0 = sin( 8.0/60.0 * RadPerDeg );
18
19
20 UT = 1.0+TimeZone;
21 *UTRise = -999.0;
22 *UTSet = -999.0;
23 Rise = Set = 0;
24 ym = SinH(year, month, day, UT-1.0) - SinH0;
25
26 while ( (UT <= 24.0+TimeZone) ) {
27
28 y0 = SinH(year, month, day, UT) - SinH0;
29 yp = SinH(year, month, day, UT+1.0) - SinH0;
30
31 Interp(ym, y0, yp, &xe, &ye, &z1, &z2, &nz);
32
33 switch(nz){
34
35 case 0:
36 break;
37 case 1:
38 if (ym < 0.0){
39 *UTRise = UT + z1;
40 Rise = 1;
41 } else {
42 *UTSet = UT + z1;
43 Set = 1;
44 }
45 break;
46 case 2:
47 if (ye < 0.0){
48 *UTRise = UT + z2;
49 *UTSet = UT + z1;
50 } else {
51 *UTRise = UT + z1;
52 *UTSet = UT + z2;
53 }
54 Rise = 1;
55 Set = 1;
56 break;
57 }
58 ym = yp;
59 UT += 2.0;
60
61 }
62
63 if (Rise){
64 *UTRise -= TimeZone;
65 *UTRise = hour24(*UTRise);
66 } else {
67 *UTRise = -999.0;
68 }
69
70 if (Set){
71 *UTSet -= TimeZone;
72 *UTSet = hour24(*UTSet);
73 } else {
74 *UTSet = -999.0;
75 }
76
77 }
78
79
UTTohhmm(double UT,int * h,int * m)80 void UTTohhmm(double UT, int *h, int *m){
81
82
83 if (UT < 0.0) {
84 *h = -1.0;
85 *m = -1.0;
86 } else {
87 *h = (int)UT;
88 *m = (int)((UT-(double)(*h))*60.0+0.5);
89 if (*m == 60) {
90 /*
91 * If it was 23:60 this should become 24:00
92 * I prefer this designation to 00:00. So dont reset h to 0 when it goes above 24.
93 */
94 *h += 1;
95 *m = 0;
96 }
97 }
98
99 }
100
101
102
103
104
105
Interp(double ym,double y0,double yp,double * xe,double * ye,double * z1,double * z2,int * nz)106 void Interp(double ym, double y0, double yp, double *xe, double *ye, double *z1, double *z2, int *nz){
107
108 double a, b, c, d, dx;
109
110 *nz = 0;
111 a = 0.5*(ym+yp)-y0;
112 b = 0.5*(yp-ym);
113 c = y0;
114 *xe = -b/(2.0*a);
115 *ye = (a*(*xe) + b) * (*xe) + c;
116 d = b*b - 4.0*a*c;
117
118 if (d >= 0){
119 dx = 0.5*sqrt(d)/fabs(a);
120 *z1 = *xe - dx;
121 *z2 = *xe+dx;
122 if (fabs(*z1) <= 1.0) *nz += 1;
123 if (fabs(*z2) <= 1.0) *nz += 1;
124 if (*z1 < -1.0) *z1 = *z2;
125 }
126
127
128
129 }
130
131
132
133
SinH(int year,int month,int day,double UT)134 double SinH(int year, int month, int day, double UT){
135
136 double TU, frac(), jd();
137 double RA_Moon, DEC_Moon, gmst, lmst, Tau, Moon();
138 double angle2pi();
139
140 TU = (jd(year, month, day, UT) - 2451545.0)/36525.0;
141
142 /* this is more accurate, but wasteful for this -- use low res approx.
143 TU2 = TU*TU;
144 TU3 = TU2*TU;
145 Moon(TU, &LambdaMoon, &BetaMoon, &R, &AGE);
146 LambdaMoon *= RadPerDeg;
147 BetaMoon *= RadPerDeg;
148 epsilon = (23.43929167 - 0.013004166*TU - 1.6666667e-7*TU2 - 5.0277777778e-7*TU3)*RadPerDeg;
149 RA_Moon = angle2pi(atan2(sin(LambdaMoon)*cos(epsilon)-tan(BetaMoon)*sin(epsilon), cos(LambdaMoon)));
150 DEC_Moon = asin( sin(BetaMoon)*cos(epsilon) + cos(BetaMoon)*sin(epsilon)*sin(LambdaMoon));
151 */
152
153 MiniMoon(TU, &RA_Moon, &DEC_Moon);
154 RA_Moon *= 15.0*RadPerDeg;
155 DEC_Moon *= RadPerDeg;
156
157 /*
158 * Compute Greenwich Mean Sidereal Time (gmst)
159 */
160 UT = 24.0*frac( UT/24.0 );
161 /* this is for the ephemeris meridian???
162 gmst = 6.697374558 + 1.0027379093*UT + (8640184.812866+(0.093104-6.2e-6*TU)*TU)*TU/3600.0;
163 */
164 gmst = UT + 6.697374558 + (8640184.812866+(0.093104-6.2e-6*TU)*TU)*TU/3600.0;
165 lmst = 24.0*frac( (gmst-Glon/15.0) / 24.0 );
166
167 Tau = 15.0*lmst*RadPerDeg - RA_Moon;
168 return( SinGlat*sin(DEC_Moon) + CosGlat*cos(DEC_Moon)*cos(Tau) );
169
170
171 }
172
173