/*******************************************************************************
+*  Fractal set generation for quadratic maps.
+*  Copyright (C) 2011, Raúl Durán Díaz, raul.duran@uah.es
+*
+*  This program is free software: you can redistribute it and/or modify
+*  it under the terms of the GNU General Public License as published by
+*  the Free Software Foundation, either version 3 of the License, or
+*  (at your option) any later version.
+*
+*  This program is distributed in the hope that it will be useful,
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*  GNU General Public License for more details.
+*
+*  You should have received a copy of the GNU General Public License
+*  along with this program.  If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
/*******************************************************************************
+*
+* 9.12.2005
+*
+* $Id: genfractal.c,v 1.2 2011/10/18 09:02:34 cvsmgr Exp $
+* $Name:  $
+*
+* Falta por hacer el cálculo del mínimo A para Gc. Ahora se usa la
+* iteración especial con valor de 50.0
+*
*******************************************************************************/
# include <stdlib.h>
# include <stdio.h>
# include <string.h>
# include <math.h>
# include <stdarg.h>
# include "mypng.h"
# include "genfractal.h"
# include "mathutil.h"

enum Zona {ZONA_Fbc_MAS, ZONA_Fbc_MENOS, ZONA_Fc_MAS, ZONA_Fc_MENOS, ZONA_Gc, ZONA_00};

static double minA = 1.0;
static double a1   = 1.0;
static double b1   = 1.0;
static double c1   = 1.0;
static double a2   = 1.0;
static double b2   = 1.0;
static double c2   = 1.0;

enum Zona zona = ZONA_Fbc_MAS;

/*******************************************************************************
+* El siguiente if a 1 para imagen en color, 0 para blanco y negro.
*******************************************************************************/
# if 1
static void PNGColorAssign(png_colorp p, int iter, int itermax)
{
double ratio = 767.0*((double) iter / (double) itermax);
int    c     = lround(ratio);

# if 0
if (c != 0 && c != 767) printf("iter: %d, ratio=%lf, c=%d\n", iter, ratio, c);
# endif
   switch (c/256)
   {
      case 0:
        p->blue = c%256;
# if 0
if (c!=0 && c!=255) printf("Caso 0\n");
# endif
      break;

      case 1:
        p->blue  = 255 - c%256;
        p->green = c%256;
# if 0
if (c!=0 && c!=255) printf("Caso 1\n");
# endif
      break;

      default:
        p->green = 255 - c%256;
        p->red   = c%256;
# if 0
if (c!=0 && c!=255) printf("Caso default\n");
# endif
   }
}
# endif

/*******************************************************************************
+* El siguiente if a 0 para imagen en color, 1 para blanco y negro.
*******************************************************************************/
# if 0
static void PNGColorAssign(png_colorp p, int iter, int itermax)
{
double ratio = 255.0*((double) iter / (double) itermax);

/*******************************************************************************
+* El siguiente if a 1 para fondo blanco, fractal negro.
*******************************************************************************/
# if 1
   p->blue = p->green = p->red = 255 - lround(ratio);
# endif
/*******************************************************************************
+* El siguiente if a 1 para fondo negro, fractal blanco.
*******************************************************************************/
# if 0
   p->blue = p->green = p->red = lround(ratio);
# endif
}
# endif

static double ComputeZR(double zr, double zi)
{
   return a1*zr*zr + 2.0*b1*zr*zi + c1*zi*zi;
}

static double ComputeZI(double zr, double zi)
{
   return a2*zr*zr + 2.0*b2*zr*zi + c2*zi*zi;
}

static double maximo(double a, double b)
{
   if (a > b)
      return a;
   else
      return b;
}

static double minimo(double a, double b)
{
   if (a < b)
      return a;
   else
      return b;
}

/*******************************************************************************
+* Función A_t para la parte F^+_{bc}
*******************************************************************************/
static double Ap_t(double t)
{
   return pow(1. - pow(t,2),2)/pow(1. + pow(t,2),2) +
          pow(t,2)/pow(1. + pow(t,2),2)*
          (
           pow(c1,2)*pow(t,2)  +
           4.*pow(b1,2)        +
           4.*b1*c1*t
          );
}


static void ComputeMinAplus00(void)
{
double t = sqrt(2./(2.+pow(c1,2)));

   minA = minimo(minA, Ap_t(t));
   minA = minimo(minA, Ap_t(-t));

   return;
}

