Sunset/sunrise calculation

Just starting out? Need help? Post your questions and find answers here.
Post Reply
l1marik
User
User
Posts: 52
Joined: Tue Jun 30, 2020 6:22 am

Sunset/sunrise calculation

Post 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
Top
User avatar
NicTheQuick
Addict
Addict
Posts: 1539
Joined: Sun Jun 22, 2003 7:43 pm
Location: Germany, Saarbrücken
Contact:
Contact NicTheQuick

Re: Sunset/sunrise calculation

Post by NicTheQuick »

It's a bit complicated but I found this formula for you: https://www.esrl.noaa.gov/gmd/grad/solc ... areqns.PDF
The english grammar is freeware, you can use it freely - But it's not Open Source, i.e. you can not change it or publish it in altered way.
Top
infratec
Always Here
Always Here
Posts: 7712
Joined: Sun Sep 07, 2008 12:45 pm
Location: Germany

Re: Sunset/sunrise calculation

Post by infratec »

If you simply enter sunrise in the search field, you will find:

viewtopic.php?f=12&t=64652&hilit=sunrise
Top
Mesa
Enthusiast
Enthusiast
Posts: 460
Joined: Fri Feb 24, 2012 10:19 am

Re: Sunset/sunrise calculation

Post 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.
Top
CalamityJames
User
User
Posts: 82
Joined: Sat Mar 13, 2010 4:50 pm

Re: Sunset/sunrise calculation

Post 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).
Top
ebs
Enthusiast
Enthusiast
Posts: 564
Joined: Fri Apr 25, 2003 11:08 pm

Re: Sunset/sunrise calculation

Post 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)
Top
dige
Addict
Addict
Posts: 1430
Joined: Wed Apr 30, 2003 8:15 am
Location: Germany
Contact:
Contact dige

Re: Sunset/sunrise calculation

Post 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.
"Daddy, I'll run faster, then it is not so far..."
Top
Post Reply

Return to “Coding Questions”