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

/*#include "simulation.vogl.h"*/

/********************************************************************************/
/*										*/
/*			Global constants					*/
/*										*/
/********************************************************************************/

#define organism 1			/*Species tag*/

#define graphic_output 0		/*0: No graphics*/
					/*1: Grahical display*/

#define location_output	1

#define spp_no 2

#define quality 1.0			/*Used to compute movement rate*/

#define error 0.0000000001		/*Rounding error term*/
#define Pi 3.1415927

#define invasion 0

#define binwidth 0.05

#define bin_max	11

#define file_out "output.results.dat"	/*Results file*/

#define interval 51

#define movement_kernel 4	/***
				1: Exponential kernel
				2: Normal kernel
				3: Uniform kernel
				4: Bivariate Normal
				***/
#define wfoo_kernel 1 		/***
				Conspecific weighting function 
				1:Normal
				2:Flat
				3:Non -JAE paper model
				***/


#define radial_max 21			/**How many bins for co-varaince functions**/
#define realisations 1			/**Number of realisations**/

#define time_lag 0.1
#define C2lag 10000.1

/********************************************************************************/
/*										*/
/*		Global variables, parameters and arrays				*/
/*										*/
/********************************************************************************/

double xmax, xmin;
double ymax, ymin;

double M1[spp_no], M2[spp_no];		/***Movement rate parameters***/

int population_size[spp_no];
int init_pop_size[spp_no];

FILE *realisation_out, *realisation_No;

double birth_sum, death_sum;

double extent = 0.0;	/***Tracks size of numberline***/
double dt=0.0;
double tick = 0.0;	
double C2tick=0.0;
double elapsed = 0.0;	/***Time passed***/
double tmax;		

double mean_move[spp_no];			/***Mean for movement kernel***/
double stdev[spp_no];
double move_max[spp_no];
double conversion[spp_no];

double B1[spp_no];			/***Environment independent and dpedent birth rates***/

double D1[spp_no];
double D2[spp_no][spp_no];			/***Density independent and dependent death rates***/
double alpha[spp_no][spp_no];

double wfoostdv[spp_no][spp_no];		/***Neighbourhood kernel standard deviation***/
double wfoormax[spp_no][spp_no];		/***Neighbourhood kernel max radius of effect***/

double **wtfoo;

int tag=0;	/***Used to tag individuals***/

FILE *C2_file;

/********************************************************************************/
/*										*/
/*			Function declaration					*/
/*										*/
/********************************************************************************/

void birth();
double check_environment();
double drand48();
double get_distance();
void neighbourhood_weights();
void movement();
void recompute_probs();
void death();
void predate();

/********************************************************************************/
/*										*/
/*		Doubly linked list to hold individuals				*/
/*										*/
/********************************************************************************/

struct individual
{
double xpos, ypos;	/*Location in the landscape*/
int species;		/*Species tag*/
double number_line;	/*Flag on the numberline*/
double event_prob;	/*Probability of an event*/

double *n_effects;	/*Array to hold neighbourhood effects*/
double env_effect;

double soil_depth;

double move_event;
double birth_event;
double death_event;
double attack_event;
double eaten_event;

int times_moved;
int serial_no;

double time_in_wood;

struct individual *prior, *next;
};

typedef struct individual element;

	/***Pointers to the start, current, and end individual in list***/
struct individual *start, *end, *neighbour;

#include"births.h"
#include"deaths.h"

/********************************************************************************/
/*										*/
/*			Birth (called to initialise population)			*/
/*										*/
/********************************************************************************/

void birth(spp, xpos, ypos, depth)
int spp;
double xpos, ypos, depth;
{
struct individual *current;
double x1, y1;

	/***If first element in the list***/
if(!start)
	{
	start = (struct individual*) malloc(sizeof(element));
	current = end = start;

	current->prior = NULL;
	current->next = NULL;
	}

	/***Otherwise add new element to end of the list***/
else
	{
	current = end;
	current->next = (struct individual*) malloc(sizeof(element));

	current = current->next;
	current->next = NULL;
	current->prior = end;
	end = current;

	}

	current->species = spp;
	current->xpos = xpos;
	current->ypos = ypos;
	current->time_in_wood=0.0;
	
		tag = population_size[spp];

			/***Calculate the event probability***/
				x1 = xpos;
				y1 = ypos;

				current->serial_no=tag;

				current->birth_event=current->death_event=current->move_event=0.0;
		   
				current->times_moved = 0; /***Records no of time individual moves***/

				current->time_in_wood=0.0;

				current->event_prob=0.0;
		/***Now move it away from it's 'mother'***/
		/***Its neighbourhood effects are computed from here***/
			movement(current);
			birth_weights(current);
}

