I have added the 2005er Version from this page that also use the year for calculation for a bit more accuracy and used the code above for it:
https://lexikon.astronomie.info/zeitgleichung/neu.html
But I'm not 100% sure if all is correct, don't understand all calculations there. The results seams to be correct, but the example C++ Code on the end of the page have more formulas in it.
For example: the page says for the 01.01.2000 that the noon time is used, so Date(2000,1,1,12,0,0). But should I use for the date for calculation also the noon time? Was not sure so I used noon because sounds more logical for me, but must not be right.
Edit: Ok, there is a error in it, if I try for example 30. May 2019, I get 16:54AM and 21:39PM, 12h too much, correct would be 4:54AM and 9:39PM. Don't know where the error is.
Edit2: must add 2 If Lines at the Code, hope this will fix it.
Edit3: forget it, only fixed another error. XD
Code: Select all
EnableExplicit
ImportC ""
time(*tloc = #Null) ; easiest method to get the wordtime without timezone
EndImport
#PI2 = #PI * 2 ; 2PI
Procedure.s SunCalc(lat.f, lon.f, ti.i, tz.d)
Define GetUp$, GetOf$
Define.d T, B, M, L, e, DK, RA, RAm, h = -0.0145, zeitglg, zeitdiff, AMOZ, UMOZ, sunrise, sunset
; lat = latitude (N/S), lon=longitude (E/W), T=Unix Timestamp Worldtime, 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 )
;
; https://lexikon.astronomie.info/zeitgleichung/neu.html
If ti = 0
ti = time()
EndIf
T = ((Date(Year(ti),Month(ti),Day(ti),12,0,0)-946728000)/86400)/36525 ; 946728000 = Date(2000,1,1,12,0,0)
B = Radian(lat)
M = #PI2 * (0.993133 + 99.997361 *T) ; radians
L = #PI2 * (0.7859453 + M/#PI2 + (6893*Sin(M)+72*Sin(2*M)+6191.2*T) / 1296000) ; radians
e = #PI2 * (23.43929111 + (-46.8150*T - 0.00059*T*T + 0.001813*T*T*T)/3600)/360 ; radians
DK = ASin(Sin(e)*Sin(L)) ; radians
RA = ATan(Tan(L)*Cos(e)) ; radians
If RA < 0 : RA + #PI : EndIf
If L > #PI : RA + #PI : EndIf
RAm = 18.71506921 + 2400.0513369*T +(0.25862 - 0.000000172*T)*T*T ; Stunden
RAm - (Round(RAm/24,#PB_Round_Down)*24)
zeitglg = 1.0027379 * (RAm - (RA*3.81972)) ; Stunden
zeitdiff = 12*ACos((Sin(h) - Sin(B)*Sin(DK)) / (Cos(B)*Cos(DK)))/#PI ; Stunden
;~~~~~ Sunrise
AMOZ=12-zeitdiff-zeitglg;
sunrise=AMOZ-lon/15.0+tz;
;~~~~~ Sunset
UMOZ=12+zeitdiff-zeitglg;
sunset=UMOZ-lon/15.0+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")
ProcedureReturn "Sunrise = " + GetUp$ + "AM - Sunset = " + GetOf$ + "PM"
EndProcedure
;~~~~~ Los Angeles , CA (34.052659,-118.218384,0,-8)
; SunCalc(Latitude,Longitude, Unix Timestamp = Worldtime without timezone in seconds since 01.01.1970, Coordinated Universal Time = UTC or Time Zone Offset from GMT = Greenwich Mean Time [in England] )
Debug SunCalc(54.3, 10.1, 0, 1)
