Precision calculation (a few minutes of time) of sunrise and sunset

Mesa · Post by **Mesa** » Tue Sep 24, 2024 10:45 am

FYI:  it's a bug free of my own version of this code somewhere in the forum.
Comments in french and english.


;COMPUTES SUNRISE SUNSET TIMES

;+========================================================+
;|                                                        |
;| Calcul de l'heure TU du lever et du coucher du Soleil  |
;|                                                        |
;| Précision : La minute de temps                         |
;| Plateforme: Multiplateforme                            |
;| ToDo : Reste à optimiser le code (réduction du code)   |
;| Auteur : Mesa 30 Janvier 2013                          |
;|                                                        |
;+========================================================+


;Grenoble +45° 11' 16" N (Nord = positif, Sud = négatif)
; -05° 43' 37" E (Est = négatif, Ouest = Positif)
;Altitude Min. 204 m Max. 475 m


; Bureau des longitudes : heure de lever du soleil
; http://www.imcce.fr/fr/ephemerides/phenomenes/rts/rts.php

; Soleil le 01/01/2013
; Lieu : 05°43’37" E / 45°11’16" N
; Date Lever Passage au méridien Coucher
; Temps Universel heure azimut heure hauteur heure azimut
; 2013-01-01 07:18 302.90 11:41 21.85 16:04 57.14
; Date Aube Crépuscule
; astronom. nautique civil civil nautique astronom.
; Temps Universel heure heure heure heure heure heure
; 2013-01-01 05:29 06:05 06:42 16:39 17:17 17:52







;ENGLISH
;+============================================================+
;|                                                            |
;| Calculation of the UT time of sunrise and sunset           |
;|                                                            |
;| Precision: The minute of time                              |
;| Platform: Multiplatform                                    |
;| ToDo: Remains to optimize the code (reduction of the code) |
;| Author: Mesa January 30, 2013                              |
;|                                                            |
;+============================================================+

;Grenoble +45° 11' 16" N (North = positive, South = negative)
; -05° 43' 37" E (East = negative, West = Positive)
;Altitude Min. 204 m Max. 475 m

; Bureau des longitudes: sunrise time
; http://www.imcce.fr/fr/ephemerides/phenomenes/rts/rts.php

; Sun on 01/01/2013
; Location: 05°43’37" E / 45°11’16" N
; Date Sunrise Meridian crossing Sunset
; Universal Time time azimuth time altitude time azimuth
; 2013-01-01 07:18 302.90 11:41 21.85 16:04 57.14
; Date Dawn Dusk
; astronom. nautical civil civil nautical astronom.
; Universal Time time hour hour hour hour hour
; 2013-01-01 05:29 06:05 06:42 16:39 17:17 17:52


EnableExplicit

Procedure.d Frac(Xfrac.d)
  ProcedureReturn Xfrac - Int(Xfrac)
EndProcedure

Procedure.d dms(Xdms.d)
  ;conversion decimales en degres minutes secondes
  ;conversion decimals to degrees minutes seconds
  ProcedureReturn Int(Xdms) + Frac(Int(Frac(Xdms) * 60) / 100) + Frac(Frac(Frac(Xdms) * 60) * 0.006)
EndProcedure