/********************************************************************************/
/*										*/
/*			Initialise the population				*/
/*										*/
/********************************************************************************/

initialise_population()
{
double xpos, ypos;
int spp;
int i,j, k;
struct individual *current;

population_size[0]=population_size[1]=0;

for(spp=0; spp<2; spp++)
for(i=1; i<=init_pop_size[spp]; i++)
	{
		population_size[spp]++;
		xpos = drand48()*xmax;
		ypos = drand48()*ymax;
		fflush(stdout);	

		birth(spp, xpos, ypos);
	}

current=start;

/***This prints the list to screen***/
/***Useful for checks***/
/*	for(;;)
		{
		printf("spp: %d\n", current->species);

		if(current->next)
			current=current->next;
		else
			break;
		}*/

	/***Now compute the envent probabilities***/
		recompute_probs();

		printf("Initialised the population\n");
		fflush(stdout);
}

/********************************************************************************/
/*										*/
/*			Recompute event probabilities				*/
/*										*/
/********************************************************************************/


	#include "checks.h"

void recompute_probs()
{
struct individual *current;
int census = 0;

	/***Get new neighbourhood weightings***/
		neighbourhood_weights();

	if(end->next != NULL)
		{
		printf("END POINTER IS WRONG\n");
		exit(1);
		}

}


/********************************************************************************/
/*										*/
/*			Birth event						*/
/*										*/
/********************************************************************************/

void newborn(current)
struct individual *current;
{
	/***When a new individual is born it must be placed 'near' to its
		'mother'***/
double xpos, ypos;
int spp;

		spp = current->species;

/*		printf("BIRTH: %d\n", population_size[spp]+1);*/

	/***Give it its mother's position in space***/
		xpos = current->xpos;
		ypos = current->ypos;

	/***Update the population size***/
		population_size[spp]++;

		/***New individual added to the end of the list***/
			birth(spp, xpos, ypos);
	
/*			printf("Mother: x: %lf\ty: %lf\t", xpos, ypos);*/



}

/********************************************************************************/
/*										*/
/*			Read in parameters					*/
/*										*/
/********************************************************************************/

get_params()
{
/*zzzz*/
tmax = 500.0;

xmin=ymin=0.0;
xmax=ymax=sqrt(4.0);

init_pop_size[0] = rint(100.0*xmax*ymax);
init_pop_size[1] = rint(100.0*xmax*ymax);

B1[0]=0.4;
B1[1]=0.0;

D1[0] = 0.1;
D1[1] = 0.05;

D2[0][0]=0.001;	/***Competition***/
D2[0][1]=0.00;	D2[1][1]=0.00;	D2[1][0]=0.00;	

		/***Predation***/
alpha[0][0]=alpha[0][1]=alpha[1][1]=0.0;

alpha[1][0]=0.002;

conversion[0]=0.0;
conversion[1]=0.3;

M2[0]=M2[1]=M1[0]=M1[1]=0.0;

M2[1]=0.00;

M1[0]=0.0;
M1[1]=0.0;

mean_move[0]=0.0;
stdev[0]=0.05;
move_max[0] = 0.45;

mean_move[1]=0.0;
stdev[1]=0.05;
move_max[1] = 0.45;

wfoostdv[0][0] = 0.05;
wfoormax[0][0] = 0.45;

wfoostdv[0][1] = 0.05;
wfoormax[0][1] = 0.45;

wfoostdv[1][1] = 0.05;
wfoormax[1][1] = 0.45;

wfoostdv[1][0] = 0.05;
wfoormax[1][0] = 0.45;
}


/********************************************************************************/
/*										*/
/*			Normal distribution					*/
/*										*/
/********************************************************************************/
double normal()
{
double r1, r2, x1, y, z;

/*do{
	r1 = drand48();
	r2 = drand48();
	y = pow(((2.0 * r1)-1),2.0) + pow(((2.0*r2)-1), 2.0);
}while(y>=1.0);

z = sqrt(-2.0*log(y)/(y));
x1 = ((2.0*r1)-1)*z;*/

/*x1 = fabs(x1);*/

	y = drand48();
	y = -2.0 * log(y);
	z = drand48();
	z = z * 2.0 * Pi;

	x1 = sqrt(y) * cos(z);

return(x1);
}

/********************************************************************************/
/*										*/
/*			Determine the event					*/
/*										*/
/********************************************************************************/

