Heures de passage au méridien / lever / coucher du soleil

rem ============================================================================
rem   Calcul des heures de passage au méridien / lever / coucher du soleil
rem ============================================================================
Init()
GUI()
end
rem ============================================================================
SUB Init()
' Label
    label clic
' Constantes
    dim k  : k  = .0172024 : ' Vitesse angulaire moyenne en radian/jour de la terre autour du soleil
    dim jm : jm = 308.67   : ' Jour du périhélie
    dim jl : jl = 21.55    : ' Jour de la longitude moyenne nulle
    dim e  : e  = .0167    : ' Excentricité de l'orbite de la terre
    dim ob : ob = .4091    : ' Obliquité de l'écliptique en radian
    dim dr : dr = pi/180
    dim hr : hr = pi/12
' Hauteur du soleil au lever et au coucher    
    dim ht : ht = -50/60 : ht = ht * dr  : ' <---- prendre aussi ht = -40/60
' Fuseau horaire et coordonnées géographiques
    dim fh : ' Fuseau horaire
    dim lo : ' lo = lo * dr : Longitude en degré, positive à l'est
    dim la : ' la = la * dr : Latitude  en degré
' Date
    dim jo : ' jour
    dim mo : ' IF mo < 30 then mo = mo + 12 : ' mois
' Heure TU du milieu de la journée
    dim h  : ' h = 12 - lo / hr
' Nombre de jours écoulés depuis le 1 Mars O h TU
    dim j  : ' j = INT(30.61 * (mo + 1)) + jo + h / 24 - 123
' Anomalie et longitude moyenne
    dim m : ' m = k * (j - jm)
    dim l : ' l = k * (j - jl)
'  Longitude vraie
    dim s : ' s = l + 2 * e * SIN(m) + 1.25 * e * e * SIN(2 * m)
' Coordonnées rectangulaires du soleil dans le repère équatorial
    dim x : ' x = cos(s)
    dim y : ' y = cos(ob) * sin(s)
    dim z : ' z = sin(ob) * sin(s)
' Equation du temps et déclinaison
    dim r  :  r = l : ' Rotation de l'angle r ...
    dim rx,ry       : ' ... autour de l'axe Z
    dim et : ' et = atn(y/x)
    dim dc : ' dc = atn(z/sqr(1-z*z))
' Heure de passage au méridien
    dim pm : ' pm = h + fh + et / hr
    dim hs : ' hs = INT(pm) : pm = 60 * (pm - hs)
    dim ms : ' ms = INT(pm) : pm = 60 * (pm - ms)
' Angle horaire au lever et au coucher
    dim ah
    dim cs : ' cs = (SIN(ht) - SIN(la) * SIN(dc)) / COS(la) / COS(dc)
             ' IF cs >  1  "Le soleil ne se lève pas"   : end
             ' IF cs < -1  "Le soleil ne se couche pas" : end
             ' IF cs = 0 then ah = pi / 2 ELSE ah = ATN(SQR(1 - cs * cs) / cs)
             ' IF cs < 0 THEN ah = ah + pi
' ------------------------------------------------------------------------------
END_SUB
rem ============================================================================
SUB GUI()
    alpha 10 : top 10,50  : left 10,10 : font_bold 10 : font_size 10,12
    font_name 10,"tahoma" : color 10,0,0,0 : font_color 10,0,255,0
    caption 10," Calcul des heures de passage au méridien / lever / coucher du soleil "
    caption 00," Calcul des heures de passage au méridien / lever / coucher du soleil "
    alpha 20 : top 20,100 : left 20,50 : font_bold 20 : font_size 20,12
    font_name 20,"tahoma" : font_color 20,0,0,255
    caption 20,"Fuseau horaire : "
    spin 30 : top 30,100  : left 30,200  : min 30,-12 : max 30,12 : font_bold 30
    hint 30,"Valeurs possibles ± 12" : position 30,1  : ' <--- Fuseau horaire de la ville de Moknine
' Pour Paris : position 30,2 : ' <--- GMT + 2
    font_color 30,255,0,0 : font_size 30,12
    alpha 40 : top 40,150 : left 40,90 : font_bold 40 : font_size 40,12
    font_name 40,"tahoma" : font_color 40,0,0,255
    caption 40,"Longitude :"
    edit 50 : top 50,150  : left 50,200 : font_bold 50 : font_color 50,255,0,0
    font_size 50,12 : hint 50,"En degrés décimaux (± 180°) positive à l'est"
    text 50,"10.9" : ' <--- C'est la longitude de la ville de Moknine (Tunisie)
' Pour Paris : text 50,"2.348"    
    alpha 60 : top 60,190 : left 60,100 : font_bold 60 : font_size 60,12
    font_name 60,"tahoma" : font_color 60,0,0,255
    caption 60,"Latitude  :"
    edit 70 : top 70,190  : left 70,200 : font_bold 70 : font_color 70,255,0,0
    font_size 70,12 : hint 70,"En degrés décimaux (± 90°) positive au Nord"
    text 70,"35.6333" : ' <--- C'est la latitude de la ville de Moknine (Tunisie)
