/*************************************************************************************************************************************/
/*                                                                                                                                   */
/*        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  :                                                                     */
/*                                                                                                                                   */
/*                                                                                                                                   */
/*        Author of '$xtc/EquKepler.01$c' :                                                                                          */
/*                                                                                                                                   */
/*                    Jean-Francois COLONNA (LACTAMME, 20051102093745).                                                              */
/*                                                                                                                                   */
/*************************************************************************************************************************************/

#include  "INCLUDES.01.I"
                                        /* Introduit le 20051109112929...                                                            */

#define   ITERATIONS          140
#define   EPSILON             0.000001

#define   TEMPS_0             0.0
#define   PAS_TEMPS           0.1

#define   THETA_CIRCULAIRE_0  0.0
#define   PAS_THETA_CIRCULAIRE 0.1

#define   DEMI_GRAND_AXE      1.0
#define   DEMI_PETIT_AXE      0.5
#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();

main()
          {
          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.                              */

          printf("a=%f\n",DEMI_GRAND_AXE);
          printf("b=%f\n",DEMI_PETIT_AXE);
          printf("e=%f\n",EXCENTRICITE);
          printf("p=%f\n",PARAMETRE_FOCAL);
          printf("\n");

          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)...                                                                                */

                                        numerateur=(rho_circulaire*cos(theta_circulaire))-DISTANCE_FOCALE;
                                        denominateur=((PARAMETRE_FOCAL+(DISTANCE_FOCALE*EXCENTRICITE))
                                                     -(rho_circulaire*EXCENTRICITE*cos(theta_circulaire))
                                                      );
                                        /* 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))+DISTANCE_FOCALE;
                                        y=(rho_elliptique*sin(theta_elliptique));
                                        /* Calcul du point courant de l'ellipse...                                                   */

                                        printf("temps=%f (formule=%f)     thetaC=%f    thetaE=%f  point={%f,%f}\n"
                                              ,temps
                                              ,formule
                                              ,theta_circulaire
                                              ,theta_elliptique
                                              ,x,y
                                               );

                                        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;
                    }
          }