Code : Tout sélectionner
Je reposte au propre, cette version [u]non boguée[/u], de mon code qui est déjà quelque part sur le forum.
;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)