void ComputeMinAplus(void)
{
double k;
double p;
double q;
double d;
double rho;
double theta;
double l1q;
double l2q;
int    i;

   k = (pow(c1,2)-2.*pow(b1,2)+2.)/(3.*b1*c1);
   p = (-4.*k-3.*b1*c1-2.*pow(c1,2)*k+3.*b1*c1*pow(k,2)+4.*pow(b1,2)*k)/(b1*c1);
   q = (-2.*pow(k,2)-3.*b1*c1*k+2+b1*c1*pow(k,3)-2.*pow(b1,2)+2.*pow(b1,2)*pow(k,2)-pow(c1,2)*pow(k,2))/(b1*c1);

   if ((d = 81.*pow(q,2)+12.*pow(p,3)) < 0.0)
   {
      l1q = q*q/4.;
      l2q = -d/324.;
      rho = sqrt(l1q + l2q);
      theta = atan(sqrt(l2q/l1q));

      for (i = 0; i < 3; i++)
         minA = minimo(minA, Ap_t(2.*pow(rho, 1./3.)*cos(theta/3. + 2.*i*M_PI/3.) + k));
   }
   else
   {
      fprintf(stderr, "#WARNING: Discriminante positivo, sólo hay una solución.\n");
      minA = minimo(minA, Ap_t(pow(-q/2. + sqrt(d)/18.,1./3.) + pow(-q/2. - sqrt(d)/18.,1./3.) + k));
   }

   return;
}

/*******************************************************************************
+* Función A_t para la parte F^-_{bc}
*******************************************************************************/
static double Am_t(double t)
{
   return 1 + pow(t,2)/pow(1. + pow(t,2),2)*
          (
           pow(c1,2)*pow(t,2)  +
           4.*pow(b1,2)        +
           4.*b1*c1*t
          );
}

static void ComputeMinAminus00(void)
{
   if (b1 != 0.0)
   {
      minA = minimo(minA, Am_t(1.));
      minA = minimo(minA, Am_t(-1.));
   }

   return;
}

static void ComputeMinAminus(void)
{
double k;
double p;
double q;
double d;
double rho;
double theta;
double l1q;
double l2q;
int    i;

   k = (pow(c1,2)-2.*pow(b1,2))/(3.*b1*c1);
   p = (4.*pow(b1,2)*k-2.*pow(c1,2)*k+3.*b1*c1*pow(k,2)-3.*b1*c1)/(b1*c1);
   q = (-2.*pow(b1,2)+2.*pow(b1,2)*pow(k,2)-pow(c1,2)*pow(k,2)+b1*c1*pow(k,3)-3.*b1*c1*k)/(b1*c1);

   if ((d = 81.*pow(q,2)+12.*pow(p,3)) < 0.0)
   {
      l1q = q*q/4.;
      l2q = -d/324.;
      rho = sqrt(l1q + l2q);
      theta = atan(sqrt(l2q/l1q));

      for (i = 0; i < 3; i++)
         minA = minimo(minA, Am_t(2.*pow(rho, 1./3.)*cos(theta/3. + 2.*i*M_PI/3.) + k));
   }
   else
   {
      fprintf(stderr, "#WARNING: Discriminante positivo, sólo hay una solución.\n");
      minA = minimo(minA, Am_t(pow(-q/2. + sqrt(d)/18.,1./3.) + pow(-q/2. - sqrt(d)/18.,1./3.) + k));
   }

   return;
}

/*******************************************************************************
+* Función A_t y sus derivadas para la parte F^+_{c}
*******************************************************************************/
static double AFcMas_t(double t)
{
   return 1/pow(1. + pow(t,2),2)*
          (
           pow(c1,2)*pow(t,4)  +
           4.*c1*pow(t,3)      +
           4.*pow(t,2)         +
           1.
          );
}

static void ComputeMinAFcMas0(void)
{

   minA = minimo(minA, AFcMas_t(sqrt(2.)/.2));
   minA = minimo(minA, AFcMas_t(-sqrt(2.)/.2));

   return;
}

static void ComputeMinAFcMas(void)
{
double k;
double p;
double q;
double d;
double rho;
double theta;
double l1q;
double l2q;
int    i;

   k = 1./3.*(pow(c1,2)-2.)/c1;
   p = (-2.*k*pow(c1,2) + 4.*k + 3.*c1*pow(k,2) - 3.*c1)/c1;
   q = (-1. - pow(c1,2)*pow(k,2) + 2.*pow(k,2) + c1*pow(k,3) - 3.*c1*k)/c1;

   if ((d = 81.*pow(q,2)+12.*pow(p,3)) < 0.0)
   {
      l1q = q*q/4.;
      l2q = -d/324.;
      rho = sqrt(l1q + l2q);
      theta = atan(sqrt(l2q/l1q));

      for (i = 0; i < 3; i++)
         minA = minimo(minA, AFcMas_t(2.*pow(rho, 1./3.)*cos(theta/3. + 2.*i*M_PI/3.) + k));
   }
   else
   {
      fprintf(stderr, "#WARNING: Discriminante positivo, sólo hay una solución.\n");
      minA = minimo(minA, AFcMas_t(pow(-q/2. + sqrt(d)/18.,1./3.) + pow(-q/2. - sqrt(d)/18.,1./3.) + k));
   }

   return;
}

