/*************************************************************************************************************************************/
/*                                                                                                                                   */
/*        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.02$c' :                                                                                          */
/*                                                                                                                                   */
/*                    Jean-Francois Colonna (LACTAMME, 20051107105318).                                                              */
/*                                                                                                                                   */
/*************************************************************************************************************************************/

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

#define   ITERATIONS          64
#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      300
#define   DEMI_PETIT_AXE      50
#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.                                                                 */

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

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

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

                                        store(BLANC
                                             ,(int)(x)+(dimX/2)
                                             ,(int)(y)+(dimY/2)
                                              );
                                        /* Trace de l'ellipse...                                                                     */

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



Copyright © Jean-François Colonna, 2021-2021.
Copyright © CMAP (Centre de Mathématiques APpliquées) UMR CNRS 7641 / Ecole Polytechnique, 2021-2021.