#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define PI 3.1415926535897932384626433

// QUICKSORT INCLUDED!!
// REAL DATA INPUT!!
// INFINITE CELLS!!
// INNER DYNAMICS!!

/*************************************************************************************/
/***************************** Declare New Structures ********************************/
/*************************************************************************************/
      

typedef struct{                     // Create a new structure called 'coord' to contain cell coordinates                                                                    
	float x;                       // Create 'x' member of coord 
	float y;                       // Create 'y' member of coord 
	float z;                       // Create 'z' member of coord 
	float hsl;
	float iptg;              
    int lysis;                     // Lysis state  (on or off)
	int condition;                 // 0 - not , 1= qs, 2 = lysis, 3 = dead
	float randx;                   // Create 'x' member of coord 
	float randy;                   // Create 'y' member of coord 
	float randz;                   // Create 'z' member of coord 
	float radii;
	int tumble;          
	} data;
	
typedef struct{
        float x;
        float y;
        float z;
        }coord;


/*************************************************************************************/
/***************************** Declare Subfunctions **********************************/
/*************************************************************************************/

void initial_cell_number(int *cell_number, float density, float observed_radius, FILE *outfile10);

void initial_cell_prop(int cell_number, float observed_radius, data *cell, FILE *outfile10);

void movement(int cell_number, data *cell, coord source);

void radii_function(int cell_number, data *cell, coord source, FILE *outfile10);

void sort(int i,int cell_number,float *radii,int *order,data *cell, FILE *outfile10);

void quorum_sensing(int i, int cell_number, float *effective_radius, float qs_density, data *cell, float *radii, int *order, FILE *outfile10);	

void inner_dynamics(int cell_number, data *cell, int *number_death, int *number_alive, int *number_qs, FILE *outfile6);

void cell_increase(int cell_number, int *increase_cell_number, float observed_radius, float density, float *radii, FILE *outfile10);

void increase_cell_prop(int cell_number, int increase_cell_number, float observed_radius, data *cell);

void cell_decrease(int cell_number, int *decrease_cell_number, data *cell);

void density_function(int i, int cell_number, float observed_radius, float *radii, FILE *outfile2);

////////////////////////////////

float vsphere(float r);

void quicksort(float *radii,int *order,int first,int last);

int partition(float *radii,int *order,int first,int last);

const int PoissonRandomNumber(const double lambda);


int main(){


	/*************************************************************************************/
	/****************************** Seed random number ***********************************/
	/*************************************************************************************/

	srand(time(NULL));

	/*************************************************************************************/
	/********************************* Variables Known ***********************************/
	/*************************************************************************************/

	int i;                             // Used as a counter for time                                                                                
	int j;                             // Used as a counter for cells 
	int k;                             // Spare counter

    int qs_number;
    int counter = 0;

	//////////////////////////////
	int time_step = 20000;              // Number of time iterations (second)
	
	float observed_radius = 0.005;      // Radius within which 'observe' - determine properties (meter)

	float density = 1.0*pow(10,12);	// TEMP: Initial density of cells (meter^3) 
	
	float qs_density = 7*pow(10,14);   // Density which causes quorum sensing (meter^3) $$Cellular COntrol Of Synthesis$$
	    

	/*************************************************************************************/
	/********************************* Variables Unknown *********************************/
	/*************************************************************************************/

	int cell_number = 0;                   	// Number of cells at each time step

	int increase_cell_number = 0;        	// Increase in cells each time step

	int decrease_cell_number = 0;			// Decrease in cells each time step

	float effective_radius = 0;            	// Minimum quroum sensing radius

	int number_death = 0;

	int number_alive = 0;

	int number_qs = 0;

	/*************************************************************************************/
	/************************************ Scale ******************************************/
	/*************************************************************************************/       

	float scale = 5000;

	density /= scale;
	
	qs_density /= scale;


	/*************************************************************************************/
	/****************************** Initialise Variables *********************************/
	/*************************************************************************************/                                                      

	coord source;
	source.x = 0.0;                
	source.y = 0.0;
	source.z = 0.0;

	/*************************************************************************************/
	/************************************* Outfiles **************************************/
	/*************************************************************************************/


	// Outfile 1 lists results as a single column
	FILE *outfile1;
	outfile1=fopen("output/chemotaxis_single.dat","w"); 

	// Outfile 2 lists density aginst time step and position
	FILE *outfile2;
	outfile2=fopen("output/density.dat","w");

	// Outfile 3 lists cell number aginst time step
	FILE *outfile3;
	outfile3=fopen("output/cell_number.dat","w");

	// Outfile 4 lists effective radius aginst time step
	FILE *outfile4;
	outfile4=fopen("output/effective_radius.dat","w");
	
	// Outfile 5 lists hsl and iptg aginst time step
	FILE *outfile5;
	outfile5=fopen("output/environment.dat","w");

	// Outfile 6 lists condition 
	FILE *outfile6;
	outfile6=fopen("output/condition.dat","w");


    // Outfile 10 Check
	FILE *outfile10;
	outfile10=fopen("output/check.dat","w");


	//////////////////////////////////////////////////////////
	printf("main\n");

	initial_cell_number(&cell_number, density, observed_radius, outfile10);						                     // FUNCTION: Calculate the initial cell number

	data *cell;                      									               
	cell = (data *) malloc (cell_number * sizeof(data));     			                                     // Initialise memory 

	float *radii;
	radii = (float *) malloc (cell_number * sizeof(float));   							                     // Initialise memory 

    int *order;
	order = (int *) malloc (cell_number * sizeof(int));   
    	  
	initial_cell_prop(cell_number, observed_radius, cell, outfile10);		                                             // FUNCTION: Assign properties to initial cells (position, steady state, etc)
	
	printf("Cell Number = %d\n",cell_number);

    	for(i=0; i<time_step; ++i){  										              

        movement(cell_number, cell, source);                                                                 // FUNCTION: Movement through chemotaxis

        radii_function(cell_number, cell, source, outfile10);									                         // FUNCTION: Calculates the radial distance between each cell and the source

        sort(i, cell_number, radii, order, cell, outfile10);

        //density_function(i, cell_number, observed_radius, radii, outfile2);				                 // FUNCTION: Calculates the density of cells
        
		quorum_sensing(i, cell_number, &effective_radius, qs_density, cell, radii, order, outfile10);	             // FUNCTION: Calculates the effective density    

		inner_dynamics(cell_number, cell, &number_death, &number_alive, &number_qs, outfile6);                                                                    // FUNCTION: Run inner dynamics

		cell_decrease(cell_number, &decrease_cell_number, cell);					                         // FUNCTION: Calculates decrease cell number from lysis and removes cells by reordering

          cell_number -= decrease_cell_number;   

		cell_increase(cell_number, &increase_cell_number, observed_radius, density, radii, outfile10); 	 // FUNCTION: Calculates the influx of cells

		/////////////////////////////////////////////////////////////////////////
		cell = (data*)realloc(cell,(cell_number+increase_cell_number)*sizeof(data));        	             // Vary Memory Allocation as required
		radii = (float*)realloc(radii,(cell_number+increase_cell_number)*sizeof(float));                     // Initialise memory 
          order = (int*)realloc(order,(cell_number+increase_cell_number)*sizeof(int));                     // Initialise memory 
		/////////////////////////////////////////////////////////////////////////

		increase_cell_prop(cell_number, increase_cell_number, observed_radius, cell);	                     // FUNCTION: Assign properties to increased cells (position, steady state, etc)

		cell_number += increase_cell_number;                                                               // Update cell number


        ///////////////////  OUTPUT  //////////////////////
		//printf("output\n");

		// Cell Positions
		for(j=0; j<cell_number; ++j){                  
			 fprintf(outfile1,"%d\t%d\t%lf\t%lf\t%lf\n",i,j,cell[j].x,cell[j].y,cell[j].z);
             }

		// Cell Number
		fprintf(outfile3,"%d\t%d\t%d\t%d\n",i,cell_number,increase_cell_number,decrease_cell_number);

		// Effective radius
		fprintf(outfile4,"%d\t%f\n",i,effective_radius);
		
		// Cell Condition
	     fprintf(outfile6,"%f\t%f\t%f\t%d\n",((float)number_alive/cell_number),((float)number_qs/cell_number),((float)number_death/cell_number),cell_number);
       
	
  
		if(counter == 100){
			printf("%d\n",i);
			printf("Cell Number = %d\n",cell_number);
			counter = 0;
			}
		counter +=1;
        
        
        }

	fclose(outfile1);
	fclose(outfile2);
	fclose(outfile3);
	fclose(outfile4);
	fclose(outfile5);
	fclose(outfile6);

	fclose(outfile10);
	
	free(cell);

	printf("Success\n");
	

	}


