MoonRise MoonSet calculations

[LwarenceD]

Greetings,

Finally started my Holiday leave from work. (YAY!)

This is my holiday project to try to learn VBA and Excel a little better and stop asking those inane questions =)

The BASIC code is below.

IF anyone has done this before, I'd like to take a look at the xla and see how everything meshes.

If anyone else wants to take a crack at this, perhaps we can share notes;

10 REM        MOONRISE-MOONSET
15 GOSUB 170
20 INPUT "LAT, LONG (DEG)";B5,L5
25 INPUT "TIME ZONE (HRS)";H
30 L5=L5/360: Z0=H/24
35 GOSUB 760: T=(J-2451545)+F
40 GOSUB 245: T=T+Z0
45 REM
50 REM      POSITION LOOP
55 FOR I=1 TO 3
60 GOSUB 495: M(I,1)=A5
65 M(I,2)=D5: M(I,3)=R5: T=T+0.5
70 NEXT
75 IF M(2,1)>M(1,1) THEN 85
80 M(2,1)=M(2,1)+P2
85 IF M(3,1)>M(2,1) THEN 95
90 M(3,1)=M(3,1)+P2
95 Z1=R1*(90.567-41.685/M(2,3))
100 S=SIN(B5*R1): C=COS(B5*R1)
105 Z=COS(Z1): M8=0: W8=0: PRINT
110 A0=M(1,1): D0=M(1,2)
115 FOR C0=0 TO 23
120 P=(C0+1)/24
125 F0=M(1,1):F1=M(2,1):F2=M(3,1)
130 GOSUB 225: A2=F
135 F0=M(1,2):F1=M(2,2):F2=M(3,2)
140 GOSUB 225: D2=F
145 GOSUB 285: A0=A2:D0=D2:V0=V2
150 NEXT
155 GOSUB 450: REM  SPECIAL MSG?
160 END
165 REM
170 REM        CONSTANTS
175 DIM M(3,3)
180 P1=3.14159265: P2=2*P1
185 R1=P1/180: K1=15*R1*1.0027379
190 S$="MOONSET AT  "
195 R$="MOONRISE AT "
200 M1$="NO MOONRISE THIS DATE"
205 M2$="NO MOONSET THIS DATE"
210 M3$="MOON DOWN ALL DAY"
215 M4$="MOON UP ALL DAY"
220 RETURN
225 REM    3-POINT INTERPOLATION
230 A=F1-F0: B=F2-F1-A
235 F=F0+P*(2*A+B*(2*P-1))
240 RETURN
245 REM     LST AT 0H ZONE TIME
250 T0=T/36525
255 S=24110.5+8640184.813*T0
260 S=S+86636.6*Z0+86400*L5
265 S=S/86400: S=S-INT(S)
270 T0=S*360*R1
275 RETURN
280 REM
285 REM  TEST AN HOUR FOR AN EVENT
290 L0=T0+C0*K1: L2=L0+K1
295 IF A2<A0 THEN A2=A2+2*P1
300 H0=L0-A0: H2=L2-A2
305 H1=(H2+H0)/2: REM  HOUR ANGLE
310 D1=(D2+D0)/2: REM  DEC
315 IF C0>0 THEN 325
320 V0=S*SIN(D0)+C*COS(D0)*COS(H0)-Z
325 V2=S*SIN(D2)+C*COS(D2)*COS(H2)-Z
330 IF SGN(V0)=SGN(V2) THEN 440 
335 V1=S*SIN(D1)+C*COS(D1)*COS(H1)-Z
340 A=2*V2-4*V1+2*V0: B=4*V1-3*V0-V2
345 D=B*B-4*A*V0: IF D<0 THEN 440
350 D=SQR(D)
355 IF V0<0 AND V2>0 THEN PRINT R$;
360 IF V0<0 AND V2>0 THEN M8=1
365 IF V0>0 AND V2<0 THEN PRINT S$;
370 IF V0>0 AND V2<0 THEN W8=1
375 E=(-B+D)/(2*A)
380 IF E>1 OR E<0 THEN E=(-B-D)/(2*A)
385 T3=C0+E+1/120: REM ROUND OFF
390 H3=INT(T3): M3=INT((T3-H3)*60)
395 PRINT USING "##:##";H3;M3;
400 H7=H0+E*(H2-H0)
405 N7=-COS(D1)*SIN(H7)
410 D7=C*SIN(D1)-S*COS(D1)*COS(H7)
415 A7=ATN(N7/D7)/R1
420 IF D7<0 THEN A7=A7+180
425 IF A7<0 THEN A7=A7+360
430 IF A7>360 THEN A7=A7-360
435 PRINT USING ",  AZ ###.#";A7
440 RETURN
445 REM
450 REM   SPECIAL MESSAGE ROUTINE
455 IF M8=0 AND W8=0 THEN 475
460 IF M8=0 THEN PRINT M1$
465 IF W8=0 THEN PRINT M2$
470 GOTO 485
475 IF V2<0 THEN PRINT M3$
480 IF V2>0 THEN PRINT M4$
485 RETURN
490 REM
495 REM   FUNDAMENTAL ARGUMENTS
500 L=0.606434+0.03660110129*T
505 M=0.374897+0.03629164709*T
510 F=0.259091+0.03674819520*T
515 D=0.827362+0.03386319198*T
520 N=0.347343-0.00014709391*T
525 G=0.993126+0.00273777850*T
530 L=L-INT(L): M=M-INT(M)
535 F=F-INT(F): D=D-INT(D)
540 N=N-INT(N): G=G-INT(G)
545 L=L*P2: M=M*P2: F=F*P2
550 D=D*P2: N=N*P2: G=G*P2
555 V=0.39558*SIN(F+N)
560 V=V+0.08200*SIN(F)
565 V=V+0.03257*SIN(M-F-N)
570 V=V+0.01092*SIN(M+F+N)
575 V=V+0.00666*SIN(M-F)
580 V=V-0.00644*SIN(M+F-2*D+N)
585 V=V-0.00331*SIN(F-2*D+N)
590 V=V-0.00304*SIN(F-2*D)
595 V=V-0.00240*SIN(M-F-2*D-N)
600 V=V+0.00226*SIN(M+F)
605 V=V-0.00108*SIN(M+F-2*D)
610 V=V-0.00079*SIN(F-N)
615 V=V+0.00078*SIN(F+2*D+N)
620 U=1-0.10828*COS(M)
625 U=U-0.01880*COS(M-2*D)
630 U=U-0.01479*COS(2*D)
635 U=U+0.00181*COS(2*M-2*D)
640 U=U-0.00147*COS(2*M)
645 U=U-0.00105*COS(2*D-G)
650 U=U-0.00075*COS(M-2*D+G)
655 W=0.10478*SIN(M)
660 W=W-0.04105*SIN(2*F+2*N)
665 W=W-0.02130*SIN(M-2*D)
670 W=W-0.01779*SIN(2*F+N)
675 W=W+0.01774*SIN(N)
680 W=W+0.00987*SIN(2*D)
685 W=W-0.00338*SIN(M-2*F-2*N)
690 W=W-0.00309*SIN(G)
695 W=W-0.00190*SIN(2*F)
700 W=W-0.00144*SIN(M+N)
705 W=W-0.00144*SIN(M-2*F-N)
710 W=W-0.00113*SIN(M+2*F+2*N)
715 W=W-0.00094*SIN(M-2*D+G)
720 W=W-0.00092*SIN(2*M-2*D)
725 REM
730 REM    COMPUTE RA, DEC, DIST
735 S=W/SQR(U-V*V)
740 A5=L+ATN(S/SQR(1-S*S))
745 S=V/SQR(U):D5=ATN(S/SQR(1-S*S))
750 R5=60.40974*SQR(U)
755 RETURN
760 REM     CALENDAR --> JD
765 INPUT "Y,M,D ";Y,M,D
770 G=1: IF Y<1582 THEN G=0
775 D1=INT(D): F=D-D1-0.5
780 J=-INT(7*(INT((M+9)/12)+Y)/4)
785 IF G=0 THEN 805
790 S=SGN(M-9): A=ABS(M-9)
795 J3=INT(Y+S*INT(A/7))
800 J3=-INT((INT(J3/100)+1)*3/4)
805 J=J+INT(275*M/9)+D1+G*J3
810 J=J+1721027+2*G+367*Y
815 IF F>=0 THEN 825
820 F=F+1: J=J-1
825 RETURN
900 REM  ***************************
910 REM  THIS PROGRAM COMPUTES THE
920 REM  TIMES OF MOONRISE AND MOON-
930 REM  SET ANYWHERE IN THE WORLD.
940 REM  FROM SKY & TELESCOPE, JULY,
950 REM  1989, PAGE 78.
960 REM  ***************************

