#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/// MPI include file here

/* global objects: */

int* iters;
double XYinc[2] = {0.0,0.0}; // Distance between pixels
double SWV[2] = {0.0,0.0}; //  Beginning coordinates: real and imaginery
int XYdots  =  0;
int Niter  = 1000;
int Procs = 1;



int MINVAL(int s, int*a)
{
   int v, e;

   v = a[0];

   for ( e = 0; e < s; e++ )
   {
       if ( v > a[e] ) v = a[e];
   
   }
   
   return(v);
}

int MAXVAL(int s, int*a)
{
   int v, e;

   v = a[0];

   for ( e = 0; e < s; e++ )
   {
       if ( v < a[e] ) v = a[e];
   
   }
   
   return(v);
}
   
     

void calcola( int X0,  int X1,  int Y0,  int Y1)
{

   int pos;
   int       IX, IY, IZ;
   double     CA, CB;
   double     ZA, ZB;
   double     ZSIZE ;
   double     ZNEWA, ZNEWB;

   for ( IX = X0; IX < X1; IX++ ) {
     /* New column  */
      CA = XYinc[0] * IX + SWV[0] ;
      for ( IY = Y0; IY < Y1; IY++ ) {
         ZA = 0.0  ;
         ZB = 0.0  ;
         CB = XYinc[1] * IY + SWV[1]  ;
         ZSIZE = (CA*CA + CB*CB )  ;
         /* Iteration loop for each ( Z**2 + C )   */
         for ( IZ = 0; IZ <= Niter ; IZ++ ) {
            if ( ZSIZE >= 4 ) break; 
            ZNEWA = ( ZA*ZA ) - ( ZB*ZB ) + CA  ;
            ZNEWB = ( 2*ZA*ZB ) + CB  ;
            ZSIZE = ( ZNEWA*ZNEWA + ZNEWB*ZNEWB )  ;
            ZA   = ZNEWA  ;
            ZB   = ZNEWB  ;
         }  
         /* Results  */
         iters[IY*XYdots+IX] = IZ;
      } 
   }  
}

void IntData2pgm(int s1, int s2, int* idata, char* name) 
{
   /*
      Simple subroutine to dump integer data in a PGM format
   */
   int i, j;
   FILE* ouni;
   int  vp, vs;
   float  rmin, rmax;
   char  fname[80];
   /*
      Write on unit 700 with  PGM format
   */
   strcpy(fname,name);
   strcat(fname,".pgm\0");
   ouni = fopen(fname,"w");
   if ( ! ouni ) {
      fprintf(stderr,"!!!! Error write access to file %s\n",fname);
   }
   /*  Magic code */
   fprintf(ouni,"P2\n");
   /*  Dimensions */
   fprintf(ouni,"%d %d\n",s1, s2);
   /*  Maximum value */
   fprintf(ouni,"255\n");
   /*  Values from 0 to 255 */
   rmin = MINVAL(s1*s2,idata); rmax = MAXVAL(s1*s2,idata);
   vs = 0;
   for ( j = 0; j < s2; j++ ) {
      for ( i = 0; i < s1; i++ ) {
         vp = (int) ( (idata[i+(j*s1)] - rmin) * 255.0 / (rmax - rmin) );
         vs++;
         fprintf(ouni,"%4d ", vp);
         if ( vs >= 10 ) {
            fprintf(ouni," \n");
            vs = 0;
         }
      }
      fprintf(ouni," ");
      vs = 0;
   }
   fclose(ouni);
   return;
}

void IntData2ppm(int s1, int s2, int* idata, char* name) 
{
   /*
      Simple subroutine to dump integer data in a PPM format
   */
   int i, j;
   FILE* ouni;
   int  vp, vs, c, r, g, b;
   float  fvp, rmin, rmax;
   char  fname[80];

   strcpy(fname,name);
   strcat(fname,".ppm\0");
   ouni = fopen(fname,"w");
   if ( ! ouni ) {
      fprintf(stderr,"!!!! Error write access to file %s\n",fname);
   }
   /*  Magic code */
   fprintf(ouni,"P3\n");
   /*  Dimensions */
   fprintf(ouni,"%d %d\n",s1, s2);
   /*  Maximum value */
   fprintf(ouni,"255\n");
   /*  Values from 0 to 255 */
   rmin = MINVAL(s1*s2,idata); rmax = MAXVAL(s1*s2,idata);
   vs = 0;
   for ( j = 0; j < s2; j++ ) {
      for ( i = 0; i < s1; i++ ) {
         fvp = ( (idata[i+(j*s1)] - rmin) * 255.0 / (rmax - rmin) );
         vp = (int)fvp % 256;
         c = vp / 42;   /* colour */
         switch (c) {
            case 0:
               vp = (int) (fvp - 1.0) * 255.0 / 42.0;
               r = vp; g = 0; b = 0;
               break;
            case 1:
               vp = (int) (fvp - 42.0) * 255.0 / 42.0;
               r = vp; g = vp; b = 0;
               break;
            case 2:
               vp = (int) (fvp - 84.0) * 255.0 / 42.0;
               r = 0; g = vp; b = 0;
               break;
            case 3:
               vp = (int) (fvp - 126.0) * 255.0 / 42.0;
               r = 0; g = vp; b = vp;
               break;
            case 4:
               vp = (int) (fvp - 168.0) * 255.0 / 42.0;
               r = 0; g = 0; b = vp;
               break;
            case 5:
               vp = (int) (fvp - 210.0) * 255.0 / 42.0;
               r = vp; g = 0; b = vp;
               break;
            case 6:
               r = 0; g = 0; b = fvp/10.0;
               break;
            default:
               r = 0; g = 0; b = 0;
         }
         vs++;
         fprintf(ouni,"%4d %4d %4d ", r,g,b);
         if ( vs >= 10 ) {
            fprintf(ouni," \n");
            vs = 0;
         }
      }
      fprintf(ouni," ");
      vs = 0;
   }
   fclose(ouni);
   return;
}

