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

/***************************************************************************
GetDepth()
	This function calculates the depth value of a given (x,y) pair in the
	basin matrix given the current water level. Negative values indicate
	depth to standing water.
Arguments: 		struct Basin *basin			pointer to basin parameters
					double *water					pointer to the current water elevation
					double *depth					pointer to the depth of a cell
					int x								x location
					int y								y location

Return Value:	void
****************************************************************************/
void GetDepth(struct Basin *basin, struct Run *run, double **topo, double *water, double *depth, int x, int y)
{
	double d;

	GetElev(basin, run, topo, &d, x, y);
	if(d!=-9999) {					/*check for blank cells*/
		if(*water > d) {				/*points that have standing water*/
			*depth = (*water-d);
		} else {							/*points that don't have standing water*/
			*depth = -(d-*water);   /*returns depth to water table*/
		}
	} else {
		*depth = d;
	}
}
/***************************************************************************
GetElev()
	This function calculates a z value for a given (x,y) pair. 
	
	WARNING: In GetTopo() elevations are 
	standardized relative to the maximum flooding level in the basins 
	(248.5 AMSL) so that positive values are never flooded. Normal level is 
	1 m (247.5 AMSL). 
Arguments: 		struct Basin *basin			pointer to basin parameters
					double *zmin					pointer to the zmin of a cell
					int x								x location
					int y								y location
					struct Run *run      pointer to run parameters

Return Value:	void
****************************************************************************/
void	GetElev(struct Basin *basin, struct Run *run, double **topo, double *z, int x, int y)
{
	double a, c, d, r, k, i, j;

	if(run->basinfile[0]) {
		*z = topo[x][y];
	} else {
		a = (double)basin->area[0];
		r = basin->xyratio[0];
		d = -basin->zmin[0];
		c = basin->curve[0];
		i = (double)(y - basin->ymax[0]/2);
		j = (double)(x - basin->xmax[0]/2);
		k = (pow((d)*PI/a*((j*j)/r+(i*i)*r),c)-d);
		if (k > basin->zmax[0]) {
			*z = (double)(basin->zmax[0]);
		} else {
			*z = k;
		}
	}
}
/***************************************************************************
GetSize()
	This function gets the dimensions of the basin matrix based on the
	parameters of the basin.
Arguments: 		struct Basin *basin		pointer to basin parameters

Return Value:	void
****************************************************************************/
void GetSize(struct Basin *basin)
{
	double a, b, c, d, r, zmax, cell;

	a = (double)basin->area[0];
	b = (double)basin->buffer[0];
	r = basin->xyratio[0];
	d = -basin->zmin[0];
	c = basin->curve[0];
	zmax = basin->zmax[0];
	cell = basin->cellsize[0];


	basin->ymax[0] = (int) ((2*(b+(sqrt(pow(zmax+d,1/c)*a/(d*PI*r)))))/cell);
	basin->xmax[0] = (int) ((2*(b+(sqrt(pow(zmax+d,1/c)*a*r/(d*PI)))))/cell);
	if(basin->ymax[0]%2 == 0) {
		basin->ymax[0]++;
	}
	if(basin->xmax[0]%2 == 0) {
		basin->xmax[0]++;
	}
}
/***************************************************************************
GetWaterElev()
	This function calculates the elevation of standing water at 4 seasons
	using a linear decrease in elevation.

Arguments: 		double *water			pointer to the water depth array
					struct Basin *basin	pointer to basin parameters
					struct Run *run      pointer to run parameters

Return Value:	void
****************************************************************************/
void GetWaterElev(struct Basin *basin, double *water, int step)
{
	double amp, freq, slope, season, year;
	int yearint;

	slope = -(double)basin->yearint[0]/3;
	freq = 2*(double)basin->drawfreq[0];
	amp = (double)basin->drawint[0];
	season = (double)(step % 4);
	year=(((double)step-(double)season)/4);
	yearint = (int)year;

	if(basin->watercos[0]) {
#ifdef STARTLOW
		*water = ((-amp/2)*(1-(cos(2*PI*(year-(freq/2))/freq))))+season*slope;  /*starts run at minwater*/
#else
		*water = ((-amp/2)*(1-(cos(2*PI*year/freq))))+season*slope; 				/*starts at water ele = 0*/
#endif
	} else {
		if(yearint % basin->drawfreq[0] != 0 || year == 0) {
			*water = season*slope;
		} else {
			if(season == 0) {
				printf("\n###### DRAWDOWN YEAR  ########");
			}
			*water = -basin->drawint[0]+(season*slope);
		}
	}
	if(*water < basin->zmin[0]) {
		*water = basin->zmin[0];
	}
#ifdef DEBUG
		printf("\n");
		printf("\nLevel %i = %2.3f", season, *water);
		if(season == 3) {
			ScreenPause();
		}
#endif
}