Scripts PARI/GP

J2000-sol.gp

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);
}
    
[Source]

saisons.gp

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);
    
[Source]

julien.gp

\\ 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]);
}
[Source]

trigo.gp

\\ 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;
    
[Source]

format.gp

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]));
[Source]

taylor.gp

\\ originally contributed by Ilya Zakharevich

\\ sample function
f(x) = sin(x)

\\ plot Taylor polynomials of f of index first + i*step <= ordlim
\\ for x in [xmin,xmax].
plot_taylor(xmin=-5, xmax=5, ordlim=16, first=1, step=1) = 
{
  local(T,Taylor_array,s,t,w,h,dw,dh,cw,ch,gh,h1, extrasize = 0.6);

  default(seriesprecision,ordlim+1);
  ordlim -= first; ordlim \= step; ordlim += first;
  T = f(tt); Taylor_array = vector(ordlim+1);
  forstep(i=ordlim+1, 1, -1,
    T += O(tt^(1 + first + (i-1)*step));
    Taylor_array[i] = truncate(T)
  );

  t = plothsizes();
  w=floor(t[1]*0.9)-2; dw=floor(t[1]*0.05)+1; cw=t[5];
  h=floor(t[2]*0.9)-2; dh=floor(t[2]*0.05)+1; ch=t[6];
  h1=floor(h/1.2);

  plotinit(2, w+2*dw, h+2*dh);
  plotinit(3, w, h1);
  s = plotrecth(3, x=xmin,xmax, f(x), 2+8+16+32);
  gh=s[4]-s[3];

  plotinit(3, w, h);
  plotscale(3,s[1],s[2],s[3]-gh*extrasize/2,s[4]+gh*extrasize/2);
  plotrecth(3, x=xmin,xmax,
    concat(f(x), subst(Taylor_array, tt, x)), 4);
  plotclip(3);
  plotcopy(3, 2, dw, dh);

  plotmove(2,floor(dw+w/2-15*cw),floor(dh/2));
  plotstring(2,"Multiple Taylor Approximations");
  plotdraw([2, 0, 0])
}
[Source]

factorisation.gp

\\ JMS - 20 juillet 2007
\\ Mise en forme d'un facteur.
facteur(f) = {
    local(s);
    if(polcoeff(f[1],0),s=Str("("Strtex(f[1])")"),s=Strtex(f[1]));
    if(f[2]-1,s=Str(s"^{"f[2]"}"));
    return(s);
}
\\ Factorisation proprement dite.
factorisation(P) = {
    local(f,n,s,i,c);
    f = factor(P); 	\\ matrice des facteurs de P
    t = matsize(f)[1];	\\ nombre de facteurs de P
    \\ Concaténation des facteurs.
    s=facteur(f[1,]);
    for(i=2,t,s=Str(s" "facteur(f[i,])));
    \\ Recherche du coefficient résultant.
    c=pollead(f[1,1]);
    for(i=2,t,c=c*pollead(f[i,1]));
    c=pollead(P)/c;
    \\ Ajout du coefficient latexifié...
    if(c-1,s=Str(Strtex(c)" "s));
    return(s);
}
[Source]

sumpolyk.gp

sumpolyk(P,a_,b_) = {
    local(p,m,Q,y);
    p = poldegree(P,k);
    y = [];
    for(m=1,p+2,y=concat(y,sum(k=1,m,eval(P))));
    Q = polinterpolate(y);
    Q = Q - subst(Q,x,a_-1);
    Q = subst(Q,x,b_);
    return(factorisation(Q));
}
[Source]