/*	#include "graphics.h"*/
	#include "stochastic.correlations.h"
	#include "predation.h"

void next_event()
{
double pick;
struct individual *current;
double sum=0.0;
double rn;
int spp;
		/***If population goes extinct, write to file then exit***/
/*		for(spp=0;spp<2;spp++)*/
		if(population_size[0] == 0 && population_size[1] == 0)
			{
			fprintf(realisation_No, "%lf\t%lf\n", elapsed, 0.0, 0.0);
			fclose(realisation_No);
			exit(1);
			}

/*		recompute_probs();*/
			
		if(elapsed > tick)
			{
				fprintf(realisation_No, "%lf\t", elapsed);
			for(spp=0;spp<2;spp++)
				{
				fprintf(realisation_No, "%lf\t", (double)population_size[spp]/(xmax*xmax));
				}
				fprintf(realisation_No, "\n");
				fflush(realisation_No);

			if(elapsed>C2tick)
				{
/*				compute_C2();*/

				C2tick+=C2lag;
				}

			tick+=time_lag;
			}

/*	printf("EXTENT:%lf\t", extent);*/

	sum=0.0;
	current=start;
/*	while(current!=NULL)
		{
		current->event_prob =	current->birth_event + current->death_event
					+ current->move_event + current->attack_event;

		sum += current->event_prob;
			
		current=current->next;
		}*/
/*	printf("SUM:%lf\n", sum);*/

	/***Pick a random number from within the number line***/
		pick = drand48()*extent;

		current=start;
		sum=0.0;

	/***Time to next event***/
		elapsed += -log(1.0-drand48())/extent;

	current=start;
	while(current!=NULL)
		{
		current->event_prob =	current->birth_event + current->death_event
					+ current->move_event + current->attack_event;

		sum += current->event_prob;
			
		if(sum>pick)
			break;

		current=current->next;
		}

	/***Now determine which type of event will occur***/

	rn = drand48();

/*	printf("Time: %lf\t", elapsed);*/
/*	printf("D: %lf\tB: %lf\tE: %lf\trn: %lf\t", death_sum, birth_sum, extent, rn);

	fflush(stdout);
*/
/*	check1();*/

	if(current->species==0){
	if(rn < current->birth_event/(current->event_prob))
		{
		newborn(current);

		}
	else if(rn > current->birth_event/(current->event_prob) &&
		rn < (current->birth_event+current->death_event)/(current->event_prob))
		{
		death(current);
		}
	else
		{
		death_weights(current);
		movement(current);
		birth_weights(current);
		}
	}
	else{

	if(rn < current->birth_event/(current->event_prob))
		{
		newborn(current);
		}
	else if(rn > current->birth_event/(current->event_prob) &&
		rn < (current->birth_event+current->death_event)/(current->event_prob))
		{
		death(current);
		}
	else if(rn > (current->birth_event+current->death_event)/(current->event_prob)
		&& rn< (current->birth_event+current->death_event+current->move_event)/(current->event_prob))
		{
		death_weights(current);
		movement(current);
		birth_weights(current);
		}
	else
		{
		predate(current);
		}

	}

/*	printf("No[1]: %d\tNo[2]: %d\n", population_size[0], population_size[1]);*/

/*	if(graphic_output == 1)
		draw_plants();*/

}

/********************************************************************************/
/*										*/
/*			Movement event						*/
/*										*/
/********************************************************************************/

void movement(current)
struct individual *current;
{
double radius, angle;
double x1,y1;
int wrong_in=0;
int wrong_out=0;
int spp;

	spp = current->species;

if(movement_kernel != 4)
do{
if (movement_kernel == 1)
	{
	radius = -log(drand48())*mean_move[spp];	
	}
else if (movement_kernel == 2)
	{
	do{
	radius = mean_move[spp] + fabs(normal())*stdev[spp];
	}while(radius>move_max[spp]);
	}
else if(movement_kernel == 3)
	{
	radius = drand48()*mean_move[spp];
	}
}while(radius>(double)xmax);

/***Move the individual***/

