#include <iostream>
#include "gr4j.h"
#include "math.h"

#include <fstream>
using namespace std;
//---------------------------------------------------------------------

Cgr4j::Cgr4j()
{
	bSnowModule = false;
	//_InitParms();
	//_CalcUH();
}
//---------------------------------------------------------------------

Cgr4j::Cgr4j(InputData Input)
{
	bSnowModule = false;
	// Simulation(Input);
}
//---------------------------------------------------------------------

Cgr4j::Cgr4j(InputData Input, double* temperature, CSnowModule *SnowModule)
{
	bSnowModule = true;
	this->temperature = temperature;
	this->SnowModule = SnowModule;
}
//---------------------------------------------------------------------
// using the model with an active SNOW MODULE
Cgr4j::Cgr4j(InputData Input, double* temperature, double T_Rain, double T_Melt, double DD_FAC)
{
	bSnowModule = true;
	this->temperature = temperature;
	if (bSnowModule)
		SnowModule = new CSnowModule(temperature, Input.precipitation, T_Rain, T_Melt, DD_FAC, Input.size);
}
//---------------------------------------------------------------------

Cgr4j::~Cgr4j()
{
	// delete all double*
}
//---------------------------------------------------------------------


///////////////////////////////////////////////////////////////////////
//
//			PUBLIC FUNCTION SIMULATION
//
///////////////////////////////////////////////////////////////////////

//---------------------------------------------------------------------
// call this function to make the calculations automatically
void Cgr4j::Simulation(InputData& Input, GR4J_Parameters& Settings)
{
	//-----------------------------------------
	// assign time series data to parameter Input
	this->Input = Input;
	//-----------------------------------------

	// initialize pointers (arrays)
	// calculate time series averages (mean pcp, Q_obs, pET)
	// assign settings to Parameter X
	_InitParms(Settings);
	//-----------------------------------------

	// initialize model for first time step
	// intitialize storages
	_InitStorages(X, LevelProdStore[0], LevelRouteStore[0]);
	//-----------------------------------------

	// Calculate the discrete ordinates of unit hydrographs
	// UH1 and UH2 respectively
	_CalcUH(UH1, UH2, nUH, X.X4);
	//-----------------------------------------

	// simulate streamflow for given period in Input
	for (int i = 0; i < Input.size; i++)
	{
		// !!!! assign values to: ES, WS, PS, and PR
		// !!!! assign LevelProdStore[i]
		if (bSnowModule)
		{
			if (i == 0) // if first time step
				LevelProdStore[i] = _SimProductionStore(X, Input.precipitation[i], Input.pET[i],
													LevelProdStore[i],ES[i], WS[i], PS[i], PR[i]);
			else
				LevelProdStore[i] = _SimProductionStore_Snow(X, SnowModule, i, Input.precipitation[i],
															Input.pET[i], temperature[i],
															LevelProdStore[i-1], ES[i], WS[i], PS[i], PR[i]);

		} else // if snow module is not used
		{
			if (i == 0) // if first time step
				LevelProdStore[i] = _SimProductionStore(X, Input.precipitation[i], Input.pET[i],
													LevelProdStore[i],ES[i], WS[i], PS[i], PR[i]);
			else
				LevelProdStore[i] = _SimProductionStore(X, Input.precipitation[i], Input.pET[i],
													LevelProdStore[i-1],ES[i], WS[i], PS[i], PR[i]);
		}
		//-----------------------------------------

		// !!!! updates current production store values (LevelProdStore[i])
		// !!!! updates PR value (Quantity of water that reaches the routing function)
		percolation[i] = _SimPercolation(LevelProdStore[i], PR[i]);
		//-----------------------------------------

		// simulate contribution from unit hydrographs
		UH1_Contrib = _Sim_UH_Contrib(UH1_Contrib, PR[i], UH1, nUH);
		UH2_Contrib = _Sim_UH_Contrib(UH2_Contrib, PR[i], UH2, nUH*2);
		//-----------------------------------------

		// simulate groundwater exchange
		if (i == 0)
			GWExchange[i] = _SimWaterExchange(X, LevelRouteStore[i]);
		else
			GWExchange[i] = _SimWaterExchange(X, LevelRouteStore[i-1]);
		//-----------------------------------------

		// calculate contribution from routing store (Q_r[i])
		// !!!! updates current routing store value (LevelRouteStore[i])
		if (i == 0)
			Q_r[i] = _SimReservoirOutflow(LevelRouteStore[i], LevelRouteStore[i], UH1_Contrib[0], GWExchange[i]);
		else
			Q_r[i] = _SimReservoirOutflow(LevelRouteStore[i], LevelRouteStore[i-1], UH1_Contrib[0], GWExchange[i]);
		//-----------------------------------------

		// simulate contribution from second unit hydrograph 2
		Q_d[i] = _Sim_UH2_ComponentOutflow(UH2_Contrib[0], GWExchange[i]);
		//-----------------------------------------

		// simulate total streamflow
		Q_sim[i] = Q_r[i] + Q_d[i];
	}
}