/*

   MAND - Compute the points outside the Mandelbrot set
 
    Computes the iteration number for each pixel in a given squared area
 
    Syntax: 
    mandel Sreal Simag XYdots Range procs Niter 
               
  INPUT:            
     REAL(8) Sreal = Real coord. of lower left vertex 
     REAL(8) Simag = Imaginary coord. of lower left vertex 
     INTEGER XYdots = Horizontal and vertical pixels 
     INTEGER Range = Square area side length 
     INTEGER procs = Independent computing units 
     INTEGER Niter = Maximum iterations
               
 ************************************************************************** */

int main( int argc, char *argv[])    /* Mandel */

{

   /*  Mandel variables */
   double Sreal, Simag, Range = 0.0;  // Beginning coords and zoom factor
   double tempo2, tempo1, tempro;
   char  fname[40];

   int dX, dY, t, X0, X1, Y0, Y1;
   int st;
   FILE* inp_unit;
   
   /*  MPI variables */
   int ierr, my_rank, numprocs;
   int LetsGo, rc;

   /// MPI initialization: define my_rank and numprocs
   my_rank = 0; numprocs = 1;
      
   LetsGo = 1;
   if ( my_rank == 0 ) {
   	   /* Master process looks for input data */

      strcpy(fname,"mandel.inp");
      inp_unit = fopen(fname,"r");
      if ( ! inp_unit ) {
         fprintf(stderr,"!!!! Error read access to file %s\n",fname);
      }

      ierr = fscanf(inp_unit,"%lf, %lf, %d, %lf, %d, %d\n",&Sreal, &Simag, &XYdots, &Range, &Procs, &Niter);
      if ( ierr != 6 ) {
   	     fprintf(stderr, "Error reading input data. Examples are:\n");
   	     fprintf(stderr, "-2.000, -2.000, 1000, 4.000, 1, 1000\n" );
   	     fprintf(stderr, "-1.500, -0.500, 1000, 0.500, 1, 1000\n" );
   	     fprintf(stderr, "-1.250, -0.200, 1000, 0.100, 1, 1000\n" );
   	     fprintf(stderr, "-1.250, -0.180, 1000, 0.030, 1, 1000\n" );
   	     fprintf(stderr, "-1.235, -0.171, 1000, 0.006, 1, 1000\n" );
   	     fprintf(stderr, "-1.233, -0.169, 1000, 0.001, 1, 1000\n" );
   	     fprintf(stderr, "-1.233, -0.169, 1000, 0.0001, 1, 1000\n" );
   	     LetsGo = 0;
      }
      SWV[0] = Sreal; SWV[1] = Simag;
      fprintf(stdout,"%lf, %lf, %d, %lf, %d, %d\n",SWV[0], SWV[1], XYdots, Range, Procs, Niter);
      fprintf(stdout," ierr = %d\n",ierr);

      if ( LetsGo ) {
         /*  Master allocates entire matrix */
         iters = (int*) malloc(4*XYdots*XYdots);
         if ( iters == NULL ) {
   	        fprintf(stderr, "Error allocating ITERS(%d,%d)\n",XYdots,XYdots);
      	     LetsGo = 0;
         }
      }
   }
   /// Send start flag to all processes
   if ( ! LetsGo ) return(-1);
      /// Send data to processes: SWV, XYdots, Range, Niter */

   fprintf(stderr, "Input data are: %lf %lf %d %lf %d %d\n", 
              SWV[0], SWV[1], XYdots, Range, Procs, Niter);

   /* Every process executes */
   Procs = numprocs;   /* procs data required for compatibility reasons */

   /* ---> Compute complex spacing between pixels <---  */
   XYinc[0] = Range / (float)XYdots;
   XYinc[1] = Range / (float)XYdots;
  
     dY = ( XYdots + (Procs - 1) ) / Procs;
       
       X0 = 0; X1 = XYdots;

       Y0 = dY * my_rank;
       Y1 = dY * ( my_rank + 1);
       if ( Y1 > XYdots ) Y1 = XYdots;
        
       if ( my_rank != 0 ) {   	   
         /*  Processes allocate local matrix */
         iters = (int*) malloc(4*XYdots*XYdots);
         if ( iters == NULL ) {
   	        fprintf(stderr, "Error allocating ITERS(%d,%d)\n",XYdots,XYdots);
         }
       }
       calcola(X0, X1, Y0, Y1);
	
   if ( my_rank == 0 ) {
      /* Master process gathers results */

      IntData2pgm(XYdots,XYdots,iters,"mandel");
      IntData2ppm(XYdots,XYdots,iters,"mandel");
   	   	   
   }
   /// Terminate MPI processes
   return(0);
}