/*******************************************************************************
+* Función A_t y sus derivadas para la parte F^-_{c}
*******************************************************************************/
static double AFcMenos_t(double t)
{
   return 1/pow(1. + pow(t,2),2)*
          (
           pow(c1,2)*pow(t,4)  -
           4.*c1*pow(t,3)      +
           4.*pow(t,2)         +
           1.
          );
}

static void ComputeMinAFcMenos0(void)
{

   minA = minimo(minA, AFcMenos_t(sqrt(2.)/.2));
   minA = minimo(minA, AFcMenos_t(-sqrt(2.)/.2));

   return;
}

static void ComputeMinAFcMenos(void)
{
double k;
double p;
double q;
double d;
double rho;
double theta;
double l1q;
double l2q;
int    i;

   k = -1./3.*(pow(c1,2)-2.)/c1;
   p = (2.*k*pow(c1,2) - 4.*k + 3.*c1*pow(k,2) - 3.*c1)/c1;
   q = (1. + pow(c1,2)*pow(k,2) - 2.*pow(k,2) + c1*pow(k,3) - 3.*c1*k)/c1;

   if ((d = 81.*pow(q,2)+12.*pow(p,3)) < 0.0)
   {
      l1q = q*q/4.;
      l2q = -d/324.;
      rho = sqrt(l1q + l2q);
      theta = atan(sqrt(l2q/l1q));

      for (i = 0; i < 3; i++)
         minA = minimo(minA, AFcMenos_t(2.*pow(rho, 1./3.)*cos(theta/3. + 2.*i*M_PI/3.) + k));
   }
   else
   {
      fprintf(stderr, "#WARNING: Discriminante positivo, sólo hay una solución.\n");
      minA = minimo(minA, AFcMenos_t(pow(-q/2. + sqrt(d)/18.,1./3.) + pow(-q/2. - sqrt(d)/18.,1./3.) + k));
   }

   return;
}

/*******************************************************************************
+* Función para la parte G_{c}
*******************************************************************************/
static void ComputeMinAGc(void)
{
double ypsilon;
double Ypsilon;
double tau1;

   ypsilon = 3*pow(c1,9)    - 9*pow(c1,8)    + 36*pow(c1,7)   -
             111*pow(c1,6)  + 540*pow(c1,5)  - 1107*pow(c1,4) +
             2592*pow(c1,3) - 3321*pow(c1,2) + 6237*c1-10044;

   if (ypsilon > 0.0)
   {
      Ypsilon = pow((c1-1.)*pow(c1+1.,2) + sqrt(ypsilon)/9., 1./3.);
      tau1 = (3*pow(Ypsilon,2) + 3*(1-c1)*Ypsilon -
             (pow(c1,3)-pow(c1,2)+3*c1-15))/(3*c1*Ypsilon);
   }
}

/*******************************************************************************
+* Cálculo del mínimo de la B y computación del Modúlo máximo.
*******************************************************************************/

static double ComputeMod2Max(double alpha, double beta)
{
double minB = 0.0;

   if (zona == ZONA_Fbc_MAS)
      minB = c1*alpha - sqrt(pow(2.*b1*alpha,2)+pow(2.*beta-c1*alpha,2));
   else if (zona == ZONA_Fbc_MENOS)
      minB = c1*alpha+2.*beta - sqrt(pow(2.*b1*alpha,2)+pow(c1*alpha,2));
   else if (zona == ZONA_Fc_MAS)
      minB = beta + alpha*c1 - sqrt(pow(beta-alpha*c1,2) + 4.*pow(alpha,2));
   else if (zona == ZONA_Fc_MENOS)
      minB = beta + alpha*c1 - sqrt(pow(beta-alpha*c1,2) + 4.*pow(alpha,2));
   else if (zona == ZONA_Gc)
      minB = alpha*(c1 + 1.) - sqrt(pow(-c1+1.,2)*pow(alpha,2)+4.*pow(alpha-beta,2));
   else if (zona == ZONA_00)
      minB = -2.*sqrt(pow(alpha,2) + pow(beta,2));
   else
   {
      fprintf(stderr, "#ERROR: Zona no implementada en el cálculo de min B.\n");
      return 0.;
   }

   return maximo(0., (1.-minB)/minA);
}