/******************************************************************************************/
/******************************************************************************************/
/******************************************************************************************/
/********************************** SUBFUNCTIONS ******************************************/
/******************************************************************************************/
/******************************************************************************************/
/******************************************************************************************/

/******************************************************************************************/
/******************************** initial_cell_number *************************************/
/******************************************************************************************/

void initial_cell_number(int *cell_number, float density, float observed_radius, FILE *outfile10){
//printf("intitial_cell_number\n");

	float volume;
	volume = vsphere(observed_radius);
	*cell_number =  (int)density * volume;
	if(*cell_number == 0){printf("Error: Initial cell number is zero.\n"); exit(1);}
	
fprintf(outfile10,"initial cell number = %d\n",*cell_number);
	}

/******************************************************************************************/
/******************************** initial_cell_prop ***************************************/
/******************************************************************************************/

void initial_cell_prop(int cell_number, float observed_radius, data *cell, FILE *outfile10){
//printf("initial_cell_prop\n");

	int j;

    j=0;
	while(j<cell_number){                               				                         

		cell[j].x = ((float)(rand() - RAND_MAX/2) / (RAND_MAX/2))*observed_radius;      // Cartesian X-coordinate
		cell[j].y = ((float)(rand() - RAND_MAX/2) / (RAND_MAX/2))*observed_radius;    	// Cartesian Y-coordinate
		cell[j].z = ((float)(rand() - RAND_MAX/2) / (RAND_MAX/2))*observed_radius;      // Cartesian Z-coordinate
		
		cell[j].radii =  sqrt( pow((cell[j].x),2) + pow((cell[j].y ),2) + pow((cell[j].z),2));
		
		cell[j].tumble = 0;													// Set tumbles to zero	
		cell[j].condition = 0;												// Set condition to zero (alive)
		cell[j].lysis = 3000;	                                               		// Analysis
		
		if(cell[j].radii <= observed_radius){
//fprintf(outfile10,"initial_positions\t%f\t%f\t%f\n",cell[j].x,cell[j].y,cell[j].z);
            ++j;
            }
         }
	}


/******************************************************************************************/
/************************************** movement ******************************************/
/******************************************************************************************/

void movement(int cell_number, data *cell, coord source){
//printf("movement\n");

	int j; 
	float distance;

	float N = 10.0;                                                                          // Surface concentration of aspartate                                                                                                                  
	float D = 0.9 * pow(10,-5);                                                              // Diffusion coefficient of aspartate from $Response Tuning In Bacterial Chemotaxis$
	float cell_speed = 20 * pow(10,-6);                                                      // Cell speed from $Computerized Analysis Of Chemotaxis....$ 


	for(j=0; j<cell_number; ++j){    									                     // Initialise the CELL loop
		if(cell[j].tumble == 0){     									                     // If tumbles = 0 then enter random tumble generator

			cell[j].randx = ((float)(rand() - RAND_MAX/2) / (RAND_MAX/2));  	             // Calculate a random number in [-1,1]
			cell[j].randy = ((float)(rand() - RAND_MAX/2) / (RAND_MAX/2));  	             // Calculate a random number in [-1,1]
			cell[j].randz = ((float)(rand() - RAND_MAX/2) / (RAND_MAX/2));  	             // Calculate a random number in [-1,1]

			//fprintf(outfile10,"Random number between [-1,1] = %lf\n",random[j].x); // Check
			//fprintf(outfile10,"Random number between [-1,1] = %lf\n",random[j].y); // Check
			//fprintf(outfile10,"Random number between [-1,1] = %lf\n",random[j].z); // Check

			double norm = sqrt( pow(cell[j].randx,2) + pow(cell[j].randy,2) + pow(cell[j].randz,2));    // Calculate the length of the randomly generated vector 

			//fprintf(outfile10,"Norm = %lf\n",norm); // Check

			cell[j].randx = ( cell[j].randx / norm ) * cell_speed; 		                       // Recalculate random components in terms of cell speed
			cell[j].randy = ( cell[j].randy / norm ) * cell_speed; 		                       // Recalculate random components in terms of cell speed
			cell[j].randz = ( cell[j].randz / norm ) * cell_speed; 		                       // Recalculate random components in terms of cell speed

			//fprintf(outfile10,"Cell %d random x-component = %e\n",j,random[j].x); // Check    
			//fprintf(outfile10,"Cell %d random y-component = %e\n",j,random[j].y); // Check    
			//fprintf(outfile10,"Cell %d random z-component = %e\n",j,random[j].z); // Check                   

			cell[j].x += cell[j].randx; 								                       // New position in cartesian coordinates
			cell[j].y += cell[j].randy; 								                       // New position in cartesian coordinates
			cell[j].z += cell[j].randz; 								                       // New position in cartesian coordinates

			// Calculate the distance from source AFTER random tumble
			distance = sqrt( pow((cell[j].x - source.x),2) + pow((cell[j].y - source.y),2) + pow((cell[j].z - source.z),2));

			//fprintf(outfile10,"Distance before %f \n",movement[j].dist_before);  // Check    
			//fprintf(outfile10,"Distance after  %f \n",movement[j].dist_after);  // Check     

			////////////////////////// Modify tumbles (Concentration) /////////////////////////////

			if(distance < cell[j].radii){
				cell[j].tumble = 20; 								// Value calculated from 0.05 per sec $Computerized Analysis Of Chemotaxis....$
				}
			else{
				cell[j].tumble = 3;  								// Value calculated from 0.3 per sec $Computerized Analysis Of Chemotaxis....$
				}
			}  

			else{ 
				cell[j].x += cell[j].randx;
				cell[j].y += cell[j].randy;
				cell[j].z += cell[j].randz;
				} 

			cell[j].tumble -= 1; 

		}  

	} 
	
	
	
