/***************************************************************
								WINDISPER:
	A SPATIOTEMPORAL MODEL FOR TREE SEED DISPERSAL BY WIND
						 Version 1.0.10.1999
							by Ran Nathan
****************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <malloc.h>
#include <time.h>
#include "getdata.h"
#include "parse.h"
#include "rand.h"

/* GENERAL CONSTANTS */
#define PI 	                                3.1415927
#define KARMAN                              0.4
/*  Von Karmam constant of the logarithmic wind profile equation.
	In most source = 0.40, but in some others = 0.41 */

site_t *site;

int nodata=-9999;
unsigned long int xllcorner=0, yllcorner=0;
/*  Default variables required to the ARC/Info input/output files */
float Wind_Array[WIND_SECTORS][WIND_ARRAY_PARAMETERS];
/*  An array that stores directional wind data */

int Cell_Width;
float Seeds_Inside_This_Period=0;
float Seeds_Inside_All_Periods=0;
float Seeds_Outside_This_Period=0;
float Seeds_Outside_All_Periods=0;
float Max_DisDis_This_Period=0;
float Max_DisDis_All_Periods=0;

int Rows;
int Columns;

BOOL GetGrid;
BOOL PeriodOut;
BOOL WindRandom;

float Mean_Tree_Height;
float Std_Tree_Height;
float Mean_Portion_Height;
float Std_Portion_Height;
float Max_Portion_Height;
float Min_Portion_Height;
float Mean_Terminal_Velocity;
float Std_Terminal_Velocity;
float Portion_Wingless;
float Threshold_Windspeed;
float Mean_Vertical_Windspeed;
float Std_Vertical_Windspeed;
float Site_Roughness;
float Site_Displacement;
float WS_Roughness;
float WS_Displacement;
float Measurement_Height;

char Period_Input_File[1000]=DEFAULT_PERIOD_INPUT_FILE;
/* Period Parameter File Name */
char Period_Output_File[1000]=DEFAULT_PERIOD_OUTPUT_FILE;
/* Period Output File Name */

void Start(void);
void Reset_Landscape(void);
void Dispersal(char*);
float Wind(float*, float*, float*, float*);
float Dispersal_Distance(float*, float*, float*, float*, float*);
float Release_Height(void);
float Terminal_Velocity(void);
float Vertical_Windspeed(void);
void Deposition_Cell(int, int, float, float);
void This_Period_Output(char *);
void All_Periods_Output(void);
void Add_and_Reset(void);


void main()
{
	int period, number_of_periods;
	char period_filename[1000]; /*period file name */

	Start();
	number_of_periods = Get_All_Periods_Data();

	if (GetGrid)
		Get_Trees_Grid();
	else
		Get_Trees_XY();

	for (period=1; period<=number_of_periods; period++) {
		printf("\n\nPeriod %d of %d:", period, number_of_periods);
		sprintf(period_filename,"%s%d",Period_Input_File,period);
		/* generates file name for each period */
		Dispersal(period_filename);
		if (PeriodOut) {
			This_Period_Output(period_filename);
			Add_and_Reset();
		}
	}
	All_Periods_Output();

	printf("\n\nDone!\n");

	free(site);
}

void Start(void)
{
	 printf("\n                            WINDISPER\n");
	 printf("         A SPATIOTEMPORAL MODEL FOR TREE SEED DISPERSAL BY WIND\n");
	 printf("                          Version 1.0.10.1999\n");
	 printf("                            by Ran Nathan\n");
}

void Dispersal(char *period)
{
	int   i=0, j=0, k=0;
	float SeeCro=0.0, DisDis=0.0, DaysPeriod;
	float MeanSeeRel, StdSeeRel, MeanLnWin, StdLnWin, MaxLnWin, MinLnWin;
	float WinDir;
	time_t t;

	idum = -time(&t);

	DaysPeriod = Get_This_Period_Data(period, &MeanSeeRel, &StdSeeRel,
								&MeanLnWin, &StdLnWin, &MaxLnWin, &MinLnWin);
	printf("\n");
	for (i = 0; i < Rows; i++) {
		printf("\rWait... %3.0f%% completed\r", 100.0*(i+1)/Rows );
		for (j = 0; j < Columns; j++) {
			if (site[i*Columns + j].source) {
			/* for each source cell... */
				SeeCro = floor(DaysPeriod*Normal_Random(MeanSeeRel, StdSeeRel)+0.5);
				if (SeeCro < 0.0)
					SeeCro = 0.0;
				for (k = 0; k < SeeCro; k++) {
				/* for each dispersed seed... */
					DisDis = Dispersal_Distance(&MeanLnWin, &StdLnWin, &MaxLnWin,
									&MinLnWin, &WinDir);
					if ( DisDis > Max_DisDis_This_Period )
						Max_DisDis_This_Period = DisDis;
					Deposition_Cell(i, j, DisDis, WinDir);
				}
			}
		}
	}
}

