/***

	Predator-prey model in continuous space

	Option to have prey either density dependent
	or density dependent population growth
	in the absence of the predator.

***/

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

/****************************************************************/
/*                                                              */
/*                GLOBAL VARIABLES + CONSTANTS                  */
/*                                                              */
/****************************************************************/
#define max_bin 20
	/*This is the radius of the second moment we
			are tracking*/

#define out_file "numerical2.rcfs.dat"
#define No_file	 "numerical2.No.dat"
#define density_file "results.dat"

#define o 0
#define e 1

	/***These are the 3 corner weightings for the power 2 closure***/
		#define alpha	4.0	/***i corner weighting***/
		#define beta	1.0	/***j corner weighting***/
		#define Gamma	1.0	/***k corner weighting***/

#define resolution 0

#define check_moments	1
#define nspecies	2

#define closure_method 2	/*2: Power 2, symmetric*/
				/*3: Power 3*/
				/*1: Power1, symmetric*/
			
#define movement_kernel 4	/*1:Normal distribution*/
				/*2:Exponential distribution*/
				/*3:Uniform distribution*/
				/*4:Bivariate Normal*/

#define weight_kernel 1 	/*1: Normal kernel*/
			
#define xbin 20
			/*These are the number of bins within the*/
#define ybin 20		/*second moment grid*/

#define Pi 3.1415927

#define tmax 25000.0
#define bin_width 0.05	/*Equal to max_bin/xbin*/

double move_stdv[nspecies] = {0.05, 0.05};
double rmax[nspecies] = {0.15, 0.15};

double W_stdv[nspecies][nspecies];

#define radial_max	4
#define rad	4

double N[nspecies] = {100.0, 100.0};		/*Stores first moment*/

FILE *data, *out10, *first;

double C2_bin[nspecies][nspecies][2*(xbin-1)+1] [2*(ybin-1)+1];	/*Tracks the cross-correlation function at each distance*/
double C2_rhs[nspecies][nspecies][2*(xbin-1)+1] [2*(ybin-1)+1];	/*Used to compute the rhs at each time step*/

double prob_move [nspecies] [2*(radial_max-1)+1] [2*(radial_max-1)+1];		/***Movement/dispersal kernel***/
double Wc [nspecies] [nspecies] [2*(radial_max-1)+1] [2*(radial_max-1)+1];	/***Neighbour weighting function***/

double B1[nspecies] = {0.4, 0.0};			
double D1[nspecies] = {0.1, 0.05};	/***Density independent death rate***/
double D2[nspecies][nspecies];		/***Density dependent death rate***/

double conv[nspecies][nspecies];
double attack[nspecies][nspecies];

double M1[nspecies] = {0.0, 0.0};

double mean_move = 0.0;		/*This is the mean length per movement event 
				where the distribution is normal*/

#define dt 1.0/20.0			/*This is the time increment for each step
				of the integration*/

double equilibrium= 0.001*dt;

double elapsed;
double tock;
FILE *temp_out, *detail, *eq_m;
int C2flag=0;

double moment_closure();
void write_moments();
double dndt[nspecies];
void write_eq();

/****************************************************************/
/*								*/
/*		Initialise Second moment bins			*/
/*								*/
/****************************************************************/

/*This needs to be called once at the beginning of the program*/

void init_C2()
{
int i, j, spp1, spp2;

attack[0][0]=0.0;
attack[0][1]=0.002;
attack[1][0]=0.00;
attack[1][1]=0.0;

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

D2[0][0]=0.001;
D2[0][1]=0.00;	D2[1][1]=0.00;	D2[1][0]=0.00;

W_stdv[0][0]=0.05;
W_stdv[0][1]=0.05;	W_stdv[1][1]=0.05;	W_stdv[1][0]=0.05;

for(i=0; i<2*(xbin-1)+1; i++)
for(j=0; j<2*(ybin-1)+1; j++)
for(spp1=0;spp1<nspecies;spp1++)
for(spp2=0;spp2<nspecies;spp2++)
	{
	C2_bin[spp1][spp2][i][j] = N[spp1]*N[spp2];
	}
}

