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