Home
JAQForum Ver 24.01
Log In or Join  
Active Topics
Local Time 06:50 02 Aug 2025 Privacy Policy
Jump to

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 : Sun Rise and Set including Twilight program for the CMM2 Gen 1

Author Message
RicM
Regular Member

Joined: 05/02/2022
Location: Australia
Posts: 52
Posted: 01:12am 22 Feb 2023
Copy link to clipboard 
Print this post

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: Australia
Posts: 6283
Posted: 01:23am 22 Feb 2023
Copy link to clipboard 
Print this post

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: Australia
Posts: 3292
Posted: 03:11am 22 Feb 2023
Copy link to clipboard 
Print this post

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: Australia
Posts: 1642
Posted: 07:26am 22 Feb 2023
Copy link to clipboard 
Print this post

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 States
Posts: 265
Posted: 09:32am 22 Feb 2023
Copy link to clipboard 
Print this post

Here's a not so quick and dirty version.


sun_riseset.zip
 
Geoffg

Guru

Joined: 06/06/2011
Location: Australia
Posts: 3292
Posted: 11:45am 22 Feb 2023
Copy link to clipboard 
Print this post

  cdeagle said  Here's a not so quick and dirty version.

You are right... 47KB.  Wow!  It is also big!

Geoff
Geoff Graham - http://geoffg.net
 
Print this page


To reply to this topic, you need to log in.

The Back Shed's forum code is written, and hosted, in Australia.
© JAQ Software 2025