	if(movement_kernel != 4)
		{
		angle = drand48() * 2.0 * Pi;
		current->xpos = current->xpos + (cos(angle) * radius);
		current->ypos = current->ypos + (sin(angle) * radius);
		}
	if(movement_kernel == 4)
		{
		/***Bivariate Normal dispersal kernel***/

			do{
			x1 = mean_move[spp] + normal()*stdev[spp];
			y1 = mean_move[spp] + normal()*stdev[spp];

			radius = sqrt( x1*x1 + y1*y1 );
			}while(radius>move_max[spp]);

			current->xpos += x1;
			current->ypos += y1;
		}

/***Allow for periodic boundaries***/
if(current->xpos >= (double)xmax)
	current->xpos = current->xpos - (double)(xmax + xmin);	/***xpos > xmax***/
if(current->xpos < (double)xmin)
	current->xpos = current->xpos + (double)(xmax - xmin);	/***xpos < xmin***/
if(current->ypos >= (double)ymax)
	current->ypos = current->ypos - (double)(ymax + ymin);	/***ypos > ymax***/
if(current->ypos < (double)ymin)
	current->ypos = current->ypos + (double)(ymax - ymin);	/***ypos < ymin***/

/*printf("nx: %lf\tny: %lf\tdistance: %lf\n", current->xpos, current->ypos, radius);
fflush(stdout);*/

		/***Check individual is still within the arena***/
if(current->xpos > (double)xmax || current->xpos < (double)xmin || current->ypos > (double)ymax || current->ypos < (double)ymin)
	{
	printf("INDIVIDUAL OUT OF ARENA. EXITING.....\n");
	exit(1);
	}
}

/********************************************************************************/
/*										*/
/*				Die die die					*/
/*										*/
/********************************************************************************/

void death(current)
struct individual *current;
{
int num=0;
int spp;
struct individual *temp;

	death_weights(current);

	temp=current;
	spp = current->species;
	population_size[spp] --;

	/***If last element left in list***/
	if(current == end && current == start)
		{
		end = start = NULL;
		}
	/***If end individual is too die***/
	else if(current == end)
		{
		end = current->prior;
		current->prior->next = NULL;
		}
	/***If first in list***/
	else if(current == start)
		{
		current->next->prior == NULL;
		start = current->next;
		}
	/***If mid-list individual to die***/
	else
		{
		current->prior->next = current->next;
		current->next->prior = current->prior;
		}

	free(temp);
	current = start;

/*	while(current!=NULL)
		{
		num++;
		current=current->next;
		}*/
/*	if(num!=population_size[spp])
		{
		printf("EEEEK\n");
		exit(1);
		}*/
}

/********************************************************************************/
/*										*/
/*			Get pair distance					*/
/*										*/
/********************************************************************************/

	/***This function returns the distance between two elements***/
	/***It allows for periodic boundaries			   ***/

double get_distance(x1, y1, x2, y2, radius)
double x1, y1, x2, y2, radius;
{
double distance, d1,d2,d3;
double dx, dy;
double nx1, ny1;

	distance = sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );

	/***Now check to see whether the pair is 'near' the edge of the arena***/

		nx1 = x2;
		ny1 = y2;

	/***If both individuals are near the edge, then move the neighbour over***/
	if(x1-radius<(double)xmin && x2+radius>(double)xmax)
		{
		/***If neighbour near xmax and target near xmin***/
			nx1 = (double)x2-xmax;
		}
	if(x1+radius>(double)xmax && x2-radius<(double)xmin)
		{
		/***If neighbour near xmin and target near xmax***/
			nx1 = (double)xmax+x2;
		}
	if(y1-radius<(double)ymin && y2+radius>(double)ymax)
		{
		/***If neighbour near ymax and target near ymin***/
			ny1 = (double)y2-ymax;
		}
	if(y1+radius>(double)ymax && y2-radius<(double)ymin)
		{
		/***If neighbour near ymin and target near ymax***/
			ny1 = (double)ymax+y2;
		}

		/***Now compute the smallest distance between the two***/
	d1 = sqrt( (x1-x2)*(x1-x2) + (y1-ny1)*(y1-ny1) );
	if(d1 < distance)
		distance = d1;

	d2 = sqrt( (x1-nx1)*(x1-nx1) + (y1-y2)*(y1-y2) );
	if(d2 < distance)
		distance = d2;

	d3 = sqrt( (x1-nx1)*(x1-nx1) + (y1-ny1)*(y1-ny1) );
	if(d3 < distance)
		distance = d3;

	return(distance);
}

/********************************************************************************/
/*										*/
/*	This function computes the neighbourhood weights for events		*/
/*										*/
/********************************************************************************/