Cheers

 

[Oorang]

I started to work my way through this... But I can't remember BASIC syntax in detail... Does anyone remember what this bit is saying:

61 M(I, 1) = A5

 

[Oorang]

I circled back around to this for fun and I now have a version that will execute (see below) but I'm pretty sure it's got some logical errors in yet. If you are still interested in this, if you could post a "truth table of correct input/output" it would be helpful. Also if you happen to know which BASIC this was targeting it would be helpful as I suspect some of the issues are related to Scope.

Public Sub MoonRiseMoonSet()
    '***************************
    'THIS PROGRAM COMPUTES THE
    'TIMES OF MOONRISE AND MOON-
    'SET ANYWHERE IN THE WORLD.
    'FROM SKY & TELESCOPE, JULY,
    '1989, PAGE 78.
    '***************************
    Dim YY(3, 3)
    Dim S$
    'MOONRISE-MOONSET
    GoSub SetUpConstants
    B5 = InputBox("LAT, LONG (DEG)", , "38,85")
    B5 = Split(B5, ",")
    L5 = B5(1)
    B5 = B5(0)
    H = InputBox("TIME ZONE (HRS)", , -5)
    L5 = L5 / 360
    Z0 = H / 24
    GoSub Calendar
    T = (J - 2451545) + F
    GoSub LstAt0hZoneTime
    T = T + Z0
    Rem
    'POSITION LOOP
    For I = 1 To 3
        GoSub FundamentalArguments
        YY(I, 1) = A5
        YY(I, 2) = D5
        YY(I, 3) = R5
        T = T + 0.5
    Next
    If Not (YY(2, 1) > YY(1, 1)) Then
        YY(2, 1) = YY(2, 1) + P2
    End If
    If Not (YY(3, 1) > YY(2, 1)) Then
        YY(3, 1) = YY(3, 1) + P2
    End If
    Z1 = R1 * (90.567 - 41.685 / YY(2, 3))
    S = Sin(B5 * R1)
    C = Cos(B5 * R1)
    Z = Cos(Z1)
    M8 = 0
    W8 = 0
