/*
			fractales.c
			===========

     *....................................................*
     .        warning, this is a 'work in progress'       .
     .                                                    .
     .        don't trust prototypes for building         .
     .        serious applications...                     .
     *....................................................*

*/

#include	<stdio.h>
#include	<stdlib.h>
#include	<math.h>
#include	<time.h>

#include        <tthimage.h>

#include	"fractales.h"

/*::------------------------------------------------------------------::*/
void Image_fractales_print_version(int flag)
{
printf("*** Fractales (%s) v 0.0.12\n", __FILE__);
if (flag)
    {
    printf("\tcompiled %s , %s\n", __DATE__, __TIME__);
    }
}
/*::------------------------------------------------------------------::*/
/*
 *	fonction de comparaison de deux variables 'double'
 */
int	compare_double(double a, double b)
{

return 0;
}
/*::------------------------------------------------------------------::*/
/*
 *	cette fonction a été repompé (sans honte :) du livre de
 *	Roger T. Stevens, "Fractal Programming in C" éditée par
 *	M&T Books.
 *	Ceci dit, elle est un peu trop monolithique, puisque l'équation
 *	à résoudre est codée en dur.
 */
int Newton_0(Image_Desc *im, int maxiter)
{
int	col, row, i;
double	dcol, drow;
double	Xmax, Xmin, Ymax, Ymin;
double  deltaX, deltaY, denom;
double	X, Y, Xsquare, Ysquare, Xold, Yold;

Xmax = Ymax = 3.5;	Xmin = Ymin = -3.5;
deltaX = (Xmax-Xmin)/(double)im->width;
deltaY = (Ymax-Ymin)/(double)im->height;

for (col=0; col<im->width; col++)
	{
	dcol = (double)col;
	for (row=0; row<im->height; row++)
		{
		drow = (double)row;
		X = Xmin + dcol * deltaX;
		Y = Ymin + drow * deltaY;
		Xsquare = 0.0;
		Ysquare = 0.0;
		Xold = Yold = 42.0;
		for (i=0; i<maxiter; i++)
			{
			Xsquare = X * X;
			Ysquare = Y * Y;
			denom = 3*((Xsquare-Ysquare)*(Xsquare-Ysquare) +
				4*Xsquare*Ysquare);
			if (denom < 0.0000001)
				denom = 0.0000001;

			X = .6666667*X + (Xsquare - Ysquare)/denom;
			Y = .6666667*Y + 2*X*Y/denom;
			if ((Xold==X) && (Yold==Y))
				break;
			}		/* finbo sur maxiter */

		if (X > 0)
			{
			Image_plotRGB(im, col, row, 0, 0, 0);
			}
		else
			{
			if (Y > 0)
				Image_plotRGB(im, col, row, 255, 0, 0);
			else
				Image_plotRGB(im, col, row, 0, 0, 255);
			}
		}
	if (0 == col%5)
		fprintf(stderr, "newton 0: col %5d / %5d\r", col, im->width);
	}

return 0;
}
/*::------------------------------------------------------------------::*/
/*
 * C'est pas encore fini ce truc ?
 */
int GingerBreadMan(Image_Desc *im, int maxiter)
{
double	x1, y1, x2, y2;
int	iter;

x1 = 10.0;
y1 = 0.333;

for (iter=0; iter<maxiter; iter++)
	{
	x2 = 1.0 - y1 + abs(x1);
	y2 = x1;

	printf("%9g  %9g\n", x2, y2);

	x1 = x2;
	y1 = y2;
	}

return 0;
}
/*::------------------------------------------------------------------::*/
/*
	En fait, il faudrait travailler sur deux images...
*/
int Image_Random_Walker_0(Image_Desc *im, Image_Desc *dst,
					int nbpass, int nbsteps)
{
int	pass, step, foo;
int	x, y, ra, ga, ba, rb, gb, bb;

if ( (foo=Image_compare_desc(im, dst)) )
    {
    fprintf(stderr, "Random_Walker_0: images are differents %d\n", foo);
    return foo;
    }

for (pass=0; pass<nbpass; pass++)
	{
	x = rand() % im->width;
	y = rand() % im->height;

	ra = Image_R_pixel(im, x, y);
	ga = Image_G_pixel(im, x, y);
	ba = Image_B_pixel(im, x, y);

	/* printf("\tP%d\t\tX%d\tY%d\t%d,%d,%d\n", pass, x, y, ra, ga, ba); */

	for (step=0; step<nbsteps; step++)
		{
		x += (rand()%3)-1;
		if (x<0)		x += im->width;
		if (x>=im->width)	x -= im->width;
		y += (rand()%3)-1;
		if (y<0)		y += im->height;
		if (y>=im->height)	y -= im->height;

		rb = Image_R_pixel(im, x, y);
		gb = Image_G_pixel(im, x, y);
		bb = Image_B_pixel(im, x, y);

		Image_plotRGB(dst, x, y,
			(ra+rb+rb)/3, (ga+gb+gb)/3, (ba+bb+bb)/3);
		}
	}

return 0;
}

/*::------------------------------------------------------------------::*/

