Versión 1.0 - julio de 2005

fractal.c

/*
 * Compilar con la opcion c99:
 * gcc -std=c99 -c filename.c -o filename.o
 *
 * Enlazar con la opcion de librerías matematicas y con allegro:
 * gcc  -lm -lX11 -o filename filename.o `allegro-config --libs` 
 *
 * En Anjuta, ejecutar en terminal allegro-config --libs,
 * y su salida son las librerías que tenemos que pasar 
 * al gcc. Así que vamos a Opciones-> Comandos -> Construir 
 * Archivo y le ponemos ahí todas las opciones.
 */

/*
fractal.c 	
Este fractal es la portada del libro:
Análisis Numérico Las matemáticas del Cálculo científico, de Kincaid, y que
por desgracia está descatalogado. En él se explica cómo programarlo y aquí 
también claro, ahí va:

Dibuja los fractales de las cuencas de atracción de las 
raíces cúbicas de la unidad. El polinomio es p(z)=z**3-1. Tiene una raíz 
real (1) y dos complejas: (-1+sqrt(3)*I)/2 y (-1-sqrt(3)*I)/2 donde I es 
sqrt(-1).
Mediante el método de Newton realizamos 40 iteraciones y si queda a menos 
de 0.25 de una raíz, consideramos que converge a ella, y le asignamos un 
color. 
Hay puntos que no convergen a ninguna raíz y se llaman conjunto de Julia 
(1918). Estos puntos forman un fractal. Así que el límite entre cuencas es 
donde se muestra la geometría fractal.
*/
#include 
#include "allegro.h"
#include 
#include 
int
main (void)
{
  PALETTE palette;
  int i, j, n;
  int ancho_pantalla, alto_pantalla;
  int amarillo, verde, rojo, azul;
  double x, y, a, id, jd;
  _Complex double z, r2, r3;


/* Obtenemos el tamaño de la pantalla */
  Display *disp = XOpenDisplay (NULL);
  Screen *scrn = DefaultScreenOfDisplay (disp);
  alto_pantalla = scrn->height;
  ancho_pantalla = scrn->width;
/*Iniciar modo grafico*/
  if (allegro_init () != 0)
    return 1;
  install_keyboard ();
  set_color_depth (16);
  if (set_gfx_mode (GFX_AUTODETECT, ancho_pantalla, alto_pantalla, 0, 0) != 0)
    {
      if (set_gfx_mode (GFX_AUTODETECT, 640, 480, 0, 0) != 0)
	{
	  allegro_message ("Error setting graphics mode\n%s\n",
			   allegro_error);
	  return 1;
	}
    }
  set_palette (black_palette);
  set_palette (palette);
  amarillo = makecol (255, 255, 50);
  verde = makecol (0, 255, 0);
  azul = makecol (0, 0, 255);
  rojo = makecol (255, 0, 0);

/* Código del fractal*/

  r2 = (-1 + sqrt (3) * I) / 2;
  r3 = (-1 - sqrt (3) * I) / 2;

  for (i = 0; i <= ancho_pantalla; i++)
    {
      for (j = 0; j <= alto_pantalla; j++)
	{

/* Le restamos ancho_pantalla/2 y alto_pantalla/2 para que el punto (0,0) 
del plano complejo quede en la mitad de la pantalla.
*/
	  id = i - ancho_pantalla / 2;
	  jd = j - alto_pantalla / 2;
	  z = (id + jd * I );
	  for (n = 1; n <= 40; n++)
	    z = z - (z / 3) + (1 / (3 * z * z));
	  x = cabs (z - r2);
	  /* Aunque haga los cálculos para todas las raíces no merece 
	     la pena poner else en los if porque tarda lo mismo. */
	  if (x <= 0.25)
	    putpixel (screen, i, j, rojo);
	  x = cabs (z - r3);
	  if (x <= 0.25)
	    putpixel (screen, i, j, amarillo);
	  x = cabs (z - 1);
	  if (x <= 0.25)
	    putpixel (screen, i, j, azul);
	}
    }
  do
    {
    }
  while (!keypressed ());
  return 0;
}

END_OF_MAIN ();