Debug.Print
    A0 = YY(1, 1)
    D0 = YY(1, 2)
    For C0 = 0 To 23
        P = (C0 + 1) / 24
        F0 = YY(1, 1)
        F1 = YY(2, 1)
        F2 = YY(3, 1)
        GoSub ThreePointInterpolation
        A2 = F
        F0 = YY(1, 2)
        F1 = YY(2, 2)
        F2 = YY(3, 2)
        GoSub ThreePointInterpolation
        D2 = F
        GoSub TestAnHourForAnEvent
        A0 = A2
        D0 = D2
        V0 = V2
    Next
    GoSub SpecialMessage    'SPECIAL MSG?
    End
    Rem
    'CONSTANTS
    '175 Dim YY(3, 3)
SetUpConstants:
    P1 = 3.14159265
    P2 = 2 * P1
    R1 = P1 / 180
    K1 = 15 * R1 * 1.0027379
    S$ = "MOONSET AT  "
    R$ = "MOONRISE AT "
    M1$ = "NO MOONRISE THIS DATE"
    M2$ = "NO MOONSET THIS DATE"
    M3$ = "MOON DOWN ALL DAY"
    M4$ = "MOON UP ALL DAY"
    Return
    '3-POINT INTERPOLATION
ThreePointInterpolation:
    A = F1 - F0
    B = F2 - F1 - A
    F = F0 + P * (2 * A + B * (2 * P - 1))
    Return
    'LST AT 0H ZONE TIME