#define  XMIN	(-1.5)
#define  XMAX	(1.5)
#define  YMIN	(-1.5)
#define  YMAX	(1.5)

int	Image_Julia_0_0(Image_Desc *im, RGB_map *pmap, int fuzz,
						double cx, double cy)
{
int	col, row;
double	deltaX, deltaY, dx, dy, x2, y2;
long	debut, fin;
int	iter, r, g, b, maxiter;

debut = time(NULL);

fprintf(stderr, "Julia_0_0: %dx%d\n", im->width, im->height);

deltaX = (XMAX - XMIN) / im->width;
deltaY = (YMAX - YMIN) / im->height;

maxiter = 0;

/*
 *	si pas de palette définie, en créer une ? 
 */
for (row=0; row<im->height; row++)
	{
	dy = YMIN + (double)row * deltaY;

	for (col=0; col<im->width; col++)
		{
#if DEBUG_LEVEL
		printf("%5d  %5d     ", row, col);
#endif

		dx = XMIN + (double)col * deltaX;

		iter = 0;
		x2 = y2 = 0.0;

		while (	(iter < fuzz) &&
			((x2 + y2) < 4.0) )
			{
			x2 = dx * dx;
			y2 = dy * dy;
			dx = 2 * dx * dy + cx;
			dy = x2 - y2 + cy;
			iter++;
			}

		if (iter > maxiter)
			maxiter = iter;

		if (iter >= fuzz)
			{
			r = g = b = 255;
			}
		else
			{
			r = pmap->red[iter%pmap->nbre];
			g = pmap->green[iter%pmap->nbre];
			b = pmap->blue[iter%pmap->nbre];
			}

		Image_plotRGB(im, col, row, r, g, b);
		
#if DEBUG_LEVEL
		printf("  %3d \r", iter);
#endif
		}
	}

fin = time(NULL);

printf("duree du calcul: %ld secondes, maxiter: %d\n", fin-debut, maxiter);

return 0;
}
/*::------------------------------------------------------------------::*/

#define  XA    -2.0
#define  XB     2.0
#define  YA    -2.0
#define  YB     1.3

int Image_Mandelbrot_0(Image_Desc *image, int maxiter)
{
int         Y, X;
int         r, g, b;
double      zr, zi, fx, fy, real, imag, fiter, dist2;
long        debut, fin, temps;         /* pour chronometrer */
int         iter, out, foo;
double      tmpr, tmpi;

debut = time(NULL);

out = 0;

for (Y=0; Y<image->height; Y++)
    {
    fprintf(stderr, "ligne %4d sur %d\r", Y, image->height);
    fy = (double)Y;
    real = ( fy / (double)image->height * (YB-YA) ) + YA;

    for (X=0; X<image->width; X++)
        {
        fx = (double)X;
        imag = ( fx / (double)image->width * (XB-XA) ) + XA;

        iter = 0;
        zr = zi = 0.0;

        do
            {
            iter++;

            dist2 = zr*zr + zi*zi;
            if (dist2 > 4.0)
                {
                out++;
                break;
                }

	    /*
             *  calculer le mandelbrot
             *  ----------------------
             *      la formule est:
             *      z0 = 0 + 0i
             *      z1 = z0^2 + c;
             *      ...
	     */
            tmpr = zr*zr - zi*zi;
            tmpi = 2.0 * zr * zi;

            zr = tmpr + real;
            zi = tmpi + imag;

            } while (iter < maxiter);

        fiter = (double)iter / (double)maxiter;

        Image_color_x(zr, zi, &r, &g, &b);
        Image_plotRGB(image, X, Y, r, g, b);

        }

    fflush(stderr);
    }

fin = time(NULL);

temps = fin - debut;
if (temps < 500)
    {
    foo = 's';
    }
else
    {
    foo = 'm';
    temps /= 60;
    }

printf("\nMandelbrot_0: temps d'execution = %ld %c\n", temps, foo);
printf("Mandelbrot_0: nombre de points dehors = %d\n", out);

return 0;
}
/*::------------------------------------------------------------------::*/
int Lorenz_Orbits(Image_Desc *img, int iters,
				double a, double b, double c, double dt)
{
int		pass;
double		dx, dy, dz, dxp, dyp, dzp;
FILE		*fp;

if (NULL == (fp=fopen("lorenz.dat", "w"))) {
	perror("lorenz datafile");
	exit(1);
	}

fprintf(stderr, "%s ( %p %d / %f %f %f %f )\n", __func__, img, iters,
				a, b, c, dt);

dx = dy = dz = 0.0001;

for (pass=0; pass<iters; pass++) {

	fprintf(fp, "%6d  %8f %8f %8f\n", pass, dx, dy, dz);

	dxp = dx - (a*dx + a*dy) * dt;
	dyp = dy + (b*dx - dy - dz*dz) * dt;
	dzp = dz - (c*dz + dx*dy) * dt;

	dx = dxp, dy = dyp, dz = dzp;
	}

fclose(fp);

return 42;
}
/*::------------------------------------------------------------------::*/
			/*
			 * et un PopCorn ?
			 */
/*::------------------------------------------------------------------::*/