Calcul des positions heliocentriques

[blotfib]

Bonsoir

J'ai trouvé les formules pour calculer la position du soleil par rapport à la terre (ci dessous Calchaz() ), mais comment les modifier pour calculer la position de la terre par rapport au soleil et les autres planetes aussi % au soleil

Merci

//---------------------------------------------------------------------------
// Julian century depuis 01/01/2000 à 12h UT
double JulianCentury (int date, int month, int year, double UT)
{
	if (month<=2) {month=month+12; year=year-1;}
  return ((int)(365.25*year) + (int)(30.6001*(month+1)) - 15 + 1720996.5 + date + UT/24.0 - 2451545.0)/36525;
} 
//---------------------------------------------------------------------------
//Solar Coordinates (according to: Jean Meeus: Astronomical Algorithms), accuracy of 0.01 degree 
double SoleilLongitude( double T)   // L
{
double M, L0, DL, L, Lm;
	M = 357.52910 + 35999.05030*T - 0.0001559*T*T - 0.00000048*T*T*T; // mean anomaly, degree
	L0 = 280.46645 + 36000.76983*T + 0.0003032*T*T; // mean longitude, degree 
	DL = (1.914600 - 0.004817*T - 0.000014*T*T)*sin(deg2rd*M)
				+ (0.019993 - 0.000101*T)*sin(deg2rd*2*M) + 0.000290*sin(deg2rd*3*M); 
	L=L0 + DL; // true longitude, deg sans modulo
	Lm=L-360*(int)(L/360);
	if (Lm<0)	Lm+=360;
	return Lm*deg2rd; // true longitude, rd 
}
//---------------------------------------------------------------------------
// convert ecliptic longitude L (rd) to right ascension RA and declination delta
void eclipticL2RA_delta( double T, double L, double * RA, double * delta)
{
double eps, X, Y, Z, R;
// old approxim  eps = deg2rd*23.43999; // obliquity of ecliptic
	eps = deg2rd*(23.4393 -0.013*T);	// obliquity of ecliptic
	X = cos(L); Y = cos(eps)*sin(L); Z = sin(eps)*sin(L); R = sqrt(1.0-Z*Z); 
	*delta = atan2(Z, R); // in rd 
	*RA = 2*atan2(Y, X+R); //  in rd 
}
//---------------------------------------------------------------------------
//compute sidereal time at Greenwich (according to: Jean Meeus: Astronomical Algorithms ) 
double TempsSideral( double T)   // theta0
{
	double theta0, theta0m;
//	T siecle julien depuis 2000
	theta0=280.46061837 + 360.98564736629*T*36525 + 0.000387933*T*T - T*T*T/38710000.0;
	theta0m=theta0-360*(int)(theta0/360);
	if (theta0m<0)	theta0m+=360;
	return theta0m*deg2rd;    // TempsSideralen rd
}
//---------------------------------------------------------------------------
//convert tau, delta to horizon coordinates of the observer (altitude h, azimuth az) , needs latitude : beta
void equatorial2horizontal( double beta, double tau, double delta, double * h, double * az)
{
	*h=asin( sin(beta)*sin(delta) + cos(beta)*cos(delta)*cos(tau));	// h en rd
	*az =atan2( -sin(tau), cos(beta)*tan(delta) - sin(beta)*cos(tau))-M_PI;	// az en rd du SUD
	if (*az<0)	*az+=2*M_PI;
}
//---------------------------------------------------------------------------
void Calchaz()
{
bool DSTApplies, Verbose=0;
double TimeZone;
double T, theta0, theta;  // T : Julien Century , theta0 : temps sideral à Greenwich
char SS[SSlenght];

	DSTApplies=(DST !=0) && CalcTimeZone(date.a, date.m, date.j, DST, TimeZoneTown);
	if (DSTApplies) TimeZone=1+TimeZoneTown;
						else  TimeZone=TimeZoneTown;
	date.h=heurelocale-TimeZone;	// heure GMT
	if (date.h<0) {date.h+=24; date.j--;}
	 if (date.j<0) printf("### erreur j<0 en rattrapant h<0");//{date.j+=???; date.m--}
	T=JulianCentury (date.j, date.m, date.a, date.h);
	if (Verbose)	printf("T=%g JC\n",T);
	Soleil.L=SoleilLongitude(T);
	if (Verbose) printf("Soleil Longitude L=%g° ie : %s\n", Soleil.L*rd2deg, rd2hminsec(Soleil.L, SS));
  eclipticL2RA_delta( T, Soleil.L, & Soleil.RA, & Soleil.delta);
  theta0=TempsSideral( T);
	theta=theta0+longitude;
	if (Verbose) printf("theta0=%g° ie : %s\n", theta0*rd2deg, rd2hminsec(theta0, SS));
	if (Verbose) printf("theta=%g°  ie heure sidérale %s\n", theta*rd2deg, rd2hminsec(theta, SS));
	Soleil.tau=theta-Soleil.RA;
	if (Verbose) printf("RA=%g° ie : %s\n", Soleil.RA*rd2deg, rd2hminsec(Soleil.RA, SS));
	if (Verbose) printf("delta=%g°  ie DE %s\n", Soleil.delta*rd2deg, rd2degminsec(Soleil.delta, SS));
	if (Verbose) printf("tau=%g°  ie angle horaire %s\n", Soleil.tau*rd2deg, rd2hminsec(Soleil.tau, SS));
	equatorial2horizontal( latitude, Soleil.tau, Soleil.delta, & Soleil.h, & Soleil.az);
	if (Verbose) printf("az=%g° (orgine Sud) ie : azimut %s\n", Soleil.az*rd2deg, rd2degminsec(Soleil.az, SS));
	if (Verbose) printf("------az=%g° (orgine Nord) ie -180 : azimut %s\n", Soleil.az*rd2deg-180, rd2degminsec(Soleil.az-M_PI, SS));
	if (Verbose) printf("h=%g°  ie altitude %s\n", Soleil.h*rd2deg, rd2degminsec(Soleil.h, SS));
}

à l'aide de :
http://www.stjarnhimlen.se/comp/tutorial.html
http://stjarnhimlen.se/comp/ppcomp.html#5
http://www.geoastro.de/elevaz/basics/index.htm
-----