## Moonrise, Moonset

### July 1, 2014

I translated the program to Python, but never got it to work. And I don’t have a BASIC interpreter that can run the Sky & Telescope program. I would appreciate it if anyone could tell me what I have done wrong, or if there is a problem in the original program:

# moonrise / moonset from __future__ import division from math import sqrt, pi, sin, cos, atan def signum(x): if x < 0: return -1 if x > 0: return 1 return 0 def moon(B5, L5, H, Mo, D, Y): # lat, long, hours west of GMT, month, day, year # positive latitude is north, negative latitude is south # positive longitude is east, negative longitude is west # lat and long given in degrees and decimals of a degree # moon(38.6272, -90.1978, 6, 6, 30, 2014) # st louis, today # GOSUB 170 constants Ma = [[0 for x in range(4)] for y in range(4)] # we dimension to 4 instead of 3 because BASIC # arrays are 1-based and Python arrays are 0-based P2 = 2 * pi R1 = pi / 180 K1 = 15 * R1 * 1.0027379 L5 = L5 / 360 Z0 = H / 24 # GOSUB 760 calendar to julian date G = 1 if Y < 1582: G = 0 D1 = int(D); F = D - D1 - 0.5 J = -1 * int(7 * (int((Mo+9)/12) + Y) / 4) J3 = 0 # not initialized in BASIC program if G <> 0: S = signum(Mo - 9) A = abs(Mo - 9) J3 = int(Y + S * int(A / 7)) J3 = -1 * int((int(J3 / 100) + 1) * 3/4) J = J + int(275 * Mo / 9) + D1 + G * J3 J = J + 1721027 + 2*G + 367*Y if F < 0: F = F + 1 J = J - 1 T = (J - 2451545) + F # GOSUB 245 lunar sidereal time at GMT time zone T0 = T / 36525 S = 24110.5 + 8640184.813 * T0 S = S + 86636.6 * Z0 + 86400 * L5 S = S / 86400 S = S - int(S) T0 = S * 360 * R1 T = T + Z0 # POSITION LOOP for I in range(1, 4): # GOSUB 495 fundamental arguments L = 0.606434 + 0.03660110129 * T M = 0.374897 + 0.03629164709 * T F = 0.259091 + 0.03674819520 * T D = 0.827362 + 0.03386319198 * T N = 0.347343 - 0.00014709391 * T G = 0.993126 + 0.00273777850 * T L = L - int(L); M = M - int(M) F = F - int(F); D = D - int(D) N = N - int(N); G = G - int(G) L = L * P2; M = M * P2; F = F * P2 D = D * P2; N = N * P2; G = G * P2 V = 0.39558 * sin(F + N) V = V + 0.08200 * sin(F) V = V + 0.03257 * sin(M - F - N) V = V + 0.01092 * sin(M + F + N) V = V + 0.00666 * sin(M - F) V = V - 0.00644 * sin(M + F - 2*D + N) V = V - 0.00331 * sin(F - 2*D + N) V = V - 0.00304 * sin(F - 2*D) V = V - 0.00240 * sin(M - F - 2*D - N) V = V + 0.00226 * sin(M + F) V = V - 0.00108 * sin(M + F - 2*D) V = V - 0.00079 * sin(F - N) V = V + 0.00078 * sin(F + 2*D + N) U = 1 - 0.10828 * cos(M) U = U - 0.01880 * cos(M - 2*D) U = U - 0.01479 * cos(2*D) U = U + 0.00181 * cos(2*M - 2*D) U = U - 0.00147 * cos(2*M) U = U - 0.00105 * cos(2*D - G) U = U - 0.00075 * cos(M - 2*D + G) W = 0.10478 * sin(M) W = W - 0.04105 * sin(2*F + 2*N) W = W - 0.02130 * sin(M - 2*D) W = W - 0.01779 * sin(2*F + N) W = W + 0.01774 * sin(N) W = W + 0.00987 * sin(2*D) W = W - 0.00338 * sin(M - 2*F - 2*N) W = W - 0.00309 * sin(G) W = W - 0.00190 * sin(2*F) W = W - 0.00144 * sin(M + N) W = W - 0.00144 * sin(M - 2*F - N) W = W - 0.00113 * sin(M + 2*F + 2*N) W = W - 0.00094 * sin(M - 2*D + G) W = W - 0.00092 * sin(2*M - 2*D) # compute right ascension, declination, distance S = W / sqrt(U - V*V) A5 = L + atan(S / sqrt(1 - S*S)) S = V / sqrt(U); D7 = atan(S / sqrt(1 - S*S)) R5 = 60.40974 * sqrt(U) Ma[I][1] = A5 Ma[I][2] = D7 Ma[I][3] = R5 T = T + 0.5 if Ma[2][1] <= Ma[1][1]: Ma[2][1] = Ma[2][1] + P2 if Ma[3][1] <= Ma[2][1]: Ma[3][1] = Ma[3][1] + P2 Z1 = R1 * (90.567 - 41.685/Ma[2][3]) S = sin(B5 * R1); C = cos(B5 * R1) Z = cos(Z1); M8 = 0; W8 = 0 A0 = Ma[1][1]; D0 = Ma[1][2] V0 = 0 # not initialized in BASIC program for C0 in range(0, 24): P = (C0 + 1) / 24 F0 = Ma[1][1]; F1 = Ma[2][1]; F2 = Ma[3][1] A = F1 - F0; B = F2 - F1 - A F = F0 + P * (2*A + B*(2*P-1)) A2 = F F0 = Ma[1][2]; F1 = Ma[2][2]; F2 = Ma[3][2] A = F1 - F0; B = F2 - F1 - A F = F0 + P * (2*A + B*(2*P-1)) D2 = F # GOSUB 285 test an hour for an event L0 = T0 + C0 * K1; L2 = L0 + K1 if A2 <> 0: A2 = A2 + 2*pi H0 = L0 - A0; H2 = L2 - A2 H1 = (H2 + H0) / 2 # hour angle D1 = (D2 + D0) / 2 # declination if C0 <= 0: V0 = S * sin(D0) + C * cos(D0) * cos(H0) - Z V2 = S * sin(D2) + C * cos(D2) * cos(H2) - Z if signum(V0) <> signum(V2): V1 = S * sin(D1) + C * cos(D1) * cos(H1) - Z A = 2*V2 - 4*V1 + 2*V0; B = 4*V1 - 3*V0 - V2 D = B*B - 4*A*V0 if D >= 0: D = sqrt(D) if V0 0: print "MOONRISE AT", if V0 0: M8 = 1 if V0 > 0 and V2 0 and V2 1 or E < 0: E = (-1*B - D) / (2 * A) T3 = C0 + E + 1/120 H3 = int(T3); M3 = int((T3 - H3) * 60) print H3, ":", M3, H7 = H0 + E * (H2 - H0) N7 = -1 * cos(D1) * sin(H7) D7 = C * sin(D1) - S * cos(D1) * cos(H7) A7 = atan(N7 / D7) / R1 if D7 < 0: A7 = A7 + 180 if A7 360: A7 = A7 - 360 print ", AZ", A7 A0 = A2; D0 = D2; V0 = V2 # GOSUB 450 special message if M8 0 or W8 0: if M8 == 0: print "NO MOONRISE THIS DATE" if W8 == 0: print "NO MOONSET THIS DATE" else: if V2 0: print "MOON UP ALL DAY"

You can run the program at http://programmingpraxis.codepad.org/98sCt5at.

Advertisements

I think WordPress ate some less than signs. See line numbers 355, 360, 365, and 425.

Fixed. Thank you.

I see 2 problems so far: division and a variable F3 is introduced, that is never used. F3 should be F2. If I change all division by floating point divivion and change varaible F3 into F2, the program gives a sunrise and sunset, which is one hour earlier than the right answer. Code is on Codepad. I do not have time now to make it run in codepad, which does not run Python 2.7. I have to quit for the Soccer Worldcup.

Forget my last remark. It runs in codepad.

Thanks for finding that.

Regarding the World Cup: My statistics were low this morning, but zoomed up after Argentina scored its goal. The effect of the World Cup was very obvious. There must be a lot of my readers who like soccer.

Like you, I will be watching the USA team in a few minutes.

On Tue, Jul 1, 2014 at 2:41 PM, Programming Praxis wrote:

>

Parts of the BASIC code appear to be missing. For example, BASIC version of the Julian date subroutine (GOSUB 760) is completely different than your Python code.

Well, it’s the same arithmetic. The BASIC program assumes that you want moonrise and moonset for today, and calculates that first, then asks you to input a different day. My program just takes the month, day and year as input parameters.