#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

const double pi=3.14159265358979323846;

typedef int planete;

struct Position {
	double Dist;
	double Dec;
	double Asc;
	double t;
	double k ;
	double JJ;
};
struct Date {
	int jour;
	int mois;
	int annee;
};
typedef struct Position position;
typedef struct Date date;

void affhms (double t); // affiche t° en heures, minutes, secondes et décimales.
void affds (double t); // affiche t° en degrés, secondes et décimales.
void aff (); // affiche le tableau de résultats

date dat;
planete k;

void aff_nom_planete (planete k)
{ switch (k)
{ case 0:cout<<"Soleil ";break;
case 1:cout<<"Lune ";break;
case 2:cout<<"Mercure";break;
case 3:cout<<"Venus ";break;
case 4:cout<<"Mars ";break;
case 5:cout<<"Jupiter";break;
case 6:cout<<"Saturne";break;
case 7:cout<<"Uranus ";break;
case 8:cout<<"Neptune";break;

}
}

double JJ(date d)
{ int y,m,a,b;
double jj;
if (d.mois>2) {y=d.annee;m=d.mois;}
if ((d.mois==1)||(d.mois==2)) {y=d.annee-1;m=d.mois+12;}
if (10000*d.annee+100*d.mois+d.jour>=15821015) {a=floor(y/100);b=2-a+floor(a/4);} else b=0;
jj=floor(365.25*y)+floor(30.6001*(m+1))+d.jour+1720994.5+b;
return jj;
}

position CalculPosition(planete k,date dat)
{

position p;
p.k=0;
double T; // temps sidéral (en siecles juliens de 36525 jours à partir du 0,5 janvier 1900 TE)
double t; // temps sidéral (en radian)
double o; // obliguité de l'écliptique
double L; // Longitude moyenne
double M; // anomalie moyenne
double E; // anomalie excentrique
double e; // exentricité
double a; // demi grand axe (constant)
double v; // anomalie vrai
double R; // Rayon vecteur en UA
double Lg; // Longituusing namespace std;de vrai du soleil

p.JJ=JJ(dat);
double J=floor(p.JJ+0.5)-0.5; // jour julien à 0H00
double j=p.JJ-J; // portion du jour julien restant
T=(J-2415020.0)/36525;
t=0.276919398+100.0021359*T+0.000001075*T*T; // temps sidéral à 0H00 GMT exprimé en révolutions

t=t+j*1.002737908; // temps sidéral à l'heure correspondant
// à la portion de jour jj
t=fmod(t,1);
t=t*2*pi; // en radian
p.t=t;

T=(JJ(dat)-2415020.0)/36525;
o=23.452294-0.0130125*T-0.00000164*T*T+0.000000503*T*T*T; o=fmod(o,360); o=o*pi/180;
L=279.69668+36000.76892*T+0.0003025*T*T; L=fmod(L,360); L=L*pi/180;
M=358.47583+35999.04975*T-0.000150*T*T-0.0000033*T*T*T; M=fmod(M,360); M=M*pi/180;
e=0.01675104-0.0000418*T-0.000000126*T*T;

E=M;
double erreur=1;
//while (erreur>1e-7 /*1e-4*/) {E1=M+e*sin(E);erreur=fabs(E1-E);E=E1;}
while (erreur>1e-9) {erreur = (M+e*sin(E)-E)/(1-e*cos(E));E += erreur;}
v=2*atan(sqrt((1+e)/(1-e))*tan(E/2));
a=1.0000002;
R=a*(1-e*cos(E));
Lg=L+v-M;

if (k!=1) // différent de Lune
{

double A; // Ascention droite
double D; // Déclinaison
double i,w,W;
double r,u,l,b;
// l,b,r coordonnées écliptique héliocentriques; longitude l, latitude b, rayon r
double lambda,delta,beta;

if (!k) { p.Asc=atan2(cos(o)*sin(Lg),cos(Lg));
p.Dec=asin(sin(o)*sin(Lg));
p.Dist=R; }
else
{

switch (k) {
case 2: // mercure:
L= 178.179078 + 149474.07078 *T + 0.0003011 *T*T ;
a= 0.3870986 ;
e= 0.20561421 + 0.00002046 *T - 0.000000030 *T*T ;
i= 7.002881 + 0.0018608 *T - 0.0000183 *T*T ;
w= 28.753753 + 0.3702806 *T + 0.0001208 *T*T ;
W= 47.145944 + 1.1852083 *T + 0.0001739 *T*T ;
break;
case 3: // venus:
L= 342.767053 + 58519.21191 *T + 0.0003097 *T*T ;
a= 0.7233316 ;
e= 0.00682069 - 0.00004774 *T + 0.000000091 *T*T ;
i= 3.393631 + 0.0010058 *T - 0.0000010 *T*T ;
w= 54.384186 + 0.5081861 *T - 0.0013864 *T*T ;
W= 75.779647 + 0.8998500 *T + 0.0004100 *T*T ;
break;
case 4: // mars:
L= 293.737334 + 19141.69551 *T + 0.0003107 *T*T ;
a= 1.5236883 ;
e= 0.09331290 + 0.000092064 *T - 0.000000077 *T*T ;
i= 1.850333 - 0.0006750 *T + 0.0000126 *T*T ;
w= 285.431761 + 1.0697667 *T + 0.0001313 *T*T + 0.00000414 *T*T*T ;
W= 48.786442 + 0.7709917 *T - 0.0000014 *T*T - 0.00000533 *T*T*T ;
break;
case 5: // jupiter:
L= 238.049257 + 3036.301986 *T + 0.0003347 *T*T - 0.00000165 *T*T*T ;
a= 5.202561 ;
e= 0.04833475 + 0.000164180 *T - 0.0000004676 *T*T - 0.0000000017 *T*T*T;
i= 1.308736 - 0.0056961 *T + 0.0000039 *T*T ;
w= 273.277558 + 0.5994317 *T + 0.00070405 *T*T + 0.00000508 *T*T*T ;
W= 99.443414 + 1.0105300 *T + 0.00035222 *T*T - 0.00000851 *T*T*T ;
break;
case 6: // saturne:
L= 266.564377 + 1223.509884 *T + 0.0003245 *T*T - 0.0000058 *T*T*T ;
a= 9.554747 ;
e= 0.05589232 - 0.00034550 *T - 0.000000728 *T*T + 0.00000000074 *T*T*T;
i= 2.492519 - 0.0039189 *T - 0.00001549 *T*T + 0.00000004 *T*T*T ;
w= 338.307800 + 1.0852207 *T + 0.00097854 *T*T + 0.00000992 *T*T*T ;
W= 112.790414 + 0.8731951 *T - 0.00015218 *T*T - 0.00000531 *T*T*T ;
break;
case 7: // uranus:
L= 244.197470 + 429.863546 *T + 0.0003160 *T*T - 0.00000060 *T*T*T ;
a= 19.21814 ;
e= 0.0463444 - 0.00002658 *T + 0.000000077 *T*T ;
i= 0.772464 + 0.0006253 *T + 0.0000395 *T*T ;
w= 98.071581 + 0.9857650 *T - 0.0010745 *T*T - 0.00000061 *T*T*T ;
W= 73.477111 + 0.4986678 *T + 0.0013117 *T*T ;
break;
case 8: // neptune:
L= 84.457994 + 219.885914 *T + 0.0003205 *T*T - 0.00000060 *T*T*T ;
a= 30.10957 ;
e= 0.00899704 + 0.000006330 *T - 0.000000002 *T*T ;
i= 1.779242 - 0.0095436 *T - 0.0000091 *T*T ;
w= 276.045975 + 0.3256394 *T + 0.00014095 *T*T + 0.000004113 *T*T*T;
W= 130.681389 + 1.0989350 *T + 0.00024987 *T*T - 0.000004718 *T*T*T;
break;
}

M=L-w-W;

double cu, cP, cQ, cV, cZ;
double cH, cS, cN, cT;
double cA, cB, cC, cD;

//----------- perturbation par attractions mutuelles ------------

switch (k) {
case 5: // jupiter:
cu = 0.1 + T/5;
cP = 237.475 + 3034.9061 *T; // angle en degrés
cQ = 265.916 + 1222.1139 *T; // angle en degrés
cV = 5*cQ - 2*cP; // angle en degrés
cZ = cQ - cP; // angle en degrés

cP=fmod(cP,360);cP=cP*pi/180;
cQ=fmod(cQ,360);cQ=cQ*pi/180; // conversion degré -> radian
cV=fmod(cV,360);cV=cV*pi/180;
cZ=fmod(cZ,360);cZ=cZ*pi/180;

// perturbation dans la longitude moyenne, en degrés:
cA = + (0.3314 - 0.0103*cu - 0.0047*cu*cu)*sin(cV)
+ (0.0032 - 0.0644*cu + 0.0021*cu*cu)*cos(cV)
+ 0.0136*sin(cZ) + 0.0185*sin(2*cZ) + 0.0067*sin(3*cZ)
+ (0.0073*sin(cZ) + 0.0064*sin(2*cZ) - 0.0338*cos(cZ))*sin(cQ)
- (0.0357*sin(cZ) + 0.0063*cos(cZ) + 0.0067*cos(2*cZ))*cos(cQ);
// perturbation dans l'excentricité, en unités de la sixième décimale:
cB = + (361 + 13*cu)*sin(cV) + (129 - 58*cu)*cos(cV)
+ (128*cos(cZ) - 676*sin(cZ) - 111*sin(2*cZ))*sin(cQ)
+ (146*sin(cZ) - 82 + 607*cos(cZ) + 99*cos(2*cZ) + 51*cos(3*cZ))*cos(cQ)
- (96*sin(cZ) + 100*cos(cZ))*sin(2*cQ)
- (96*sin(cZ) - 102*cos(cZ))*cos(2*cQ);
// perturbation dans la longitude du périhélie, en degrés:
cC = + (0.0072 - 0.0031*cu)*sin(cV) - 0.0204*cos(cV)
+ (0.0073*sin(cZ) + 0.0340*cos(cZ) + 0.0056*cos(2*cZ))*sin(cQ)
+ (0.0378*sin(cZ) + 0.0062*sin(2*cZ) - 0.0066*cos(cZ))*cos(cQ)
- 0.0054*sin(cZ)*sin(2*cQ) + 0.0055*cos(cZ)*cos(2*cQ);
// perturbation dans le demie-grand axe, en unités de la sixième décimale:
cD = - 263*cos(cV) + 205*cos(cZ) + 693*cos(2*cZ) + 312*cos(3*cZ)
+ 299*sin(cZ)*sin(cQ)
+ (204*sin(2*cZ)-337*cos(cZ))*cos(cQ);

L += cA;
M += cA - cC/e;
e += cB/1000000;
a += cD/1000000;
break;

case 6: // saturne:
cu = 0.1 + T/5;
cP = 237.475 + 3034.9061 *T; // angle en degrés
cQ = 265.916 + 1222.1139 *T; // angle en degrés
cV = 5*cQ - 2*cP; // angle en degrés
cZ = cQ - cP; // angle en degrés

cP=fmod(cP,360);cP=cP*pi/180;
cQ=fmod(cQ,360);cQ=cQ*pi/180; // conversion degré -> radian
cV=fmod(cV,360);cV=cV*pi/180;
cZ=fmod(cZ,360);cZ=cZ*pi/180;

// perturbation dans la longitude moyenne, en degrés:
cA = + (- 0.8142 + 0.0181*cu + 0.0167*cu*cu)*sin(cV)
+ (- 0.0105 + 0.1609*cu - 0.0041*cu*cu)*cos(cV)
- 0.1488*sin(cZ) - 0.0408*sin(2*cZ) - 0.0152*sin(3*cZ)
+ (0.0089*sin(cZ) - 0.0165*sin(2*cZ))*sin(cQ)
+ (0.0813*cos(cZ) + 0.0150*cos(2*cZ))*sin(cQ)
+ (0.0856*sin(cZ) + 0.0253*cos(cZ) + 0.0144*cos(2*cZ))*cos(cQ)
+ 0.0092*sin(2*cZ)*sin(2*cQ);
// perturbation dans l'excentricité, en unités de la sixième décimale:
cB = + (- 793 + 255*cu)*sin(cV) + (1338 + 123*cu)*cos(cV)
+ 1241*sin(cQ) + (39-62*cu)*sin(cZ)*sin(cQ)
+ (2660*cos(cZ) - 469*cos(2*cZ) - 187*cos(3*cZ) - 82*cos(4*cZ))*sin(cQ)
- (1270*sin(cZ) + 420*sin(2*cZ) +150*sin(3*cZ))*cos(cQ)
- 62*sin(4*cZ)*cos(cQ)
+ (221*sin(cZ) - 221*sin(2*cZ) - 57*sin(3*cZ))*sin(2*cQ)
- (278*cos(cZ) - 202*cos(2*cZ))*sin(2*cQ)
- (284*sin(cZ) + 159*cos(cZ))*cos(2*cQ)
+ (216*cos(2*cZ) + 56*cos(3*cZ))*cos(2*cQ);
// perturbation dans la longitude du périhélie, en degrés:
cC = + (0.0771 + 0.0072*cu)*sin(cV)
+ (0.0458 - 0.0148*cu)*cos(cV)
- (0.0758*sin(cZ) + 0.0248*sin(2*cZ) + 0.0086*sin(3*cZ))*sin(cQ)
- (0.0726 + 0.1504*cos(cZ) - 0.0269*cos(2*cZ) - 0.0101*cos(3*cZ))*cos(cQ)
- (0.0136*sin(cZ) - 0.0136*cos(2*cZ))*sin(2*cQ)
- (0.0137*sin(cZ) - 0.0120*sin(2*cZ))*cos(2*cQ)
+ (0.0149*cos(cZ) - 0.0131*cos(2*cZ))*cos(2*cQ);
// perturbation dans le demie-grand axe, en unités de la sixième décimale:
cD = + 2933*cos(cV) + 33629*cos(cZ) - 3081*cos(2*cZ)
- 1423*cos(3*cZ) - 671*cos(4*cZ)
+ (1098 - 2812*sin(cZ) + 688*sin(2*cZ))*sin(cQ)
+ (2138*cos(cZ) - 999*cos(2*cZ) - 642*cos(3*cZ))*sin(cQ)
- 890*cos(cQ)
+ (2206*sin(cZ) - 1590*sin(2*cZ) - 647*sin(3*cZ))*cos(cQ)
+ (2885*cos(cZ) + 2172*cos(2*cZ))*cos(cQ)
- 778*cos(cZ)*sin(2*cQ) - 856*sin(cZ)*cos(2*cQ);
L += cA;
M += cA - cC/e;
e += cB/1000000;
a += cD/1000000;
break;

case 7: // uranus:

cu = 0.1 + T/5;
cP = 237.475 + 3034.9061 *T; // angle en degrés
cQ = 265.916 + 1222.1139 *T; // angle en degrés

cH = 284.02 + 8.51*T; // angle en degrés
cS = 243.52 + 428.47*T; // angle en degrés
cZ = cS - cP; // angle en degrés
cN = cS - cQ; // angle en degrés
cT = 200.25 - 209.98*T; // angle en degrés

cP=fmod(cP,360);cP=cP*pi/180;
cQ=fmod(cQ,360);cQ=cQ*pi/180; // conversion degré -> radian
cH=fmod(cH,360);cH=cH*pi/180;
cT=fmod(cT,360);cT=cT*pi/180;
cS=fmod(cS,360);cS=cS*pi/180;
cZ=fmod(cZ,360);cZ=cZ*pi/180;
cN=fmod(cN,360);cN=cN*pi/180;

cA = 0.8643*sin(cH) + 0.0360*sin(2*cH) + (0.0822 - 0.0068*cu)*cos(cH);
cB = 2098*cos(cH) - 335*sin(cH) + 131*cos(2*cH);
cC = 0.1203*sin(cH) + 0.0195*cos(cH) + 0.0062*sin(2*cH);
cD = - 3825*cos(cH);

L += cA;
M += cA - cC/e;
e += cB/1000000;
a += cD/1000000;
break;

case 8: // neptune:

cH = 284.02 + 8.51*T; // angle en degrés
cT = 200.25 - 209.98*T; // angle en degrés
cZ = 153.71 + 2816.42*T; // angle en degrés
cN = 182.15 + 1003.62*T; // angle en degrés

cH=fmod(cH,360);cH=cH*pi/180;
cT=fmod(cT,360);cT=cT*pi/180; // conversion degré -> radian
cZ=fmod(cZ,360);cZ=cZ*pi/180;
cN=fmod(cN,360);cN=cN*pi/180;

cA = - 0.5898*sin(cH) - 0.0561*cos(cH) - 0.0243*sin(2*cH);
cB = + 439*sin(cH) + 426*cos(cH) + 113*sin(2*cH) + 109*cos(2*cH);
cC = + 0.0240*sin(cH) - 0.0253*cos(cH);
cD = - 817*sin(cH) + 8189*cos(cH) + 781*cos(2*cH);

L += cA;
M += cA - cC/e;
e += cB/1000000;
a += cD/1000000;
break;
}

L=fmod(L,360);L=L*pi/180;
i=fmod(i,360);i=i*pi/180; // conversion degré -> radian
w=fmod(w,360);w=w*pi/180;
W=fmod(W,360);W=W*pi/180;
M=fmod(M,360);M=M*pi/180;

E=M,erreur=1;
//while (erreur>1e-8 ) {ni++;E1=M+e*sin(E);erreur=fabs(E1-E);E=E1;}

while (erreur>1e-9) {erreur = (M+e*sin(E)-E)/(1-e*cos(E));E += erreur;}
v=2*atan(sqrt((1+e)/(1-e))*tan(E/2));
r=a*(1-e*cos(E));
u=L+v-M-W;
l = W + atan2(cos(i)*sin(u),cos(u));
b=asin(sin(u)*sin(i));

if (k==7) // égal à uranus
{ l = l + pi/180 *
( + (0.0101 - 0.0010*cu)*sin(cS+cN)
- (0.0386 - 0.0020*cu)*cos(cS+cN)
+ (0.0350 - 0.0010*cu)*cos(2*cS+cN)
- 0.0148*sin(cZ) + 0.0099*sin(cT) + 0.0088*sin(2*cT) );
r = r - 0.02595 + 0.00498*cos(cZ)
- 0.00123*cos(cS) + 0.00335*cos(cN)
+ (0.00579*cos(cS) - 0.00116*sin(cS) + 0.00139*cos(2*cS))*sin(cN)
+ (0.00135*cos(cS) + 0.00570*sin(cS) + 0.00139*sin(2*cS))*cos(cN)
+ 0.00090*cos(2*cT) + 0.00089*(cos(cT) - cos(3*cT));
}
if (k==8) // égal à neptune
{ l = l + pi/180 * ( 0.0096*sin(cZ) + 0.0052*sin(cN) );
r = r - 0.04060 + 0.00499*cos(cZ) + 0.00274*cos(cN)
+ 0.00204*cos(cT) + 0.00105*cos(2*cT);
}

//Lg et R et o

lambda=Lg+atan2(r*cos(b)*sin(l-Lg), r*cos(b)*cos(l-Lg) + R);
delta=sqrt(R*R+r*r+2*r*R*cos(b)*cos(l-Lg));
beta=asin(r/delta*sin(b));

A=atan2(sin(lambda)*cos(o)-tan(beta)*sin(o), cos(lambda));
D=asin(sin(beta)*cos(o)+cos(beta)*sin(o)*sin(lambda));
p.Asc=A; p.Dec=D; p.Dist=delta;

}
}
else
{
double Lp,M,Mp,D,F,W,e,lambda,B,w1,w2,beta,g,t,d,i;

Lp= 270.434164 + 481267.8831 *T - 0.001133 *T*T + 0.0000019 *T*T*T;
M= 358.475833 + 35999.0498 *T - 0.000150 *T*T - 0.0000033 *T*T*T;
Mp= 296.104608 + 477198.8491 *T + 0.009192 *T*T + 0.0000144 *T*T*T;
D= 350.737486 + 445267.1142 *T - 0.001436 *T*T + 0.0000019 *T*T*T;
F= 11.250889 + 483202.0251 *T - 0.003211 *T*T - 0.0000003 *T*T*T;
W= 259.183275 - 1934.1420 *T + 0.002078 *T*T + 0.0000022 *T*T*T;

Lp=fmod(Lp,360);
M=fmod(M,360);
Mp=fmod(Mp,360);
D=fmod(D,360);
F=fmod(F,360);
W=fmod(W,360);

W=W*pi/180; // conversion en radian de W !!!

Lp += + 0.000233 * sin((51.2+20.2*T)*pi/180);
M += - 0.001778 * sin((51.2+20.2*T)*pi/180);
Mp += + 0.000817 * sin((51.2+20.2*T)*pi/180);
D += + 0.002011 * sin((51.2+20.2*T)*pi/180);

t= 0.003964 * sin((346.560+132.870*T-0.0091731*T*T)*pi/180); Lp+=t; Mp+=t; D+=t; F+=t;

Lp += + 0.001964 *sin(W);
Mp += + 0.002541 *sin(W);
D += + 0.001964 *sin(W);
F += - 0.024691 *sin(W);
F += - 0.004328 *sin(W+(275.05-2.30*T)*pi/180);

Mp=Mp*pi/180;
M=M*pi/180;
D=D*pi/180;
F=F*pi/180;

e= 1 - 0.002495*T - 0.00000752*T*T;

lambda = Lp + 6.288750 *sin(Mp)
+ 1.274018 *sin(2*D-Mp)
+ 0.658309 *sin(2*D)
+ 0.213616 *sin(2*Mp)
- 0.185596 *sin(M) *e
- 0.114336 *sin(2*F)
+ 0.058793 *sin(2*D-2*Mp)
+ 0.057212 *sin(2*D-M-Mp) *e
+ 0.053320 *sin(2*D+Mp)
+ 0.045874 *sin(2*D-M) *e
+ 0.041024 *sin(Mp-M) *e
- 0.034718 *sin(D)
- 0.030465 *sin(M+Mp) *e
+ 0.015326 *sin(2*D-2*F)
- 0.012528 *sin(2*F+Mp)
- 0.010980 *sin(2*F-Mp)
+ 0.010674 *sin(4*D-Mp)
+ 0.010034 *sin(3*Mp)
+ 0.008548 *sin(4*D-2*Mp)
- 0.007910 *sin(M-Mp+2*D) *e
- 0.006783 *sin(2*D+M) *e
+ 0.005162 *sin(Mp-D)
+ 0.005000 *sin(M+D) *e
+ 0.004049 *sin(Mp-M+2*D) *e
+ 0.003996 *sin(2*Mp+2*D)
+ 0.003862 *sin(4*D)
+ 0.003665 *sin(2*D-3*Mp)
+ 0.002695 *sin(2*Mp-M) *e
+ 0.002602 *sin(Mp-2*F-2*D)
+ 0.002396 *sin(2*D-M-2*Mp) *e
- 0.002349 *sin(Mp+D)
+ 0.002249 *sin(2*D-2*M) *e*e
- 0.002125 *sin(2*Mp+M) *e
- 0.002079 *sin(2*M) *e*e
+ 0.002059 *sin(2*D-Mp-2*M) *e*e
- 0.001773 *sin(Mp+2*D-2*F)
- 0.001595 *sin(2*F+2*D)
+ 0.001220 *sin(4*D-M-Mp) *e
- 0.001110 *sin(2*Mp+2*F)
+ 0.000892 *sin(Mp-3*D)
- 0.000811 *sin(M+Mp+2*D) *e
+ 0.000761 *sin(4*D-M-2*Mp) *e
+ 0.000717 *sin(Mp-2*M) *e*e
+ 0.000704 *sin(Mp-2*M-2*D) *e*e
+ 0.000693 *sin(M-2*Mp+2*D) *e
+ 0.000598 *sin(2*D-M-2*F) *e
+ 0.000550 *sin(Mp+4*D)
+ 0.000538 *sin(4*Mp)
+ 0.000521 *sin(4*D-M) *e
+ 0.000486 *sin(2*Mp-D);

lambda=lambda*pi/180;

B = + 5.128189 *sin(F)
+ 0.280606 *sin(Mp+F)
+ 0.277693 *sin(Mp-F)
+ 0.173238 *sin(2*D-F)
+ 0.055413 *sin(2*D+F-Mp)
+ 0.046272 *sin(2*D-F-Mp)
+ 0.032573 *sin(2*D+F)
+ 0.017198 *sin(2*Mp+F)
+ 0.009267 *sin(2*D+Mp-F)
+ 0.008823 *sin(2*Mp-F)
+ 0.008247 *sin(2*D-M-F) *e
+ 0.004323 *sin(2*D-F-2*Mp)
+ 0.004200 *sin(2*D+F+Mp)
+ 0.003372 *sin(F-M-2*D) *e
+ 0.002472 *sin(2*D+F-M-Mp) *e
+ 0.002222 *sin(2*D+F-M) *e
+ 0.002072 *sin(2*D-F-M-Mp) *e
+ 0.001877 *sin(F-M+Mp) *e
+ 0.001828 *sin(4*D-F-Mp)
- 0.001803 *sin(F+M) *e
- 0.001750 *sin(3*F)
+ 0.001570 *sin(Mp-M-F) *e
- 0.001487 *sin(F+D)
- 0.001481 *sin(F+M+Mp) *e
+ 0.001417 *sin(F-M-Mp) *e
+ 0.001350 *sin(F-M) *e
+ 0.001330 *sin(F-D)
+ 0.001106 *sin(F+3*Mp)
+ 0.001020 *sin(4*D-F)
+ 0.000833 *sin(F+4*D-Mp)
+ 0.000781 *sin(Mp-3*F)
+ 0.000670 *sin(F+4*D-2*Mp)
+ 0.000606 *sin(2*D-3*F)
+ 0.000597 *sin(2*D+2*Mp-F)
+ 0.000492 *sin(2*D+Mp-M-F) *e
+ 0.000450 *sin(2*Mp-F-2*D)
+ 0.000439 *sin(3*Mp-F)
+ 0.000423 *sin(F+2*D+2*Mp)
+ 0.000422 *sin(2*D-F-3*Mp)
- 0.000367 *sin(M+F+2*D-Mp) *e
- 0.000353 *sin(M+F+2*D) *e
+ 0.000331 *sin(F+4*D)
+ 0.000317 *sin(2*D+F-M+Mp) *e
+ 0.000306 *sin(2*D-2*M-F) *e*e
- 0.000283 *sin(Mp+3*F);

w1=0.0004664*cos(W);
w2=0.0000754*cos(W+(275.05-2.30*T)*pi/180);
beta=B*(1-w1-w2);
beta=beta*pi/180;

g = + 0.950724
+ 0.051818 *cos(Mp)
+ 0.009531 *cos(2*D-Mp)
+ 0.007843 *cos(2*D)
+ 0.002824 *cos(2*Mp)
+ 0.000857 *cos(2*D+Mp)
+ 0.000533 *cos(2*D-M) *e
+ 0.000401 *cos(2*D-M-Mp) *e
+ 0.000320 *cos(Mp-M) *e
- 0.000271 *cos(D)
- 0.000264 *cos(M+Mp) *e
- 0.000198 *cos(2*F-Mp)
+ 0.000173 *cos(3*Mp)
+ 0.000167 *cos(4*D-Mp)
- 0.000111 *cos(M) *e
+ 0.000103 *cos(4*D-2*Mp)
- 0.000084 *cos(2*Mp-2*D)
- 0.000083 *cos(2*D+M) *e
+ 0.000079 *cos(2*D+2*Mp)
+ 0.000072 *cos(4*D)
+ 0.000064 *cos(2*D-M+Mp) *e
- 0.000063 *cos(2*D+M-Mp) *e
+ 0.000041 *cos(M+D) *e
+ 0.000035 *cos(2*Mp-M) *e
- 0.000033 *cos(3*Mp-2*D)
- 0.000030 *cos(Mp+D)
- 0.000029 *cos(2*F-2*D)
- 0.000029 *cos(2*Mp+M) *e
+ 0.000026 *cos(2*D-2*M) *e*e
- 0.000023 *cos(2*F-2*D+Mp)
+ 0.000019 *cos(4*D-M-Mp) *e;

g=g*pi/180;

p.Dist=6378.14/sin(g); // en kilomètres.
p.Dist/=1000; // en milier de kilomètres.
p.Asc=atan2(sin(lambda)*cos(o)-tan(beta)*sin(o), cos(lambda));
p.Dec=asin(sin(beta)*cos(o)+cos(beta)*sin(o)*sin(lambda));

d=acos(cos(lambda-Lg)*cos(beta))*180/pi;
i=180-d-0.1468*sin(d*pi/180)*(1-0.0549*sin(Mp))/(1-0.0167*sin(M));
p.k=(1+cos(i*pi/180))/2;

}
return p;
}

//--------------------------------------------------
int main(int argc, char* argv[])

{
cout<<"Version 2.2"<<endl<<endl;

char c;
const int Lignemax=20;
char ch[Lignemax+1];
//RgnHandle cursorRgn;
//EventRecord event;

cout<<setprecision(5);
cout<<"Jour : "<<flush;cin>>dat.jour;//if (! cin) {cin.clear();cin.putback('1');}cin.getline(ch, Lignemax);
cout<<"Mois : "<<flush;cin>>dat.mois;//if (! cin) {cin.clear();cin.putback('1');}cin.getline(ch, Lignemax);
cout<<"Annee : "<<flush;cin>>dat.annee;//if (! cin) {cin.clear();cin.putback('1');}cin.getline(ch, Lignemax);
aff();

do { cout<<"Jour : "<<flush;cin>>dat.jour;//if (! cin) {cin.clear();cin.putback('1');}cin.getline(ch, Lignemax);
cout<<"Mois : "<<flush;cin>>dat.mois;//if (! cin) {cin.clear();cin.putback('1');}cin.getline(ch, Lignemax);
cout<<"Annee : "<<flush;cin>>dat.annee;//if (! cin) {cin.clear();cin.putback('1');}cin.getline(ch, Lignemax);
aff();

} while (1);
return 0;
}

