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.

About these ads

Pages: 1 2 3

7 Responses to “Moonrise, Moonset”

  1. Mike said

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

  2. programmingpraxis said

    Fixed. Thank you.

  3. Paul said

    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.

  4. Paul said

    Forget my last remark. It runs in codepad.

  5. programmingpraxis said

    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:

    >

  6. Mike said

    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.

  7. programmingpraxis said

    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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 643 other followers

%d bloggers like this: