/*----------------------------------------------------------*/
/*                                                          */
/*         AUTHOR : Eric VIOLARD                            */
/*         E-MAIL : violard@icps.u-strasbg.fr               */
/*         ORGANISM : Université Louis Pasteur (Strasbourg) */
/*         CREATION : 23/11/02                              */
/*                                                          */
/* ---------------------------------------------------------*/
  
/*  

     Programme permettant de dessiner l'ensemble de Mandelbrot.
 
     Le programme demande à l'utilisateur les dimensions de la fenêtre graphique.

*/

#include <stdio.h> /* pour utiliser les fonctions printf et scanf */
#include <graph.h> /* pour utiliser les fonctions graphiques */

/* --- fonctions utiles */

float min(float a, float b)
     /* minimum de 2 réels */
{
  float r;

  r=a;
  if(b<a) r=b;
  return r;
}

float max(float a, float b)    
    /* maximum de 2 réels */
{
  float r;

  r=b;
  if(a>b) r=a;
  return r;
}

/* --- */

/* --- Définition du type des nombres complexes --- */

/* munis des opérations : 

  - partie réelle
  - partie imaginaire
  - construction d'un nombre complexe
  - addition 
  - multiplication
  - module
*/
 
typedef struct { float reel; float imag; } complexe; 
/*
   type des nombre complexes définis par ses parties réelles et imaginaires 
*/


float partie_reelle(complexe z)
     /* partie réelle d'un nombre complexe */
{
  return z.reel;
}

float partie_imaginaire(complexe z)
     /* partie imaginaire d'un nombre complexe */
{
  return z.imag;
}

complexe cons_complexe(float reel, float imag)
     /* construction d'un nombre complexe à partir de sa partie réelle et sa partie imaginaire */
{
  complexe z;
  
  z.reel = reel;
  z.imag = imag;
 
  return z;
}


complexe somme_complexe(complexe z1, complexe z2)      
     /* somme de 2 nombres complexes */
{
    return cons_complexe(
      partie_reelle(z1) + partie_reelle(z2),
      partie_imaginaire(z1) + partie_imaginaire(z2)
    );
}

complexe produit_complexe(complexe z1, complexe z2)     
     /* produit de 2 nombres complexes */
{
    return cons_complexe(
      partie_reelle(z1)*partie_reelle(z2) - 
      partie_imaginaire(z1)*partie_imaginaire(z2),
      partie_imaginaire(z1)*partie_reelle(z2) + 
      partie_reelle(z1)*partie_imaginaire(z2)
    );
}

float module_complexe(complexe z)  
    /* module d'un nombre complexe */
{
    return 
      partie_reelle(z)*partie_reelle(z) + 
      partie_imaginaire(z)*partie_imaginaire(z);
}

void ecriture_complexe(complexe z)  
     /* affiche un nombre complexe */
{
    printf("%f+i*%f\n",partie_reelle(z),partie_imaginaire(z));
}


/* --- */

/* --- Définition du type des intervalles fermés de nombres réels --- */

/* munis des opérations : 

  - construction d'un intervalle, 
  - borne inf
  - borne sup
  - thalès
*/
 
typedef struct { float a; float b; } intervalle; 
/* 
   type des intervalles fermés définis par 2 réels a et b 
        [a,b] si a<=b ou
        [b,a] si b<=a  
*/


intervalle cons_intervalle(float a,float b)
     /* construction d'un intervalle à partir de 2 réels */
{
  intervalle I;
 
  I.a = a;
  I.b = b;

  return I;
}


float borne_inf(intervalle I)
     /* borne inférieure d'un intervalle */
{
  return min(I.a,I.b);
}

float borne_sup(intervalle I)
     /* borne supérieure d'un intervalle */
{
  return max(I.a,I.b);
}


float thales(float x1, intervalle I1, intervalle I2)
     /* détermine le réel x1 de l'intervalle I1 
        correspondant au réel x2 de l'intervalle I2 
        en utilisant Thalès */ 
{
  float x2;

  x2 = (x1-I1.a) * ((I2.b-I2.a)/(I1.b-I1.a)) + I2.a;

  return x2;
}

void ecriture_intervalle(intervalle I) 
     /* affiche un intervalle */
{
    printf("[%f..%f]\n",borne_inf(I),borne_sup(I));
}

/* --- */

#define N 200 /* borne utilisée pour déterminer la convergence des suites */



main()
{
  

  /* déclarations */
  
  int largeur,hauteur; /* dimensions de la fenêtre */
  complexe zmin=cons_complexe(-2.05,-1.3),
           zmax=cons_complexe(0.55,1.3); /* défini l'intervalle des nombres complexes considérés */
  int X,Y; /* les coordonnées d'un pixel de la fenêtre graphique */
  complexe c; /* le nombre complexe dans l'intervalle [zmin,zmax] et qui correpond au pixel (X,Y)  */
  complexe z;  /* termes de la suite (zn) */
  int n; /* indice des termes de la suite (zn) */


  /* --- initialisation de la fenêtre graphique */

  printf("Entrer les dimensions de la fenêtre graphique (largeur x hauteur) en nombre de pixels : "); 
  scanf("%d %d",&largeur,&hauteur); 
  initgraph(largeur,hauteur);
  cleargraph();


  /* --- tracé de l'ensemble --- */ 

  setcolor(0.,1.,0.); /* en vert */


  /* on parcours tous les pixels */

  for(X=0;X<largeur;X++)
  for(Y=0;Y<hauteur;Y++) 
    {

      /* --- construction de c, le nombre complexe correspondant au pixel (X,Y) --- */
      { 
      float x,y; /* partie réelle et partie imaginaire de c (variables temporaires) */

      x = thales((float)X,cons_intervalle(0.,(float)largeur),
                   cons_intervalle(partie_reelle(zmin),partie_reelle(zmax)));
      y = thales((float)Y,cons_intervalle(0.,(float)hauteur),
                   cons_intervalle(partie_imaginaire(zmin),partie_imaginaire(zmax)));
      c = cons_complexe(x,y);
      }
 
      /* --- calcul de termes de la suite (zn) --- */
      
      z = cons_complexe(0.,0.);      
      n = 0;

      while(module_complexe(z)<2 && n<=N)
        {
	  z = somme_complexe(produit_complexe(z,z),c);   
	  n = n+1;
        }

      
      /* --- illumination éventuelle du pixel (X,Y) --- */

      if(module_complexe(z)<2) putpixel(X,Y);
    }

  flushgraph(); /* pour forcer l'affichage */


  /* --- terminaisons --- */

  waitgraph(); /* attente d'un frappe au clavier */
  closegraph(); /* fermeture de la fenêtre graphique */

}