/*************************************************************************************************************************************/
/*                                                                                                                                   */
/*        T E S T   D E   L A   C O N S T R U C T I O N   D ' U N E   C O U R B E   D E   V O N   K O C H                            */
/*        D A N S   L E   C A D R E   D U   M O O C   S U R   L A   G E O M E T R I E   F R A C T A L E  :                           */
/*                                                                                                                                   */
/*                                                                                                                                   */
/*        Author of '$xtc/VonKoch.11$c' :                                                                                            */
/*                                                                                                                                   */
/*                    Jean-Francois Colonna (LACTAMME, AAAAMMJJhhmmss).                                                              */
/*                                                                                                                                   */
/*************************************************************************************************************************************/

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

extern    void      *malloc();

extern    double    sqrt();

extern    double    cos();
extern    double    sin();

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   Xnormalise(x)                                                                                                                 \
                    ((int)((x-XG)/(XD-XG)*dimX))
#define   Ynormalise(y)                                                                                                                 \
                    ((int)((y-YB)/(YH-YB)*dimY))
#define   IMAGE(x,y)                                                                                                                    \
                    (*(image + ((int)(((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.                                                  */

#define   XG        (0.00)
#define   XD        (1.00)
                                        /* Extremites Gauche et Droite sur l'axe des 'X'.                                            */
#define   YB        (0.00)
#define   YH        (1.00)
                                        /* Extremites Bas et Haut sur l'axe des 'Y'.                                                 */

typedef   struct    complexe  {
                              double    reelle;
                              double    imaginaire;
                              }         complexe;
                                        /* 'struct' de definition d'un nombre complexe...                                            */
#define   Reelle(z)           (z.reelle)
#define   Imaginaire(z)       (z.imaginaire)

#define   Cinit(z,x1,y1)                                                                                                                \
                    {                                                                                                                   \
                    Reelle(z) = x1;                                                                                                     \
                    Imaginaire(z) = y1;                                                                                                 \
                    }
#define   Cegal(z,z1)                                                                                                                   \
                    {                                                                                                                   \
                    Cinit(z,Reelle(z1),Imaginaire(z1));                                                                                 \
                    }
#define   Csomme(z,z1,z2)                                                                                                               \
                    {                                                                                                                   \
                    complexe  zT;                                                                                                       \
                    Reelle(zT) = Reelle(z1) + Reelle(z2);                                                                               \
                    Imaginaire(zT) = Imaginaire(z1) + Imaginaire(z2);                                                                   \
                    Cegal(z,zT);                                                                                                        \
                    }
#define   Cdifference(z,z1,z2)                                                                                                          \
                    {                                                                                                                   \
                    complexe  zT;                                                                                                       \
                    Reelle(zT) = Reelle(z1) - Reelle(z2);                                                                               \
                    Imaginaire(zT) = Imaginaire(z1) - Imaginaire(z2);                                                                   \
                    Cegal(z,zT);                                                                                                        \
                    }
#define   Cproduit(z,z1,z2)                                                                                                             \
                    {                                                                                                                   \
                    complexe  zT;                                                                                                       \
                    Reelle(zT) = (Reelle(z1)*Reelle(z2)) - (Imaginaire(z1)*Imaginaire(z2));                                             \
                    Imaginaire(zT) = (Reelle(z2)*Imaginaire(z1)) + (Reelle(z1)*Imaginaire(z2));                                         \
                    Cegal(z,zT);                                                                                                        \
                    }
#define   Crotation(z,z1,rho,theta)                                                                                                     \
                    {                                                                                                                   \
                    complexe  zR;                                                                                                       \
                    Cinit(zR,rho*cos(theta),rho*sin(theta));                                                                            \
                    Cproduit(z,z1,zR);                                                                                                  \
                    }
#define   Cmodule2(z)                                                                                                                   \
                    ((Reelle(z)*Reelle(z)) + (Imaginaire(z)*Imaginaire(z)))

unsigned  char      *image;
                                        /* Definition de l'image a generer...                                                        */

double    seuil=0.002;
                                        /* Definition du seuil d'arret.                                                              */
double    angle=PI/3.0;
double    rapport=3.0;
#define   HAUTEUR                                                                                                                       \
                    sin(angle)/rapport
#define   MILIEU                                                                                                                        \
                    2
                                        /* Definition de la transformation...                                                        */

construction(Az,Ez)
complexe  Az,Ez;
          {
          double    Ax=Reelle(Az),Ay=Imaginaire(Az);
          double    Ex=Reelle(Ez),Ey=Imaginaire(Ez);

          complexe  difference;

          Cdifference(difference,Ez,Az);

          if        (sqrt(Cmodule2(difference)) < seuil)
                    {
                    store(BLANC,Xnormalise(Reelle(Az)),Ynormalise(Imaginaire(Az)));
                    store(BLANC,Xnormalise(Reelle(Ez)),Ynormalise(Imaginaire(Ez)));
                                        /* Marquage et arret de la recursivite...                                                    */
                    }
          else
                    {
                    double    Bx,By;
                    complexe  Bz;

                    double    Cx,Cy;
                    complexe  Cz;

                    double    Hx=(Ex+Ax)/MILIEU,Hy=(Ey+Ay)/MILIEU;
                    double    AEx_reduit=(Ex-Ax)/rapport,AEy_reduit=(Ey-Ay)/rapport;
                    double    HCx=(Ey-Ay)*HAUTEUR,HCy=(Ex-Ax)*HAUTEUR;

                    double    Dx,Dy;
                    complexe  Dz;

                    Bx = Ax + AEx_reduit;
                    By = Ay + AEy_reduit;

                    Cx = Hx - HCx;
                    Cy = Hy + HCy;

                    Dx = Ex - AEx_reduit;
                    Dy = Ey - AEy_reduit;

                    Cinit(Bz,Bx,By);
                    Cinit(Cz,Cx,Cy);
                    Cinit(Dz,Dx,Dy);

                    construction(Az,Bz);
                    construction(Bz,Cz);
                    construction(Cz,Dz);
                    construction(Dz,Ez);
                    }

          return(OK);
          }

main()
          {
          complexe  Az,Ez;
                                        /* Definition des extremites de la courbe.                                                   */
          int       x,y;
                                        /* Definition des coordonnees d'initialisation de l'image.                                   */

          Get(dimX,"dimX");
          Get(dimY,"dimY");
                                        /* Recuperation des dimensions en 'X' et en 'Y' de l'image a generer.                        */

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

          Cinit(Az,XG,(YB+YH)/2);
          Cinit(Ez,XD,(YB+YH)/2);
                                        /* Initialisation des extremites de la courbe.                                               */

          construction(Az,Ez);
                                        /* Construction de la courbe.                                                                */

          write(1,image,dimX*dimY);
                                        /* Sortie de l'image...                                                                      */
          }