/*
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 .
*/
#include
#include
#include
#include
#include
// 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= 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 && overdiffd = PRAD && overdiffc = PRAD && overdiffc xmax){xmax=cx[i][j];}
if(cx[i][j]ymax){ymax=cy[i][j];}
if(cy[i][j]zmax){zmax=cz[i][j];}
if(cz[i][j]