/****************************************************************/
/*								*/
/* julia.c:	a program that calculates the julia		*/
/*		set (fractal).					*/
/*		This is supposed to be an example of		*/
/*		how the create a simple interface to		*/
/*		XDim						*/
/*								*/
/*		to compile do:					*/
/*			cc -o julia julia.c			*/
/*								*/
/*		Design: Walter Benzing 1995			*/
/*		NO RIGHTS RESERVED				*/
/*								*/
/****************************************************************/

/*
 * standard includes
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

/*
 * Define Constants
 */

#define NUM_ITER 30		/* Number of iterations */
#define MAX_VAL 9.0		/* Outside begins infinite */

/*
 * Define Complex Sructure
 */

typedef struct
  {
   double real;
   double imag;
  } myComplex;

/*
 * Complex math as macros
 * WARNING:	SIDE EFFECT => c0,c1 and c2 MUST be DIFFERENT
 */

#define CADD(c0,c1,c2)\
  {\
   c0.real = c1.real+c2.real;\
   c0.imag = c1.imag+c2.imag;\
  }

#define CSUB(c0,c1,c2)\
  {\
   c0.real = c1.real-c2.real;\
   c0.imag = c1.imag-c2.imag;\
  }

#define CMUL(c0,c1,c2)\
  {\
   c0.real = c1.real*c2.real-c1.imag*c2.imag;\
   c0.imag = c1.real*c2.imag+c1.imag*c2.real;\
  }

#define CDIV(c0,c1,c2)\
  {\
   c0.imag = c2.real*c2.real+c2.imag*c2.imag;\
   c0.real = (c1.real*c2.real+c1.imag*c2.imag);\
   c0.imag = (c1.imag*c2.real-c1.real*c2.imag);\
  }
  
#define CABS(c0) (c0.real*c0.real+c0.imag*c0.imag)

/************************************************************************/
/*									*/
/*	Save a Data Field in a suitable Way for XDim			*/
/*									*/
/*	fileName:	Name of the Data file				*/
/*	width:		width of the Data Field				*/
/*	height:		height of the data field			*/
/*	field:		pointer to a (double) array			*/
/*	factX,offX:		X-Value = factX*i + offX; 0<=i<width	*/
/*	factY,offY:		Y-Value = factY*j + offY; 0<=j<height	*/
/*									*/
/************************************************************************/


void EX_SaveData(char *fileName, int width, int height, double *field,
	double factX, double offX, double factY, double offY)
  {
   int i, size;
   double *dp;
   FILE *out;
   
   if(field)
     {
      out = fopen(fileName, "w");
      if(out)
	{
	 fprintf(out, "XDim data file:\n");	/* tell it's an XDim data File*/
	 fprintf(out, "Max X: %d\n", width);	/* number of points in x,y */
	 fprintf(out, "Max Y: %d\n", height);
	 fprintf(out, "Dim X: %s\n", "real");	/* name on axis */
	 fprintf(out, "Dim Y: %s\n", "imag");
	 fprintf(out, "Dim Z: %s\n", "fun");
	 fprintf(out, "Offset X: %e\n", offX);	/* X-Offset */
	 fprintf(out, "Form Factor X: %e\n", factX);	/* X-Factor */
	 fprintf(out, "Offset Y: %e\n", offY);	/* Y-Offset */
	 fprintf(out, "Form Factor Y: %e\n", factY);	/* Y-Factor */
	 fprintf(out, "Data:\n");		/* Start Data output */
	 size = width*height;
	 dp = field;
	 for(i=0; i<size; i++)
	    {
             fprintf(out,"%e\n",*dp++);
	    }
	}
      else
      	printf("Couldn't open %s for write !!", fileName);      
     }
  }

/************************************************************************/
/*									*/
/*	Caculate Value							*/
/*	returns z=f(x,y)						*/
/*									*/
/*	with two parameters param1 and param2				*/
/*									*/
/*	This example calculates the julia set				*/
/*	this is not intended to be fast or clever,			*/
/*	just a basic example						*/
/*									*/
/* PUT YOUR FAVORIT FUNCTION HERE !!!!					*/
/************************************************************************/

double EX_myFunc(double x, double y, double param1, double param2)
  {
   int i;
   myComplex z0, z1, z2;
   double val;
   
   z0.real = x;
   z0.imag = y;
   z1.real = param1;
   z1.imag = param2;
   for(i=0; i<NUM_ITER; i++)
      {
       CMUL(z2,z0,z0);
       CADD(z0,z2,z1);
       val =CABS(z0);
       if(val > MAX_VAL)
         break;
      }
   if(val > MAX_VAL)
     val = MAX_VAL;
   return(MAX_VAL-val);
  }

/************************************************************************/
/*									*/
/*	Main Program							*/
/*									*/
/*	Read in parameters etc.						*/
/*									*/
/************************************************************************/

int main(int argc, char **argv)
   {
    char fileName[256];
    double minX, maxX, factX, minY, maxY, factY, x, y, z, par1, par2;
    double *field, *dp, fact;
    int width, height,i,j;
    
    /*** Second argument is filename if available ***/

    printf("This Program calculates the julia set of z*z+c\n");
    printf("Data format is suitable for XDim\n\n");
    if(argc > 1)
      strncpy(fileName, argv[1],256);
    else
      {
       printf("file name: ");
       scanf("%s", fileName);
      }
    width = 200;
    height = 200;
    minX = -2.0;
    maxX = 2.0;
    minY = -2.0;
    maxY = 2.0;
    par1 = 1.0;
    par2 = 0.0;
    printf("width (200): ");
    scanf("%d", &width);
    printf("height (200): ");
    scanf("%d", &height);
    printf("Min X [-2]: ");
    scanf("%le", &minX);
    printf("Max X [+2]: ");
    scanf("%le", &maxX);
    printf("Min Y [-2]:");
    scanf("%le", &minY);
    printf("Max Y [+2]: ");
    scanf("%le", &maxY);
    printf("c-real [-1]: ");
    scanf("%le", &par1);
    printf("c-imag [0]: ");
    scanf("%le", &par2);
    if(width <=1)
      width = 10;
    if(height <=1)
      height = 10;

    /*
     * note (width-1), (height-1)
     */

    factX = (maxX-minX)/(width-1);
    factY = (maxY-minY)/(height-1);
    /*
     * allocate memory
     */
    field = malloc(width*height*sizeof(double));
    if(field)
      {
       /*
        * calculate data
        */
       dp = field;
       y = minY;
       fact = 100.0/(double)(height-1);
       for(i=0; i<height; i++)
          {
           y += factY;
           x = minX;
           for(j=0; j< width; j++)
              {
               x += factX;
               *dp++ = EX_myFunc(x,y,par1,par2);
              }
           printf("\r%5.2f", (double)i*fact);
          }
     printf("\nReady!\n");
     /*
      * Save Data
      */
     EX_SaveData(fileName, width, height, field,
	factX, minX, factY, minY);
    }
    else
      printf("Couldn't allocate %d bytes!!\n", width*height*sizeof(double));
   }