Procedure.s dmsS(Xdms.d)
  ;conversion decimales en degres minutes secondes
  ;conversion decimals to degrees minutes seconds
  Protected h$ = Str(Int(Xdms)) + "h "
  Protected m.d = 100*Frac(Xdms)
  Protected m$ = Str(Int(m)) + "m "
  Protected s.d = 100*Frac(m)
  Protected s$ = Str(Int(Round(s, #PB_Round_Nearest ))) + "s"
  
  ProcedureReturn h$ + m$ + s$
EndProcedure

Procedure.d dec(xdec.d)
  ;conversion degres minutes secondes en decimales
  ;conversion degrees minutes seconds to decimals
  ProcedureReturn Int(Xdec) + Int(frac(Xdec) * 100) / 60 + (frac(Xdec) * 100 - Int(frac(Xdec) * 100)) * 100 / 3600
EndProcedure

Procedure.d jd2(Annee, Mois, jour, Heure, Minute, Seconde)
  ; Jour Julien à 0h TU (Temps Universel )
  ; --Ne fonctionne pas avant le 15 Octobre 1582
  
  ; Julian Day at 0h UT (Universal Time)
  ; --Does not work before October 15, 1582
  
  Protected j.d
  Protected m
  Protected a
  Protected jj.d
  
  a = Annee
  m = Mois
  j = jour + (Heure + Minute / 60 + Seconde / 3600) / 24
  
  jj = 367*a - Int(1.75 * (Int((m + 9) / 12) + a))
  jj - Int(0.75 * (1 + Int((Int((m - 9) / 7) + a) * 0.01)))
  jj + Int((275*m) / 9) + j + 1721028.5
  
  ProcedureReturn jj
EndProcedure
Procedure.d LeverCoucherSoleil(LongitudeT.d, LatitudeT.d, Jour, Mois, An, Lever = #True, Option.f =  - 0.61)
  
  ; Option
  ; -18<= Option <-12 : Crépuscule astronomique
  ; -12<= Option <-6 : Crépuscule nautique
  ; -6 <= Option <réfraction : Crépuscule civil
  ;
  ; réfraction = -36.6'= -0.61° : Ephémérides du Bureau des longitudes
  ; = -34' : Ephéméride nautique
  ; = centre du Soleil : Ephémérides astronomiques
  ; = bord inférieur du Soleil : Ephémérides pour navigateurs
  ;
  
  ; Lever
  ; Lever=#True : Calcul du lever du Soleil
  ; Lever=#False : Calcul du coucher du Soleil
  
  ; Vérifiez que les longituse EST soient négatives et OUEST positive
  
  
  ;=====================================================================
  ; Option
  ; -18<= Option <-12 : Astronomical twilight
  ; -12<= Option <-6 : Nautical twilight
  ; -6 <= Option <refraction : Civil twilight
  ;
  ; refraction = -36.6'= -0.61° : Ephemeris of the Bureau of Longitudes
  ; = -34' : Nautical ephemeris
  ; = center of the Sun : Astronomical ephemeris
  ; = lower edge of the Sun : Ephemeris for navigators
  ;
  
  ; Sunrise
  ; Sunrise=#True : Calculation of sunrise
  ; Sunrise=#False : Calculation of sunset
  
  ; Check that the EAST longituse are negative and WEST positive
  
  Protected.d FA, H0, FI, LO, JD, T0, S0, TU, TS, T
  Protected.d LL, M1, C, TL, EX, AL, X, DE, AH, TM
  Protected IN0, i
  
  FA = 180 / #PI
  H0 = 0
  
  ; Radian
  FI = Radian(dec(LatitudeT))
  LO = Radian(dec(LongitudeT))
  
  
  ; Jour julien de la date à 00h:00'00" (ne pas changer l'heure)
  ; Julian day of the date at 00h:00'00" (do not change the time)
  JD = jd2(An, Mois, Jour, 0, 0, 0)
  ;Debug JD
  
  
  ;Temps Sideral à Greenwitch ;Sidereal Time in Greenwitch
  T0 = (JD - 2415020.0) / 36525
  S0 = 6.6460656 + 2400.051262*T0 + 0.00002591*T0*T0
  S0 = S0 - 24*Int(S0 / 24)
  TU = 12
  TS = S0 + TU * 1.002737908
  ;Debug TS
  
  
  ;Position de la Terre sur son orbite (ou soleil);Position of the Earth in its orbit (or sun)
  T  = T0 + TU / 24 / 36525
  LL = Radian(279.69668 + T * 36000.76892 + T * T * 0.0003025)
  M1 = Radian(358.47583 + T * 35999.04975 - T * T * 0.000151 - T * T * T * 0.0000033)
  C  = Radian((1.91946 - 0.004789*T - T * T * 0.000014) * Sin(M1) + (0.020094 - T * 0.0001) * Sin(2*M1) + 0.000293*Sin(3*M1))
  
  
  ;Longitue ecliptique ;Ecliptic longitude
  TL = LL + C
  ;Coordonnées écliptique ;Ecliptic coordinates
  EX = Radian(23.452294 - 0.0130125*T - 0.00000164*T * T + T * T * T * 0.000000503)
  
  AL = ATan(Cos(EX) * Tan(TL))
  If Cos(TL) < 0
    AL + #PI
  EndIf
  X  = Sin(EX) * Sin(TL)
  DE = ATan(X / Sqr( - X * X + 1))
  ;Angle Horaire ;hour angle
  X = (Sin(H0) - Sin(FI) * Sin(DE)) / Cos(FI) / Cos(DE)
  If Abs(X) > 1
    MessageRequester("impossible", StrD(Degree(H0)))
    End
  EndIf
  AH =  - ATan(X / Sqr( - X * X + 1)) + #PI / 2
  
  TM = TS - Degree((LO + AL)) / 15
  TM = 36 - TM
  TM = TM - 24*Int(TM / 24)
  
  ;Lever ;sunrise
  If Lever = #True
    IN0 =  - 1
  Else
    IN0 = 1
  EndIf
  
  H0 = Radian(Option)
  
  For i = 1 To 3
    ;Position de la Terre sur son orbite (ou soleil) ;Position of the Earth in its orbit (or sun)
    T  = T0 + TU / 24 / 36525
    LL = Radian(279.69668 + T * 36000.76892 + T * T * 0.0003025)
    M1 = Radian(358.47583 + T * 35999.04975 - T * T * 0.000151 - T * T * T * 0.0000033)
    C  = Radian((1.91946 - 0.004789*T - T * T * 0.000014) * Sin(M1) + (0.020094 - T * 0.0001) * Sin(2*M1) + 0.000293*Sin(3*M1))
    
    ;Longitue ecliptique
    TL = LL + C
    ;Coordonnées écliptique
    EX = Radian(23.452294 - 0.0130125*T - 0.00000164*T * T + T * T * T * 0.000000503)
    
    AL = ATan(Cos(ex) * Tan(TL))
    If Cos(TL) < 0
      AL + #PI
    EndIf
    X  = Sin(EX) * Sin(TL)
    DE = ATan(X / Sqr( - X * X + 1))
    ;Angle Horaire ;hour angle
    X = (Sin(H0) - Sin(FI) * Sin(DE)) / Cos(FI) / Cos(DE)
    If Abs(X) > 1
      MessageRequester("impossible", StrD(Degree(H0)))
      End
    EndIf
    AH =  - ATan(X / Sqr( - X * X + 1)) + #PI / 2
    
    TS = Degree((IN0*AH + AL + LO)) / 15
    While TS < S0
      TS = TS + 24
    Wend
    TU = (TS - S0) / 1.002737908
    TU = TU - 24*Int(TU / 24)
  Next i
  ProcedureReturn TU
EndProcedure

Procedure isLeapYear(Year)
  If (Year % 4 = 0 And Year % 100) Or Year % 400 = 0
    ProcedureReturn #True
  Else
    ProcedureReturn #False
  EndIf
EndProcedure

Procedure SunSetSunRise_By_Noaa(Y, M, D, longitude, lat)
  ;https://www.esrl.noaa.gov/gmd/grad/solcalc/solareqns.PDF
  ;====> very poor accuracy but quicker
  
  If Y < 1970 Or Y > 2037:ProcedureReturn - 1:EndIf
  
  Protected y0 = 365, G.d, eqtime.d, decl.d, time_offset.d, tst.d, ha.d, phi.d, theta.d
  Protected sunset.d, sunrise.d, snoon.d
  Protected H.d = 0.0, Hour, mn, sc
  Protected Cos_phi.d, Cos_180moinstheta.d
  
  If isLeapYear(Y):y0 + 1:EndIf
  
  ;the fractional Year (G) is calculated, in radians.	
  G = 2 * #PI / y0 * (DayOfYear(Date(Y, M, D, 0, 0, 0)) - 1 + (H - 12) / 24)
  
  ;equation of time (in minutes) And the solar declination angle (in radians).
  eqtime = 229.18 * (0.000075 + 0.001868*Cos(G) - 0.032077*Sin(G) - 0.014615*Cos(2*G) - 0.040849*Sin(2*G))
  decl   = 0.006918 - 0.399912*Cos(G) + 0.070257*Sin(G) - 0.006758*Cos(2*G) + 0.000907*Sin(2*G) - 0.002697*Cos(3*G) + 0.00148*Sin(3*G)
  
  ; ;the true solar time is calculated. First the time offset is found, in minutes, And then the true solar time, in minutes.
  ; time_offset = eqtime + 4*longitude - 60*timezone
  ; ; where eqtime is in minutes, longitude is in degrees (positive To the east of the Prime Meridian), timezone is in hours from UTC (U.S. Mountain Standard Time = -7 hours).
  ; tst = hr*60 + mn + sc/60 + time_offset
  ; ; where hr is the Hour (0 - 23), mn is the Minute (0 - 59), sc is the Second (0 - 59).
  ; ; The solar hour angle, in degrees, is:
  ; ha = (tst / 4) - 180
  ; ; The solar zenith angle (phi) can then be found from the hour angle (ha), latitude (lat) And solar declination (decl) using the following equation:
  ; Cos_phi = Sin(Radian(lat))*Sin(decl) + Cos(Radian(lat))*Cos(decl)*Cos(Radian(ha))
  ; ; And the solar azimuth (theta, degrees clockwise from north) is found from:
  ; Cos_180moinstheta= -((Sin(lat)*Cos(phi)-Sin(decl))/Cos(lat)*Sin(phi))
  
  
  
  ; Sunrise/Sunset Calculations
  ;.............................
  ; For the special Case of sunrise Or sunset, the zenith is set To 90.833 (the approximate correction For atmospheric refraction at sunrise And sunset, And the size of the solar disk), And the hour angle becomes:
  ha = ACos(Cos (Radian(90.833)) / (Cos(Radian(lat)) * Cos(decl)) - Tan(Radian(lat)) * Tan(decl))
  
  ; where the positive number corresponds To sunrise, negative To sunset.
  ; Then the UTC time of sunrise (Or sunset) in minutes is:
  ;BEWARE longitude in degrees, positive to the east
  sunrise = 720 - 4 * (-longitude + Degree(ha)) - eqtime
  Protected ah.d =  - ha
  
  sunset = 720 - 4 * (-longitude + Degree(ah)) - eqtime
  ; where longitude And hour angle are in degrees And the equation of time is in minutes.
  
  ; Solar noon For a given location is found from the longitude (in degrees, positive To the east of the Prime Meridian) And the equation of time (in minutes):
  snoon = 720 + 4*longitude - eqtime
  
  
  Macro formatme(x)
    If x < 0:x + 360:EndIf
    If x >= 360:x = x - 360*Int(x / 360):EndIf
  EndMacro
  
  
  formatme(sunset)
  Debug "sunset " + StrD(sunset / 15) + " " + dmsS(dms((sunset / 15)))
  formatme(sunrise)
  Debug "sunriset " + StrD(sunrise / 15) + " "  + dmsS(dms(sunrise / 15))
  formatme(snoon)
  Debug "noon " + StrD(snoon / 15) + " " + dmsS(dms((snoon / 15)))
EndProcedure


;Affichage
; 	Debug "Le 01/01/2013 à Grenoble"
; 	Debug ""
; 	Debug "Lever/Coucher Civil: "
; 	Debug dms(LeverCoucherSoleil(-5.4337,45.1116,1,1,2013))          ;07h17'37" ~ 07h 18'
; 	Debug dms(LeverCoucherSoleil(-5.4337,45.1116,1,1,2013,#False))
; 	Debug ""
; 	Debug "Lever/CoucherNautique: "
; 	Debug dms(LeverCoucherSoleil(-5.4337,45.1116,1,1,2013,#True,-12))
; 	Debug dms(LeverCoucherSoleil(-5.4337,45.1116,1,1,2013,#False,-12))
; 	Debug ""
; 	Debug "Lever/CoucherAstronomique (Nuit noire): "
; 	Debug dms(LeverCoucherSoleil(-5.4337,45.1116,1,1,2013,#True,-18))
; 	Debug dms(LeverCoucherSoleil(-5.4337,45.1116,1,1,2013,#False,-18))
;   Debug ""

Debug "The 12 august 2020 in Paris, FRANCE"
;année bissextile ;leap year
;https://promenade.imcce.fr/fr/pages5/585.html


Define Long.d = -2.20 ;02°20' E ;négatif à l'EST
Define Lat.d  = 48.50 ;48°50' N

Define D = 12
Define M = 8
Define Y = 2020

Debug ""
Debug "Civilian Sunrise / Sunset: "
; 	Debug dms(LeverCoucherSoleil(Long,Lat,D,M,Y))          ;Result 4h 41m 11s UTC, sunrise should be 04h 38' UTC
; 	Debug dms(LeverCoucherSoleil(Long,Lat,D,M,Y,#False))   ;Result 19h 9m 11s UTC, sunset should be 19h 09' UTC
Debug dmsS(dms(LeverCoucherSoleil(Long, Lat, D, M, Y))) + " UTC"
Debug dmsS(dms(LeverCoucherSoleil(Long, Lat, D, M, Y, #False))) + " UTC"

Debug ""
Debug "Nautical Sunrise / Sunset: "
Debug dmsS(dms(LeverCoucherSoleil(Long, Lat, D, M, Y, #True,  - 12))) + " UTC"
Debug dmsS(dms(LeverCoucherSoleil(Long, Lat, D, M, Y, #False,  - 12))) + " UTC"
Debug ""
Debug "Astronomical Sunrise / Sunset (Dark Night): "
Debug dmsS(dms(LeverCoucherSoleil(Long, Lat, D, M, Y, #True,  - 18))) + " UTC"
Debug dmsS(dms(LeverCoucherSoleil(Long, Lat, D, M, Y, #False,  - 18))) + " UTC"







Debug ""
Debug "******************************************************"
Debug "SunSetSunRise_By_Noaa (very poor accuracy)"
Debug "*********************************************"
SunSetSunRise_By_Noaa(Y, M, D, long, lat)

jacdelad · Post by **jacdelad** » Tue Sep 24, 2024 2:04 pm

Hi Mesa,
interesting code. I would change it into a library, but that's a matter of taste.
Serious question: Longitude and latitude are given as float. Your example:

Code: Select all

Define Long.d = -2.20 ;02°20' E ;négatif à l'EST
Define Lat.d  = 48.50 ;48°50' N

If I'm not in error, 50' <> 0.5°, so, is your calculation right??? Sorry, if I'm wrong.

threedslider · Post by **threedslider** » Tue Sep 24, 2024 5:06 pm

Nice ! Very technical and they work well too

Mesa · Post by **Mesa** » Wed Sep 25, 2024 9:22 am

@jacdelad : Thanks. Inside the code, a procedure change longitudes and lat in deg min sec into decimal degrees, so they are no mistakes (because it's easier (for me) to use long and lat in deg min sec

).

M.

jacdelad · Post by **jacdelad** » Wed Sep 25, 2024 12:32 pm

Ah ok. I still find it a bit confusing, entering a decimal should be xx.yy degress. I would prefer to schemes: xx.yy° and xx°yy'zz", this seems more logical to me. But I can change this myself, if I need to.

Thanks for sharing.

PureBasic Forums - English

Precision calculation (a few minutes of time) of sunrise and sunset

Precision calculation (a few minutes of time) of sunrise and sunset

Re: Precision calculation (a few minutes of time) of sunrise and sunset

Re: Precision calculation (a few minutes of time) of sunrise and sunset

Re: Precision calculation (a few minutes of time) of sunrise and sunset

Re: Precision calculation (a few minutes of time) of sunrise and sunset