/*******************************************************************************
+* Parte fija del cómputo.
*******************************************************************************/
static void IterComputation(png_colorpp image, int itermax)
{
int         row, col, iter;
double      mod2;
double      mod2max;
double      zrtemp, zitemp;
double      zr, zi;
double      cr, ci;
# if 0
png_color   p;
# endif

   for (row = 0; row < NUM_ROWS; row++)
   {
      for (col = 0; col < NUM_COLS; col++)
      {
         cr = PixToCoordX(col);
         ci = PixToCoordY(row);

         mod2max = ComputeMod2Max(cr, ci);

         zr = 0.0;
         zi = 0.0;

         for (iter = 0; iter < itermax; iter++)
         {
            zrtemp = ComputeZR(zr, zi);
            zitemp = ComputeZI(zr, zi);

            zr = zrtemp + cr;
            zi = zitemp + ci;

            mod2 = zr*zr + zi*zi;
            if (mod2 > mod2max)
            {
# if 0
               printf("Nos salimos en la: %d\n", iter);
# endif
               break;
            }
         }
# if 0
         if (iter >= itermax)
            printf("Nos salimos por exceso, en (%lf, %lf)\n", cr, ci);
# endif

         PNGColorAssign(&image[row][col], iter, itermax);
      }
   }

# if 0
   p.red = 0; p.green = 255; p.blue = 0;
   DrawHorizLine(image, 0.0, p);
   DrawHorizLine(image, 1.0, p);
   DrawHorizLine(image, -1.0, p);

   DrawVertiLine(image, 0.0, p);
   DrawVertiLine(image, 1.0, p);
   DrawVertiLine(image, -1.0, p);
# endif
}

/*******************************************************************************
+* Parte fija del cómputo.
*******************************************************************************/
static void IterComputationEspecial(png_colorpp image, int itermax, double modmax)
{
int         row, col, iter;
double      mod2;
double      mod2max;
double      zrtemp, zitemp;
double      zr, zi;
double      cr, ci;
int         flag_changed = 0,
            mandel_zone  = 0,
            is_break     = 0;
# if 0
png_color   p;
# endif

   for (row = 0; row < NUM_ROWS; row++)
   {
      for (col = 0; col < NUM_COLS; col++)
      {

         cr = PixToCoordX(col);
         ci = PixToCoordY(row);

         mod2max = modmax;

         zr = 0.0;
         zi = 0.0;

         is_break = 0;
         for (iter = 0; iter < itermax; iter++)
         {
            zrtemp = ComputeZR(zr, zi);
            zitemp = ComputeZI(zr, zi);

            zr = zrtemp + cr;
            zi = zitemp + ci;

            mod2 = zr*zr + zi*zi;
            if (mod2 > mod2max)
            {
               if (mandel_zone)
               {
                  mandel_zone = 0;
                  flag_changed = 1;
               }
               is_break = 1;
               break;
            }
         }
         if (is_break == 0)
         {
            if (mandel_zone == 0)
            {
               mandel_zone = 1;
               flag_changed = 1;
            }
         }
         else
         {
            if (mandel_zone == 1)
            {
               mandel_zone = 0;
               flag_changed = 1;
            }
         }

         if (flag_changed)
         {
# if 0
            printf("%15.10lf %15.10lf\n", cr, ci);
# endif
            flag_changed = 0;
         }

         PNGColorAssign(&image[row][col], iter, itermax);
      }
# if 0
      printf("\n");
# endif
   }

# if 0
   p.red = 0; p.green = 255; p.blue = 0;
   DrawHorizLine(image, 0.0, p);
   DrawHorizLine(image, 1.0, p);
   DrawHorizLine(image, -1.0, p);

   DrawVertiLine(image, 0.0, p);
   DrawVertiLine(image, 1.0, p);
   DrawVertiLine(image, -1.0, p);
# endif
}

/*******************************************************************************
+* GenFractal
+* Genera la imagen fractal para las zonas \Delta^+ y \Delta^-
*******************************************************************************/
int GenFractal(char *img_name, double alpha,
                               double beta,
                               int    itermax,
                               double side_len,
                               double offset_x,
                               double offset_y)
{
png_colorpp image = NULL;
png_text    info_text[NUM_TEXT];
int         status;
double      c_temp;

/*******************************************************************************
+* Parámetros de tamaño para la imagen. Reserva de memoria para ella.
*******************************************************************************/
   SetSide(side_len);
   SetOffsetX(offset_x);
   SetOffsetY(offset_y);

   if (!(image = AllocatePNG(NUM_ROWS, NUM_COLS)))
      return ERROR;

   BuildTextPtr(info_text, alpha, beta, itermax, side_len, offset_x, offset_y);

/*******************************************************************************
+* Definición de la forma cuadrática.
+* Parte fija:
*******************************************************************************/
   a1 = 0.0;
   a2 = 1.0;
   b2 = 0.0;

/*******************************************************************************
+* Definición de la forma cuadrática.
+* Chequeo de condiciones sobre alpha y beta.
*******************************************************************************/
   if (beta == 0.0)   /* Beta ha de ser diferente de cero */
      return ERROR;

   if ((alpha + 2.*beta) == 0.0) /* alpha + 2*beta != 0.0 */
      return ERROR;

/*******************************************************************************
+* Definición de la forma cuadrática.
+* Parte variable en función de alpha y beta:
*******************************************************************************/
   c_temp = (4.*pow(alpha,3) - 27.*pow(beta,2))/pow(alpha + 2.*beta,3);

   if (c_temp < 0.0)
   {
   /* Zona F^-_{bc} */
      c2 = 1.0;
      zona = ZONA_Fbc_MENOS;
   }
   else
   {
   /* Zona F^+_{bc} */
      c2 = -1.0;
      zona = ZONA_Fbc_MAS;
   }

   b1 = -c2*(alpha - beta)/(alpha + 2.*beta);
   c1 = sqrt((-c2)*c_temp);

   if (zona == ZONA_Fbc_MAS)
   {
      if (c1 == 2.*b1  ||  c1 == -2.*b1)
      {
         fprintf(stderr, "#ERROR: Valores de b1 y c1 inválidos.\n");
         return ERROR;
      }
   }

/*******************************************************************************
+* Cálculo del mínimo de A, según valores de b1, c1 y zonas.
*******************************************************************************/
   minA = 1.;
   if (zona == ZONA_Fbc_MAS)
   {
      if (b1 == 0.0  ||  c1 == 0.0)
         ComputeMinAplus00();
      else
         ComputeMinAplus();
   }
   else if (zona == ZONA_Fbc_MENOS)
   {
      if (b1 == 0.0  ||  c1 == 0.0)
         ComputeMinAminus00();
      else
         ComputeMinAminus();
   }
   else
   {
      fprintf(stderr, "#ERROR: Zona no implementada.\n");
      return ERROR;
   }

/*******************************************************************************
+* Computación y generación de la imagen.
*******************************************************************************/
   IterComputation(image, itermax);
   status = WritePNG(img_name, image, info_text, NUM_TEXT);
   DeallocatePNG(image, NUM_ROWS, NUM_COLS);
   return status;
}

/*******************************************************************************
+* GenFractal
+* Genera la imagen fractal para las zonas C^+ y C^-
*******************************************************************************/
int GenFractalLinear(char *img_name, double alpha,
                                     double beta,
                                     int    itermax,
                                     double side_len,
                                     double offset_x,
                                     double offset_y)
{
png_colorpp image = NULL;
png_text    info_text[NUM_TEXT];
int         status;

/*******************************************************************************
+* Parámetros de tamaño para la imagen. Reserva de memoria para ella.
*******************************************************************************/
   SetSide(side_len);
   SetOffsetX(offset_x);
   SetOffsetY(offset_y);

   if (!(image = AllocatePNG(NUM_ROWS, NUM_COLS)))
      return ERROR;

   BuildTextPtr(info_text, alpha, beta, itermax, side_len, offset_x, offset_y);

/*******************************************************************************
+* Definición de la forma cuadrática.
+* Parte fija:
*******************************************************************************/
   a1 = 0.0;
   a2 = 1.0;
   b2 = 0.0;
   c2 = 0.0;

/*******************************************************************************
+* Definición de la forma cuadrática.
+* Chequeo de condiciones sobre beta.
*******************************************************************************/
   if (beta == 0.0)   /* Beta ha de ser diferente de cero */
      return ERROR;

   if (alpha + 2.*beta != 0.0)
      return ERROR;   /* alpha + 2*beta ha de ser cero */

/*******************************************************************************
+* Definición de la forma cuadrática.
+* Parte variable en función de beta: setting de zonas.
*******************************************************************************/
   if (beta*(27.0+32.0*beta) < 0.0)
   {
      b1 = -1.0;
      zona = ZONA_Fc_MENOS;
   }
   else
   {
      b1 = 1.0;
      zona = ZONA_Fc_MAS;
   }

   c1 = b1/9.0*sqrt(fabs(3.0*(27.0+32.0*beta)/beta));

/*******************************************************************************
+* Cálculo del mínimo de A, según valores de b1, c1 y zonas.
*******************************************************************************/
   minA = 1.;
   if (zona == ZONA_Fc_MAS)
   {
      if (c1 == 0.0)
         ComputeMinAFcMas0();
      else
         ComputeMinAFcMas();
   }
   else if (zona == ZONA_Fc_MENOS)
   {
      if (c1 == 0.0)
         ComputeMinAFcMenos0();
      else
         ComputeMinAFcMenos();
   }
   else
   {
      fprintf(stderr, "#ERROR: Zona no implementada.\n");
      return ERROR;
   }

/*******************************************************************************
+* Computación y generación de la imagen.
*******************************************************************************/
   IterComputation(image, itermax);
   status = WritePNG(img_name, image, info_text, NUM_TEXT);
   DeallocatePNG(image, NUM_ROWS, NUM_COLS);
   return status;
}

/*******************************************************************************
+* GenFractal
+* Genera la imagen fractal para el eje X (menos el origen).
*******************************************************************************/
int GenFractalAlpha(char *img_name, double alpha,
                                    double beta,
                                    int    itermax,
                                    double side_len,
                                    double offset_x,
                                    double offset_y)
{
png_colorpp image = NULL;
png_text    info_text[NUM_TEXT];
int         status;

/*******************************************************************************
+* Parámetros de tamaño para la imagen. Reserva de memoria para ella.
*******************************************************************************/
   SetSide(side_len);
   SetOffsetX(offset_x);
   SetOffsetY(offset_y);

   if (!(image = AllocatePNG(NUM_ROWS, NUM_COLS)))
      return ERROR;

   BuildTextPtr(info_text, alpha, beta, itermax, side_len, offset_x, offset_y);

/*******************************************************************************
+* Definición de la forma cuadrática.
+* Parte fija:
*******************************************************************************/
   a1 = 1.0;
   a2 = 0.0;
   b1 = 1.0;
   b2 = -1.0;
   c2 = 0.0;

/*******************************************************************************
+* Definición de la forma cuadrática.
+* Chequeo de condiciones sobre alpha y beta, y setting de zona.
*******************************************************************************/
   if (alpha == 0.0  ||  beta != 0.0)   /* Alpha ha de ser diferente de cero */
      return ERROR;                     /* y beta igual a cero.              */

   zona = ZONA_Gc;

/*******************************************************************************
+* Definición de la forma cuadrática.
+* Parte variable en función de alpha:
*******************************************************************************/

   c1 = (4.0*alpha - 9.0)/alpha/12.0;

/*******************************************************************************
+* Cálculo del mínimo de A, según valores de b1, c1 y zonas.
*******************************************************************************/
   ComputeMinAGc();
   minA = 1.;

# if 0
   if (zona == ZONA_Gc)
   {
      fprintf(stderr, "#ERROR: Zona no implementada.\n");
      return ERROR;
   }
   else
   {
      fprintf(stderr, "#ERROR: Zona no implementada.\n");
      return ERROR;
   }
# endif

/*******************************************************************************
+* Computación y generación de la imagen.
*******************************************************************************/
# if 0
   IterComputation(image, itermax);
# endif
   IterComputationEspecial(image, itermax, 50.0);
   status = WritePNG(img_name, image, info_text, NUM_TEXT);
   DeallocatePNG(image, NUM_ROWS, NUM_COLS);
   return status;
}