/******************************************************************************************/
/*********************************** radii_function ***************************************/
/******************************************************************************************/

void radii_function(int cell_number, data *cell, coord source, FILE *outfile10){
//printf("radii_function\n");

	int j;

	for(j=0; j<cell_number; ++j){
		cell[j].radii =  sqrt( pow((cell[j].x - source.x),2) + pow((cell[j].y - source.y),2) + pow((cell[j].z - source.z),2));
//fprintf(outfile10,"radii = %f\n",cell[j].radii);
        }
	}
	
	
/******************************************************************************************/
/*********************************** sort_function ***************************************/
/******************************************************************************************/
	
void sort(int i,int cell_number,float *radii,int *order,data *cell, FILE *outfile10){
//printf("sort\n");

     int j;

    for(j=0; j<cell_number; ++j){
        radii[j] = cell[j].radii;
        order[j] = j;
        }

    quicksort(radii,order,0,cell_number-1);								                              // Sort the array radii
    
    if(i==0){                                                                 	// CHECK
        for(j=0; j<cell_number; ++j){                                         	
            fprintf(outfile10,"cell[%d].radii = %f\t",j,cell[j].radii);       	
            fprintf(outfile10,"radii[%d] = %f\t",j,radii[j]);                 	
            fprintf(outfile10,"order[%d] = %d\n",j,order[j]);                 	
            }                                                                   
        }                                                                       
    }
	
	
/******************************************************************************************/
/*********************************** quorum_sensing ***************************************/
/******************************************************************************************/

void quorum_sensing(int i, int cell_number, float *effective_radius, float qs_density, data *cell, float *radii, int *order, FILE *outfile10){
//printf("quorum_sensing\n");

	int j;
	int k;
	int l;
	float volume;
	float density;
	
    *effective_radius = 0.0;                                		// Set effective radii to zero

	for(j=cell_number-1; j>=20; --j){  								// Countdown from cell_number to 2 (largest radii to smallest radii)
		
        volume = vsphere(radii[j]);  								// Volume of a sphere     
		density = j / volume;  										// Calculate the density                                                        
		
        if(density > qs_density){  									// If the density is greater than the density required for quorum sensing
            *effective_radius = radii[j];                 			// Record the radius at which density fails to be high enough
        
            for(k=j; k>0; --k){    
                l = order[k];
			 if(cell[l].condition == 0){
		           cell[l].condition = 1;                              // Quorum sensed
				 }
			 }
            
            break;
            }
        }
	}

/******************************************************************************************/
/************************************ inner_dynamics **************************************/
/******************************************************************************************/
       