float Dispersal_Distance(float* MeanLnWin, float* StdLnWin, float* MaxLnWin,
	float* MinLnWin, float *WinDir)
{
	float H=0.0, F=0.0, W=0.0, Uz=0.0, DisDis=0.0;
	float DisHei, RouLen, FriVel, temp;

	DisHei = Site_Displacement;
	RouLen = Site_Roughness;

	if (Random2(&idum) <= Portion_Wingless)
		DisDis = 0.0;
	else {
		H = Release_Height();
		if ( H <= (DisHei+RouLen) )
	/* 	A limitation of the logarithmic wind profile */
			DisDis = 0.0;
		else {
			do {
				do {
					F = Terminal_Velocity();
					W = Vertical_Windspeed();
				} while( (F - W) <= 0.0 );
				/* 	This is the W < F constraint to assure finite DisDis */
				*WinDir = Wind(MeanLnWin, StdLnWin, MaxLnWin, MinLnWin);
				do {
					temp = Normal_Random(*MeanLnWin, *StdLnWin);
				} while (temp < *MinLnWin || temp > *MaxLnWin);
				Uz = exp(temp);
				/*  Selects a random windspeed, from a lognormal distribution */
				FriVel = Uz * KARMAN / (log((Measurement_Height - WS_Displacement)/
									WS_Roughness));
				/*  Claculates the friction velocity in the reference station,
					assuming	it is equal to that of the simulated site (Fields &
					Sharpe 1980; p. 18).*/

				if (Threshold_Windspeed < exp(*MinLnWin))
					Threshold_Windspeed = exp(*MinLnWin);
				Uz = ( FriVel/KARMAN ) * log( (H - DisHei)/RouLen );
			} while ( Uz < Threshold_Windspeed && Uz > exp(*MaxLnWin) );

			DisDis =	(FriVel/(KARMAN*(F-W)))*( (H-DisHei)*log((H-DisHei)/
				(exp(1)*RouLen) ) + RouLen );
		}
	}

	return (DisDis);
}

float Wind(float *MeanLnWin, float *StdLnWin, float *MaxLnWin, float *MinLnWin)
{
	float lottery=0.0, temp=0.0, WinDir=0;
	int i=0;

	lottery = Random2(&idum);

	if (WindRandom)
	/*	random wind direction, assuming equal probability for all directions */
		while ( i < WIND_SECTORS && lottery >= (i+0.5)/WIND_SECTORS ) i++;
	else {
	/* for specified wind data */
		while ( i < WIND_SECTORS && lottery >= Wind_Array[i][0] ) i++;
		*MeanLnWin = Wind_Array[i][1];
		*StdLnWin = Wind_Array[i][2];
		*MaxLnWin = Wind_Array[i][3];
		*MinLnWin = Wind_Array[i][4];
	}

	temp = (2*i + Random2(&idum) - 0.5)*PI/WIND_SECTORS;
	if (temp < 0 )
		WinDir = 2*PI - abs(temp);
	else
		if (temp > 2*PI)
			WinDir = temp - 2*PI;
		else
			WinDir = temp;

	return (WinDir);
}

float Release_Height(void)
{
	float RelHei, tree_height, portion;

	do {
		tree_height = Normal_Random(Mean_Tree_Height, Std_Tree_Height);
	} while (tree_height < 0.0);

	do {
		portion = Normal_Random(Mean_Portion_Height, Std_Portion_Height);
	} while (portion > Max_Portion_Height || portion < Min_Portion_Height);

	RelHei = tree_height * portion;

	return (RelHei);
}

float Terminal_Velocity(void)
{
	float TerVel;

	do
		TerVel = Normal_Random(Mean_Terminal_Velocity, Std_Terminal_Velocity);
	while (TerVel <= 0.0);

	return (TerVel);
}

float Vertical_Windspeed(void)
{
	float VerWin;

	VerWin = Normal_Random(Mean_Vertical_Windspeed, Std_Vertical_Windspeed);

	return VerWin;
}

void Deposition_Cell(int i, int j, float DisDis, float WinDir)
{
	int  dep_i=0, dep_j=0;

	dep_i = floor(0.5 + i - DisDis*cos(WinDir)/Cell_Width);
	dep_j = floor(0.5 + j + DisDis*sin(WinDir)/Cell_Width);

	if ( ( (dep_i>=0) && (dep_j>=0) ) && ( (dep_i<Rows) && (dep_j<Columns ) ) ) {
		site[dep_i*Columns + dep_j].this_period_seeds ++;
		Seeds_Inside_This_Period++;
	}
	else
		Seeds_Outside_This_Period++;
}