// call free_memory() after Monte Carlo Simulations

//---------------------------------------------------------------------

///////////////////////////////////////////////////////////////////////
//
//		SIMULATE STREAMFLOW
//
///////////////////////////////////////////////////////////////////////

//---------------------------------------------------------------------

// simulate water level in production store
double Cgr4j::_SimProductionStore(GR4J_Parameters X, double pcp, double pet,
								  double LevelProdStore,
								  double& ES, double& WS,
								  double& PS, double& PR)
{
	ES = 0.0;		// quantity of water that will evaporate from production store
	WS = 0.0;		// netto pet / X.X1 (capacity of production store)
	PS = 0.0;		// part of netto precipitation that fills the production store
	PR = 0.0;		// netto precipitation - PS

	if (pcp >= pet)
	{
		WS = (pcp - pet) / X.X1;
		if (WS > 13.0) WS = 13.0;
		PS = X.X1 * (1.0 - pow(LevelProdStore / X.X1, 2.0)) * tanh(WS)
			 / (1.0 + LevelProdStore / X.X1 * tanh(WS));
		PR = pcp - pet - PS;
	} else {
		WS = (pet - pcp) / X.X1;
		if (WS > 13.0) WS = 13.0;
		ES = LevelProdStore * (2.0 - LevelProdStore / X.X1) * tanh(WS) 
			 / (1.0 + LevelProdStore / X.X1 * tanh(WS));
		PS = 0.0;
		PR = 0.0;
	}
	return(LevelProdStore - ES + PS);
}
//---------------------------------------------------------------------

// simulate water level in production store
// taking snow melt into account using the SNOW MODULE
double Cgr4j::_SimProductionStore_Snow(GR4J_Parameters X, 
										CSnowModule* SnowModule, int current_day,
										double pcp, double pet, double temperature,
										double LevelProdStore,
										double& ES, double& WS,
										double& PS, double& PR)
{
	ES = 0.0;		// quantity of water that will evaporate from production store
	WS = 0.0;		// netto pet / X.X1 (capacity of production store)
	PS = 0.0;		// part of netto precipitation that fills the production store
	PR = 0.0;		// netto precipitation - PS

	if (pcp >= pet)
	{
		WS = (pcp - pet) / X.X1;
		if (WS > 13.0) WS = 13.0;
		PS = X.X1 * (1.0 - pow(LevelProdStore / X.X1, 2.0)) * tanh(WS)
			 / (1.0 + LevelProdStore / X.X1 * tanh(WS));
		PR = pcp - pet - PS;
	} else {
		WS = (pet - pcp) / X.X1;
		if (WS > 13.0) WS = 13.0;
		ES = LevelProdStore * (2.0 - LevelProdStore / X.X1) * tanh(WS) 
			 / (1.0 + LevelProdStore / X.X1 * tanh(WS));
		PS = 0.0;
		PR = 0.0;
	}

	if (temperature < SnowModule->Get_T_Rain())
	{
		PS = 0.0;
		PR = 0.0;
	}

	if (temperature > SnowModule->Get_T_Melt())
	{
		if (pcp + SnowModule->Get_MeltRate(current_day) >= pet)
		{
			WS = (pcp + SnowModule->Get_MeltRate(current_day) - pet) / X.X1;
			if (WS > 13.0) WS = 13.0;
			PS = X.X1 * (1.0 - pow(LevelProdStore / X.X1, 2.0)) * tanh(WS)
				 / (1.0 + LevelProdStore / X.X1 * tanh(WS));
			PR = pcp + SnowModule->Get_MeltRate(current_day) - pet - PS;
		} else {
			WS = (pet - pcp + SnowModule->Get_MeltRate(current_day)) / X.X1;
			if (WS > 13.0) WS = 13.0;
			ES = LevelProdStore * (2.0 - LevelProdStore / X.X1) * tanh(WS) 
				 / (1.0 + LevelProdStore / X.X1 * tanh(WS));
			PS = 0.0;
			PR = 0.0;
		}
	}

	// test
	if (temperature < SnowModule->Get_T_Melt() && temperature > SnowModule->Get_T_Rain())
	{
		if (pcp + SnowModule->Get_MeltRate(current_day) >= pet)
		{
			WS = (pcp + SnowModule->Get_MeltRate(current_day) - pet) / X.X1;
			if (WS > 13.0) WS = 13.0;
			PS = X.X1 * (1.0 - pow(LevelProdStore / X.X1, 2.0)) * tanh(WS)
				 / (1.0 + LevelProdStore / X.X1 * tanh(WS));
			PR = pcp + SnowModule->Get_MeltRate(current_day) - pet - PS;
		} else {
			WS = (pet - pcp + SnowModule->Get_MeltRate(current_day)) / X.X1;
			if (WS > 13.0) WS = 13.0;
			ES = LevelProdStore * (2.0 - LevelProdStore / X.X1) * tanh(WS) 
				 / (1.0 + LevelProdStore / X.X1 * tanh(WS));
			PS = 0.0;
			PR = 0.0;
		}
	}
	/*
	if (temperature > SnowModule->Get_T_Melt())
		PR += SnowModule->Get_MeltRate(current_day);
	if (temperature < SnowModule->Get_T_Melt() && temperature > SnowModule->Get_T_Rain())
		PR += SnowModule->Get_MeltRate(current_day);
	*/

	return(LevelProdStore - ES + PS);


}