void inner_dynamics(int cell_number, data *cell, int *number_death, int *number_alive, int *number_qs, FILE *outfile6){
//printf("inner_dynamics\n");

    int j;
    
    *number_alive = 0;
    *number_qs = 0;
   
	for(j=0; j<cell_number; ++j){
        
        if(cell[j].condition == 0){
            ++*number_alive;
            }    
            
        else if(cell[j].condition == 1){
		  --cell[j].lysis; 
            if(cell[j].lysis == 0){
                cell[j].condition = 2;
			 ++*number_death;
                } 
		  else{
			 ++*number_qs;
			 } 
            }    
        }
	}


/******************************************************************************************/
/*********************************** cell_decrease ****************************************/
/******************************************************************************************/

void cell_decrease(int cell_number, int *decrease_cell_number, data *cell){
//printf("cell_decrease\n");

	int j;
     int k;
     int l;
    
    *decrease_cell_number = 0;
	for(j=0; j<cell_number; ++j){
		if(cell[j].condition == 2){        		          // condition 3 means dead
			*decrease_cell_number += 1;						// Increase the number of dead cells
			for(k=j; k<cell_number-1; ++k){
				cell[k] = cell[k+1];					  // Remove cell from memory
				} 
			}
		}
	}


/******************************************************************************************/
/********************************** cell_increase *****************************************/
/******************************************************************************************/

void cell_increase(int cell_number, int *increase_cell_number, float observed_radius, float density, float *radii, FILE *outfile10){
//printf("cell_increase\n");

	int j;	

	int critical_number;													// Cell number at 'critical_radius'
	int critical_diff;														// Number of cells in the critical zone

	float critical_radius;													// Radius outside which the density must be the predefined radius
	float critical_volume;													// Volume of the critical zone
	float critical_density;													// Density of the critical zone
	float density_diff;														// Density difference between the predefined density and density in the critical zone

	critical_radius = 0.8 * observed_radius;								// Calculate 'critical_radius' as a fraction of predefined 'observed_radius'

	for(j=cell_number-1; j>=0; --j){	                                    // NOTE: cell_number-1 used!!!										

		if(radii[j] < critical_radius){										// When the radius falls below the 'critical_radius' set 'critical_number' to j
			critical_number = j;										
			break;
			}
		}

	critical_diff = cell_number - critical_number;								// Calculate number of cells within critical zone

	critical_volume = vsphere(radii[cell_number-1]) - vsphere(radii[critical_number]);	// // NOTE: cell_number-1 used!! Calculate volume of critical zone

	critical_density = critical_diff / critical_volume;							// Calculate density within critical zone

fprintf(outfile10,"cell_number = %d\n",cell_number);
fprintf(outfile10,"critical_number = %d\n",critical_number);
fprintf(outfile10,"critical_diff = %d\n",critical_diff);
fprintf(outfile10,"critical_volume = %f\n",critical_volume);
fprintf(outfile10,"density = %f\n",density);
fprintf(outfile10,"critical_density = %f\n",critical_density);

	density_diff = density - critical_density;

	if(density_diff > 0){ 
		*increase_cell_number = (int)(density_diff * critical_volume);
		}
	else{
		*increase_cell_number = 0;
		}
			
fprintf(outfile10,"increase_cell_number = %d\n",*increase_cell_number);

	}


/******************************************************************************************/
/******************************** cell_increase_prop **************************************/
/******************************************************************************************/

void increase_cell_prop(int cell_number, int increase_cell_number, float observed_radius, data *cell){
//printf("increase_cell_prop\n");

	int j;

	for(j=cell_number; j<(cell_number+increase_cell_number); ++j){				

		float random_radius = (0.9 * observed_radius) + (0.1 * observed_radius * ((float)(rand() - RAND_MAX/2) / (RAND_MAX/2)));   // Radial distance - enter them in critical zone
		float random_theta  = ((double)rand() / RAND_MAX) * 180;                  		// Theta angle - random
		float random_alpha  = ((double)rand() / RAND_MAX) * 360;                  		// Alpha angle - random

		cell[j].x = random_radius * cos(random_alpha) * sin(random_theta);       		// Convert to cartesian x-coordinate
		cell[j].y = random_radius * sin(random_alpha) * sin(random_theta);       		// Convert to cartesian y-coordinate
		cell[j].z = random_radius * cos(random_theta);                           		// Convert to cartesian z-coordinate

		cell[j].tumble = 0;													// Set tumbles to zero
		cell[j].condition = 0;												// Set condition to zero (alive)
        	cell[j].lysis = 3000;												// Analysis
        
		}     
	}


/******************************************************************************************/
/********************************* density_function ***************************************/
/******************************************************************************************/