void neighbourhood_weights()
{
struct individual *focal, *neighbour;
double x1,y1,x2,y2;
int i,j;
int num=0;
double effect, norm;
double distance;
double weighting=0.0;
double sum=0.0;
int m,n;
FILE *weight;
int spp1, spp2;

	weight = fopen("weight.test.dat", "w");

	/***Set pointer to beginning of list***/
		focal = start;

		birth_sum = death_sum = extent = 0.0;
		num = 0;
	while(focal!=NULL)
		{
		num++;
		weighting=0.0;
		x1 = focal->xpos;
		y1 = focal->ypos;

		spp1 = focal->species;

		focal->death_event = D1[spp1];
		focal->birth_event = B1[spp1];
		focal->move_event = M1[spp1];
		focal->attack_event=0.0;

		neighbour = start;

			{
			while(neighbour!=NULL)
				{
				x2 = neighbour->xpos;
				y2 = neighbour->ypos;
				spp2 = neighbour->species;

				/***Compute the pair distance allowing for periodic boundaries***/
					distance = get_distance(x1,y1, x2,y2, wfoormax[spp1][spp2]);

			/***Compute neighbourhood effect on death***/
				if(neighbour != focal)
				if(distance<=wfoormax[spp1][spp2])
					{
					num++;
				/***Bivariate normal competition kernel***/
					if(wfoo_kernel == 1)
						{
						norm = 2.0 * Pi * wfoostdv[spp1][spp2] * wfoostdv[spp1][spp2];

/*			norm = 2.0 * Pi * wfoostdv[spp1][spp2] * wfoostdv[spp1][spp2] * (1.0 - exp(-wfoormax[spp1][spp2]*wfoormax[spp1][spp2]/(2.0*wfoostdv[spp1][spp2]*wfoostdv[spp1][spp2])));*/

						effect = distance/wfoostdv[spp1][spp2];
						effect = effect * effect;
						
						effect = exp(-0.5 * effect);

						effect = effect/norm;

						focal->attack_event+= effect*alpha[spp1][spp2];
						}
				/***Flat kernel***/
					else if(wfoo_kernel == 2)
						{
						effect = 1.0/(Pi * wfoormax[spp1][spp2] * wfoormax[spp1][spp2]);
						}

					effect *=D2[spp1][spp2];
					
					weighting += effect;
					sum+=effect;

/*					fprintf(weight,"x1: %lf y1: %lf x2: %lf y2: %lf distance: %lf effect: %lf\n", x1,y1, x2,y2, distance, effect);
					fflush(weight);*/
					}
				
				neighbour = neighbour->next;
				}

				focal->death_event += weighting;
				/*printf("%lf\n", focal->death_event);*/
				death_sum += focal->death_event;
				extent += focal->death_event;

				birth_sum += focal->birth_event;
				extent += focal->birth_event;

				extent+=focal->attack_event;
				extent+=focal->move_event;
			}
		focal = focal->next;

		}

/*	printf("Sum: %lf\n", death_sum/(double)population_size[0]);
	exit(1);*/
}



/********************************************************************************/
/*										*/
/*			Main (nice and small as usual)				*/
/*										*/
/********************************************************************************/
	#include "spatial.h"

int main(argv, argc)
int argv;
char *argc[];
{
int seed;
int i, k;
FILE *test;
struct individual *current;
double mean =0.0;
int number;
char filename[80];
int new_pop;
double xpos, ypos;
int flag =0;

	test = fopen ("test.dat", "w");

	seed = atoi(argc[2]);

	get_params();

	M1[0]= atof(argc[3]);

	srand48(seed);

	initialise_population();
	recompute_probs();

/*	if(graphic_output==1)
		{
		initialize_vogl();
		doublebuffer();
		}*/

	i=atoi(argc[1]);

	/***Reset the clock***/
		elapsed = 0.0;

	k = 0;

	sprintf(filename, "No.output.%d.results.dat", i);
	realisation_No = fopen(filename, "w");

	sprintf(filename, "C2.output.%d.results.dat", i);
	C2_file = fopen(filename, "w");

/*	check2();*/

do{	
	number = 0;

		/***Get individual from number line***/
			next_event();
/***Invasion of second species, when first species at stochastic equilibrium***/
if(invasion==1)
	if(flag==0)
		{
		if(floor(elapsed)==100)
			{
			printf("INVASION OF NEW SPECIES\n");
			sleep(2);
		/***Now make 15 new individuals***/
		/***Newborn are randomly positioned in the arena***/
			for(new_pop=0;new_pop<10;new_pop++)
				{
				xpos=drand48()*xmax;
				ypos=drand48()*ymax;

				birth(1, xpos, ypos);
				population_size[1]++;
				}
			flag=1;
			}
		}
		
}while(elapsed<=(tmax+0.1));


/***Write final locations to file***/	
	if(location_output == 1)
		write_locations(i);

fprintf(test, "%d\n", number);
fflush(test);

fflush(stdout);
/*sleep(10);*/

return(1);

}