/*************************************************************************************************************************************/
/* */
/* R E S O L U T I O N D E L ' E Q U A T I O N D E K E P L E R */
/* P E R M E T T A N T A U N C O R P S C E L E S T E D E P A R C O U R I R */
/* U N E T R A J E C T O I R E E L L I P T I Q U E */
/* A V E C A V A N C E D E P E R I H E L I E :E */
/* */
/* */
/* Author of '$xtc/EquKepler.12$c' : */
/* */
/* Jean-Francois COLONNA (LACTAMME, 20260204184424). */
/* */
/*************************************************************************************************************************************/
#include "INCLUDES.01.I"
#define ITERATIONS 18000
#define EPSILON 0.000001
#define TEMPS_0 0.0
#define PAS_TEMPS 0.02
#define THETA_CIRCULAIRE_0 0.0
#define PAS_THETA_CIRCULAIRE 0.1
#define AVANCE_PERIHELIE_0 0
#define PAS_AVANCE_PERIHELIE 0.0002
#define DEMI_GRAND_AXE 200.0
#define DEMI_PETIT_AXE 100.0
#define EXCENTRICITE (DISTANCE_FOCALE/DEMI_GRAND_AXE)
#define DISTANCE_FOCALE sqrt(EXP2(DEMI_GRAND_AXE)-EXP2(DEMI_PETIT_AXE))
#define PARAMETRE_FOCAL EXP2(DEMI_PETIT_AXE)/DEMI_GRAND_AXE
extern double cos();
extern double sin();
extern double acos();
extern double sqrt();
extern void *malloc();
static int dimX=0;
#define Xmin 0
#define Xmax (Xmin + (dimX-1))
/* Definition des abscisses. */
static int dimY=0;
#define Ymin 0
#define Ymax (Ymin + (dimY-1))
/* Definition des ordonnees. */
#define IMAGE(x,y) \
(*(image + (((y-Ymin)*dimX) + (x-Xmin)))) \
/* Acces a un point de l'image. */
#define store(n,x,y) \
{ \
if (((x>=Xmin) && (x<=Xmax)) && ((y>=Ymin) && (y<=Ymax))) \
{ \
IMAGE(x,y) = n; \
} \
else \
{ \
} \
} \
/* Rangement d'un point valide d'une image. */
main()
{
int x,y;
/* Definition des coordonnees. */
unsigned char *image;
/* Definition de l'image a generer... */
int n;
double temps=TEMPS_0;
double rho_circulaire=DEMI_GRAND_AXE;
double theta_circulaire=THETA_CIRCULAIRE_0;
/* Angle 'theta' de parcours du cercle circonscrit a l'ellipse. */
double avance_perihelie=AVANCE_PERIHELIE_0;
/* Avance de perihelie initiale. */
Get(dimX,"dimX");
Get(dimY,"dimY");
/* Recuperation des dimensions en 'X' et en 'Y' de l'image a generer. */
image=(unsigned char*)malloc(dimY*dimX);
/* Definition de l'image a generer... */
for (y=Ymin ; y<=Ymax ; y++)
{
for (x=Xmin ; x<=Xmax ; x++)
{
store(NOIR,x,y);
/* Initialisation de l'image finale... */
}
}
store(BLANC
,(int)(dimX/2)
,(int)(dimY/2)
);
/* Trace de l'un des 2 foyers qui est le centre de rotation... */
for (n=1 ; n<=ITERATIONS ; n++)
{
double pas_theta=PAS_THETA_CIRCULAIRE;
int iterer=0;
while (iterer == 0)
{
double formule=theta_circulaire - (EXCENTRICITE*sin(theta_circulaire));
/* Formule d'equivalence entre le temps (de parcours de l'ellipse) et le 'theta' du cercle : */
/* */
/* thetaC - e.sin(thetaC) = temps */
/* */
/* d'ou : */
/* */
/* thetaC = fonction1(temps). */
/* */
if (ABSO(formule-temps) < EPSILON)
{
double numerateur,denominateur;
double rho_elliptique,theta_elliptique;
/* {rhoE,thetaE} de parcours de l'ellipse, l'origine etant alors le foyer de droite de */
/* l'ellipse. */
int correction;
/* Afin que 'thetaE' puisse aller au-dela de pi... */
double x,y;
/* {x,y} du point courant de l'ellipse, l'origine etant le centre de l'ellipse (et du */
/* cercle)... */
double xA,yA;
/* {xA,yA} du point courant de l'ellipse avec avance de perihelie... */
numerateur=(rho_circulaire*cos(theta_circulaire))-DISTANCE_FOCALE;
denominateur=((PARAMETRE_FOCAL+(DISTANCE_FOCALE*EXCENTRICITE))
-(rho_circulaire*EXCENTRICITE*cos(theta_circulaire))
);
/* On notera le 20260205214030 que dans 'v $xrr/AvanceDePerihelie.11$K denominateur' on */
/* a retrouve, lors de sa mise au point (contenant a la fois ce 'c' et le '$K' de */
/* 'v $xrk/lorenz.11$K'), cette expression de la forme : */
/* */
/* denominateur=(A */
/* -B */
/* ); */
/* */
/* sous cette forme, on trouve : */
/* */
/* denominateur == A */
/* */
/* alors que la forme : */
/* */
/* denominateur=(A-B); */
/* */
/* donne correctement : */
/* */
/* denominateur == A-B */
/* */
/* Mystere... */
/* On a : */
/* */
/* p */
/* rhoE = ------------------- */
/* 1 + e.cos(thetaE) */
/* */
/* (l'origine etant le foyer de droite de l'ellipse) d'ou : */
/* */
/* x = a.cos(thetaC) = rhoE.cos(thetaE) + d */
/* */
/* (l'origine etant le centre de l'ellipse -et du cercle-) d'ou : */
/* */
/* thetaE = fonction2(thetaC). */
/* */
/* et enfin, les coordonnees {x,y} du point courant de l'ellipse : */
/* */
/* x = rhoE.cos(thetaE) + d */
/* y = rhoE.sin(thetaE) */
/* */
/* (l'origine etant le centre de l'ellipse -et du cercle-) avec : */
/* */
/* a = DEMI_GRAND_AXE */
/* b = DEMI_PETIT_AXE */
/* d = DISTANCE_FOCALE */
/* e = EXCENTRICITE */
/* p = PARAMETRE_FOCAL */
/* */
if (ABSO(numerateur/denominateur) <= 1.0)
{
theta_elliptique = acos(numerateur/denominateur);
}
else
{
if ((numerateur/denominateur) >= 0.0)
{
theta_elliptique = acos(+1);
}
else
{
theta_elliptique = acos(-1);
}
}
correction=theta_circulaire/PI;
/* En effet, 'acos(...)' renvoie des valeurs dans [0,PI]. Il convient donc de corriger */
/* 'thetaE' en fonction de 'thetaC'... */
if (correction >= 0)
{
if (EST_PAIR(correction))
{
theta_elliptique=((correction+0)*PI) + theta_elliptique;
}
else
{
theta_elliptique=((correction+1)*PI) - theta_elliptique;
}
}
else
{
}
rho_elliptique=PARAMETRE_FOCAL/(1+(EXCENTRICITE*cos(theta_elliptique)));
x=(rho_elliptique*cos(theta_elliptique));
y=(rho_elliptique*sin(theta_elliptique));
/* Calcul du point courant de l'ellipse... */
xA=((x*cos(avance_perihelie))-(y*sin(avance_perihelie)));
yA=((x*sin(avance_perihelie))+(y*cos(avance_perihelie)));
/* Calcul du point courant de l'ellipse avec avance de perihelie... */
store(BLANC
,(int)(xA)+(dimX/2)
,(int)(yA)+(dimY/2)
);
/* Trace de l'ellipse... */
avance_perihelie=avance_perihelie+PAS_AVANCE_PERIHELIE;
iterer++;
}
else
{
if (formule > temps)
{
pas_theta = pas_theta/2.0;
theta_circulaire = theta_circulaire-pas_theta;
}
else
{
theta_circulaire = theta_circulaire+pas_theta;
}
}
}
temps = temps + PAS_TEMPS;
}
write(1,image,dimX*dimY);
/* Sortie de l'image... */
}