Calculation of Sunrise & Sunset time

Share your advanced PureBasic knowledge/code with the community.
Post Reply
User avatar
jacdelad
Addict
Addict
Posts: 2014
Joined: Wed Feb 03, 2021 12:46 pm
Location: Riesa

Re: Calculation of Sunrise & Sunset time

Post by jacdelad »

Hi Mesa,
thanks very much. Unfortunately my french stops after "Baguette", so what exactly is the option-parameter for? And what does the output mean (like "7.1" means "7.1 hours" -> "07:06")?

Also, I think it's better to use the pb-internal date-format in one parameter instead of three parameters:

Code: Select all

;+==========================================================+
;|                                                          |
;| Calcul de l'heure TU du lever et du coucher du Soleil    |
;|                                                          |
;| Précision : La minute de temps                           |
;| Plateforme: Multiplateforme                              |
;| ToDo      :                                              |
;| Auteur    : Mesa  30 Janvier 2013                        |
;|                                                          |
;+==========================================================+



;Ville
;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


EnableExplicit

Procedure.d Frac(Xfrac.d)   
  ProcedureReturn Xfrac-Int(Xfrac)    
EndProcedure   
Procedure.d dms(Xdms.d)   
  ;conversion decimales en degres minutes secondes  
  ProcedureReturn Int(Xdms)+Frac(Int(Frac(Xdms)*60)/100)+Frac(Frac(Frac(Xdms)*60)*0.006)   
EndProcedure   
Procedure.d dec(xdec.d)   
  ;conversion degres minutes secondes en decimales  
  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 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, InputDate, Lever=#True, Option.f=-0.61)
  
  Protected Jour=Day(InputDate), Mois=Month(InputDate), An=Year(InputDate)
  ; 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
  

  
  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=dec(LatitudeT)/FA
  LO=dec(LongitudeT)/FA
  
  
  ; Jour julien de la date à 00h:00'00" (ne pas changer l'heure)
  JD=jd2(An,Mois,Jour,0,0,0)
  ;Debug JD
  
  
  ;Temps Sideral à 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)
  T= T0+TU/24/36525
  LL=279.69668+T*36000.76892+T*T*0.0003025
  LL / FA
  M1 = 358.47583+T*35999.04975-T*T*0.000151-T*T*T*0.0000033
  M1 / FA
  C = (1.91946-0.004789*T-T*T*0.000014)*Sin(M1)+(0.020094-T*0.0001)*Sin(2*M1)+0.000293*Sin(3*M1)
  C / FA
  ;Longitue ecliptique
  TL = LL +C
  ;Coordonnées écliptique
  EX = 23.452294-0.0130125*T-0.00000164*T*T+T*T*T*0.000000503
  EX / FA
  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
  X=(Sin(H0)-Sin(FI)*Sin(DE))/Cos(FI)/Cos(DE)
  If Abs(X)>1 
    MessageRequester("impossible",StrD(H0*FA))
    End
  EndIf
  AH=-ATan(X/Sqr(-X*X+1))+ #PI/2
    
  TM=TS-(LO+AL)*FA/15
  TM=36-TM
  TM=TM-24*Int(TM/24)
  
  ;Lever
  If Lever=#True
    IN0=-1
  Else
    IN0=1
  EndIf
  
  H0=Option/FA
  
  For i= 1 To 3
    ;Position de la Terre sur son orbite (ou soleil)
    T= T0+TU/24/36525
    LL=279.69668+T*36000.76892+T*T*0.0003025
    LL / FA
    M1 = 358.47583+T*35999.04975-T*T*0.000151-T*T*T*0.0000033
    M1 / FA
    C = (1.91946-0.004789*T-T*T*0.000014)*Sin(M1)+(0.020094-T*0.0001)*Sin(2*M1)+0.000293*Sin(3*M1)
    C / FA
    ;Longitue ecliptique
    TL = LL +C
    ;Coordonnées écliptique
    EX = 23.452294-0.0130125*T-0.00000164*T*T+T*T*T*0.000000503
    EX / FA
    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
    X=(Sin(H0)-Sin(FI)*Sin(DE))/Cos(FI)/Cos(DE)
    If Abs(X)>1 
      MessageRequester("impossible",StrD(H0*FA))
      End
    EndIf
    AH=-ATan(X/Sqr(-X*X+1))+ #PI/2
     
    TS=(IN0*AH+AL+LO)*FA/15
    While TS<S0
      TS=TS+24
    Wend
    TU=(TS-S0)/1.002737908
    TU=TU-24*Int(TU/24)
  Next i
  ProcedureReturn TU
