Updated: 09/10/03

Geo. Anderson
(please remove the letters "NO_SPAM" from my reply-to address)
*******************************************************
Remember, the people who paint lane stripes on the road
of life probably have different objectives than you do.
-------------------------------------------------------
Like sports cars? Check out my book:
"WINNING! A Race Driver's Handbook
at http://durequip.com/winning/
*******************************************************

'**********************************************************************
'*************** BASIC SUNSET CALCULATION PROGRAM *******************
'** By: Thomas Laureanno 10/02/92 **
'** **
'** LAT = 41.25 LONG = 71.5 ' For Newport, Rhode Island **
'** **
'** One note regarding azimuth angles: In the SOUTHERN hemisphere **
'** this program assumes that the South Pole is zero degrees, and **
'** the azimuth angle is measured COUNTER-clockwise through East. **
'** Therefore, an azimuth angle of 108 degrees is NORTH of East. **
'** When Lat and Long are input, use a negative value for Southern **
'** Latitudes and for Eastern Longitudes. **
'********************************************************************
'********************************************************************

'----------------- Program Begins ---------------------------
CLS : KEY OFF: CLEAR
PRINT "This program finds the declination of the sun, the equation"
PRINT "of time, the azimuth angles of sunrise and sunset, and the"
PRINT "times of sunrise and sunset for any point on earth."
PRINT
PRINT "Input eastern longitudes and southern latitudes as NEGATIVE."
PRINT
PRINT
REM
DIM N(12)
PL = 3.14159 / 26: J = 57.2958
INPUT "ENTER LATITUDE (FORMAT DD.MM)"; D1
INPUT "ENTER LONGITUDE (FORMAT DD.MM)"; D2
'This subroutine converts DD.MM input to DD.DD
DEGTMP = (ABS(D1) - ABS(FIX(D1))) * 100 / 60
D1 = (FIX(ABS(D1)) + DEGTMP) * SGN(D1)
DEGTMP = (ABS(D2) - ABS(FIX(D2))) * 100 / 60
D2 = (FIX(ABS(D2)) + DEGTMP) * SGN(D2)
LA = D1
IF LA < 0 THEN LA = LA + 180
IF D2 < 0 THEN D2 = D2 + 360
LO = FIX(D2 / 15) * 15: REM finds time zone beginning
TD = (D2 - LO) / 15
INPUT "ENTER MONTH,DAY (11,25)"; M, DA
IF M = 0 AND DA = 0 THEN
   M = VAL(MID$(DATE$, 1, 2))
   DA = VAL(MID$(DATE$, 4, 2))
END IF
FOR I = 1 TO 12: READ N(I): NEXT I
DATA 0,31,59,90,120,151
DATA 181,212,243,273,304,334
X = (N(M) + DA) / 7
D = .456 - 22.195 * COS(PL * X) - .43 * COS(2 * PL * X) - .156 * COS(3 * PL * X) + 3.83 * SIN(PL * X) + .06 * SIN(2 * PL * X) - .082 * SIN(3 * PL * X)
PRINT
PRINT "DECLINATION OF SUN:";
PRINT USING "###.#"; D;
PRINT " DEGREES"
E = 8.000001E-03 + .51 * COS(PL * X) - 3.197 * COS(2 * PL * X) - .106 * COS(3 * PL * X) - .15 * COS(4 * PL * X) - 7.317001 * SIN(PL * X) - 9.471001 * SIN(2 * PL * X) - .391 * SIN(3 * PL * X) - .242 * SIN(4 * PL * X)
PRINT "EQUATION OF TIME:";
PRINT USING "###.#"; E;
PRINT " MINUTES"
CL = COS(LA / J): SD = SIN(D / J): CD = COS(D / J): Y = SD / CL
IF ABS(Y) >= 1 THEN PRINT "NO SUNRISE OR SUNSET": END
Z = 90 - J * ATN(Y / SQR(1 - Y * Y))
PRINT "AZIMUTH OF SUNRISE:";
PRINT USING "####.#"; ABS(Z);
PRINT " DEGREES"
PRINT "AZIMUTH OF SUNSET: ";
PRINT USING "####.#"; 360 - ABS(Z);
PRINT " DEGREES"
ST = SIN(Z / J) / CD
IF ABS(ST) >= 1 THEN
   T = 6
   TT = 6
   'GOTO 590
ELSE
   CT = SQR(1 - ST * ST)
   T = J / 15 * ATN(ST / CT)
   TT = T
