Cette page fait suite à la détermination des saisons, les scripts présentés ici mettent en œuvre des calculs à partir d'éléments simplifiés de la théorie (Variations Séculaires des Orbites Planétaires) de Bretagnon et Simon. Ces éléments sont tirés du livre :
Les calculs sont effectués avec pari-gp. Un formulaire permet de tester en ligne le code décrit ici.
Les calculs portent ici sur les coordonnées du soleil dans un repère géocentrique. La variable t que l'on retrouve dans les fonctions exprime le temps en siècles à compter du 1 janvier 2000 à 12 heures. La signification des différentes quantités définies ici peut être laissée de côté dans un premier temps, pour l'instant elle vont servir à déterminer l'instant du passage du soleil à son périgée et ceci dans le but de déterminer, à l'aide de l'équation de Kepler seulement, les dates des saisons de l'année qui suit.
read("julien.gp");
read("trigo.gp");
read("format.gp");
\\ Excentricité de l'orbite solaire
ExcentriciteSoleil(t) = 0.016708617-0.000042037*t-0.0000001236*t*t;
\\ Longitude moyenne du soleil rapportée à l'équinoxe moyen de la date
LmEmSoleil(t) = vpr(d2r(280.4665)+d2r(36000.76983)*t+d2r(0.0003032)*t*t);
\\ Anomalie moyenne du soleil
AmEmSoleil(t) = vpr(d2r(357.52910)+d2r(35999.05030)*t-d2r(0.0001559)*t*t);
\\ Equation du centre
EcEmSoleil(t) =
{
local(M,C);
M = AmEmSoleil(t);
C = (d2r(1.9146)-d2r(0.004817)*t-d2r(0.000014)*t*t)*sin(M);
C += (d2r(0.019993)-d2r(0.000101)*t)*sin(2*M);
C += d2r(0.00029)*sin(3*M);
return(vpr(C));
}
\\ Longitude vraie du soleil
LvEmSoleil(t) = vpr(LmEmSoleil(t)+EcEmSoleil(t));
\\ Longitude apparente du soleil (rapportée à l'équinoxe vrai de la date)
LaEvSoleil(t) =
{
local(Omega);
Omega = d2r(125.04)-d2r(1934.136)*t;
return(vpr(LvEmSoleil(t)-d2r(0.0059)-d2r(0.00478)*sin(Omega)));
}
\\ Recherche du périgée à partir de l'équation M = 0
TempsPerigeeSoleilM(A) =
{
local(jd,i,M,MT,t);
i = 1; t = 0; jd = JD2000(1,1,A,0,0,0); M = AmEmSoleil((jd+t)/36525);
while (i>10^(-7) ,
t += i;
MT = AmEmSoleil((jd+t)/36525);
if (MT<M,
t -= i; i /= 2,
M = MT
);
);
return(floor(t*10^7)/10^7+0.0);
}
La fonction TempsPerigeeSoleilM procède ainsi : en partant de l'anomalie moyenne au premier janvier de l'année considérée, on incrémente d'une journée jusqu'à obtenir une anomalie inférieure (comme cette anomalie est toujours évaluée entre 0 et 2Pi, c'est le basculement à 0 en passant par 2Pi), ce qui correspond à un instant voisin de celui recherché. On affine ensuite à l'aide d'une dichotomie. Le résulat retourné est donc le nombre de jour entre le 1 janvier à 0 heure et l'instant du périgée.
Ce script contient la définition des fonctions calculant les dates des saisons, ces définitions sont précédées de la résolution de l'équation de Kepler (tpourl).
read("J2000-sol.gp");
\\ v - u où v est l'anomale vraie et u l'anomalie excentrique
VmoinsU(u,e) = asin((e-(1-sqrt(1-e*e))*cos(u))*sin(u)/(1-e*cos(u)));
\\ Calcul de t connaissant l pour un passage au périgée caractérisé
\\ par tO et lO.
tpourl(l,e,lO,tO) =
{
local(u,v,i);
v = l-lO;
u = v;
for(i=1,20,u = v-VmoinsU(u,e));
return((u-e*sin(u))*365.25/2/Pi+tO);
}
\\ Instant du début de la saison, A l'année et l (un multiple de Pi/2)
\\ la longitude caractérisant la saison.
SAISON(A,l) =
{
local(jd,tO,lO,e);
jd = JD2000(1,1,A,0,0,0);
tO = TempsPerigeeSoleilM(A);
lO = LaEvSoleil((jd+tO)/36525);
e = ExcentriciteSoleil((jd+tO)/36525);
return(CD(tpourl(l,e,lO,tO)+JD(1,1,A,0,0,0)));
}
printemps(A) = SAISON(A,2*Pi);
ete(A) = SAISON(A,2.5*Pi);
automne(A) = SAISON(A,3*Pi);
hiver(A) = SAISON(A,3.5*Pi);
Dans les scripts précédents des données angulaires sont manipulées, les fonctions suivantes assurent quelques transformations ou réductions.
\\ Conversion de degrés en radians
d2r(d) = d/180*Pi;
\\ Conversion de radians en degrés
r2d(r) = r/Pi*180;
\\ Conversion de degrés décimaux en degrés, minutes, secondes décimales
d2dms(d) =
{
local(m,s);
m = (d-floor(d))*60;
s = (m-floor(m))*60;
return([floor(d),floor(m),s]);
}
\\ Réduction d'un angle à sa valeur entre 0 et 2*Pi (valeur principale)
vpr(r) = r-floor(r/2/Pi)*2*Pi;
Les dates sont calculées en jours juliens, la conversion en dates ordinaires est assurée par les procédures suivantes.
\\ Julien est en l'honneur de Jules Scaliger, le père du mathématicien
\\ qui a introduit cette échelle de temps en 1583. Le calendrier
\\ grégorien date de 1582.
\\ Calcul de la date en jours juliens
JD(J,M,A,h,m,s) =
{
local(AT,BT,JDT) ;
if (M<3 , A = A-1; M = M+12);
AT = floor(A/100);
BT = 2-AT+floor(AT/4);
JDT = floor(365.25*(A+4716))
+floor(30.6001*(M+1))
+J
+(h+(m+s/60)/60)/24
+BT
-1524.5;
return(JDT);
}
\\ Calcul de la date en jours juliens rapportée à l'époque 2000
\\ Utile pour l'utilisation des théories planétaires actuelles
JD2000(J,M,A,h,m,s) = JD(J,M,A,h,m,s)-JD(1,1,2000,12,0,0);
\\ Calcul de la date calendaire à partir du JD
\\ La valeur retounée est le vecteur [J,M,A,h,m,s,JS]
\\ (JS est le jour de la semaine)
CD(jd) =
{
local(a,A,B,C,D,E,F,G,h,m,s,Z);
Z = floor(jd+0.5); F = jd+0.5-Z;
if (Z>=2299161,
a = floor((Z-1867216.25)/36524.25);
A = Z+1+a-floor(a/4),
A = Z;
);
B = A+1524;
C = floor((B-122.1)/365.25);
D = floor(365.25*C);
E = floor((B-D)/30.6001);
G = B-D-floor(30.6001*E)+F;
h = (G-floor(G))*24;
m = (h-floor(h))*60;
s = (m-floor(m))*60;
if (E<14,E=E-1,E=E-13);
if (floor(E)>2,C=C-4716,C=C-4715);
return([floor(G),E,C,floor(h),floor(m),floor(s),(Z+1)%7]);
}
MOIS =
["janvier","février","mars","avril","mai","juin",\
"juillet","août","septembre","octobre","novembre","décembre"];
JOURS =
["Dimanche","Lundi","Mardi","Mercredi","Jeudi","Vendredi","Samedi"];
\\ format degrés, minutes, secondes
format_dms(d) = Str(d[1]"°"d[2]"'"floor(d[3])"''");
\\ format j mois AAAA, h:m:s
format_date(d) = Str(JOURS[d[7]+1]" "d[1]" "MOIS[d[2]]" "d[3]\
", "d[4]":"d[5]":"floor(d[6]));
Le script suivant ordonne les calculs.
#!/bin/sh
cat << *** | gp -q
read("saisons.gp");
read("format.gp");
read("trigo.gp");
a = TempsPerigeeSoleilM($1);
b = JD(1,1,$1,0,0,0);
print(format_date(CD(a+b)));
print(format_date(printemps($1)));
print(format_date(ete($1)));
print(format_date(automne($1)));
print(format_date(hiver($1)));
***
Pour l'invoquer sur l'année 2002, il suffit d'entrer en console:
Les résultats retournés sont alors:
Jeudi 3 janvier 2002, 12:37:49 Mercredi 20 mars 2002, 19:13:18 Vendredi 21 juin 2002, 13:26:19 Lundi 23 septembre 2002, 5:4:56 Dimanche 22 décembre 2002, 1:21:15
Pour comparer, voici les dates des saisons telles qu'on peut les trouver dans l'annuaire du Bureau des Longitudes (éphémérides astronomiques 2002):
L'écart maximal est donc de 10 mn, ce qui est peu. Il provient du fait que nous avons considéré le mouvement de la terre autour du soleil (ou du soleil autour de la terre) comme étant purement elliptique, ce qui revient, en premier lieu, à négliger toute influence des autres planètes du système solaire.