PureBasic Forums - English

http://www.purebasic.com
https://www.purebasic.fr/english/

Sunset/sunrise calculation

https://www.purebasic.fr/english/viewtopic.php?t=75824

Page 1 of 1

Sunset/sunrise calculation

Posted: Tue Aug 11, 2020 5:20 pm
by l1marik
Has somebody here written procedure for sunset/sunrise (and the same for moon) calculation?

I would like to apply it to my gauges, of course that I can write my own, but if somebody has it, it should be better :-)

Image

Thank you.
Lukas

Re: Sunset/sunrise calculation

Posted: Tue Aug 11, 2020 6:08 pm
by NicTheQuick
It's a bit complicated but I found this formula for you: https://www.esrl.noaa.gov/gmd/grad/solc ... areqns.PDF

Re: Sunset/sunrise calculation

Posted: Tue Aug 11, 2020 7:23 pm
by infratec
If you simply enter sunrise in the search field, you will find:

viewtopic.php?f=12&t=64652&hilit=sunrise

Re: Sunset/sunrise calculation

Posted: Wed Aug 12, 2020 5:38 pm
by Mesa

Code: Select all

  
        ;+========================================================+ 
	;|                                                                                                      | 
	;| 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 
	 
	 
; 	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, 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 
		 
	 
		 
		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,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 "The 12 august 2020 in Paris, FRANCE" 
	;https://promenade.imcce.fr/fr/pages5/585.html
	
	Long.d =-2.33333333333 ;02°20' E 
	Lat.d  =48.83333333333 ;48°50' N
	
	D=12
	M=8
	Y=2020
	
	Debug "" 
	Debug "Civilian Sunrise / Sunset: " 
	Debug dms(LeverCoucherSoleil(Long,Lat,D,M,Y))          ;sunrise should be 04h38' UTC
	Debug dms(LeverCoucherSoleil(Long,Lat,D,M,Y,#False))   ;sunset should be 19h09' UTC
	Debug "" 
	Debug "Nautical Sunrise / Sunset: " 
	Debug dms(LeverCoucherSoleil(Long,Lat,D,M,Y,#True,-12)) 
	Debug dms(LeverCoucherSoleil(Long,Lat,D,M,Y,#False,-12)) 
	Debug "" 
	Debug "Astronomical Sunrise / Sunset (Dark Night): " 
	Debug dms(LeverCoucherSoleil(Long,Lat,D,M,Y,#True,-18)) 
	Debug dms(LeverCoucherSoleil(Long,Lat,D,M,Y,#False,-18))
Mesa.

Re: Sunset/sunrise calculation

Posted: Mon Dec 14, 2020 5:10 pm
by CalamityJames
Davido sent me a message ages ago, drawing my attention to this post, but I have only just seen it! Coincidentally I have just written a post about moon phases (here:viewtopic.php?f=13&t=61671 You will find a link which points to some Applesoft Basic programs which include Sunset and Sunrise. They are easy to translate into Pure Basic (I show an example in the post).

Re: Sunset/sunrise calculation

Posted: Tue Dec 15, 2020 5:07 pm
by ebs
Here is a simple routine that I found (and modified slightly). It seems fairly accurate:

Code: Select all

Procedure.s SunCalc (lat.f, lon.f, T.i, tz.d)
  Define.d B, D, zeitdiff, zeitglg, AMOZ, UMOZ, sunrise, pi = 3.1415927, h = -0.0145, sunset
  ; lat = latitude (N/S), lon=longitude (E/W), T=Day of Year Number, tz=Time Zone Offset from GMT    pi = 3.1415927 <-- use PureBasic's #PI
  ;         Latitude : ( Horizontal Map Lines )                                    Longitude : ( Vertical Map Lines )
  ;[Q] = Equator = 0.0000000 degrees Latitude                         [P] = Prime Meridian = 0.0000000 degrees Longitude @ Greenwich, England
  ;[N] = Northern Hemisphere ( Positive value , North of Equator )    [E] = Eastern Hemisphere ( Positive value , East of Prime Meridian )
  ;[S] = Southern Hemisphere ( - Negative value , South of Equator )  [W] = Western Hemisphere ( - Negative value , West of Prime Meridian )
  ;
  ; http://lexikon.astronomie.info/zeitgleichung
  
  If T = 0
    T = DayOfYear(Date())
  EndIf
  
  B = Radian(lat)
  D = 0.40954*Sin(0.0172*(T-79.349740))
  zeitdiff=12*ACos((Sin(h)-Sin(B)*Sin(D))/(Cos(B)*Cos(D)))/#PI
  zeitglg=-0.1752*Sin(0.033430*T+0.5474)-0.1340*Sin(0.018234*T-0.1939)
  
  ;~~~~~ Sunrise
  AMOZ=12-zeitdiff-zeitglg;
  sunrise=AMOZ-lon/15+tz;
  
  ;~~~~~ Sunset
  UMOZ=12+zeitdiff-zeitglg;
  sunset=UMOZ-lon/15+tz-12;
  GetUp$=RSet(Str(Int(sunrise)),2,"0") + ":" + RSet(Str(Round(Mod(sunrise,1) * 60,#PB_Round_Down)),2,"0")
  GetOf$=RSet(Str(Int(sunset)),2,"0") + ":" + RSet(Str(Round(Mod(sunset,1) * 60,#PB_Round_Down)),2,"0")
  ; GetUp$=RSet(Str(Int(sunrise)),2,"0") + ":" + RSet(Str(Round(Mod(sunrise,1) * 60,#PB_Round_Up)),2,"0")
  ; GetOf$=RSet(Str(Int(sunset)),2,"0") + ":" + RSet(Str(Round(Mod(sunset,1) * 60,#PB_Round_Up)),2,"0")
  
  ProcedureReturn "Sunrise = " + GetUp$ + "AM - Sunset = " + GetOf$ + "PM"
EndProcedure

;~~~~~ Los Angeles , CA  (34.052659,-118.218384,0,-8)
;     SunCalc(Latitude,Longitude, Day of Year Number, Coordinated Universal Time = UTC or Time Zone Offset from GMT = Greenwich Mean Time [in England] )
;Debug SunCalc(34.052659,-118.218384, 0, -8)
; Buffalo, NY
Debug SunCalc(42.8864, -78.8784, 0, -4)

Re: Sunset/sunrise calculation

Posted: Thu Dec 17, 2020 8:38 am
by dige
There is a nice site and verify the results:

Sun: https://www.suncalc.org/#/42.8864,-78.8 ... 35/324.0/2
Moon: https://www.mooncalc.org/#/42.8864,-78. ... 35/324.0/2

With only 2 minutes deviation, it seems to be quite accurate.
All times are UTC+01:00
Page 1 of 1

Powered by phpBB® Forum Software © phpBB Limited