END IF

IF D < 0 AND LA < 90 THEN T = 12 - T: TT = T

IF D > 0 AND LA > 90 THEN T = 12 - T: TT = T
T = T + TD - E / 60 - .04

T1 = INT(T): T2 = T - T1: T1$ = STR$(T1): T2 = INT((T2 * 600 + 5) / 10)
T2$ = STR$(T2): T2$ = RIGHT$(T2$, LEN(T2$) - 1)
IF INT(T2) < 10 THEN T2$ = "0" + T2$
GM = FIX(D2 / 15): REM calculate difference between GM and local time
IF CNT = 0 THEN GM = VAL(T1$) + GM: REM GMT for sunrise
IF CNT > 0 THEN GM = VAL(T1$) + 12 + GM: REM GMT for sunset
IF GM + (VAL(T2$) / 60) > 24 THEN GM = GM - 24
GM$ = STR$(GM): GM$ = RIGHT$("0" + GM$, 2)
PRINT
PRINT "TIME OF SUNRISE:"; T1$; ":"; T2$; " "; T$; "L.T. "; GM$; ":"; T2$; " GM"
T = 12 - TT: T = T + TD - E / 60 + .04
CNT = 1
T1 = INT(T): T2 = T - T1: T1$ = STR$(T1): T2 = INT((T2 * 600 + 5) / 10)
T2$ = STR$(T2): T2$ = RIGHT$(T2$, LEN(T2$) - 1)
IF INT(T2) < 10 THEN T2$ = "0" + T2$
GM = FIX(D2 / 15): REM calculate difference between GM and local time
IF CNT = 0 THEN GM = VAL(T1$) + GM: REM GMT for sunrise
IF CNT > 0 THEN GM = VAL(T1$) + 12 + GM: REM GMT for sunset
IF GM + (VAL(T2$) / 60) > 24 THEN GM = GM - 24
GM$ = STR$(GM): GM$ = RIGHT$("0" + GM$, 2)
PRINT "TIME OF SUNSET: "; T1$; ":"; T2$; " "; T$; "L.T. "; GM$; ":"; T2$; " GM"
END

