unit UHeureSoleil; // version 11-11-04
// unite pour calcul du lever et coucher du soleil à PARIS par défaut
// appel par fonction : Get_SunTime
// 11-11-04 mise en forme et commentaires

// auteur : alantell - novembre 2004

// origine : livres sur l'astronomie

interface

uses SysUtils;

type
THSunRecord = record
SAubeReel: TDateTime; // reel
SCrepReel: TDateTime;
SAubeCivi: TDateTime; // civil
SCrepCivi: TDateTime;
SAubeNaut: TDateTime; // nautique
SCrepNaut: TDateTime;
SAubeAstr: TDateTime; // astronomique
SCrepAstr: TDateTime;
end;

var
HSun: THSunRecord;

// retourne horaire civil, nautique, astronomique du lever et coucher du soleil
function Get_SunTime(gDate: TDateTime; gLatLieuDc:Double; gLonLieuDc:Double; gGMTLieu:Double): THSunRecord;

implementation

function Get_SunTime(gDate: TDateTime; gLatLieuDc:Double; gLonLieuDc:Double; gGMTLieu:Double): THSunRecord;
// AHvo : angle horaire du soleil à 0h par rapport à Greenwich
// AHMLi : angle horaire du méridien du lieu d'observation
// AHMSo : angle horaire du méridien du soleil
// LVSo : longitude vraie du soleil
// ASDSo : ascension droite du soleil
// DclSo : déclinaison du soleil
// AziSo : azimut du soleil
// HSo : hauteur du soleil au dessus de l'horizon
// AMSo : anomalie moyenne du soleil

/////////////////////// fonction internes //////////////////////////

function Interval2PiRd(X:double):double;
label Test2pi;
begin
Test2pi:
if X < 0 then begin X:= X + Pi * 2; goto Test2pi; end;
if (X = 2 * Pi) or (X > 2 * Pi) then begin X:= X - Pi * 2; goto Test2pi; end;
result:= X;
end;
// fonction ArcSinus
function ArcSin(X:double):double;
begin result:= ArcTan( X / Sqrt(-X * X +1)); end;
// fonction ArcCosinus
function ArcCos(X:double):double;
begin result:= ArcTan(-X / Sqrt(-X * X +1)) + ( Pi/2 ); end;
// fonction de conversion heure décimale en heure, minute
function HrDcToHrMn(X:double):double;
var H,M: Byte;
begin H:= Trunc(X); M:= Trunc(frac(X) * 60); result:= H + ( M/100 ); end;
// fonction conversion valeur sexagécimale en décimale
function SexToDec(X:double):double;
var Sign: ShortInt; H,M,S: Integer; Z: Double;
begin
if X < 0 then Sign:= -1 else Sign:= 1;
X:= Abs(X + 0.00005);
H:= Trunc(X); Z:= Frac(X) * 100; // heure ou degrés
M:= Trunc(Z); S:= Trunc(Frac(Z) * 100); // minute seconde
result:= (H + (M/60) + (S/3600)) * Sign;
end;

//////////////////////////////////////////////////////////////////////

var NumJrB, N, Crepuscule: Byte; Sign: ShortInt;
// gLonLieuDc, gLatLieuDc, gGMTLieu: Double;
gLonLieuRd, gLatLieuRd, gGMTLieuHd: Double;
AHMLiRd, AHvgRd, AMSo, LVSoRd, ASDSoRd, DclSoRd: Double;
SinHSoRd, AziSoRd, TpsLRd, AHMSoRd, TpsTURd: Double;
PtM,NumJr, Z, SinX, cosAz, TpsHRd, TpsHorlog, SW: Double;
AAnnee, AMois, AJour, HrSol,MnSol: Word;
PJAn, HorSol: TDateTime;
begin
// nota: il suffit de changer les coordonnées pour une autre ville
// coordonnées de Paris en degré décimal
// gLatLieuDc:= 48.50; // latitude
// gLonLieuDc:= -2.20; // longitude
// gGMTLieu:= 1; // décalage horaire
// coordonnées de Paris en radian
gLonLieuRd:= SexToDec(gLonLieuDc) * Pi/180;
gLatLieuRd:= SexToDec(gLatLieuDc) * Pi/180;
gGMTLieuHd:= SexToDec(gGMTLieu);
// NumJr = n° du jour depuis le 0 janvier de l'année en cours
DecodeDate(gDate,AAnnee,AMois,Ajour);
PJAn:= EncodeDate(AAnnee,1,1);
NumJr:= gDate - PJAn + 1;
// HrFusHd:= - gGMTLieuHd; changement de signe du décalage horaire
// AHMLiRd:= HrFusHd * Pi/12; angle horaire du méridien du lieu observation en rd
// angle horaire en rad du méridien du lieu observation
AHMLiRd:= -gGMTLieuHd * Pi/12;