/****************************************************************/
/*                                                              */
/*           COMPUTE RHS OF MOMENT EQUATIONS                    */
/*                                                              */
/****************************************************************/
	/*Don't forget to convert movement/neighbourhood kernels	*/
	/*to probability from pdf when integrating round an individual	*/
void calc_cross_rhs()
{
int i,j,k,spp1,spp2,m,n, p,q;
int W_x, W_y;
int invi, invj;
int invm, invn;
int move_x, move_y;
int end = 2*(radial_max-1)+1;
int endC2 = 2*(xbin-1)+1;

for(i=0; i<endC2; i++)
	{
/*		printf(".");
		fflush(stdout);*/

	for(j=0; j<endC2; j++)
		{
		for(spp1=0;spp1<nspecies;spp1++)
		for(spp2=0;spp2<nspecies;spp2++)
			{
			C2_rhs[spp1][spp2][i][j]=0.0;
			invi = endC2-1 -i;	/**The extra -1 is required because of the indexing**/
			invj = endC2-1 -j;	/**which goes from 0...(2*bin-1)-1**/
			
		/**************************************/
		/*	Density independent death     */
		/**************************************/
			C2_rhs[spp1][spp2][i][j] -= D1[spp1] * C2_bin[spp1][spp2][i][j];
			C2_rhs[spp1][spp2][i][j] -= D1[spp2] * C2_bin[spp2][spp1][i][j];

		/*********************************************/
		/*	Density independent movement out     */
		/*********************************************/
			C2_rhs[spp1][spp2][i][j] -= M1[spp1] * C2_bin[spp1][spp2][i][j];
			C2_rhs[spp1][spp2][i][j] -= M1[spp2] * C2_bin[spp2][spp1][i][j];

		/**************************************/
		/*		Corrections	      */
		/**************************************/
			/**************************************/
			/*	Death correction terms	      */
			/**************************************/
			if(abs(i-(xbin-1))<radial_max-1 && abs(j-(xbin-1))<radial_max-1)
				{
				W_x = i - (xbin-1) + radial_max-1;
				W_y = j - (ybin-1) + radial_max-1;
				C2_rhs[spp1][spp2][i][j] -= D2[spp1][spp2] * Wc[spp1][spp2][W_x][W_y] * C2_bin[spp1][spp2][i][j];

		/*******************************/
		/***Predation correction term***/
		/*******************************/
				C2_rhs[spp1][spp2][i][j] -= attack[spp1][spp2] * Wc[spp1][spp2][W_x][W_y] * C2_bin[spp1][spp2][i][j];

/*				}
			if(abs(invi-(xbin-1))<radial_max-1 && abs(invj-(xbin-1))<radial_max-1)
				{*/
				W_x = invi - (xbin-1) + radial_max-1;
				W_y = invj - (ybin-1) + radial_max-1;
				C2_rhs[spp1][spp2][i][j] -= D2[spp2][spp1] * Wc[spp2][spp1][W_x][W_y] * C2_bin[spp2][spp1][i][j];

		/*******************************/
		/***Predation correction term***/
		/*******************************/
				C2_rhs[spp1][spp2][i][j] -= attack[spp2][spp1] * Wc[spp2][spp1][W_x][W_y] * C2_bin[spp2][spp1][i][j];
				}
			/**************************************/
			/*	Birth correction terms        */
			/**************************************/
			if(spp1 == spp2)
				{
				if(abs(i-(xbin-1))<radial_max-1 && abs(j-(xbin-1))<radial_max-1)
					{
					move_x = i - (xbin-1) + radial_max-1;			move_y = j - (ybin-1) + radial_max-1;
					C2_rhs[spp1][spp2][i][j] += B1[spp2] * prob_move[spp2][move_x][move_y]*N[spp2];
					C2_rhs[spp1][spp2][i][j] += B1[spp1] * prob_move[spp1][move_x][move_y]*N[spp1];
					}
				}

			for(m=i-(radial_max-1); m<i+radial_max-1; m++)
			for(n=j-(radial_max-1); n<j+radial_max-1; n++)
				{
				invm = (m-i) + invi;		invn = (n-j) + invj;
				
			/**************************************/
			/*	Density independent births    */
			/**************************************/
				move_x = m - i + radial_max-1;
				move_y = n - j + radial_max-1;
		
/*			if(move_x>=0 && move_y>=0 && move_x<end && move_y<end)*/
					{

					if(m>=0 && m<endC2 && n>=0 && n<endC2)
						{

						C2_rhs[spp1][spp2][i][j] += B1[spp1] * prob_move[spp1][move_x][move_y] * C2_bin[spp1][spp2][m][n]*bin_width*bin_width;
						C2_rhs[spp1][spp2][i][j] += M1[spp1] * prob_move[spp1][move_x][move_y] * C2_bin[spp1][spp2][m][n]*bin_width*bin_width;

						C2_rhs[spp1][spp2][i][j] += B1[spp2] * prob_move[spp2][move_x][move_y] * C2_bin[spp2][spp1][m][n]*bin_width*bin_width;
						C2_rhs[spp1][spp2][i][j] += M1[spp2] * prob_move[spp2][move_x][move_y] * C2_bin[spp2][spp1][m][n]*bin_width*bin_width;

						}
					else
						{

						C2_rhs[spp1][spp2][i][j] += B1[spp1] * prob_move[spp1][move_x][move_y] * N[spp1]*N[spp2]*bin_width*bin_width;
						C2_rhs[spp1][spp2][i][j] += M1[spp1] * prob_move[spp1][move_x][move_y] * N[spp1]*N[spp2]*bin_width*bin_width;

						C2_rhs[spp1][spp2][i][j] += B1[spp2] * prob_move[spp2][move_x][move_y] * N[spp2]*N[spp1]*bin_width*bin_width;
						C2_rhs[spp1][spp2][i][j] += M1[spp2] * prob_move[spp2][move_x][move_y] * N[spp2]*N[spp1]*bin_width*bin_width;

						}

			/********************************/
			/***Predation correction terms***/
			/********************************/
					W_x = m-(xbin-1)+radial_max-1;
					W_y = n-(ybin-1)+radial_max-1;

					if(fabs(m-(xbin-1))<radial_max-1)
					if(fabs(n-(ybin-1))<radial_max-1){

						if(attack[spp2][spp1]>0.000001)
						C2_rhs[spp1][spp2][i][j] += conv[spp2][spp1]*attack[spp2][spp1] * Wc[spp2][spp1][W_x][W_y] * prob_move[spp1][move_x][move_y] * C2_bin[spp2][spp1][m][n]*bin_width*bin_width;

						if(attack[spp1][spp2]>0.000001)
						C2_rhs[spp1][spp2][i][j] += conv[spp1][spp2]*attack[spp1][spp2] * Wc[spp1][spp2][W_x][W_y] * prob_move[spp2][move_x][move_y] * C2_bin[spp1][spp2][m][n]*bin_width*bin_width;
						}

					}

			/**************************************/
			/*	Density dependent deaths      */
			/**************************************/
			if(abs(m-i)<radial_max-1 && abs(n-j)<radial_max-1)
				for(k=0; k<nspecies; k++)
					{
					W_x = (m-i) + radial_max-1;
					W_y = (n-j) + radial_max-1;

			/****Competition****/
					if(D2[spp1][k]>0.000001)
						C2_rhs[spp1][spp2][i][j] -= D2[spp1][k] * Wc[spp1][k][W_x][W_y] * moment_closure(spp1, spp2, k, i, j, (m-i+(xbin-1)), (n-j+(ybin-1)))*bin_width*bin_width;
					if(D2[spp2][k]>0.000001)
						C2_rhs[spp1][spp2][i][j] -= D2[spp2][k] * Wc[spp2][k][W_x][W_y] * moment_closure(spp2, spp1, k, i, j, (m-i+(xbin-1)), (n-j+(ybin-1)))*bin_width*bin_width;

			/****Predation****/
					if(attack[spp1][k]>0.000001)
						C2_rhs[spp1][spp2][i][j] -= attack[spp1][k] * Wc[spp1][k][W_x][W_y] * moment_closure(spp1, spp2, k, i, j, (m-i+(xbin-1)), (n-j+(ybin-1)))*bin_width*bin_width;
					if(attack[spp2][k]>0.000001)
						C2_rhs[spp1][spp2][i][j] -= attack[spp2][k] * Wc[spp2][k][W_x][W_y] * moment_closure(spp2, spp1, k, i, j, (m-i+(xbin-1)), (n-j+(ybin-1)))*bin_width*bin_width;

/***Creation of predator parent-offspring pairs***/
			/***Corrections***/
					if(spp1==spp2){
					move_x = i - (xbin-1) + radial_max-1;
					move_y = j - (ybin-1) + radial_max-1;

					if(move_x>=0 && move_x<2*radial_max-1 && move_y>=0 && move_y <2*radial_max-1){
							{
							C2_rhs[spp1][spp2][i][j]+=attack[k][spp1]*conv[k][spp1]*prob_move[spp1][move_x][move_y] * Wc[k][spp1][W_x][W_y] * C2_bin[k][spp1][m-i+xbin-1][n-j+ybin-1]*bin_width*bin_width;
							C2_rhs[spp1][spp2][i][j]+=attack[k][spp2]*conv[k][spp2]*prob_move[spp2][move_x][move_y] * Wc[k][spp2][W_x][W_y] * C2_bin[k][spp2][m-i+xbin-1][n-j+ybin-1]*bin_width*bin_width;
							}

						}	
					}

				if(spp1 == 1 || spp2 == 1)
					for(p=m-rad-1; p<m+rad-1;p++)
					for(q=n-rad-1; q<n+rad-1;q++)
						{
			/****************************************/
			/*	Gains through predation		*/
			/****************************************/
						move_x = (m-i) + radial_max-1;
						W_x = m - p +radial_max-1;	/*Xi''*/

						move_y = (n-j) + radial_max-1;
						W_y = n - q +radial_max-1;

						if(W_x>=0 && W_y>=0 && W_x<2*radial_max-1 && W_y<2*radial_max-1)
							{

						if(conv[k][spp1]>0.00001)
C2_rhs[spp1][spp2][i][j] += attack[k][spp1]*conv[k][spp1]*prob_move[spp1][move_x][move_y]*Wc[k][spp1][W_x][W_y]*moment_closure(spp1, spp2, k, m, n, (m-p+(xbin-1)), (n-q+(ybin-1)))*bin_width*bin_width*bin_width*bin_width;

						if(conv[k][spp2]>0.00001)
C2_rhs[spp1][spp2][i][j] += attack[k][spp2] * conv[k][spp2] * 
prob_move[spp2][move_x][move_y] * Wc[k][spp2][W_x][W_y] * moment_closure(spp2, spp1, k, m, n, (m-p+(xbin-1)), (n-q+(ybin-1)))*bin_width*bin_width*bin_width*bin_width;

						}

							
						}

					}
				}

			}
		}

	}
}