' Pour Paris : text 70,"48.853"
    alpha 80 : top 80,230 : left 80,80 : font_bold 80 : font_size 80,12
    font_name 80,"tahoma" : font_color 80,0,0,255
    caption 80,"Date (jour) :"
    edit 90 : top 90,230 : left 90,200 : font_bold 90 : font_color 90,255,0,0
    font_size 90,12 : hint 90,"Jour"
    text 90,date_day
    alpha 100 : top 100,270 : left 100,80 : font_bold 100 : font_size 100,12
    font_name 100,"tahoma" : font_color 100,0,0,255
    caption 100,"Date (mois) :"
    edit 110 : top 110,270 : left 110,200 : font_bold 110 : font_color 110,255,0,0
    font_size 110,12 : hint 110,"Mois"
    text 110,date_month
    
    memo 120 : top 120,310 : left 120,10 : width 120,550 : height 120,100
    font_bold 120 : font_size 120,12 : font_color 120,0,255,0 : hide 120
    color 120,0,0,0 : font_name 120,"tahoma"
    button 130 : top 130,150 : left 130,400 : font_bold 130 : font_size 130,12
    caption 130,"Calculer" : on_click 130,clic
    button 140 : top 140,200 : left 140,400 : font_bold 140 : font_size 140,12
    caption 140,"Quitter" : on_click 140,clic
    
END_SUB
rem ============================================================================
' Rotation de l'angle r autour de l'axe z
SUB Rotation_Z()
    rx = 0 + COS(r) * x + SIN(r) * y
    ry = 0 - SIN(r) * x + COS(r) * y
    x  = rx : y = ry
END_SUB
rem ============================================================================
Clic:
   select number_click
       case 130 : Calculer()
       case 140 : if message_confirmation_yes_no("Vous voulez vraiment quitter ?") = 1 then terminate
   end_select
return
rem ============================================================================
SUB Calculer()
    dim_local t$,zero$
    show 120
    clear 120 : t$ = ""
    fh = val(text$(30))  : ' fuseau horaire
    lo = val(text$(50))  : lo = lo * dr : ' longitude
    la = val(text$(70))  : la = la * dr : ' latitude
    jo = val(text$(90))  : ' jour
    mo = val(text$(110)) : IF mo < 3 THEN mo = mo + 12 : ' mois
    h = 12 - (lo / hr)   : ' Heure TU du milieu de la journée
    j = int(30.61 * (mo + 1)) + jo + (h / 24) - 123 : ' Nombre de jours écoulés depuis le 1 Mars O h TU
    m = k * (j - jm)     : l = k * (j - jl) : ' Anomalie et longitude moyenne
    s = l + 2 * e * sin(m) + 1.25 * e * e * sin(2 * m) : ' Longitude vrai    
    x = COS(s) : y = COS(ob) * SIN(s) : z = SIN(ob) * SIN(s) : '  Coordonnées rectangulaires du soleil dans le repère équatorial
    r = l : Rotation_Z()
    et = ATN(y / x)                    : ' Equation du temps
    dc = ATN(z / SQR(1 - z * z))       : ' Déclinaison
    pm = h + fh + et / hr              : ' heure de passage au méridien
    hs = INT(pm) : pm = 60 * (pm - hs) : ' heure
    ms = INT(pm) : pm = 60 * (pm - ms) : ' minute
    t$ = " Heure du passage du Soleil au méridien (midi vrai)      : " + str$(hs) + ":"
    t$ = t$ + str$(ms) + ":" + str$(int(pm))
    item_add 120,t$  
    cs = (SIN(ht) - SIN(la) * SIN(dc)) / COS(la) / COS(dc) : ' Angle horaire au lever et au coucher
    if cs > 1  then item_add 120,"Le Soleil ne se lève pas !" : exit_sub
    if cs < -1 then item_add 120,"Le Soleil ne se couche pas !" : exit_sub
    if cs = 0 THEN ah = pi / 2 : ELSE : ah = ATN(SQR(1 - cs * cs) / cs)
    if cs < 0 then ah = ah + pi
    pm = h + fh + (et - ah) / hr  : ' Lever du Soleil
    IF pm < 0 THEN pm = pm + 24
    if pm > 24 then pm = pm - 24
    hs = INT(pm) : pm = 60 * (pm - hs)
    if len(str$(hs)) < 2 then zero$ = "0" : else : zero$ = ""
    t$ = " Heure du lever du Soleil ************ (heure locale) : " + zero$ + str$(hs) + ":" + str$(int(pm))
    item_add 120,t$
    pm = h + fh + (et + ah) / hr  : ' Coucher du Soleil
    IF pm > 24 THEN pm = pm - 24
    if pm < 0 then pm = pm + 24
    hs = INT(pm) : pm = 60 * (pm - hs)
    if len(str$(hs)) < 2 then zero$ = "0" : else : zero$ = ""
    t$ = " Heure du coucher du Soleil ********** (heure locale) : " + zero$ + str$(hs) + ":" + str$(int(pm))
    item_add 120,t$
        
END_SUB
rem ============================================================================