6 #define DegPerRad 57.29577951308232087680
7 #define RadPerDeg 0.01745329251994329576
9 extern double Glon
, SinGlat
, CosGlat
, TimeZone
;
12 void MoonRise(int year
, int month
, int day
, double LocalHour
, double *UTRise
, double *UTSet
){
14 double UT
, ym
, y0
, yp
, SinH0
;
15 double xe
, ye
, z1
, z2
, SinH(), hour24();
18 SinH0
= sin( 8.0/60.0 * RadPerDeg
);
25 ym
= SinH(year
, month
, day
, UT
-1.0) - SinH0
;
27 while ( (UT
<= 24.0+TimeZone
) ) {
29 y0
= SinH(year
, month
, day
, UT
) - SinH0
;
30 yp
= SinH(year
, month
, day
, UT
+1.0) - SinH0
;
32 Interp(ym
, y0
, yp
, &xe
, &ye
, &z1
, &z2
, &nz
);
66 *UTRise
= hour24(*UTRise
);
73 *UTSet
= hour24(*UTSet
);
81 void UTTohhmm(double UT
, int *h
, int *m
){
89 *m
= (int)((UT
-(double)(*h
))*60.0+0.5);
92 * If it was 23:60 this should become 24:00
93 * I prefer this designation to 00:00. So dont reset h to 0 when it goes above 24.
107 void Interp(double ym
, double y0
, double yp
, double *xe
, double *ye
, double *z1
, double *z2
, int *nz
){
109 double a
, b
, c
, d
, dx
;
116 *ye
= (a
*(*xe
) + b
) * (*xe
) + c
;
120 dx
= 0.5*sqrt(d
)/fabs(a
);
123 if (fabs(*z1
) <= 1.0) *nz
+= 1;
124 if (fabs(*z2
) <= 1.0) *nz
+= 1;
125 if (*z1
< -1.0) *z1
= *z2
;
135 double SinH(int year
, int month
, int day
, double UT
){
137 double TU
, frac(), jd();
138 double RA_Moon
, DEC_Moon
, gmst
, lmst
, Tau
, Moon();
141 TU
= (jd(year
, month
, day
, UT
) - 2451545.0)/36525.0;
143 /* this is more accurate, but wasteful for this -- use low res approx.
146 Moon(TU, &LambdaMoon, &BetaMoon, &R, &AGE);
147 LambdaMoon *= RadPerDeg;
148 BetaMoon *= RadPerDeg;
149 epsilon = (23.43929167 - 0.013004166*TU - 1.6666667e-7*TU2 - 5.0277777778e-7*TU3)*RadPerDeg;
150 RA_Moon = angle2pi(atan2(sin(LambdaMoon)*cos(epsilon)-tan(BetaMoon)*sin(epsilon), cos(LambdaMoon)));
151 DEC_Moon = asin( sin(BetaMoon)*cos(epsilon) + cos(BetaMoon)*sin(epsilon)*sin(LambdaMoon));
154 MiniMoon(TU
, &RA_Moon
, &DEC_Moon
);
155 RA_Moon
*= 15.0*RadPerDeg
;
156 DEC_Moon
*= RadPerDeg
;
159 * Compute Greenwich Mean Sidereal Time (gmst)
161 UT
= 24.0*frac( UT
/24.0 );
162 /* this is for the ephemeris meridian???
163 gmst = 6.697374558 + 1.0027379093*UT + (8640184.812866+(0.093104-6.2e-6*TU)*TU)*TU/3600.0;
165 gmst
= UT
+ 6.697374558 + (8640184.812866+(0.093104-6.2e-6*TU
)*TU
)*TU
/3600.0;
166 lmst
= 24.0*frac( (gmst
-Glon
/15.0) / 24.0 );
168 Tau
= 15.0*lmst
*RadPerDeg
- RA_Moon
;
169 return( SinGlat
*sin(DEC_Moon
) + CosGlat
*cos(DEC_Moon
)*cos(Tau
) );