;- https://gml.noaa.gov/grad/solcalc/solareqns.PDF
;-
;- https://gml.noaa.gov/grad/solcalc/sunrise.html
;- https://www.timeanddate.com/sun/australia/brisbane
this was in GWBASIC ~1982
Code: Select all
12900 '----SUNRISE / SUNSET ---------------------------------------------------
12950 'Basel 47.33 -07.35, south and east negative
13000 'Zurich 47.22 -08.32
13050 'Bern 46.56 -07.26
13060 'VALLETTA 35.7 -14.3
13090 CLS
13100 '
13150 DIM N(12)
13200 D1=46.56:D2=-7.26 'BERN
13210 '
13250 PL=3.141592654#/26:J=57.29577951#
13300 GOSUB 16750
13350 LA = D1
13400 IF LA < 0 THEN LA = LA + 180
13450 IF D2 < 0 THEN D2 = D2 + 360
13500 LO = FIX(D2/15)*15 :REM finds time zone beginning
13550 TD=(D2-LO)/15
13600 '
13650 M$=MID$(DATE$,1,2):M=VAL(M$)
13700 DA$=MID$(DATE$,4,2):DA=VAL(DA$)
13750 '
13800 FOR I=1 TO 12: READ N(I):NEXT I
13850 DATA 0,31,59,90,120,151
13900 DATA 181,212,243,273,304,334
13950 X=(N(M)+DA)/7
14000 '
14050 D=.4560001-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)
14100 '
14150 'LOCATE 3,2
14200 PRINT"DECLINATION OF SUN:";
14250 PRINT USING"###.#";D;
14300 PRINT" DEGREES"
14350 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)
14400 '
14450 PRINT"EQUATION OF TIME:";
14500 PRINT USING"###.#";E;
14550 PRINT" MINUTES"
14600 CL=COS(LA/J):SD=SIN(D/J):CD=COS(D/J):Y=SD/CL
14650 IF ABS(Y)=>1 THEN PRINT"NO SUNRISE OR SUNSET":GOTO 17150
14700 Z = 90 - J*ATN(Y/SQR(1-Y*Y))
14750 PRINT"AZIMUTH OF SUNRISE:";
14800 PRINT USING"####.#";ABS(Z);
14850 PRINT" DEGREES"
14900 PRINT"AZIMUTH OF SUNSET: ";
14950 PRINT USING"####.#";360-ABS(Z);
15000 PRINT" DEGREES"
15050 ST=SIN(Z/J)/CD
15100 IF ABS(ST)>=1 THEN T=6:TT=6:GOTO 15300
15150 CT=SQR(1-ST*ST)
15200 T=J/15*ATN(ST/CT)
15250 TT=T
15300 'x
15350 IF D<0 AND LA<90 THEN T=12-T:TT=T
15400 IF D>0 AND LA>90 THEN T=12-T:TT=T
15450 T=T+TD-E/60-.04
15500 GOSUB 16050
15550 'LOCATE 3,2
15600 'PRINT "TIME OF SUNRISE BERN:";T1$;":";T2$;" ";T$;"L.T. ";GM$;":";T2$;" GM";
15650 PRINT "SWITZERLAND BERN SUNRISE= ";T1$;":";T2$;" ";
15700 T=12-TT:T=T+TD-E/60+.04
15750 CNT=1
15800 T12=12
15850 GOSUB 16050
15900 'PRINT "TIME OF SUNSET= ";T1$;":";T2$;" ";T$;"L.T. ";GM$;":";T2$;" GM";" ";SD$
15950 PRINT " SUNSET= ";T1$;":";T2$;" ";SD$
16000 GOTO 17150
16050 'x
16100 T1=INT(T):T2=T-T1
16150 IF SD$="ST" THEN LET T1=T1+1
16200 T1$=STR$(T1+T12):T2=INT((T2*600+5)/10):IF T2=60 THEN 16250 ELSE 16300
16250 T2=0:T1=T1+1:T1$=STR$(T1+T12)
16300 T2$=STR$(T2):T2$=RIGHT$(T2$,LEN(T2$)-1)
16350 IF INT(T2)<10 THEN T2$="0"+T2$
16400 GM = FIX(D2/15) :REM calculate difference between GM and local time
16450 IF CNT = 0 THEN GM = VAL(T1$)+GM :REM GMT for sunrise
16500 IF CNT > 0 THEN GM = VAL(T1$)+12+GM :REM GMT for sunset
16550 IF GM +(VAL(T2$)/60)> 24 THEN GM = GM - 24
16600 GM$ = STR$(GM) :GM$ = RIGHT$("0"+GM$,2)
16650 RETURN
16700 '
16750 'x This subroutine converts DD.MM input to DD.DD
16800 DEGTMP = (ABS(D1)-ABS(FIX(D1))) *100/60
16850 D1 = (FIX(ABS(D1))+DEGTMP)*SGN(D1)
16900 DEGTMP = (ABS(D2)-ABS(FIX(D2))) *100/60
16950 D2 = (FIX(ABS(D2))+DEGTMP)*SGN(D2)
17000 RETURN
17050 '
17100 '
17150 'END
30030 INPUT "PUSH ENTER";nothing
30040 A$=INKEY$
30050 IF LEN (A$)=0 THEN 30040
30060 IF ASC(A$)=13 THEN 31040
30080 'A$=INKEY$:IF A$="" THEN 30080
31030 'SHELL "CD\":SYSTEM
31040 PRINT "END"