void density_function(int i, int cell_number, float observed_radius, float *radii, FILE *outfile2){
//printf("density\n");
          
 
    int j;
	int k;

	int   partition;												// The number of partitions to 'break' observed_radius 					
	float partition_size;										// The size of each partition
	float threshold_radii;						

	partition = 50;											    // TEMP

	struct density_prop{										// Define new structure to record density properties
		int   number[partition];		
		float radii[partition];
		float volume[partition];
		float density[partition];
		};				

	struct density_prop store;									// Declare a structure 'density_prop' called 'store'

	for(k=0; k<partition; ++k){									// Set the variables to zero
		store.number[k] = 0;		
		store.radii[k]  = 0;
		store.volume[k] = 0;
		store.density[k] = 0;
		}

	partition_size = observed_radius / partition; 				// Calculate the 'partition_size'

	
	
    threshold_radii = partition_size;							// Initialise the threshold radius for 'do-while' loop
    j=0;														// Initialise the variables to zero for 'do-while' loop
	k=1;														// Initialise the variables to zero for 'do-while' loop
	do{	
        if(j==cell_number){
            break;
            }
        							   									
		if(radii[j] < (k*partition_size) ){					       // If: radii is less than current threshold then count and move to next cell
            store.number[k] += 1;
			++j;                                      
			}

		else{											// Else: store current values and increase threshold
            store.radii[k]   = k*partition_size;
			store.volume[k]  = vsphere(k*partition_size) - vsphere((k-1)*partition_size);
			
            if(store.number[k] == 0){
                store.density[k] = 0;     
                }
            else{
                store.density[k] = store.number[k] / store.volume[k] ;
                }
                
			++k;
			}
			
		}while(k < partition);

	for(k=0; k<partition; ++k){
        threshold_radii = k*partition_size;
		fprintf(outfile2, "%d\t%f\t%f\t%d\t%f\n",i,store.radii[k],log(store.density[k]),store.number[k],store.volume[k]);
     	}

     }


/******************************************************************************************/
/******************************************************************************************/
/******************************************************************************************/
/********************************* SUBSUBFUNCTIONS ****************************************/
/******************************************************************************************/
/******************************************************************************************/
/******************************************************************************************/

/******************************************************************************************/
/************************************** vsphere *******************************************/
/******************************************************************************************/

float vsphere(float r){

	float volume;
	volume = (4/3) * PI * pow(r,3);
	return(volume);
}

/******************************************************************************************/
/********************************** quicksort *********************************************/
/******************************************************************************************/

void quicksort(float *radii, int *order, int first,int last){

	int pivot;

	if (first<last){
		pivot=partition(radii,order,first,last);
		quicksort(radii,order,first,pivot-1);
		quicksort(radii,order,pivot+1,last);	
		}
	}


/******************************************************************************************/
/********************************** partition *********************************************/
/******************************************************************************************/

int partition(float *radii,int *order,int first,int last){

	int pivot;
	float pivotvalue;
	float temp;
	float temp_order;
	int i;

	pivot=first;
	pivotvalue=radii[first];
	
	for(i=first; i<=last; ++i){
		if(radii[i]<pivotvalue){
			++pivot;
			if(i!=pivot){
				
                temp=radii[pivot];
				radii[pivot]=radii[i];
				radii[i]=temp;
				
                temp_order=order[pivot];
				order[pivot]=order[i];
				order[i]=temp_order;
				}
			}

		}

	temp=radii[pivot];
	radii[pivot]=radii[first];
	radii[first]=temp;
	
	temp_order=order[pivot];
	order[pivot]=order[first];
	order[first]=temp_order;

	return(pivot);
	}


/******************************************************************************************/
/********************************* PoissonRandomNumber ************************************/
/******************************************************************************************/

const int PoissonRandomNumber(const double lambda){
	int k=0;                             //Counter
	const int max_k = 2*lambda;          //k upper limit
	double p = 1.*rand()/((RAND_MAX));   //uniform random number
	double P = exp(-lambda);             //probability
	double sum=P;                        //cumulant
	if (sum>=p) return 0;                //done allready
		for (k=1; k<max_k; ++k) {            //Loop over all k:s
			P*=lambda/(double)k;               //Calc next prob
			sum+=P;                            //Increase cumulant
			if (sum>=p) break;                 //Leave loop
			}

	return k;                            //return random number
	}