//////////// calcul du lever à Pi / 2 /////////////////////
PtM:= Pi / 2;
//temps approximatif = angle position soleil sur écliptique suivant nb de jours
// = nombre de jour + ((PtM + longitude ) / 2 * Pi)
AHvgRd:= NumJr + ((PtM + gLonLieuRd) / (2 * Pi));
// anomalie moyenne du soleil 1j -> 2 * Pi / 365,25 = 0,01720242383896
// AMSo:= AMSo - 0.0574039;
AMSo:= (AHvgRd * 2 * Pi / 365.25) - 0.0574039;
// LVSo = longitude vraie du soleil
LVSoRd:= AMSo + Sin(AMSo) * 0.0334405 + Sin(2*AMSo) * 0.000349066 + 4.93289;
LVSoRd:= Interval2PiRd(LVSoRd); // intervalle 0-2Pi
Z:= frac(LVSoRd / (Pi / 2)); // déterminaison du cadran en divisant par Pi/2
// si Z = 0 on ajoute 1s d'arc soit 0.00000484814 rd
if Z = 0 Then LVSoRd:= LVSoRd + 0.00000484814;
N:= 0;
if LVSoRd > (Pi / 2) then N:= 1;
if LVSoRd > (3 * Pi / 2) then N:= 2;
//////// ascension droite du soleil /////////
// ASDSoRd = ArcTan(tan(LVSoRd) * cos(déclinaison)
// cos(déclinaison) : cos 23°26,4' = 0.91746 // tan(LVSoRd)
ASDSoRd:= sin(LVSoRd) / cos(LVSoRd);
ASDSoRd:= ArcTan(ASDSoRd * 0.91746);
// ajustement du cadran
case N of
0: ASDSoRd:= ASDSoRd;
1: ASDSoRd:= ASDSoRd + Pi;
2: ASDSoRd:= ASDSoRd + (Pi * 2);
end;
//////// déclinaison du soleil /////////
// DclSoRd = ArcSin(sin(LVSoRd) * sin(déclinaison))
// sin(déclinaison) : sin 23°26,4' = 0.39782
SinX:= sin(LVSoRd) * 0.39782;
DclSoRd:= ArcSin(SinX);
//////// conversion des coordonnées //////////
// crépuscule lorsque le centre du soleil est en dessous de l'horizon
// angle H = centre du soleil par rapport à l'horizon d'où H < 0°
// crépuscule astronomique H = 18° crépuscule nautique H = 12°
// crépuscule civil H = 6° crépuscule réel H = 0°50'
// valeur du sin(H) suivant option choisie
for Crepuscule:= 1 to 4 do begin
SinHSoRd:= -0.0; // initialisation
case Crepuscule of
1: SinHSoRd:= -0.014543897651582; // HSo = - 0° 50'
2: SinHSoRd:= -0.104528463267653; // HSo = - 6°
3: SinHSoRd:= -0.207911690817759; // HSo = - 12°
4: SinHSoRd:= -0.309016994374947; // HSo = - 18°
end;
cosAz:= SinHSoRd - (Sin(DclSoRd) * Sin(gLatLieuRd));
cosAz:= cosAz / (Cos(DclSoRd) * Cos(gLatLieuRd));
Z:= Abs(cosAz);
if Z <= 1 then begin // ajustement azimut du soleil
AziSoRd:= ArcCos(cosAz); // au coucher: Azi = arccos(cosAz)
AziSoRd:= 2 * Pi - AziSoRd; // au lever : Azi = 2 * pi - arccos(cosAz)
// heure local apparent du noeud ascendant
// année tropique entre 2 équinoxes = 365j 5h 48mn 45,975s
// en jour décimal 365.24219878473 j
Z:= AHvgRd * 2 * pi / 365.24219878473;
TpsLRd:= AziSoRd + ASDSoRd - Z - 1.73364; // angle horaire local
//AHMSoRd:= TpsLRd;
TpsTURd:= TpsLRd + gLonLieuRd; // temps universel
TpsHRd:= TpsTURd - AHMLiRd; // temps horloge
Z:= Interval2PiRd(TpsHRd) * 12/Pi; // Décimal en sexagésimal: angle*12/Pi
TpsHorlog:= HrDcToHrMn(Z); // heure décimale en heure, minute
HrSol:= Trunc(TpsHorlog); MnSol:= Round(frac(TpsHorlog)*100);
HorSol:= EncodeTime(HrSol,MnSol,0,0);
// affichage horaire
case Crepuscule of
1: HSun.SAubeReel:= HorSol;
2: HSun.SAubeCivi:= HorSol;
3: HSun.SAubeNaut:= HorSol;
4: HSun.SAubeAstr:= HorSol;
end;
end
else TpsHorlog:= 0;
end;