void This_Period_Output(char *period)
/*
	Generates two output files:
	1. General info on seed dispersal simulation for this period
	2. An ARC/info grid file for model's prediction of the periodic
		post-dispersal	distribution of seeds within the simulated landscape.
*/
{
	FILE *fp;
	int i, j;
	char FileName[1000];

	sprintf(FileName, "%s.out", period);
	fp = fopen(FileName, "w");
	printf("\nOutput file '%s' opened", FileName);

	fprintf(fp, "\nTotal #seeds released: %.0f", Seeds_Inside_This_Period +
																Seeds_Outside_This_Period);
	fprintf(fp, "\nTotal #seeds deposited within landscape: %.0f",
																Seeds_Inside_This_Period );
	fprintf(fp, "\nTotal #seeds deposited outside landscape: %.0f",
																Seeds_Outside_This_Period );
	fprintf(fp, "\nFraction of seeds deposited within landscape: %.2f%",
		Seeds_Inside_This_Period * 100.0 / (Seeds_Inside_This_Period +
																Seeds_Outside_This_Period) );
	fprintf(fp, "\nMaximum dispersal distance: %.2f", Max_DisDis_This_Period );

	fclose(fp);
	printf("\nOutput file '%s' closed\n", FileName);

	sprintf(FileName, "%s.grd", period);
	fp = fopen(FileName, "w");
	printf("\nOutput file '%s' opened", FileName);

	fprintf(fp, "ncols         %d\n", Columns);
	fprintf(fp, "nrows         %d\n", Rows);
	fprintf(fp, "xllcorner     %lu\n", xllcorner);
	fprintf(fp, "yllcorner     %lu\n", yllcorner);
	fprintf(fp, "cellsize      %d\n", Cell_Width);
	fprintf(fp, "NODATA_value  %d", nodata);
	for (i = 0; i < Rows; i++) {
		fprintf(fp, "\n", i);
		for (j = 0; j < Columns; j++)
			fprintf(fp, "%d ", site[i*Columns + j].this_period_seeds);
	}

	fclose(fp);
	printf("\nOutput file '%s' closed\n", FileName);
}

void All_Periods_Output(void)
/*
	Generates two output files:
	1. General info on seed dispersal simulation for all periods
	2. An ARC/info grid file for model's prediction of the total
		post-dispersal	distribution of seeds within the simulated landscape.
*/
{
	FILE *fp;
	int i, j;

	fp = fopen("all_runs.txt", "w");
	printf("\nOutput file 'all_runs.txt' opened");

	fprintf(fp, "\nTotal #seeds released: %.0f", Seeds_Inside_All_Periods +
																Seeds_Outside_All_Periods);
	fprintf(fp, "\nTotal #seeds falling within landscape: %.0f",
																Seeds_Inside_All_Periods );
	fprintf(fp, "\nTotal #seeds falling outside landscape: %.0f",
																Seeds_Outside_All_Periods );
	fprintf(fp, "\nFraction of seeds falling within landscape: %.2f%",
		Seeds_Inside_All_Periods * 100.0 / (Seeds_Inside_All_Periods +
																Seeds_Outside_All_Periods) );
	fprintf(fp, "\nMaximum dispersal distance: %.2f", Max_DisDis_All_Periods );

	fclose(fp);
	printf("\nOutput file 'all_runs.txt' closed\n");

	fp = fopen("all_runs.grd", "w");
	printf("\nOutput file 'all_runs.grd' opened");

	fprintf(fp, "ncols         %d\n", Columns);
	fprintf(fp, "nrows         %d\n", Rows);
	fprintf(fp, "xllcorner     %lu\n", xllcorner);
	fprintf(fp, "yllcorner     %lu\n", yllcorner);
	fprintf(fp, "cellsize      %d\n", Cell_Width);
	fprintf(fp, "NODATA_value  %d", nodata);
	for (i = 0; i < Rows; i++) {
		fprintf(fp, "\n", i);
		for (j = 0; j < Columns; j++)
			fprintf(fp, "%d ", site[i*Columns + j].all_periods_seeds);
	}

	fclose(fp);
	printf("\nOutput file 'all_runs.grd' closed\n");
}

void Add_and_Reset(void)
{
	int i;

	for (i = 0; i < Columns*Rows; i++) {
		site[i].all_periods_seeds += site[i].this_period_seeds;
		site[i].this_period_seeds = 0;
	}

	Seeds_Inside_All_Periods += Seeds_Inside_This_Period;
	Seeds_Outside_All_Periods += Seeds_Outside_This_Period;
	Seeds_Inside_This_Period = 0;
	Seeds_Outside_This_Period = 0;

	if (Max_DisDis_All_Periods < Max_DisDis_This_Period)
		Max_DisDis_All_Periods = Max_DisDis_This_Period;
	Max_DisDis_This_Period = 0;
}