void aff()
{

position p = CalculPosition((planete)0,dat); // 0 = Soleil

cout<<dat.jour<<"/"<<dat.mois<<"/"<<dat.annee<<" (Jour Julien : "<<setprecision(11)<<JJ(dat)<<") (Temps siderale : ";
cout.setf(ios::right);
affhms(p.t);
cout<<")"<<endl;
cout<<"------------------------------------------------------------------"<<endl;
cout<<setprecision(5);

cout<<"Planetes";
cout<<" "<<"Ascension droite";
cout<<" "<<"Declinaison";
cout<<" "<<"Distance";
cout<<" "<<"Eclairement"<<endl;

cout<<"------------------------------------------------------------------"<<endl;
for (int k=0;k<9;k++)
{
cout<<setw(8);aff_nom_planete (k); // k est une planete
position p = CalculPosition(k,dat); // k est une planete

cout <<" ";affhms(p.Asc);
cout <<" ";affds(p.Dec);
cout <<" ";
cout<<setw(6)<<floor(p.Dist*10000)/10000;
if (k==1) {cout<<" "<<floor(100*p.k)<<"%";}
cout<<endl;

}
cout<<endl;
}

void affhms(double t)
{ while (t<0) t=t+2*pi;
while (t>2*pi) t=t-2*pi;
t=t*12/pi;
cout<<setw(2);
cout<<floor(t)<<"h "<<flush;
t=60*fmod(t,1);
cout<<setw(2)<<floor(t)<<"m ";
t=60*fmod(t,1);
cout<<setw(2)<<floor(t)<<"s"<<flush;
}

void affds(double t)
{ while (t<0) t=t+2*pi;
while (t>pi) t=t-2*pi;
t=t*180/pi;
if (t<0) cout<<"-";else cout<<" ";
cout<<setw(2)<<floor(fabs(t))<<" deg ";
t=60*fmod(t,1);
cout<<setw(2)<<floor(fabs(t))<<"'"<<flush;
}