/****************************************************************/
/*                                                              */
/*          CARRY OUT THE EULER STEP                            */
/*                                                              */
/****************************************************************/

void step()
{
int m, n, a, b;
int i,j;
int W_x, W_y;
int spp1, spp2;

dndt[0]=N[0];
dndt[1]=N[1];

/***Euler step for the first moments***/
	for(spp1=0;spp1<nspecies;spp1++)
		{
		N[spp1] += dt * N[spp1] * (B1[spp1] - D1[spp1]);

			for(m=0; m<2*(xbin-1)+1; m++)
			for(n=0; n<2*(ybin-1)+1; n++)
				for(spp2=0;spp2<nspecies;spp2++)
					{
					W_x = m-(xbin-1) + radial_max-1;
					W_y = n-(ybin-1) + radial_max-1;
						
				if(W_x>=0 && W_x<2*radial_max-1 && W_y>=0 && W_y<2*radial_max-1)
						{
	N[spp1] -= (double)(dt * D2[spp1][spp2] * Wc[spp1][spp2][W_x][W_y] * C2_bin[spp1][spp2][m][n] * bin_width * bin_width);
	N[spp1] -= (double)(dt * attack[spp1][spp2] * Wc[spp1][spp2][W_x][W_y] * C2_bin[spp1][spp2][m][n] * bin_width * bin_width);

	N[spp1] += (double)(dt * attack[spp2][spp1] * conv[spp2][spp1] * Wc[spp2][spp1][W_x][W_y] * C2_bin[spp2][spp1][m][n] * bin_width * bin_width);

						}
					}
		if(N[spp1]<0.0)
			N[spp1]=0.0;
		}

/*	printf("dN1/dt=%lf\tdN2/dt=%lf\n", N[0]-dndt[0], N[1]-dndt[1]);
	fflush(stdout);*/

/***Euler step for the second moments and neighbour***/
	for(m=0; m<2*(xbin-1)+1; m++)
	for(n=0; n<2*(ybin-1)+1; n++)
		for(spp1=0;spp1<nspecies;spp1++)
		for(spp2=0;spp2<nspecies;spp2++)
			{
				C2_bin[spp1][spp2][m][n] += (double)(dt *(double) C2_rhs[spp1][spp2][m][n]);
				if(C2_bin[spp1][spp2][m][n]<0.0)
					C2_bin[spp1][spp2][m][n]=0.0;
			}

if(elapsed>dt)
if(fabs(dndt[0]-N[0])<equilibrium && fabs(dndt[1]-N[1])<equilibrium)
		{
		fprintf(detail, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", B1[0], M1[1], N[0], N[1], C2_bin[0][0][xbin-1][ybin-1]/(N[0]*N[0]), C2_bin[0][1][xbin-1][ybin-1]/(N[0]*N[1]), C2_bin[1][1][xbin-1][ybin-1]/(N[1]*N[1]), elapsed);
		fflush(detail);
		fclose(detail);
		printf("Equilibrium reached, exiting\n");
		write_eq();
		exit(1);
		}

}

/****************************************************************/
/*								*/
/*		Moment closure					*/
/*								*/
/****************************************************************/

double moment_closure(i,j,k, x1,y1, x2,y2)
int i,j,k, x1,x2, y1,y2;
{
int flag=0;
double approximation;
int x11,y11, x22, y22;
double part1, part2, part3, part4, part5, part6;
FILE *moment1, *moment2;
double U,V;
/***The dimensions of the grid of the second moment***/
	int end = 2*(xbin-1)+1;

	/***For the power 2, symmetric moment closure there are 6 parts.
	All other closure methods tried are based on these building blocks.
	For more detail see Chapter 21 of Dieckmann, Law & Metz 2000. Ecological
	Geometry. Cambridge Univesity Press***/


		if(x1>=0 && x1<end && y1>=0 && y1<end)
			part3 = part1 = C2_bin[i][j][x1][y1];
		else
			part3 = part1 = N[i]*N[j];

		if(x2>=0 && x2<end && y2>=0 && y2<end)
			part5 = part2 = C2_bin[i][k][x2][y2];
		else
			part5 = part2 = N[i]*N[k];

		if(xbin-1+x2-x1>=0 && xbin-1+x2-x1<end && ybin-1+y2-y1>=0 && ybin-1+y2-y1<end)
			part6 = part4 = C2_bin[j][k][(xbin-1)+x2-x1][(ybin-1)+y2-y1];
		else
			part6 = part4 = N[j]*N[k];


if(closure_method == 2)
		{
		/***Symmetric, power 2***/
		approximation = (1.0/(alpha+Gamma)) * (
				    alpha * part1 * part2/N[i]
				    	+
				    beta * part3 * part4/N[j]
				       +
				    Gamma * part5 * part6/N[k]
				       -
				    beta * N[i] * N[j] * N[k]
				       );
		}
else if(closure_method == 3)
		{
		/***Power 3***/
		approximation = part1 * part2 * part4/ (N[i]*N[j]*N[k]);
		}
 if(closure_method == 1)
		{
		/***Power 1, symmetric -same as Bolker & Pacala 1997/1999/2000***/
		approximation = part1*N[k] + part2*N[j] + part4*N[i] - 2.0*N[i]*N[j]*N[k];
		}
/*if(i==j==k==o)
	printf("MC: %lf\n", approximation);				*/

/***Check to see whether the closure goes -ve***/
if(approximation<0.000000001)
	{
	C2flag=1;
	}
/*
	if(flag==1)
		printf("-ve\t%d\t%d\t%d\n", i,j,k);*/

		return(approximation);
}

/****************************************************************/
/*                                                              */
/*	LOOK UP TABLE FOR MOVEMENT PROBABILITIES		*/
/*                                                              */
/****************************************************************/

/*This function computes the probability of moving from one bin 
to another based upon the movement kernal. The results are then 
normalized for a discrete distance of the parent distribution
beyond which the probability of movement to this bin is zero*/

void look_up()
{
int i, j;
double sum[2] = {0.0, 0.0};
double normalize[2]={0.0, 0.0};
double distance, dx, radial_d;
double z, PDF, m_normal, s_normal, trans;
double increment, x, y;
double move_sum[2]={0.0,0.0};
double move_sum1[2]={0.0, 0.0};
double move_sum2[2]={0.0, 0.0};

FILE *ptr, *out1;
int spp1;

	ptr = fopen("out.dat", "w");
	out1 = fopen("error.dat", "w");

	/***The method used here is to compute the radial p.d.f. then to 
	convert to the Cartesian form. The problem is that in doing so we have 
	to divide the p.d.f by (distance*2*pi). Hence when distance=0 we have a 
	problem (dividing by zero). To get round this we initially ignore the 
	p.d.f. at distance zero, compute the p.d.f. for all other distances.
	The p.d.f. at zero is then = 1 - sum(p.d.f. all other distances)***/

for(i=0; i<2*(radial_max-1)+1; i++)
for(j=0; j<2*(radial_max-1)+1; j++)
for(spp1=0;spp1<nspecies;spp1++)
	{
		prob_move[spp1][i][j]=0.0;

	/*Compute the relative distances of the cells to the central cell*/
		x = (double)i - (radial_max-1);
		y = (double)j - (radial_max-1);

		x = fabs(x*bin_width);
		y = fabs(y*bin_width);

		distance = x*x + y*y;
		distance = sqrt(distance);
	
	if(distance < rmax[spp1])
		if(movement_kernel == 4)
			{
			trans = 2.0 * move_stdv[spp1] * move_stdv[spp1] * Pi;

			m_normal = (distance)/move_stdv[spp1];
			m_normal = m_normal * m_normal;

			PDF = exp(-0.5 * m_normal)/trans;

			move_sum[spp1] += PDF;

			prob_move[spp1][i][j] = PDF;
		}
	}

move_sum1[0]=move_sum1[1]=1.0/(bin_width*bin_width);

/***Renormalise***/

	for(i=0;i<2*(radial_max-1)+1;i++)
		{
		for(j=0;j<2*(radial_max-1)+1;j++)
			{
			fprintf(ptr, "%d\t%d\t%d\t%d\t", i, j, i-radial_max+1, j-radial_max+1);

			for(spp1=0;spp1<nspecies;spp1++)
				{
				prob_move[spp1][i][j] *= move_sum1[spp1]/move_sum[spp1];

			/***Print movement/dispersal kernel to file***/
				fprintf(ptr, "%.10lf\t", prob_move[spp1] [i] [j]);
				move_sum2[spp1] += prob_move[spp1] [i] [j];
				}
			fprintf(ptr, "\n");
			}
		}

	for(spp1=0;spp1<nspecies;spp1++)
		{
		printf("Renormalised zeroth bin: %lf\n", prob_move[spp1][radial_max-1][radial_max-1]);
		printf("Total %.10lf\n", move_sum2[spp1]);
		fprintf(ptr, "\n#Sum = %.10lf\n", move_sum2[spp1]);
		}

	fflush(ptr);
	fclose(ptr);
}


/****************************************************************/
/*								*/
/*		Weighting functions lookup tables		*/
/*								*/
/****************************************************************/

void weighting()
{
double W_trans[nspecies][nspecies]={0.0,0.0,0.0,0.0};
double W_norm[nspecies][nspecies]={0.0,0.0,0.0,0.0};
double W_PDF[nspecies][nspecies]={0.0,0.0,0.0,0.0};
int i, j;
double m,n;
int spp1, spp2;
double increment;
double threshold;
double distance, dx, dy;
double W_sum[nspecies][nspecies]={0.0,0.0,0.0,0.0};
double W_sum1[nspecies][nspecies]={0.0,0.0,0.0,0.0};
FILE *output1;

output1 = fopen("weighting.dat", "w");

	for(spp1=0;spp1<nspecies;spp1++)
	for(spp2=0;spp2<nspecies;spp2++)
		{
		W_norm[spp1][spp2] = 2.0*Pi*W_stdv[spp1][spp2]*W_stdv[spp1][spp2];
		W_sum[spp1][spp2]=W_sum1[spp1][spp2]=0.0;
		}

for(spp1=0;spp1<nspecies;spp1++)
for(spp2=0;spp2<nspecies;spp2++)
for(i=0; i<2*(radial_max-1)+1; i++)
	{
	for(j=0; j<2*(radial_max-1)+1; j++)
		{
		Wc[spp1][spp2][i][j]=0.0;
			
	/***Compute the distance from the origin***/
		dx = (double)i - (radial_max-1);
		dx = fabs(dx * bin_width);

		dy = (double)j - (radial_max-1);
		dy = fabs(dy * bin_width);

		distance = (dx*dx) + (dy*dy);
		distance = sqrt(distance);

	/***Normal kernel***/
		threshold = (W_stdv[spp1][spp2]*3.0);

		if(distance<threshold)
				{
				W_trans[spp1][spp2] = (distance)/W_stdv[spp1][spp2];
				W_trans[spp1][spp2] = W_trans[spp1][spp2] * W_trans[spp1][spp2];

				W_PDF[spp1][spp2] = exp(-0.5 * W_trans[spp1][spp2])/W_norm[spp1][spp2];

				Wc[spp1][spp2][i][j] = W_PDF[spp1][spp2];
				W_sum[spp1][spp2]+=W_PDF[spp1][spp2];
				}
		else
			Wc[spp1][spp2][i][j]=0.0;
		}
	}


	for(spp1=0;spp1<nspecies;spp1++)
	for(spp2=0;spp2<nspecies;spp2++)
		for(i=0;i<2*(radial_max-1)+1;i++)
		for(j=0;j<2*(radial_max-1)+1;j++)
			{
			W_sum1[spp1][spp2] += Wc[spp1][spp2][i][j];
			W_sum[spp1][spp2] = 1.0/(bin_width*bin_width);
			}

/***Now renormalise so that the function sums to one***/
	for(spp1=0;spp1<nspecies;spp1++)
	for(spp2=0;spp2<nspecies;spp2++)
		for(i=0;i<2*(radial_max-1)+1;i++)
		for(j=0;j<2*(radial_max-1)+1;j++)
			{
			Wc[spp1][spp2][i][j] *= (W_sum[spp1][spp2]/W_sum1[spp1][spp2]);
			}

	/***Print detail of kernels to file***/
	for(i=0;i<2*(radial_max-1)+1;i++)
		{
		for(j=0;j<2*(radial_max-1)+1;j++)
			{
			fprintf(output1, "%d\t%d\t", i-radial_max+1, j-radial_max+1);

			for(spp1=0;spp1<nspecies;spp1++)
			for(spp2=0;spp2<nspecies;spp2++)
				{
				fprintf(output1, "%.10lf\t", Wc[spp1][spp2][i][j]);
				}
			fprintf(output1, "\n");

			}
		fprintf(output1, "\n");
		}

	fclose(output1);

}

/****************************************************************/
/*                                                              */
/*								*/
/*                                                              */
/****************************************************************/


void write_moments()
{
int m,n,spp1,spp2;

	fprintf(first, "%lf\t", elapsed);

	for(spp1=0;spp1<nspecies;spp1++)
		{
		printf("\t\t%lf\t", N[spp1]);
		fprintf(first, "%lf\t", N[spp1]);
		}

	fprintf(first, "Tijk: %d\t", C2flag);
	printf("Tijk: %d\n", C2flag);

	fprintf(first, "\n");
	printf("\n");
	fflush(first);
	fflush(stdout);

/***Write second moment to file and screen***/
	for(n=xbin-1, m=xbin-1; m<2*xbin-1; m++)
		{
		fprintf(out10, "%.2lf\t%.2lf\t", elapsed, bin_width*(double)(m-(xbin-1)));
		printf("%.2lf\t%.2lf\t", elapsed, (double)(m-(xbin-1))*bin_width);	

		for(spp1=0; spp1<nspecies; spp1++)
			{
			for(spp2=0; spp2<nspecies; spp2++)
					{
					printf("%.5lf\t", C2_bin[spp1][spp2][m][n]/(N[spp1]*N[spp2]));
					fprintf(out10, "%.5lf\t", C2_bin[spp1][spp2][m][n]/(N[spp1]*N[spp2]));

					}
			}
		printf("\n");
		fprintf(out10, "\n");

		}
		fprintf(out10, "\n");
	fflush(out10);
}

/****************************************************************/
/*                                                              */
/*		Write equilibrium moments			*/
/*                                                              */
/****************************************************************/

void write_eq()
{
int spp1, spp2, i, j;

	for(spp1=0;spp1<nspecies;spp1++)
			fprintf(eq_m, "%lf\t", N[spp1]);

	fprintf(eq_m, "\n");

	for(i=0;i<2*xbin-1;i++)
	for(j=0;j<2*ybin-1;j++)
		{
		for(spp1=0;spp1<nspecies;spp1++)
		for(spp2=0;spp2<nspecies;spp2++)
			{
			fprintf(eq_m, "%lf\t", C2_bin[spp1][spp2][i][j]);
			}
		fprintf(eq_m, "\n");
		}
fclose(eq_m);

}

/****************************************************************/
/*                                                              */
/*                                                              */
/****************************************************************/

#include "C2_checks.h"

int main(argc, argv)
int argc;
char *argv[];
{
FILE *in, *out;
int i, j;

char fn[100];
double clock, tick;

out = fopen("test.dat", "w");
data = fopen("run.dat", "w");
out10 = fopen(out_file, "w");
first = fopen(No_file, "w");
temp_out = fopen("output.temp.dat", "w");
detail = fopen(density_file, "a");

	init_C2();

	B1[0]=atof(argv[1]);
	M1[0]=M1[1]=atof(argv[2]);

	sprintf(fn, "equilibrium.%.4lf.%.2lfdat", B1[0], M1[0]);
	eq_m = fopen(fn, "w");

	look_up();
	weighting();

	elapsed = 0.0;

	printf("Read in Cee\n");
	fflush(stdout);

	tick = tock = dt;

	do{

		write_moments();
		C2flag=0;
		elapsed+=dt;
/*		printf("t = %lf\n", elapsed);
		printf("\n");
		fflush(stdout);*/

		calc_cross_rhs();
		step();
/*printf("Gets here\n");*/
fflush(stdout);

/*		if(check_moments == 1)
			check();*/

	}while(elapsed<=tmax+dt);
fprintf(out, "FINISHED> WHOOPEEEEEE\n");
fflush(out);
return 0;
}