PureBasic

Code : Tout sélectionner

;+==========================================================+
;|                                                          |
;| 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))