#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifdef _OPENMP
#include <omp.h>
#endif
/* global objects: */

int *iters, *pixel;
double XINC   = 0.0, YINC   = 0.0;
double Sreal  = 0.0, Simag  = 0.0;
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 = XINC * IX + Sreal ;
      for ( IY = Y0; IY < Y1; IY++ ) {
         ZA = 0.0  ;
         ZB = 0.0  ;
         CB = YINC * IY + Simag  ;
         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 IntData2pixel(int s1, int s2, int* idata, int* pixel) 
{
   /*
      Simple subroutine to transform integer data into normilized [0-255] int
      values
   */
   int i, j, my_id, i0, i1, di;
   FILE* ouni;
   int  vp, vs;
   float  rmin, rmax;

      
   #ifdef _OPENMP
      my_id = omp_get_thread_num();
   #else
      my_id = 0;
   #endif

   /*  Normalize values from 0 to 255 */
   rmin = MINVAL(s1*s2,idata); rmax = MAXVAL(s1*s2,idata);
   
   di = (s1 + Procs -1) / Procs;
   i0 = di * my_id;
   i1 = di * (my_id+1);
   if ( i1 > s1 ) i1 = s1;
   
   for ( i = i0; i < i1; i++ ) {
      for ( j = 0; j < s2; j++ ) {
         pixel[i+(j*s1)] = (int) ( (idata[i+(j*s1)] - rmin) * 255.0 / (rmax - rmin) );
      }
   }
   return;
}

void Pixel2pgm(int s1, int s2, int* pixel, char* name) 
{
   /*
      Simple subroutine to dump normalized [0-255] data in a PGM format
   */
   int i, j, my_id, i0, i1, di;
   FILE* ouni;
   int  vp, vs;
   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 already must be from 0 to 255 */
   vs = 0;
   for ( i = 0; i < s1; i++ ) {
      for ( j = 0; j < s2; j++ ) {
         vs++;
       	 fprintf(ouni,"%4d ", pixel[i+(j*s1)]);
       	 if ( vs >= 10 ) {
       	 	 fprintf(ouni," \n");
       	 	 vs = 0;
       	 }
      }
      fprintf(ouni," ");
      vs = 0;
   }
   fclose(ouni);
   return;
}

void Pixel2ppm(int s1, int s2, int* pixel, char* name) 
{
   /*
      Simple subroutine to dump normalized [0-255] data in a PPM format
   */
   int i, j;
   FILE* ouni;
   int  vp, vs, c, r, g, b;
   float  fvp;
   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 already must be from 0 to 255 */
   vs = 0;
   for ( i = 0; i < s1; i++ ) {
      for ( j = 0; j < s2; j++ ) {
         fvp = pixel[i+(j*s1)];
         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 Range = 0.0;  
   double tempo2=0.0, tempo1=0.0, tempro;
   char  fname[40];

   int dX, dY, t, X0, X1, Y0, Y1;
   int st, ierr, numprocs=1, my_id=0;
   FILE* inp_unit;
   
   // 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);
      fprintf(stdout,"%lf, %lf, %d, %lf, %d, %d\n",Sreal, Simag, XYdots, Range, Procs, NITER);
      if (ierr != 6 ) fprintf(stdout," ierr = %d\n",ierr);
/*
      if ( ierr != 0 ) {
   	     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;
      }
*/

         iters = (int*) malloc(4*XYdots*XYdots);
         if ( iters == NULL ) {
   	        fprintf(stderr, "Error allocating ITERS(%d,%d)\n",XYdots,XYdots);
   	     }
   	     pixel = (int*) malloc(4*XYdots*XYdots);
         if ( iters == NULL ) {
   	        fprintf(stderr, "Error allocating PIXEL(%d,%d)\n",XYdots,XYdots);
   	     }
   fprintf(stderr, "Running OpenMP Mandel. Input data are: %lf %lf %d %lf %d %d\n", 
              Sreal, Simag, XYdots, Range, Procs, NITER);
    
   
   /* Every process executes */
   // Insert parallel region, pay attention to  shared or private variables
   
   // The two following statements are executed only if openmp is defined
   //numprocs = omp_get_num_threads();
   //my_id = omp_get_thread_num();
   
   Procs = numprocs;   /* procs data required for compatibility reasons */
   
   // This statement is executed only by master thread if openmp is defined
   //    tempo1 = omp_get_wtime(); 
      

   /* ---> Compute complex spacing between pixels <---  */
   XINC = Range / (float)XYdots;
   YINC = Range / (float)XYdots;
  
/*  Let's decompose on 2nd dimension  */
/*   dX = ( XYdots + (Procs - 1) ) / Procs */
     dY = ( XYdots + (Procs - 1) ) / Procs;
       
       X0 = 0; X1 = XYdots;

       Y0 = dY * my_id;
       Y1 = dY * ( my_id + 1);
       if ( Y1 > XYdots ) Y1 = XYdots;
        
      
       calcola(X0, X1, Y0, Y1);
       // Insert a barrier to syncronize all threads
	   
	          
       // This statement is executed only by master thread if openmp is defined
       //tempo2 = omp_get_wtime(); 
       
       
       IntData2pixel(XYdots,XYdots,iters,pixel);
       // Insert a barrier to syncronize all threads
       
       // The following statements are executed only by mater thread
       	   // Pixel2pgm(XYdots,XYdots,pixel,"mandel");
       	   Pixel2ppm(XYdots,XYdots,pixel,"mandel");
       	   tempro =  (tempo2 - tempo1);
       	   printf("Elapsed time is %lf \n",tempro);
       
   return(0);
}