Notice. New forum software under development. It's going to miss a few functions and look a bit ugly for a while, but I'm working on it full time now as the old forum was too unstable. Couple days, all good. If you notice any issues, please contact me.

Forum Index : Microcontroller and PC projects : Sunrise Sunset times

Author Message
TassyJim

Guru

Joined: 07/08/2011
Location: Australia
Posts: 3176
 Posted: 10:50pm 12 Dec 2011 Copy link to clipboard Print this post

To continue with another post,
I have blown the dust of a QuickBasic suntimes program I found a few years ago.
The only changes I have made are to the print routines to format the time.
When inputting location information,
latitude is negative for South and in decimal degrees.
Longitude is Positive for East and in decimal degrees.
Time zone is -11 for Eastern Summer time.

'code
5 DIM A(2), D(2)
10' Sunrise-Sunset
20 GOSUB 300
30 INPUT "Lat-deg="; B5
32 INPUT "Long-deg="; 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

320 P1 = 3.14159265
322 P2 = 2 * P1
330 DR = P1 / 180
332 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.812999999 * 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 = (0-B + D) / (2 * A)
680 IF E > 1 OR E < 0 THEN E = (0-B - D) / (2 * A)
690 T3 = C0 + E + 1 / 120 ' Round off
700 H3 = INT(T3): M3 = INT((T3 - H3) * 60)
705 Hr\$=str\$(H3):Hr\$=right\$("0"+Hr\$,2):Mn\$=str\$(M3):Mn\$=right\$(" 0"+Mn\$,2)
710 PRINT Hr\$;":";Mn\$;
720 H7 = H0 + E * (H2 - H0)
730 N7 =0 -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
782 PRINT ", azimuth ";
790 PRINT 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 - .01 * 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 = -.0001 - .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= "; Y
1183 INPUT "Month= "; M
1186 INPUT "Day= "; D

1190 G = 1: IF Y < 1583 THEN G = 0
1200 D1 = INT(D): F = D - D1 - .5
1210 J =0 -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 =0 -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.

This program appears to agree with most on-line calculators.

Jim

It all started with the ZX81....
VK7JH
http://www.c-com.com.au/MMedit.htm

Gizmo

Joined: 05/06/2004
Location: Australia
Posts: 4736
 Posted: 11:18pm 12 Dec 2011 Copy link to clipboard Print this post

Hi Jim

Works fine.

Glenn
"Treat the earth well, it was not given to you by your parents, it was lent to you by your children"

JAQ Software

DuinoMiteMegaAn
Senior Member

Joined: 17/11/2011
Location: Australia
Posts: 231
 Posted: 12:37am 13 Dec 2011 Copy link to clipboard Print this post

Comparison Results of Sunrise / Sunset Algorithms

My last posted algorithm:

Sydney, Australia Lat: -33 deg 53' South Long: -151 deg 10 East
Latitude and Longitude --> south and east ARE negative <-------

Results

Expected Results

Date of "location test" for sunrise/sunset in Sydney = December 18,2011
============================================================ =========
Internet calculated sunrise/sunset for the date = Sunrise 5:39 AM Sunset 8:04 PM - LinkB (below)
============================================================ =========

Actual Results

DECLINATION OF SUN:-22.5491 DEGREES
EQUATION OF TIME: 3.11609 MINUTES
AZIMUTH OF SUNRISE: 62.4889 DEGREES
AZIMUTH OF SUNSET: 297.511 DEGREES
Sydney, Austalia SUNRISE= 5:5 SUNSET= 19:9
Hit CR To Continue ?

=======================
TassyJim sunrise sunset
=======================

Sydney, Australia Lat: -33 deg 53' South Long: 151 deg 10 East
Posted: 13 December 2011 at 8:50am | IP Logged

latitude is negative for South and in decimal degrees.
Longitude is Positive for East and in decimal degrees. <----------

Decimal Latitude = -33.88333 Link A: (above)
Decimal Longittude = 151.16667 Link A: (above)
Time zone -11 Eastern summer time.

Actual Results

\> run

Lat-deg=? -33.88333

Long-deg=? 151.16667

Time zone-hrs=? -11

Year= ? 2011

Month= ? 12

Day= ? 18

Sunrise at 05:39, azimuth 119.16
Sunset at 20:04, azimuth 240.814

Internet calculated sunrise/sunset for the date = Sunrise 5:39 AM Sunset 8:04 PM - LinkB (above)

The only comment I have is "Spot On" TassyJim!