#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "delta.h"
#include "kirkrand.h"

/***************************************************************************
AgeSeed()
	This function gives a linear drop in seed numbers over the seedlife
Arguments: 		struct Lifehist *lifehist	pointer to lifehistory parameters
					int *seeds						pointer to the number of seeds in
															a cell
					int i								species i.d. number

Return Value:	void
****************************************************************************/
void	AgeSeed(struct Basin *basin, struct Lifehist *lifehist, struct Plant ***plant)
{
	double d, m;
	int x, y, i;

	for(y=0; y < basin->ymax[0]; y++) {
		for(x=0; x < basin->xmax[0]; x++) {
			if(plant[x][y][0].adult[0] == -9999)
				continue;
			for(i=1; i < MAXSP; i++) {
				d = (double)plant[x][y][0].seed[i];
				m = lifehist->seedsurv[i];
				plant[x][y][0].seed[i] = (int)(d*m);
			}
		}
	}
}
/***************************************************************************
Germinate()
	This function selects a percent of each age seed for germination if water
	is shallower than 0 m and the cell is empty of adults at the begining
	of the time step. Seeds can still germiante if the cell was invaded by
	rhizomes within the current step. The germination percent is species
	specific and is determined in lifehist parameters structure.
Arguments: 		struct Plant ***plant		pointer to plant structure
					struct Lifehist *lifehist	pointer to life history parameters
					struct Run *run				pointer to run parameters
					struct Basin *basin			pointer to basin parameters

Return Value:	void
****************************************************************************/
void	Germinate(struct Basin *basin, struct Lifehist *lifehist, struct Plant ***plant, struct Run *run, double **topo, double *water)
{
	int x, y, i;
	double depth;

	for(y=0; y < basin->ymax[0]; y++) {
		for(x=0; x < basin->xmax[0]; x++) {
			if(plant[x][y][0].adult[0] == -9999)
				continue;
			GetDepth(basin, run, topo, water, &depth, x, y);
			for(i=1; i < MAXSP; i++) {
				if(depth < 0 && plant[x][y][1].adult[0] == 0) {
					plant[x][y][0].sdlng[i] += (int)(((double)plant[x][y][1].seed[i])*lifehist->germrate[i]);
					plant[x][y][0].seed[i] = plant[x][y][0].seed[i] - (int)(((double)plant[x][y][1].seed[i])*lifehist->germrate[i]);
					if(plant[x][y][0].seed[i] < 0) {
						plant[x][y][0].seed[i] = 0;
					}
				}
			}
		}
	}
}
/***************************************************************************
GetProb()
	This function calculates the probablity of each species of plant 
	occurring at a given water depth based on logistic regressions using 
	frequency of occurrence of species in the MERP wetland complex during 
	1980. The regressions for species 2,3,5,6 were calculated by de Swart et
	al. (1992). The model for unvegetated areas (case 0) was calculated by EWS.
	Cases 1 (mudflat annuals) and 4 (Scirpus validus) were very rare at the 
	start of the experiment and so were not able to be modeled.
	
	
Arguments: 		struct Plant ***plant		pointer to plant structure
					struct Basin *basin			pointer to basin parameters
					struct Lifehist *lifehist  pointer to lifehist parameters
					struct Run *run            pointer to run parameters

Return Value:	void
****************************************************************************/
void GetProb(double *probptr, int *adltptr, double *dpthptr)
{
	double elev, depth, logit, max;

	depth = *dpthptr;
	elev = (247.5 - depth);

	switch(*adltptr) {
		case 0:
			/**************************************************
			This regression model was fit by EWS and is based on 
			straight elevation numbers rather than adjusted elevations as 
			used by de Swart (1992)
			***************************************************/		
			logit = 1724.5 - 6.9767*elev;
			*probptr = exp(logit)/(1+exp(logit));
			/*Scale the probability to run from 0 to 1 */
			max = 1;
			*probptr /= max;
			break;
		case 1:
			*probptr = 0;
			break;
		case 2:
			logit = -160.386 + 86.257*(elev-244)-11.6030*(elev-244)*(elev-244);
			if(elev >= 247.8)
				logit += 33.35*(elev-247.8)*(elev-247.8);
			if(elev >= 248.1)
				logit += -79.241*(elev-248.1)*(elev-248.1);
			*probptr = exp(logit)/(1+exp(logit));
			/*Scale the probability to run from 0 to 1 */
			max = 0.853;
			*probptr /= max;
			break;
		case 3:
			logit = -581.58 + 317.26*(elev-244)-43.402*(elev-244)*(elev-244);
			*probptr = exp(logit)/(1+exp(logit));
			/*Scale the probability to run from 0 to 1 */
			max = 0.1415;
			*probptr /= max;
			break;
		case 4:
			/**************************************************
			This case is the probability curve for Scirpus validus and 
			is identical to the curve fro Scolochloa (case 3).
			***************************************************/		
			logit = -581.58 + 317.26*(elev-244)-43.402*(elev-244)*(elev-244);
			*probptr = exp(logit)/(1+exp(logit));
			/*Scale the probability to run from 0 to 1 */
			max = 0.1415;
			*probptr /= max;
			break;
		case 5:
			logit = -359.49 + 210.76*(elev-244)-30.9353*(elev-244)*(elev-244);
			if(elev >= 247.7)
				logit += 95.0456*(elev-247.7)*(elev-247.7);
			if(elev >= 247.9)
				logit += -984.93*(elev-247.9)*(elev-247.9);
			*probptr = exp(logit)/(1+exp(logit));
			/*Scale the probability to run from 0 to 1 */
			max = 0.372;
			*probptr /= max;
			break;
		case 6:
			logit = -948.25 + 567.21*(elev-244)-84.8925*(elev-244)*(elev-244);
			*probptr = exp(logit)/(1+exp(logit));
			/*Scale the probability to run from 0 to 1 */
			max = 0.313;
			*probptr /= max;
			break;
		case 7:
			*probptr = 0;
			break;
		default:
			printf("\n\nError in switch statement in GetProb()\nadult = %d", *adltptr);
			exit(EXIT_FAILURE);
			break;
	}
}
/***************************************************************************
KillAnnuals()
	This function kills annual species.
Arguments: 		struct Plant ***plant		pointer to plant structure
					struct Basin *basin			pointer to basin parameters
					struct Lifehist *lifehist  pointer to lifehist parameters
					struct Run *run            pointer to run parameters

Return Value:	void
****************************************************************************/
void KillAnnuals(struct Basin *basin, struct Lifehist *lifehist, struct Plant ***plant)
{
	int x, y;

	for(y=0; y < basin->ymax[0]; y++) {
		for(x=0; x < basin->xmax[0]; x++) {
			if(plant[x][y][0].adult[0] == -9999)
				continue;
			plant[x][y][0].adult[0] *= lifehist->peren[plant[x][y][0].adult[0]];
		}
	}
}
/***************************************************************************
KillSeedling()
	This function kills all seedlings in t = 0.
Arguments: 		struct Plant ***plant		pointer to plant structure
					struct Lifehist *lifehist	pointer to lifehistory parameters

Return Value:	void
****************************************************************************/
void	KillSeedling(struct Basin *basin, struct Plant ***plant)
{
	int x, y, i;

	for(y=0; y < basin->ymax[0]; y++) {
		for(x=0; x < basin->xmax[0]; x++) {
			if(plant[x][y][0].adult[0] == -9999)
				continue;
			for(i=1; i < MAXSP; i++) {
				plant[x][y][0].sdlng[i] = 0;
			}
		}
	}
}
/***************************************************************************
SelectAdult()
	This randomly elevates one species of seedling to adult status based on
	the number of seedlings present of each species.
Arguments: 		struct Plant ***plant		pointer to plant structure
					struct Basin *basin			pointer to basin parameters

Return Value:	void
****************************************************************************/
void	SelectAdult(struct Basin *basin, struct Plant ***plant)
{
	int x, y, i, total, select;
   double d;

	for(y=0; y < basin->ymax[0]; y++) {
		for(x=0; x < basin->xmax[0]; x++) {
			if(plant[x][y][0].adult[0] == -9999)
				continue;
			for(i=1, total = 0; i < MAXSP; i++) {
				total += plant[x][y][0].sdlng[i];
			}
			if(total) {
				d=dr250();
				select = 1 + (int)(d*total);
			} else {
				select = 0;
			}
			for(i=1; i < MAXSP; i++) { /*counts a random distance into seedling array to select seedling*/
				select -= plant[x][y][0].sdlng[i];
				if(select < 0) {
					break;
				}
			}
			if(plant[x][y][0].adult[0] == 0 && select != 0) {
					plant[x][y][0].adult[0] = i;
					plant[x][y][0].origx[0] = x;
					plant[x][y][0].origy[0] = y;
					plant[x][y][0].reproduce[0] = 0; /*set time lag for vege. expansion*/
			}
		}
	}
}
/***************************************************************************
WaterDeath()
	This function kills adults based on a normal distribution around
	their mean occurence.
Arguments: 		struct Plant ***plant		pointer to plant structure
					struct Basin *basin			pointer to basin parameters
					struct Lifehist *lifehist  pointer to lifehist parameters
					struct Run *run            pointer to run parameters

Return Value:	void
****************************************************************************/
void WaterDeath(struct Basin *basin, struct Lifehist *lifehist, struct Plant ***plant, struct Run *run, double **topo, double *water)
{
	int x, y, adult;
	double depth, prob, z, zsqr, mean, width, select;

	for(y=0; y < basin->ymax[0]; y++) {
		for(x=0; x < basin->xmax[0]; x++) {
			if(plant[x][y][0].adult[0] == -9999)
				continue;
			adult = plant[x][y][0].adult[0];
			mean = lifehist->depthmean[adult];
			width = lifehist->depthwidth[adult];
			GetDepth(basin, run, topo, water, &depth, x, y);
			GetProb(&prob, &adult, &depth);
			select = dr250();
			if(select > prob && depth > mean) {  /*determine if plant is damaged */
				if(plant[x][y][0].reproduce[0] == 0) {
					plant[x][y][0].adult[0] = 0;   /*kill nonrep. plants*/
				} else {            /*make rep. plants nonrepro.*/
					plant[x][y][0].reproduce[0] = 0;
				}
			} else {                              /*if plant was not damaged it returns to rep. state*/
				if(plant[x][y][0].reproduce[0] == 0) {
					plant[x][y][0].reproduce[0] = 1;
				}
			}
		}
	}
}