/*******************************************************************************
+* GenFractal_00
+* Genera la imagen fractal para el punto origen (0,0)
*******************************************************************************/
int GenFractal_00(char *img_name, double alpha,
                                  double beta,
                                  int    itermax,
                                  double side_len,
                                  double offset_x,
                                  double offset_y)
{
png_colorpp image = NULL;
png_text    info_text[NUM_TEXT];
int         status;

/*******************************************************************************
+* Parámetros de tamaño para la imagen. Reserva de memoria para ella.
*******************************************************************************/
   SetSide(side_len);
   SetOffsetX(offset_x);
   SetOffsetY(offset_y);

   if (!(image = AllocatePNG(NUM_ROWS, NUM_COLS)))
      return ERROR;

   BuildTextPtr(info_text, alpha, beta, itermax, side_len, offset_x, offset_y);

/*******************************************************************************
+* Definición de la forma cuadrática.
*******************************************************************************/
   a1 = 0.0;
   b1 = 1.0;
   c1 = 0.0;
   a2 = 1.0;
   b2 = 0.0;
   c2 = -1.0;

/*******************************************************************************
+* Chequeo de condiciones sobre alpha y beta. Setting de zona.
*******************************************************************************/
   if (alpha != 0.0  ||  beta != 0.0)   /* Alfa y Beta han de ser cero. */
      return ERROR;

   zona = ZONA_00;

/*******************************************************************************
+* Cálculo del mínimo de A, según valores de b1, c1 y zonas.
*******************************************************************************/
   minA = 1.;
   if (zona == ZONA_00)
   {
      minA = 1.;
   }
   else
   {
      fprintf(stderr, "#ERROR: Zona no implementada.\n");
      return ERROR;
   }

/*******************************************************************************
+* Computación y generación de la imagen.
*******************************************************************************/
   IterComputation(image, itermax);
   status = WritePNG(img_name, image, info_text, NUM_TEXT);
   DeallocatePNG(image, NUM_ROWS, NUM_COLS);
   return status;
}

/*******************************************************************************
+* GenFractalEspecial1
+* Genera la imagen para F^+_c cuando c = 0.
*******************************************************************************/
int GenFractalEspecial1(char *img_name, double alpha,
                                        double beta,
                                        int    itermax,
                                        double side_len,
                                        double offset_x,
                                        double offset_y)
{
png_colorpp image = NULL;
png_text    info_text[NUM_TEXT];
int         status;

/*******************************************************************************
+* Parámetros de tamaño para la imagen. Reserva de memoria para ella.
*******************************************************************************/
   SetSide(side_len);
   SetOffsetX(offset_x);
   SetOffsetY(offset_y);

   if (!(image = AllocatePNG(NUM_ROWS, NUM_COLS)))
      return ERROR;

   BuildTextPtr(info_text, alpha, beta, itermax, side_len, offset_x, offset_y);

/*******************************************************************************
+* Definición de la forma cuadrática.
*******************************************************************************/
   a1 = 0.0;
   b1 = 1.0;
   c1 = 0.0;
   a2 = 1.0;
   b2 = 0.0;
   c2 = 0.0;

   IterComputationEspecial(image, itermax, 20.0);
   status = WritePNG(img_name, image, info_text, NUM_TEXT);
   DeallocatePNG(image, NUM_ROWS, NUM_COLS);
   return status;
}

/*******************************************************************************
+* GenFractalEspecial2
+* Genera la imagen para F^-_c cuando c = 0.
*******************************************************************************/
int GenFractalEspecial2(char *img_name, double alpha,
                                        double beta,
                                        int    itermax,
                                        double side_len,
                                        double offset_x,
                                        double offset_y)
{
png_colorpp image = NULL;
png_text    info_text[NUM_TEXT];
int         status;

/*******************************************************************************
+* Parámetros de tamaño para la imagen. Reserva de memoria para ella.
*******************************************************************************/
   SetSide(side_len);
   SetOffsetX(offset_x);
   SetOffsetY(offset_y);

   if (!(image = AllocatePNG(NUM_ROWS, NUM_COLS)))
      return ERROR;

   BuildTextPtr(info_text, alpha, beta, itermax, side_len, offset_x, offset_y);

/*******************************************************************************
+* Definición de la forma cuadrática.
*******************************************************************************/
   a1 = 0.0;
   b1 = -1.0;
   c1 = 0.0;
   a2 = 1.0;
   b2 = 0.0;
   c2 = 0.0;

   IterComputationEspecial(image, itermax, 20.0);
   status = WritePNG(img_name, image, info_text, NUM_TEXT);
   DeallocatePNG(image, NUM_ROWS, NUM_COLS);
   return status;
}

/*******************************************************************************
+* GenFractalEspecial3
+* Genera la imagen para G^_c cuando c = 0.
*******************************************************************************/
int GenFractalEspecial3(char *img_name, double alpha,
                                        double beta,
                                        int    itermax,
                                        double side_len,
                                        double offset_x,
                                        double offset_y)
{
png_colorpp image = NULL;
png_text    info_text[NUM_TEXT];
int         status;

/*******************************************************************************
+* Parámetros de tamaño para la imagen. Reserva de memoria para ella.
*******************************************************************************/
   SetSide(side_len);
   SetOffsetX(offset_x);
   SetOffsetY(offset_y);

   if (!(image = AllocatePNG(NUM_ROWS, NUM_COLS)))
      return ERROR;

   BuildTextPtr(info_text, alpha, beta, itermax, side_len, offset_x, offset_y);

/*******************************************************************************
+* Definición de la forma cuadrática.
*******************************************************************************/
   a1 = 1.0;
   b1 = 1.0;
   c1 = 0.0;
   a2 = 0.0;
   b2 = -1.0;
   c2 = 0.0;

   IterComputationEspecial(image, itermax, 20.0);
   status = WritePNG(img_name, image, info_text, NUM_TEXT);
   DeallocatePNG(image, NUM_ROWS, NUM_COLS);
   return status;
}