LstAt0hZoneTime:
    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
    Return
    Rem
    'TEST AN HOUR FOR AN EVENT
TestAnHourForAnEvent:
    L0 = T0 + C0 * K1
    L2 = L0 + K1
    If Not A20 Then
        V0 = S * Sin(D0) + C * Cos(D0) * Cos(H0) - Z
        V2 = S * Sin(D2) + C * Cos(D2) * Cos(H2) - Z
        If Sgn(V0) = Sgn(V2) Then
            Return
        End If
    End If
    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 Then
        Return
    End If
    D = Sqr(D)
    If V0 < 0 And V2 > 0 Then
Debug.Print R$;
    End If
    If V0 < 0 And V2 > 0 Then
        M8 = 1
    End If
    If V0 > 0 And V2 < 0 Then
Debug.Print S$;
    End If
    If V0 > 0 And V2 < 0 Then
        W8 = 1
    End If
    E = (-B + D) / (2 * A)
    If E > 1 Or E < 0 Then
        E = (-B - D) / (2 * A)
    End If
    T3 = C0 + E + 1 / 120    'ROUND OFF
    H3 = Int(T3)
    M3 = Int((T3 - H3) * 60)
Debug.Print USING; "##:##"; H3; M3;
    H7 = H0 + E * (H2 - H0)
    N7 = -Cos(D1) * Sin(H7)
    D7 = C * Sin(D1) - S * Cos(D1) * Cos(H7)
    A7 = Atn(N7 / D7) / R1
    If D7 < 0 Then
        A7 = A7 + 180
    End If
    If A7 < 0 Then
        A7 = A7 + 360
    End If
    If A7 > 360 Then
        A7 = A7 - 360
    End If
Debug.Print USING; ",  AZ ###.#"; A7
    Return
    Rem
    'SPECIAL MESSAGE ROUTINE
SpecialMessage:
    If M8 = 0 And W8 = 0 Then
        If V2 < 0 Then
Debug.Print M3$
        End If
        If V2 > 0 Then
Debug.Print M4$
        End If
    Else
        If M8 = 0 Then
Debug.Print M1$
        End If
        If W8 = 0 Then
Debug.Print M2$
        End If
    End If
    Return
    Rem
    'FUNDAMENTAL ARGUMENTS
FundamentalArguments:
    L = 0.606434 + 0.03660110129 * T
    M = 0.374897 + 0.03629164709 * T
    F = 0.259091 + 0.0367481952 * T
    D = 0.827362 + 0.03386319198 * T
    N = 0.347343 - 0.00014709391 * T
    G = 0.993126 + 0.0027377785 * 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.082 * 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.0024 * 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.0188 * 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.0213 * 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.0019 * 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)
    Rem
    'COMPUTE RA, DEC, DIST
    S = W / Sqr(U - V * V)
    A5 = L + Atn(S / Sqr(1 - S * S))
    S = V / Sqr(U)
    D5 = Atn(S / Sqr(1 - S * S))
    R5 = 60.40974 * Sqr(U)
    Return
Calendar:
    'CALENDAR --> JD
    Y = InputBox("Y,M,D ", , Date)
    D = Day(Y)
    M = Month(Y)
    Y = Year(Y)
    G = 1
    If Y < 1582 Then
        G = 0
    End If
    D1 = Int(D)
    F = D - D1 - 0.5
    J = -Int(7 * (Int((M + 9) / 12) + Y) / 4)
    If G <> 0 Then
        S = Sgn(M - 9)
        A = Abs(M - 9)
        J3 = Int(Y + S * Int(A / 7))
        J3 = -Int((Int(J3 / 100) + 1) * 3 / 4)
    End If
    J = J + Int(275 * M / 9) + D1 + G * J3
    J = J + 1721027 + 2 * G + 367 * Y
    If F < 0 Then
        F = F + 1
        J = J - 1
    End If
    Return
End Sub