10 ' Sunrise-Sunset
20 GOSUB 300
30 INPUT "Lat, Long (deg)";B5,L5
40 INPUT "Time zone (hrs)";H
50 L5=L5/360: Z0=H/24
60 GOSUB 1170: T=(J-2451545)+F
70 TT=T/36525+1: ' TT = centuries
80 ' from 1900.0
90 GOSUB 410: T=T+Z0
100 '
110 ' Get Sun's Position
120 GOSUB 910: A(1)=A5: D(1)=D5
130 T=T+1
140 GOSUB 910: A(2)=A5: D(2)=D5
150 IF A(2) < A(1) THEN A(2)=A(2)+P2
160 Z1=DR*90.833: ' Zenith dist.
170 S=SIN(B5*DR): C=COS(B5*DR)
180 Z=COS(Z1): M8=0: W8=0: PRINT
190 A0=A(1): D0=D(1)
200 DA=A(2)-A(1): DD=D(2)-D(1)
210 FOR C0=0 TO 23
220 P=(C0+1)/24
230 A2=A(1)+P*DA: D2=D(1)+P*DD
240 GOSUB 490
250 A0=A2: D0=D2: V0=V2
260 NEXT
270 GOSUB 820: ' Special msg?
280 END
290 '
300 ' Constants
310 DIM A(2),D(2)
320 P1=3.14159265: P2=2*P1
330 DR=P1/180: K1=15*DR*1.0027379
340 S$="Sunset at "
350 R$="Sunrise at "
360 M1$="No sunrise this date"
370 M2$="No sunset this date"
380 M3$="Sun down all day"
390 M4$="Sun up all day"
400 RETURN
410 ' LST at 0h zone time
420 T0=T/36525
430 S=24110.5+8640184.813*T0
440 S=S+86636.6*Z0+86400*L5
450 S=S/86400: S=S-INT(S)
460 T0=S*360*DR
470 RETURN
480 '
490 ' Test an hour for an event
500 L0=T0+C0*K1: L2=L0+K1
510 H0=L0-A0: H2=L2-A2
520 H1=(H2+H0)/2: ' Hour angle,
530 D1=(D2+D0)/2: ' declination,
540 ' at half hour
550 IF C0>0 THEN 570
560 V0=S*SIN(D0)+C*COS(D0)*COS(H0)-Z
570 V2=S*SIN(D2)+C*COS(D2)*COS(H2)-Z
580 IF SGN(V0)=SGN(V2) THEN 800
590 V1=S*SIN(D1)+C*COS(D1)*COS(H1)-Z
600 A=2*V2-4*V1+2*V0: B=4*V1-3*V0-V2
610 D=B*B-4*A*V0: IF D<0 THEN 800
620 D=SQR(D)
630 IF V0<0 AND V2>0 THEN PRINT R$;
640 IF V0<0 AND V2>0 THEN M8=1
650 IF V0>0 AND V2<0 THEN PRINT S$;
660 IF V0>0 AND V2<0 THEN W8=1
670 E=(-B+D)/(2*A)
680 IF E>1 OR E<0 THEN E=(-B-D)/(2*A)
690 T3=C0+E+1/120: ' Round off
700 H3=INT(T3): M3=INT((T3-H3)*60)
710 PRINT USING "##:##";H3;M3;
720 H7=H0+E*(H2-H0)
730 N7=-COS(D1)*SIN(H7)
740 D7=C*SIN(D1)-S*COS(D1)*COS(H7)
750 AZ=ATN(N7/D7)/DR
760 IF D7<0 THEN AZ=AZ+180
770 IF AZ<0 THEN AZ=AZ+360
780 IF AZ>360 THEN AZ=AZ-360
790 PRINT USING ", azimuth ###.#";AZ
800 RETURN
810 '
820 ' Special-message routine
830 IF M8=0 AND W8=0 THEN 870
840 IF M8=0 THEN PRINT M1$
850 IF W8=0 THEN PRINT M2$
860 GOTO 890
870 IF V2<0 THEN PRINT M3$
880 IF V2>0 THEN PRINT M4$
890 RETURN
900 '
910 ' Fundamental arguments
920 ' (Van Flandern &
930 ' Pulkkinen, 1979)
940 L=.779072+.00273790931*T
950 G=.993126+.0027377785*T
960 L=L-INT(L): G=G-INT(G)
970 L=L*P2: G=G*P2
980 V=.39785*SIN(L)
990 V=V-.01000*SIN(L-G)
1000 V=V+.00333*SIN(L+G)
1010 V=V-.00021*TT*SIN(L)
1020 U=1-.03349*COS(G)
1030 U=U-.00014*COS(2*L)
1040 U=U+.00008*COS(L)
1050 W=-.00010-.04129*SIN(2*L)
1060 W=W+.03211*SIN(G)
1070 W=W+.00104*SIN(2*L-G)
1080 W=W-.00035*SIN(2*L+G)
1090 W=W-.00008*TT*SIN(G)
1100 '
1110 ' Compute Sun's RA and Dec
1120 S=W/SQR(U-V*V)
1130 A5=L+ATN(S/SQR(1-S*S))
1140 S=V/SQR(U):D5=ATN(S/SQR(1-S*S))
1150 R5=1.00021*SQR(U)
1160 RETURN
1165 '
1170 ' Calendar --> JD
1180 INPUT "Year, Month, Day";Y,M,D
1190 G=1: IF Y<1583 THEN G=0
1200 D1=INT(D): F=D-D1-.5
1210 J=-INT(7*(INT((M+9)/12)+Y)/4)
1220 IF G=0 THEN 1260
1230 S=SGN(M-9): A=ABS(M-9)
1240 J3=INT(Y+S*INT(A/7))
1250 J3=-INT((INT(J3/100)+1)*3/4)
1260 J=J+INT(275*M/9)+D1+G*J3
1270 J=J+1721027+2*G+367*Y
1280 IF F>=0 THEN 1300
1290 F=F+1: J=J-1
1300 RETURN
1310 '
1320 ' This program by Roger W. Sinnott calculates the times of sunrise
1330 ' and sunset on any date, accurate to the minute within several
1340 ' centuries of the present. It correctly describes what happens in the
1350 ' arctic and antarctic regions, where the Sun may not rise or set on
1360 ' a given date. Enter north latitudes positive, west longitudes
1370 ' negative. For the time zone, enter the number of hours west of
1380 ' Greenwich (e.g., 5 for EST, 4 for EDT). The calculation is
1390 ' discussed in Sky & Telescope for August 1994, page 84.

Go to Top of Page Return to Tom's X-10 Ideas Webpage Return to Tom's X-10 Webpage