//---------------------------------------------------------------------

double Cgr4j::_SimPercolation(double& LevelProdStore, double& PR)
{
	double perc = 0.0;
	double s_temp = 0.0;
	s_temp = LevelProdStore / (pow(1.0 + pow(LevelProdStore / 2.25 / X.X1, 4),0.25));
	perc = LevelProdStore - s_temp;
	LevelProdStore = s_temp;
	PR += perc;

	return(perc);
}
//---------------------------------------------------------------------

double* Cgr4j::_Sim_UH_Contrib(double* UHx_Contrib, double PR, double* UH, int size)
{
	for (int i = 0; i < size-1; i++)
		UHx_Contrib[i] = UHx_Contrib[i+1] + UH[i] * PR;
	UHx_Contrib[size-1] = UH[size-1] * PR;
	return(UHx_Contrib);
}
//---------------------------------------------------------------------

double Cgr4j::_SimWaterExchange(GR4J_Parameters X, double LevelRouteStore)
{
	// Calculating the groundwater exchange term
	return(X.X2 * pow(LevelRouteStore / X.X3, 3.5));
}
//---------------------------------------------------------------------

double Cgr4j::_SimReservoirOutflow(double& LevelRouteStore, double LevelRouteStore_before, double uh1_1, double gw_exchange)
{
	double Q_res = 0.0;
	double R2 = 0.0;
	// updating the level in the routing store
	LevelRouteStore = max(0.0, LevelRouteStore_before + uh1_1 * 0.9 + gw_exchange);
	R2 = LevelRouteStore / ( pow(1 + pow(LevelRouteStore / X.X3, 4.0),0.25) );
	Q_res = LevelRouteStore - R2;
	LevelRouteStore = R2;

	return(Q_res);
}
//---------------------------------------------------------------------

double Cgr4j::_Sim_UH2_ComponentOutflow(double uh2_1, double GWExchange)
{
	return(max(0.0, uh2_1 * 0.1 + GWExchange));
}
//---------------------------------------------------------------------


///////////////////////////////////////////////////////////////////////
//
//			UNIT HYDROGRAPH FUNCTIONS
//
///////////////////////////////////////////////////////////////////////
//---------------------------------------------------------------------
void Cgr4j::_CalcUH(double* UH1, double* UH2, int nUH, double X4)
{
	// Calculate the discrete ordinates of unit hydrographs
	// UH1 and UH2 respectively
	UH1 = _UH1(UH1, nUH, X4);
	UH2 = _UH2(UH2, nUH*2, X4);
}
//---------------------------------------------------------------------
double Cgr4j::_SS1(double count, double tc)
{
	// tc = unit hydrograph time constant
	if (count < 0)		return(0);
	if (count < tc)		return(pow(double(count) / tc, 2.5));
	if (count >= tc)	return(1);
	return(0);
}
//---------------------------------------------------------------------
double Cgr4j::_SS2(double count, double tc)
{
	// tc = unit hydrograph time constant
	if (count < 0)		return(0);
	if (count <= tc)	return(0.5 * pow(count / tc, 2.5));
	if (count < 2*tc)	return(1 - 0.5 * pow(2 - count / tc, 2.5));
	if (count >= 2*tc)	return(1);
	return(0);
}
//---------------------------------------------------------------------
double* Cgr4j::_UH1(double* UH, int nuh, double tc)
{
	for (int i = 0; i < nuh; i++)
	{
		UH[i] = _SS1(double(i+1), tc) - _SS1(double(i), tc);
	}
	return(UH);
}
//---------------------------------------------------------------------
double* Cgr4j::_UH2(double* UH, int nuh, double tc)
{
	for (int i = 0; i < nuh; i++)
	{
		UH[i] = _SS2(double(i+1), tc) - _SS2(double(i), tc);
	}
	return(UH);
}

