![]() |
Forum Index : Microcontroller and PC projects : Sun Rise and Set including Twilight program for the CMM2 Gen 1
Author | Message | ||||
RicM Regular Member ![]() Joined: 05/02/2022 Location: AustraliaPosts: 52 |
Dear Members, Dose any one have a program for calculating the Sun Rise & Set including Twilight Times anywhere in the world. Kind Regards, RicM |
||||
TassyJim![]() Guru ![]() Joined: 07/08/2011 Location: AustraliaPosts: 6283 |
This is from 10 years ago: ' Sunrise-Sunset ' adapted by TassyJim on TBS OPTION DEFAULT FLOAT DIM A(2), D(2) P1 = 3.14159265 P2 = 2 * P1 DR = P1 / 180 K1 = 15 * DR * 1.0027379 Ss$ = "Sunset at " Sr$ = "Sunrise at " M1t$ = "No sunrise this date" M2t$ = "No sunset this date" M3t$ = "Sun down all day" M4t$ = "Sun up all day" INPUT "Lat-deg="; B5 INPUT "Long-deg="; L5 INPUT "Time zone-hrs="; H L5 = L5 / 360: Z0 = -H / 24 GOSUB getJD T = (J - 2451545) + F TT = T / 36525 + 1 ' TT = centuries from 1900.0 ' LST at 0h zone time T0 = T / 36525 S = 24110.5 + 8640184.8 * T0 S = S + 86636.6 * Z0 + 86400 * L5 S = S / 86400: S = S - INT(S) T0 = S * 360 * DR T = T + Z0 ' Get Sun's Position GOSUB _Lbl910 A(1) = A5: D(1) = D5 T = T + 1 GOSUB _Lbl910 A(2) = A5: D(2) = D5 IF A(2) < A(1) THEN A(2) = A(2) + P2 Z1 = DR * 90.833 ' Zenith dist. S = SIN(B5 * DR): C = COS(B5 * DR) Z = COS(Z1): M8 = 0: W8 = 0: PRINT A0 = A(1): D0 = D(1) DA = A(2) - A(1): DD = D(2) - D(1) FOR C0 = 0 TO 23 P = (C0 + 1) / 24 A2 = A(1) + P * DA: D2 = D(1) + P * DD GOSUB _Lbl490 A0 = A2: D0 = D2: V0 = V2 NEXT C0 IF M8 = 0 AND W8 = 0 THEN IF V2 < 0 THEN PRINT M3t$ IF V2 > 0 THEN PRINT M4t$ ELSE IF M8 = 0 THEN PRINT M1t$ IF W8 = 0 THEN PRINT M2t$ ENDIF END ' ' _Lbl490: ' Test an hour for an event L0 = T0 + C0 * K1: L2 = L0 + K1 H0 = L0 - A0: H2 = L2 - A2 H1 = (H2 + H0) / 2 ' Hour angle, D1 = (D2 + D0) / 2 ' declination, ' at half hour IF C0 >= 0 THEN V0 = S * SIN(D0) + C * COS(D0) * COS(H0) - Z ENDIF V2 = S * SIN(D2) + C * COS(D2) * COS(H2) - Z IF SGN(V0) <> SGN(V2) THEN V1 = S * SIN(D1) + C * COS(D1) * COS(H1) - Z Ax = 2 * V2 - 4 * V1 + 2 * V0: B = 4 * V1 - 3 * V0 - V2 Dy = B * B - 4 * Ax * V0 IF Dy >= 0 THEN Dy = SQR(Dy) IF V0 < 0 AND V2 > 0 THEN PRINT sR$; IF V0 < 0 AND V2 > 0 THEN M8 = 1 IF V0 > 0 AND V2 < 0 THEN PRINT Ss$; IF V0 > 0 AND V2 < 0 THEN W8 = 1 E = (0-B + Dy) / (2 * Ax) IF E > 1 OR E < 0 THEN E = (0-B - Dy) / (2 * Ax) T3 = C0 + E + 1 / 120 ' Round off H3 = INT(T3): M3 = INT((T3 - H3) * 60) Hr$=STR$(H3):Hr$=RIGHT$("0"+Hr$,2):Mn$=STR$(M3):Mn$=RIGHT$(" 0"+Mn$,2) PRINT Hr$;":";Mn$; H7 = H0 + E * (H2 - H0) N7 =0 -COS(D1) * SIN(H7) D7 = C * SIN(D1) - S * COS(D1) * COS(H7) AZ = ATN(N7 / D7) / DR IF D7 < 0 THEN AZ = AZ + 180 IF AZ < 0 THEN AZ = AZ + 360 IF AZ > 360 THEN AZ = AZ - 360 PRINT ", azimuth ";AZ ENDIF ENDIF RETURN ' _Lbl910: ' Fundamental arguments ' (Van Flandern & ' Pulkkinen, 1979) L = .779072 + .00273790931 * T G = .993126 + .0027377785 * T L = L - INT(L): G = G - INT(G) L = L * P2: G = G * P2 V = .39785 * SIN(L) V = V - .01 * SIN(L - G) V = V + .00333 * SIN(L + G) V = V - .00021 * TT * SIN(L) U = 1 - .03349 * COS(G) U = U - .00014 * COS(2 * L) U = U + .00008 * COS(L) W = -.0001 - .04129 * SIN(2 * L) W = W + .03211 * SIN(G) W = W + .00104 * SIN(2 * L - G) W = W - .00035 * SIN(2 * L + G) W = W - .00008 * TT * SIN(G) ' ' Compute Sun's RA and Dec 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 = 1.00021 * SQR(U) RETURN ' getJD: ' Calendar --> JD INPUT "Year= "; Y INPUT "Month= "; M INPUT "Day= "; Dy G = 1: IF Y < 1583 THEN G = 0 D1 = INT(Dy): F = Dy - D1 - .5 J =0 -INT(7 * (INT((M + 9) / 12) + Y) / 4) IF G <> 0 THEN S = SGN(M - 9): Ax = ABS(M - 9) J3 = INT(Y + S * INT(Ax / 7)) J3 =0 -INT((INT(J3 / 100) + 1) * 3 / 4) ENDIF 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 ENDIF RETURN ' ' This program by Roger W. Sinnott calculates the times of sunrise ' and sunset on any date, accurate to the minute within several ' centuries of the present. It correctly describes what happens in the ' arctic and antarctic regions, where the Sun may not rise or set on ' a given date. Enter north latitudes positive, west longitudes ' negative. For the time zone, enter the number of hours west of ' Greenwich (e.g., -5 for EST, -4 for EDT). The calculation is ' discussed in Sky & Telescope for August 1994, page 84. If I remember correctly, you change line 37 to suit the standard required. civil, military etc as well as start of twilight. Z1 = DR * 90.833 ' Zenith dist. Search the forum for "sunrise". I and others have made numerous posts over the years. Jim Edited 2023-02-22 11:26 by TassyJim VK7JH MMedit |
||||
Geoffg![]() Guru ![]() Joined: 06/06/2011 Location: AustraliaPosts: 3292 |
This is a great algorithm that I have successfully used in the past: Geoff Sunrise/Sunset Algorithm Source: http://williams.best.vwh.net From: Almanac for Computers, 1990 published by Nautical Almanac Office United States Naval Observatory Washington, DC 20392 NOTE: LONGTITUDE USES + FOR WEST. THIS IS THE OPPOSITE OF THE GPS MODULE VB Function: Function sun(dateinput As Date, lat, Lon, zenith, rise As Boolean) 'return the UTC of sunrise or sunset in hours 'return -1 if sun never rises, -2 if it never sets 'lat and lon are the lat/lon of the location in degrees 'zenith is the zenith in degrees '=90.8333333 for official sunset '=96 for civil twilight '=102 for nautical twilight '=108 for astronomical twilight 'rise =True for sunrise, False for sunset N1 = Int(275 * Month(dateinput) / 9) n2 = Int((Month(dateinput) + 9) / 12) n3 = (1 + Int((Year(dateinput) - 4 * Int(Year(dateinput) / 4) + 2) / 3)) nd = N1 - (n2 * n3) + Day(dateinput) - 30 lnghour = Lon / 15# If (rise) Then t = nd + ((6 - lnghour) / 24) Else t = nd + ((18 - lnghour) / 24) End If m = (0.9856 * t) - 3.289 l = md(m + 1.916 * sIn(degtorad(m)) + 0.02 * sIn(degtorad(2 * m)) + 282.634, 360) ra = radtodeg(md(Atn(0.91764 * Tan(degtorad(l))), 2 * Pi)) lquad = Int(l / 90) * 90 raquad = Int(ra / 90) * 90 ra = ra + lquad - raquad ra = ra / 15 sindec = 0.39782 * sIn(degtorad(l)) cosdec = Cos(Asn(sindec)) coshr = (Cos(degtorad(zenith)) - sindec * sIn(degtorad(lat))) / (cosdec * (Cos(degtorad(lat)) + 1E-20)) If (coshr > 1) Then sun = -1 ElseIf (coshr < -1) Then sun = -2 Else If (rise) Then h = 360 - radtodeg(Acs(coshr)) Else h = radtodeg(Acs(coshr)) End If h = h / 15 t = h + ra - 0.06571 * t - 6.622 sun = md(t - lnghour, 24) End If End Function Explanation: day, month, year: date of sunrise/sunset latitude, longitude: location for sunrise/sunset zenith: Sun's zenith for sunrise/sunset offical = 90 degrees 50' THIS IS THE CLOSEST TO THE STANDARD SUNRISE/SET civil = 96 degrees nautical = 102 degrees astronomical = 108 degrees NOTE: longitude is positive for East and negative for West NOTE: the algorithm assumes the use of a calculator with the trig functions in "degree" (rather than "radian") mode. Most programming languages assume radian arguments, requiring back and forth convertions. The factor is 180/pi. So, for instance, the equation RA = atan(0.91764 * tan(L)) would be coded as RA = (180/pi)*atan(0.91764 * tan((pi/180)*L)) to give a degree answer with a degree input for L. 1. first calculate the day of the year N1 = floor(275 * month / 9) N2 = floor((month + 9) / 12) N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3)) N = N1 - (N2 * N3) + day - 30 2. convert the longitude to hour value and calculate an approximate time lngHour = longitude / 15 if rising time is desired: t = N + ((6 - lngHour) / 24) if setting time is desired: t = N + ((18 - lngHour) / 24) 3. calculate the Sun's mean anomaly M = (0.9856 * t) - 3.289 4. calculate the Sun's true longitude L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634 NOTE: L potentially needs to be adjusted into the range [0,360) by adding/subtracting 360 5a. calculate the Sun's right ascension RA = atan(0.91764 * tan(L)) NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360 5b. right ascension value needs to be in the same quadrant as L Lquadrant = (floor( L/90)) * 90 RAquadrant = (floor(RA/90)) * 90 RA = RA + (Lquadrant - RAquadrant) 5c. right ascension value needs to be converted into hours RA = RA / 15 6. calculate the Sun's declination sinDec = 0.39782 * sin(L) cosDec = cos(asin(sinDec)) 7a. calculate the Sun's local hour angle cosH = (cos(zenith) - (sinDec * sin(latitude))) / (cosDec * cos(latitude)) if (cosH > 1) the sun never rises on this location (on the specified date) if (cosH < -1) the sun never sets on this location (on the specified date) 7b. finish calculating H and convert into hours if if rising time is desired: H = 360 - acos(cosH) if setting time is desired: H = acos(cosH) H = H / 15 8. calculate local mean time of rising/setting T = H + RA - (0.06571 * t) - 6.622 9. adjust back to UTC UT = T - lngHour NOTE: UT potentially needs to be adjusted into the range [0,24) by adding/subtracting 24 10. convert UT value to local time zone of latitude/longitude localT = UT + localOffset Geoff Graham - http://geoffg.net |
||||
Turbo46![]() Guru ![]() Joined: 24/12/2017 Location: AustraliaPosts: 1642 |
You don't say what you want if it for but Captainboing did a 'quick and dirty one' on FOTS that may suit? Bill Keep safe. Live long and prosper. |
||||
cdeagle Senior Member ![]() Joined: 22/06/2014 Location: United StatesPosts: 265 |
Here's a not so quick and dirty version. sun_riseset.zip |
||||
Geoffg![]() Guru ![]() Joined: 06/06/2011 Location: AustraliaPosts: 3292 |
You are right... 47KB. Wow! It is also big! Geoff Geoff Graham - http://geoffg.net |
||||
![]() |
![]() |
The Back Shed's forum code is written, and hosted, in Australia. | © JAQ Software 2025 |