///////////// calcul du coucher à 3/2 de Pi ////////////////
PtM:= Pi * 3 / 2;
AHvgRd:= NumJr + ((PtM + gLonLieuRd) / (2 * Pi));
AMSo:= (AHvgRd * 2 * Pi / 365.25) - 0.0574039;
LVSoRd:= AMSo + Sin(AMSo)*0.0334405 + Sin(2*AMSo)*0.000349066 + 4.93289;
LVSoRd:= Interval2PiRd(LVSoRd);
Z:= frac(LVSoRd / (Pi / 2)); if Z = 0 then LVSoRd:= LVSoRd + 0.00000484814;
N:= 0;
if LVSoRd > (Pi / 2) then N:= 1; if LVSoRd > (3 * Pi / 2) then N:= 2;
ASDSoRd:= sin(LVSoRd) / cos(LVSoRd); ASDSoRd:= ArcTan(ASDSoRd * 0.91746);
case N of
0: ASDSoRd:= ASDSoRd;
1: ASDSoRd:= ASDSoRd + Pi;
2: ASDSoRd:= ASDSoRd + (2 * Pi);
end;
SinX:= sin(LVSoRd) * 0.39782; DclSoRd:= ArcSin(SinX);
for Crepuscule:= 1 to 4 do begin
case Crepuscule of
1: SinHSoRd:= -0.014543897651582; 2: SinHSoRd:= -0.104528463267653;
3: SinHSoRd:= -0.207911690817759; 4: SinHSoRd:= -0.309016994374947;
end;
cosAz:= SinHSoRd - (Sin(DclSoRd) * Sin(gLatLieuRd));
cosAz:= cosAz / (Cos(DclSoRd) * Cos(gLatLieuRd));
Z:= Abs(cosAz);
if Z <= 1 then begin
AziSoRd:= ArcCos(cosAz);
// AziSoRd:= 2 * Pi - AziSoRd;
Z:= AHvgRd * 2 * pi / 365.24219878473;
TpsLRd:= AziSoRd + ASDSoRd - Z - 1.73364;
//AHMSoRd:= TpsLRd;
TpsTURd:= TpsLRd + gLonLieuRd;
TpsHRd:= TpsTURd - AHMLiRd;
Z:= Interval2PiRd(TpsHRd) * 12/Pi;
TpsHorlog:= HrDcToHrMn(Z);
HrSol:= Trunc(TpsHorlog);
MnSol:= Round(frac(TpsHorlog)*100);
HorSol:= EncodeTime(HrSol,MnSol,0,0);
case Crepuscule of
1: HSun.SCrepReel:= HorSol;
2: HSun.SCrepCivi:= HorSol;
3: HSun.SCrepNaut:= HorSol;
4: HSun.SCrepAstr:= HorSol;
end;
end
else TpsHorlog:= 0;
end;
result:= HSun;
end;

end.