/*******************************************************************************
+* GenFractalEspecial5
+* Genera el caso F^+_{00}
*******************************************************************************/
int GenFractalEspecial5(char *img_name, double alpha,
                                        double beta,
                                        int    itermax,
                                        double side_len,
                                        double offset_x,
                                        double offset_y)
{
png_colorpp image = NULL;
png_text    info_text[NUM_TEXT];
int         status;

/*******************************************************************************
+* Parámetros de tamaño para la imagen. Reserva de memoria para ella.
*******************************************************************************/
   SetSide(side_len);
   SetOffsetX(offset_x);
   SetOffsetY(offset_y);

   if (!(image = AllocatePNG(NUM_ROWS, NUM_COLS)))
      return ERROR;

   BuildTextPtr(info_text, alpha, beta, itermax, side_len, offset_x, offset_y);

/*******************************************************************************
+* Definición de la forma cuadrática.
*******************************************************************************/
   a1 = 0.0;
   b1 = 0.0;
   c1 = 0.0;
   a2 = 1.0;
   b2 = 0.0;
   c2 = -1.0;

   IterComputationEspecial(image, itermax, 40.0);
   status = WritePNG(img_name, image, info_text, NUM_TEXT);
   DeallocatePNG(image, NUM_ROWS, NUM_COLS);
   return status;
}

/*******************************************************************************
+* GenFractalEspecial6
+* Genera el caso F^-_{00}
*******************************************************************************/
int GenFractalEspecial6(char *img_name, double alpha,
                                        double beta,
                                        int    itermax,
                                        double side_len,
                                        double offset_x,
                                        double offset_y)
{
png_colorpp image = NULL;
png_text    info_text[NUM_TEXT];
int         status;

/*******************************************************************************
+* Parámetros de tamaño para la imagen. Reserva de memoria para ella.
*******************************************************************************/
   SetSide(side_len);
   SetOffsetX(offset_x);
   SetOffsetY(offset_y);

   if (!(image = AllocatePNG(NUM_ROWS, NUM_COLS)))
      return ERROR;

   BuildTextPtr(info_text, alpha, beta, itermax, side_len, offset_x, offset_y);

/*******************************************************************************
+* Definición de la forma cuadrática.
*******************************************************************************/
   a1 = 0.0;
   b1 = 0.0;
   c1 = 0.0;
   a2 = 1.0;
   b2 = 0.0;
   c2 = 1.0;

   IterComputationEspecial(image, itermax, 40.0);
   status = WritePNG(img_name, image, info_text, NUM_TEXT);
   DeallocatePNG(image, NUM_ROWS, NUM_COLS);
   return status;
}

/*******************************************************************************
+* GenFractalEspecial7
+* Genera el caso F^+_{bc} cuando $c=\pm 2b$
*******************************************************************************/
int GenFractalEspecial7(char *img_name, double alpha,
                                        double beta,
                                        int    itermax,
                                        double side_len,
                                        double offset_x,
                                        double offset_y)
{
png_colorpp image = NULL;
png_text    info_text[NUM_TEXT];
int         status;

/*******************************************************************************
+* Parámetros de tamaño para la imagen. Reserva de memoria para ella.
*******************************************************************************/
   SetSide(side_len);
   SetOffsetX(offset_x);
   SetOffsetY(offset_y);

   if (!(image = AllocatePNG(NUM_ROWS, NUM_COLS)))
      return ERROR;

   BuildTextPtr(info_text, alpha, beta, itermax, side_len, offset_x, offset_y);

/*******************************************************************************
+* Definición de la forma cuadrática.
+* Parte fija:
*******************************************************************************/
   a1 = 0.0;
   a2 = 1.0;
   b2 = 0.0;
   c2 = -1.0;
   b1 = beta;
   c1 = 2*beta;

/*******************************************************************************
+* Computación y generación de la imagen.
*******************************************************************************/
   IterComputationEspecial(image, itermax, 40.0);
   status = WritePNG(img_name, image, info_text, NUM_TEXT);
   DeallocatePNG(image, NUM_ROWS, NUM_COLS);
   return status;
}