/*
    willamowski_rossler.c, by Timothy Jones, models the Willamowsky-Rossler
    equations in a petri dish.

    Copyright (C) 2010 Timothy Jones

    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/>.
*/

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

#include <pngwriter.h>

// g++ willamowski_rossler.c -lpngwriter -lfreetype -lpng


int WID=480*2;
int HIT=320*2;
int PRAD=300;  //Raius of petri dish, must fit in box

double **mmatrix(int r, int c);
double **mmatrix(int r, int c)
{

int row, col;
double **matrix;
matrix = (double **)malloc((size_t)(r * sizeof(double *)));
//assert(matrix != NULL);

for (row = 0; row < r; row++) 
   {
     matrix[row] = (double *)malloc((size_t)(c * sizeof(double)));
     //  assert(matrix[row] != NULL);
     for (col = 0; col < c; col++)  matrix[row][col] = 0;  //set all to zero
    }
 return matrix;
}


int main (void)
{
  FILE *fp;
  fp=fopen("will_ross_out.dat","w");
//  double *temp;
//  temp=(double *)malloc((size_t) (3*((WID*HIT)))*sizeof(double));
  double dt=0.0005;
  double D=0.001;  //scaled D*dt/(delta_x)^2 = 10 ^ -3 === D
  double **cx;
  double **cy;
  double **cz; 
  cx = mmatrix(WID,HIT);
  cy = mmatrix(WID,HIT);
  cz = mmatrix(WID,HIT);
  double k1=31.2;
  double k_1=0.2;
  double k2=1.572;
  double k_2=0.1;
  double k3=10.8;
  double k_3=0.12;
  double k4=1.02;
  double k_4=0.01;
  double k5=16.5;
  double k_5=0.5;

  int filenum=0;
  char filename[256];
  double overdiff;
  double overdiffa;
  double overdiffb;
  double overdiffc;
  double overdiffd;
  double xmax=-10000;
  double xmin=10000;
  double ymax=-10000;
  double ymin=10000;
  double zmax=-10000;
  double zmin=10000;

  srand(time(NULL));

  

  for(int i=0; i<WID; i++){ 
    for(int j=0; j<HIT; j++){
 	overdiff=sqrt((double) (WID/2-i)*(WID/2-i) + (double) (HIT/2-j)*(HIT/2-j));
	if(( overdiff < PRAD)) {
  
         cx[i][j]=(rand()/RAND_MAX)*1;
         cy[i][j]=(rand()/RAND_MAX)*1;
         cz[i][j]=(rand()/RAND_MAX)*10;
 
	}
        else{
          cx[i][j]=0;
          cy[i][j]=0;
          cz[i][j]=0;
        }
    }}

         cx[WID/2][HIT/2]=0;
         cy[WID/2][HIT/2]=0;
         cz[WID/2][HIT/2]=0;
 
  printf("Done setup, moving on\n");
  
  for(int a=0; a<1500000; a++){

 for(int i=0; i<WID; i++){ 
    for(int j=0; j<HIT; j++){
      
  	overdiff=sqrt((double) (WID/2-i)*(WID/2-i) + (double) (HIT/2-j)*(HIT/2-j));

	if(overdiff < 110){
	  overdiffa=sqrt((double) (WID/2-(i+1))*(WID/2-(i+1)) + (double) (HIT/2-j)*(HIT/2-j));  //Upper right
	  overdiffb=sqrt((double) (WID/2-i)*(WID/2-i) + (double) (HIT/2-(j+1))*(HIT/2-(j+1)));  //Lower right
	  overdiffc=sqrt((double) (WID/2-(i-1))*(WID/2-(i-1)) + (double) (HIT/2-j)*(HIT/2-j));  //Upper left
	  overdiffd=sqrt((double) (WID/2-i)*(WID/2-i) + (double) (HIT/2-(j-1))*(HIT/2-(j-1)));  //Lower left
        }

	if( overdiff < PRAD && overdiffa < PRAD && overdiffb < PRAD && overdiffc < PRAD && overdiffd < PRAD  )  //Well in the middle of dish
     {



  cx[i][j] = cx[i][j] + dt*(k1*cx[i][j] - k_1*cx[i][j]*cx[i][j] -k2*cx[i][j]*cy[i][j] + k_2*cy[i][j]*cy[i][j] - k4*cx[i][j]*cz[i][j] + k_4)    + D*(cx[(i+1 + WID)%WID][j] + cx[(i-1 + WID)%WID][j] -2*cx[i][j] + cx[i][(j+1+HIT)%HIT] + cx[i][(j-1+HIT)%HIT] -2*cx[i][j] );




cy[i][j] = cy[i][j] + dt*( k2*cx[i][j]*cy[i][j] - k_2*cy[i][j]*cy[i][j] -k3 * cy[i][j] + k_3) + 
			  D*(cy[(i+1 + WID)%WID][j] + cy[(i-1 + WID)%WID][j] -2*cy[i][j] + cy[i][(j+1+HIT)%HIT] + cy[i][(j-1+HIT)%HIT] -2*cy[i][j] );
     

 cz[i][j] = cz[i][j] + dt*(-k4*cx[i][j]*cz[i][j] + k_4 + k5*cz[i][j] - k_5*cz[i][j]*cz[i][j]) + 
                          D*(cz[(i+1 + WID)%WID][j] + cz[(i-1 + WID)%WID][j] -2*cz[i][j] + cz[i][(j+1+HIT)%HIT] + cz[i][(j-1+HIT)%HIT] -2*cz[i][j] );



     }

	if(overdiff < PRAD && overdiffa >= PRAD && overdiffd >=PRAD)  // Upper_right
     {

  cx[i][j] = cx[i][j] + dt*(k1*cx[i][j] - k_1*cx[i][j]*cx[i][j] -k2*cx[i][j]*cy[i][j] + k_2*cy[i][j]*cy[i][j] - k4*cx[i][j]*cz[i][j] + k_4)    + D*(cx[(i+1 + WID)%WID][j] -2*cx[i][j] + cx[i][(j+1+HIT)%HIT] -2*cx[i][j] );




cy[i][j] = cy[i][j] + dt*( k2*cx[i][j]*cy[i][j] - k_2*cy[i][j]*cy[i][j] -k3 * cy[i][j] + k_3) + 
			  D*(cy[(i-1 + WID)%WID][j] -2*cy[i][j] + cy[i][(j+1+HIT)%HIT] -2*cy[i][j] );
     

 cz[i][j] = cz[i][j] + dt*(-k4*cx[i][j]*cz[i][j] + k_4 + k5*cz[i][j] - k_5*cz[i][j]*cz[i][j]) + 
                          D*(cz[(i-1 + WID)%WID][j] -2*cz[i][j] + cz[i][(j+1+HIT)%HIT] -2*cz[i][j] );


     }
 
	if(overdiff < PRAD && overdiffa >=PRAD && overdiffb >= PRAD)  // Lower_right
     {
  cx[i][j] = cx[i][j] + dt*(k1*cx[i][j] - k_1*cx[i][j]*cx[i][j] -k2*cx[i][j]*cy[i][j] + k_2*cy[i][j]*cy[i][j] - k4*cx[i][j]*cz[i][j] + k_4)    
             + D*(cx[(i-1 + WID)%WID][j] -2*cx[i][j] + cx[i][(j-1+HIT)%HIT] -2*cx[i][j] );




cy[i][j] = cy[i][j] + dt*( k2*cx[i][j]*cy[i][j] - k_2*cy[i][j]*cy[i][j] -k3 * cy[i][j] + k_3) + 
			  D*(cy[(i-1 + WID)%WID][j] -2*cy[i][j] + cy[i][(j-1+HIT)%HIT] -2*cy[i][j] );
     

 cz[i][j] = cz[i][j] + dt*(-k4*cx[i][j]*cz[i][j] + k_4 + k5*cz[i][j] - k_5*cz[i][j]*cz[i][j]) + 
                          D*(cz[(i-1 + WID)%WID][j] -2*cz[i][j] + cz[i][(j-1+HIT)%HIT] -2*cz[i][j] );


     }

	if(overdiff < PRAD && overdiffc >= PRAD && overdiffd >=PRAD) //Upper_left
     {
cx[i][j] = cx[i][j] + dt*(k1*cx[i][j] - k_1*cx[i][j]*cx[i][j] -k2*cx[i][j]*cy[i][j] + k_2*cy[i][j]*cy[i][j] - k4*cx[i][j]*cz[i][j] + k_4)    + 

  D*(cx[(i+1 + WID)%WID][j] -2*cx[i][j] + cx[i][(j+1+HIT)%HIT] -2*cx[i][j] );




cy[i][j] = cy[i][j] + dt*( k2*cx[i][j]*cy[i][j] - k_2*cy[i][j]*cy[i][j] -k3 * cy[i][j] + k_3) + 
			  D*(cy[(i+1 + WID)%WID][j] -2*cy[i][j] + cy[i][(j+1+HIT)%HIT] -2*cy[i][j] );
     

 cz[i][j] = cz[i][j] + dt*(-k4*cx[i][j]*cz[i][j] + k_4 + k5*cz[i][j] - k_5*cz[i][j]*cz[i][j]) + 
                          D*(cz[(i+1 + WID)%WID][j]  -2*cz[i][j] + cz[i][(j+1+HIT)%HIT]  -2*cz[i][j] );


     }
       
	if(overdiff < PRAD && overdiffc >= PRAD && overdiffb >=PRAD)  // Lower_left
     {

cx[i][j] = cx[i][j] + dt*(k1*cx[i][j] - k_1*cx[i][j]*cx[i][j] -k2*cx[i][j]*cy[i][j] + k_2*cy[i][j]*cy[i][j] - k4*cx[i][j]*cz[i][j] + k_4)    + D*(cx[(i+1 + WID)%WID][j] -2*cx[i][j] + cx[i][(j-1+HIT)%HIT] -2*cx[i][j] );




cy[i][j] = cy[i][j] + dt*( k2*cx[i][j]*cy[i][j] - k_2*cy[i][j]*cy[i][j] -k3 * cy[i][j] + k_3) + 
			  D*(cy[(i+1 + WID)%WID][j] -2*cy[i][j] + cy[i][(j-1+HIT)%HIT] -2*cy[i][j] );
     

 cz[i][j] = cz[i][j] + dt*(-k4*cx[i][j]*cz[i][j] + k_4 + k5*cz[i][j] - k_5*cz[i][j]*cz[i][j]) + 
                          D*(cz[(i+1 + WID)%WID][j]  -2*cz[i][j] + cz[i][(j-1+HIT)%HIT] -2*cz[i][j] );


     }


	if(overdiff < PRAD && overdiffa >= PRAD && overdiffd <PRAD && overdiffb <PRAD)  //Right center side
     {


 cx[i][j] = cx[i][j] + dt*(k1*cx[i][j] - k_1*cx[i][j]*cx[i][j] -k2*cx[i][j]*cy[i][j] + k_2*cy[i][j]*cy[i][j] - k4*cx[i][j]*cz[i][j] + k_4)    
       + D*(cx[(i-1 + WID)%WID][j] -2*cx[i][j] + cx[i][(j+1+HIT)%HIT] + cx[i][(j-1+HIT)%HIT] -2*cx[i][j] );




cy[i][j] = cy[i][j] + dt*( k2*cx[i][j]*cy[i][j] - k_2*cy[i][j]*cy[i][j] -k3 * cy[i][j] + k_3) + 
			  D*(cy[(i-1 + WID)%WID][j] -2*cy[i][j] + cy[i][(j+1+HIT)%HIT] + cy[i][(j-1+HIT)%HIT] -2*cy[i][j] );
     

 cz[i][j] = cz[i][j] + dt*(-k4*cx[i][j]*cz[i][j] + k_4 + k5*cz[i][j] - k_5*cz[i][j]*cz[i][j]) + 
                          D*(cz[(i-1 + WID)%WID][j] -2*cz[i][j] + cz[i][(j+1+HIT)%HIT] + cz[i][(j-1+HIT)%HIT] -2*cz[i][j] );



      }


	if(overdiff < PRAD && overdiffc >= PRAD && overdiffd <PRAD && overdiffb <PRAD)  //Left center side
     {

 cx[i][j] = cx[i][j] + dt*(k1*cx[i][j] - k_1*cx[i][j]*cx[i][j] -k2*cx[i][j]*cy[i][j] + k_2*cy[i][j]*cy[i][j] - k4*cx[i][j]*cz[i][j] + k_4)    + D*(cx[(i+1 + WID)%WID][j] -2*cx[i][j] + cx[i][(j+1+HIT)%HIT] + cx[i][(j-1+HIT)%HIT] -2*cx[i][j] );




cy[i][j] = cy[i][j] + dt*( k2*cx[i][j]*cy[i][j] - k_2*cy[i][j]*cy[i][j] -k3 * cy[i][j] + k_3) + 
			  D*(cy[(i+1 + WID)%WID][j]  -2*cy[i][j] + cy[i][(j+1+HIT)%HIT] + cy[i][(j-1+HIT)%HIT] -2*cy[i][j] );
     

 cz[i][j] = cz[i][j] + dt*(-k4*cx[i][j]*cz[i][j] + k_4 + k5*cz[i][j] - k_5*cz[i][j]*cz[i][j]) + 
                          D*(cz[(i+1 + WID)%WID][j] -2*cz[i][j] + cz[i][(j+1+HIT)%HIT] + cz[i][(j-1+HIT)%HIT] -2*cz[i][j] );

      }


	if(overdiff < PRAD && overdiffd >= PRAD && overdiffc <PRAD && overdiffa <PRAD)  //top center side
     {

 cx[i][j] = cx[i][j] + dt*(k1*cx[i][j] - k_1*cx[i][j]*cx[i][j] -k2*cx[i][j]*cy[i][j] + k_2*cy[i][j]*cy[i][j] - k4*cx[i][j]*cz[i][j] + k_4)    + D*(cx[(i+1 + WID)%WID][j] + cx[(i-1 + WID)%WID][j] -2*cx[i][j] + cx[i][(j+1+HIT)%HIT] -2*cx[i][j] );




cy[i][j] = cy[i][j] + dt*( k2*cx[i][j]*cy[i][j] - k_2*cy[i][j]*cy[i][j] -k3 * cy[i][j] + k_3) + 
			  D*(cy[(i+1 + WID)%WID][j] + cy[(i-1 + WID)%WID][j] -2*cy[i][j] + cy[i][(j+1+HIT)%HIT] -2*cy[i][j] );
     

 cz[i][j] = cz[i][j] + dt*(-k4*cx[i][j]*cz[i][j] + k_4 + k5*cz[i][j] - k_5*cz[i][j]*cz[i][j]) + 
                          D*(cz[(i+1 + WID)%WID][j] + cz[(i-1 + WID)%WID][j] -2*cz[i][j] + cz[i][(j+1+HIT)%HIT] -2*cz[i][j] );

      }


	if(overdiff < PRAD && overdiffb >= PRAD && overdiffc <PRAD && overdiffa <PRAD)  //bot center side
     {

 cx[i][j] = cx[i][j] + dt*(k1*cx[i][j] - k_1*cx[i][j]*cx[i][j] -k2*cx[i][j]*cy[i][j] + k_2*cy[i][j]*cy[i][j] - k4*cx[i][j]*cz[i][j] + k_4)    + D*(cx[(i+1 + WID)%WID][j] + cx[(i-1 + WID)%WID][j] -2*cx[i][j] + cx[i][(j-1+HIT)%HIT] -2*cx[i][j] );




cy[i][j] = cy[i][j] + dt*( k2*cx[i][j]*cy[i][j] - k_2*cy[i][j]*cy[i][j] -k3 * cy[i][j] + k_3) + 
			  D*(cy[(i+1 + WID)%WID][j] + cy[(i-1 + WID)%WID][j] -2*cy[i][j] + cy[i][(j-1+HIT)%HIT] -2*cy[i][j] );
     

 cz[i][j] = cz[i][j] + dt*(-k4*cx[i][j]*cz[i][j] + k_4 + k5*cz[i][j] - k_5*cz[i][j]*cz[i][j]) + 
                          D*(cz[(i+1 + WID)%WID][j] + cz[(i-1 + WID)%WID][j] -2*cz[i][j]  + cz[i][(j-1+HIT)%HIT] -2*cz[i][j] );

      }



 if(cx[i][j]>xmax){xmax=cx[i][j];}
 if(cx[i][j]<xmin){xmin=cx[i][j];}
 if(cy[i][j]>ymax){ymax=cy[i][j];}
 if(cy[i][j]<ymin){ymin=cy[i][j];}
 if(cz[i][j]>zmax){zmax=cz[i][j];}
 if(cz[i][j]<zmin){zmin=cz[i][j];}
    }}


 //}

 //  printf("xMax: %f, xMin: %f, yMax: %f, yMin: %f, zMax: %f, zMin: %f \n",xmax,xmin,ymax,ymin,zmax,zmin);

 if((a%100) == 0){
   filenum=100000000+a;
 sprintf(filename,"%i.png",filenum);
 pngwriter png(WID,HIT,0,filename);
 
 for(int i=0; i<WID; i++){ 
    for(int j=0; j<HIT; j++){

       	overdiff=sqrt((double) (WID/2-i)*(WID/2-i) + (double) (HIT/2-j)*(HIT/2-j));

	if(overdiff < PRAD)
        {
          png.plot(i,j,fabs((cx[i][j]-cx[WID/2][HIT/2])/cx[WID/2][HIT/2]),fabs((cy[i][j]-cy[WID/2][HIT/2])/cy[WID/2][HIT/2]),fabs((cz[i][j]-cz[WID/2][HIT/2])/cz[WID/2][HIT/2]));
	}
    }
    //printf("%f %f %f\n",100*(cx[i][50]-cx[WID/2][HIT/2])/cx[WID/2][HIT/2],100*(cy[i][50]-cy[WID/2][HIT/2])/cy[WID/2][HIT/2],100*(cz[i][50]-cz[WID/2][HIT/2])/cz[WID/2][HIT/2] );
    }
png.close();
 }

 }



return 0;


}