EndProcedure


;Affichage
Debug "Le 01/01/2013 à Grenoble"
Debug ""
Debug "Lever/Coucher Civil: "
Debug dms(LeverCoucherSoleil(-5.4337,45.1116,ParseDate("%dd.%mm.%yyyy","01.01.2013")))          ;07h17'37" ~ 07h 18'
Debug dms(LeverCoucherSoleil(-5.4337,45.1116,ParseDate("%dd.%mm.%yyyy","01.01.2013"),#False))
Debug ""
Debug "Lever/CoucherNautique: "
Debug dms(LeverCoucherSoleil(-5.4337,45.1116,ParseDate("%dd.%mm.%yyyy","01.01.2013"),#True,-12))
Debug dms(LeverCoucherSoleil(-5.4337,45.1116,ParseDate("%dd.%mm.%yyyy","01.01.2013"),#False,-12))
Debug ""
Debug "Lever/CoucherAstronomique (Nuit noire): "
Debug dms(LeverCoucherSoleil(-5.4337,45.1116,ParseDate("%dd.%mm.%yyyy","01.01.2013"),#True,-18))
Debug dms(LeverCoucherSoleil(-5.4337,45.1116,ParseDate("%dd.%mm.%yyyy","01.01.2013"),#False,-18))
Good morning, that's a nice tnetennba!

PureBasic 6.21/Windows 11 x64/Ryzen 7900X/32GB RAM/3TB SSD
Synology DS1821+/DX517, 130.9TB+50.8TB+2TB SSD
Top
Mesa
Enthusiast
Enthusiast
Posts: 433
Joined: Fri Feb 24, 2012 10:19 am

Re: Calculation of Sunrise & Sunset time

Post by Mesa »

If you don't use the option then the calculation is good for everyone.
The option is only used for aviators, navigators (sailors) and astronomers.

The output is in decimal hours (7.5h) but in the examples I gave, I use my dms() function which transforms decimal hours into hours, minutes, seconds and fraction of seconds: 17.5h becomes 17h 30' 00".000...

; Lever (sunrise)
; Lever=#True: Sunrise calculation
; Lever=#False: Sunset calculation





Read under only if you are an aviator, navigator (sailor) or an astronomer.
; Option (in degree)
; if -18°<= Option <-12°: Astronomical twilight. The night is really dark.
; if -12°<= Option <-6°: Nautical Twilight.The sun has set but it is still daylight, it is twilight.The sun hasn't risen yet but it's already daytime, it's dawn. Useful for amateur airplane pilots.
; if -6' <= Option <refraction: Civil Twilight (used by every one)
;
; Atmospheric refraction: The effect of the atmosphere on the time of sunrise and sunset .
; refraction = -36.6'= -0.61°: Default for civil twilight for the Ephemerides of the Bureau des longitudes (French national institut that make ephemeris)
; = -34': For Nautical ephemeris (not used in this algorithm)
; Calculation relative to the center of the sun: Astronomical ephemerides (Used by default in this algorithm)
; Calculation relative to the lower edge of the Sun: Ephemerides for navigators (Not used in this algorithm)
; (Sunrises and sunsets times are different for sailors !)


Code: Select all

Ex: Directly in h min sec with the function dms(...)
Debug "September 8th 2022 at Paris, the Eiffel tower: 48.86509418723488N, 2.292439520977688E"
Debug ""
Debug "Surise: "
Debug dms(LeverCoucherSoleil(-2.292439520977688, 48.86509418723488, ParseDate("%dd.%mm.%yyyy","08.09.2022")))  ;5.1716901104237678 ~ 05h 17'
Debug "Sunset: " 
Debug dms(LeverCoucherSoleil(-2.292439520977688, 48.86509418723488, ParseDate("%dd.%mm.%yyyy","08.09.2022"),#False));18.171414055916191 ~ 18h 17'
Debug ""
M.
Top
Post Reply

Return to “Tricks 'n' Tips”