///////////////////////////////////////////////////////////////////////
//
//		INITIALIZING PARAMETERS
//
///////////////////////////////////////////////////////////////////////
//---------------------------------------------------------------------
void Cgr4j::_InitParms(GR4J_Parameters Settings)
{
	// initialize input parameters
	X.X1 = Settings.X1;
	X.X2 = Settings.X2;
	X.X3 = Settings.X3;
	X.X4 = Settings.X4;

	// initialize unit hydrograph parameters
	nUH		= 10;
	UH1		= new double[nUH];
	UH2		= new double[nUH*2];
	UH1_Contrib = new double[nUH];
	UH2_Contrib = new double[nUH*2];
	for (int i = 0; i < nUH; i++)
		UH1_Contrib[i] = 0.0;
	for (int i = 0; i < nUH*2; i++)
		UH2_Contrib[i] = 0.0;

	// initialize arrays
	ES = new double[Input.size];
	WS = new double[Input.size];
	PS = new double[Input.size];
	PR = new double[Input.size];
	Q_r = new double[Input.size];
	Q_d = new double[Input.size];
	Q_sim = new double[Input.size];
	percolation = new double[Input.size];
	LevelProdStore = new double[Input.size];
	LevelRouteStore = new double[Input.size];
	GWExchange = new double[Input.size];
	for (int i = 0; i < Input.size; i++)
	{
		ES[i] = 0.0;
		WS[i] = 0.0;
		PS[i] = 0.0;
		PR[i] = 0.0;
		Q_r[i] = 0.0;
		Q_d[i] = 0.0;
		Q_sim[i] = 0.0;
		GWExchange[i] = 0.0;
		percolation[i] = 0.0;
		LevelProdStore[i] = 0.0;
		LevelRouteStore[i] = 0.0;
	}

	// calculate time series averages
	mean_pcp	= _SumArray(Input.precipitation, Input.size) / Input.size;
	mean_Q		= _SumArray(Input.Q_obs, Input.size) / Input.size;
	mean_pet	= _SumArray(Input.pET, Input.size) / Input.size;	
}
//---------------------------------------------------------------------

void Cgr4j::_InitStorages(GR4J_Parameters X, double& LevelProdStore, double& LevelRouteStore)
{
	double QB = 0.9 * mean_Q / X.X3;
	LevelRouteStore = 1.0 - (1.0 /(4.0 * pow(1.0 + QB,4.0) 
					  - 3.0 + pow(4.0 * QB, 0.2)
					  + pow(8.0 * QB, 0.4)));
	LevelRouteStore *= X.X3;
	LevelProdStore = X.X1 * min(1.0, sqrt(mean_Q / mean_pcp));
}
//---------------------------------------------------------------------

// If this class is used during a Monte-Carlo Simulation, this function
// must be called after each simulation. Otherwise a memory overflow will
// be caused.
void Cgr4j::free_memory()
{
	delete[] UH1;
	delete[] UH2;
	delete[] UH1_Contrib;
	delete[] UH2_Contrib;
	delete[] ES;
	delete[] WS;
	delete[] PS;
	delete[] PR;
	delete[] Q_r;
	delete[] Q_d;
	delete[] percolation;
	delete[] LevelProdStore;
	delete[] LevelRouteStore;
	delete[] GWExchange;
	delete[] Q_sim;
}

///////////////////////////////////////////////////////////////////////
//			HELPER
///////////////////////////////////////////////////////////////////////

double Cgr4j::_SumArray(double* val, int size)
{
	double sum = 0;
	for (int i = 0; i < size; i++)
		sum += val[i];
	return(sum);
}

//---------------------------------------------------------------------
//			    Assign first and last record (user selection)
//---------------------------------------------------------------------
void Cgr4j::AssignFirstLastRec(CSG_Table &pTable, int &first, int &last,
									 CSG_String date1, CSG_String date2, int dateField)
{
	int		j;
	///////////////////////////////////////////////////////////////
	// searching the first and the last record of the time range
	///////////////////////////////////////////////////////////////

	for (j = 0; j < pTable.Get_Record_Count(); j++)
	{
		if (!strcmp(pTable.Get_Record(j)->asString(dateField), date1))
		{
			first = j;
		}
		else if (!strcmp(pTable.Get_Record(j)->asString(dateField), date2))
		{
			last = j;
		}
	}
}