#include "FlowRouteCalc.h"
#include <math.h>
#include <fstream>
#include <vector>
#include <iostream>
#include <numeric>
#include <string>

// Caution: diffusion hydrograph runoff routing function has the potential to create or destroy water if not properly calibrated. Safest is use nonrouted runoff.

//Hydrologically routes runoff depths from overland impervious area, overland pervious area, and subsurface pervious area. 
//All inputs runoff depths, PAQi_m, DCIAQi_m, SubQi_m were transformed to folder area (e.g. bulk area) in FlowSumCalc.cpp or earlier. 

//2-parameter diffusion hydrograph runoff routing developed by Yang and Endreny (2013), replacing the time-area hydrograph developed by Wang et al. (2008)
//Yang, Yang, Theodore A. Endreny. "Watershed hydrograph model based on surface flow diffusion", Water Resources Research Vol 49:507-516 (2013): 
//Wang, Jun, Theodore A. Endreny, and David J. Nowak. "Mechanistic Simulation of Tree Effects in an Impervious Water Balance Model1." JAWRA Journal of the American Water Resources Association Vol(44)1: 75-85 (2008)
//Criss, and Winston (2008a). Properties of a diffusive hydrograph and the interpretation of its single parameter. Mathematical Geosciences, 40(3), 313-325. doi:10.1007/s11004-008-9145-9
//Criss andWinston (2008b). Discharge predictions of a rainfall-driven theoretical hydrograph compared to common models and observed data. Water Resources Research, 44(10). doi:10.1029/2007wr006415

//The 2-parameter diffusion model has 2 time-based parameters, alpha and beta, representing watershed scale flow diffusivity and flow celerity
//This 2-parameter diffusion model has excellent representation of natural hydrograph characteristics, including the inflection points in the rising and falling limbs, a fast rising limb and a gradual falling limb, and time to peak.
//  The 2 parameters alpha and beta represent watershed and storm influences on runoff and hydrograph features
//  If surface flow and subsurface flow have same time constants alpha and beta, the model can simulate discharge with mixed surface flow and subsurface flow.
//  If surface flow and subsurface flow have different time constants alpha and beta, a parallel model that consists of two two-parameter models can be constructed to simulate the hydrograph
// The 2-parameter surface flow diffusion hydrograph model is used for PerviousFlow (PAQ_m) and Impervious Flow (DCIAQ_m).
// The 1-parameter subsurface flow diffusion hydrograph model is used for Baseflow (SubQ_m), based on Criss and Winston (2008a).
// The time input, j, is used in calculating the discharge from the precip (pulse input) at time j in the future

// Notes on alpha and beta, from Yang and Endreny (2013):
/*
	A visual demonstration of the influence of alpha and beta on the normalized hydrograph is shown in Figure 4 of Yang and Endreny (2013).
	Setting beta = 1, and varying alpha from 0.1 and 10, alpha and beta together determine the rising limb, tmax, and falling limb when (Figure 4a);
	The smaller alpha, the faster the rising and falling limb and the smaller tmax.
	Setting beta = 1 and alpha >10, beta determines the hydrograph shape and tmax = 2(beta/3).
	Setting alpha = 1 and varying beta from 0.1 to 100, the normalized hydrographs series indicate that beta always contributes to tmax (Figure 4b).
	The larger beta, the longer the hydrograph duration.
	If alpha goes to infinity, kinematic flow celerity = 0 and the 2-parameter model takes the same form as the 1-parameter model of Criss and Winston (2008b)
*/

void FlowRouteCalc::calculateStatistical(Inputs* input, int timeStep)
{
	//cte 20200918 developing an alternative routing option to the Yang and Endreny (2013) 2-parameter diffusion hydrograph
	//	when flag = 1 use 2-parameter diffusion hydrograph, flag = 0 skip routing
	double routing_flag = 1;

	int IT = timeStep;
	int IIN = 0;	// The for-loop controller - equals to the timestep + the instance of the for-loop
	auto const stride = 3;

	//TotalTimeSteps (hr) is defined in HydroPlusConfig.xml key or redefined in Inputs::readHydroPlusConfigFile() if key = 0 or empty
	//	totalTs is time step for input RiverTimeSeries.txt, RiverSWTimeSeries.txt and RiverGWTimeSeries.txt
	int totalTs = input->SimulationNumericalParams["TotalTimeSteps"];

	//If timeStep equals 0, or first time step, then dimensionalize vectors
	if (timeStep == 0) {

		//Dimensionalize size for catchment runoff depth (m) vectors as time distributed vector values
		//Note: These four input-> vector double variables were created in Inputs.h

		//Pervious area runoff, including runon from impervious area (m)
		input->PAQi_Routed_m.insert(input->PAQi_Routed_m.end(), totalTs, 0.0);
		//Directly connected impervious area runoff (m)
		input->DCIAQi_Routed_m.insert(input->DCIAQi_Routed_m.end(), totalTs, 0.0);
		//Subsurface runoff (m)
		input->SubQi_Routed_m.insert(input->SubQi_Routed_m.end(), totalTs, 0.0);
		//Total runoff (m)
		input->Qi_Routed_m.insert(input->Qi_Routed_m.end(), totalTs, 0.0);
	}

	//Define pervious, impervious, and baseflow runoff terms using unrouted runoff depths from FLowSumCalc.cpp
	double PAQ_m = input->PAQi_m[timeStep];
	double DCIAQ_m = input->DCIAQi_m[timeStep];
	double SubQ_m = input->SubQi_m[timeStep];

	if (routing_flag == 1)
	{

		double PAQ_RT_A_h = input->RefParams["PAQ_RT_A_h"];		// Time parameter of the two parameter surface flow diffusion hydrograph model
		double PAQ_RT_B_h = input->RefParams["PAQ_RT_B_h"];		// Time parameter of the two parameter surface flow diffusion hydrograph model
		double DCIAQ_RT_A_h = input->RefParams["DCIAQ_RT_A_h"];
		double DCIAQ_RT_B_h = input->RefParams["DCIAQ_RT_B_h"];
		double SSQ_RT_h = input->RefParams["SSQ_RT_h"];

		double Time_PA_Max_h = (-0.75 * PAQ_RT_A_h) + (0.75 * PAQ_RT_A_h) * sqrt(1 + (16.0 / 9.0) * (PAQ_RT_B_h / PAQ_RT_A_h));
		double Time_DCIA_Max_h = (-0.75 * DCIAQ_RT_A_h) + (0.75 * DCIAQ_RT_A_h) * sqrt(1 + (16.0 / 9.0) * (DCIAQ_RT_B_h / DCIAQ_RT_A_h));
		double Time_Sub_Max_h = (2.0 / 3.0) * SSQ_RT_h;

		// The routine implements a for loop to delay runoff from the current timestep, (int j = 0; j < tupper; j++),
		//		in effect distributing (delaying) some of the current runoff to future timesteps 
		//		tupper is set by 10*(b * 7.435 + 4.73).  When b is 1.3, tupper becomes 143, instead of 143.9, due to being static cast as an integer
		//		tupper is the time it takes for less than X% of Qmax (due to the precip pulse) to be left
		//		The tupper coefficients were derived by Yang Yang, examining the relationship of beta and the time to X% of Qmax

		// different max distirbution times - separated for surface flows and subsurface flows
		int Time_PA_upper_h = static_cast<int>(10 * ((Time_PA_Max_h * 7.435) + 4.73));
		int Time_DCIA_upper_h = static_cast<int>(10 * ((Time_DCIA_Max_h * 7.435) + 4.73));
		int Time_Sub_upper_h = static_cast<int>(10 * ((Time_Sub_Max_h * 7.435) + 4.73));

		// for water balance tracking - Tom T
		double Delay_SubQ_TS_sum = 0.0;
		double Delay_PAQ_TS_sum = 0.0;
		double Delay_DCIAQ_TS_sum = 0.0;
		double PAQ_it = 0.1;
		double SubQ_it = 0.1;
		double DCIAQ_it = 0.1;

		//PAQ_m (pervious area runoff) computation with Yang and Endreny (2013)
		if (PAQ_m > 0.0)
		{
			// for each spot in the vector
			for (int j = 0; j < Time_PA_upper_h && IT + j < totalTs; ++j)
			{
				IIN = IT + j;    // IIN  = timestep + j

				double first = PAQ_m / sqrt(3.14) * sqrt(PAQ_RT_B_h) * pow(PAQ_it, -1.5) * exp(-PAQ_RT_B_h / PAQ_it) * exp(-PAQ_it / PAQ_RT_A_h) * exp(2 * sqrt(PAQ_RT_B_h / PAQ_RT_A_h));
				if (j == 0)
					PAQ_it = 0.0;
				PAQ_it += 1;
				double second = PAQ_m / sqrt(3.14) * sqrt(PAQ_RT_B_h) * pow(PAQ_it, -1.5) * exp(-PAQ_RT_B_h / PAQ_it) * exp(-PAQ_it / PAQ_RT_A_h) * exp(2 * sqrt(PAQ_RT_B_h / PAQ_RT_A_h));
				PAQ_it += 1;
				double third = PAQ_m / sqrt(3.14) * sqrt(PAQ_RT_B_h) * pow(PAQ_it, -1.5) * exp(-PAQ_RT_B_h / PAQ_it) * exp(-PAQ_it / PAQ_RT_A_h) * exp(2 * sqrt(PAQ_RT_B_h / PAQ_RT_A_h));
				PAQ_it += 1;

				double Delay_PAQ_TS_sum = first + second + third; //std::accumulate(PAQ_it, PAQ_it + stride, 0.0);
				// Pervious Flow (at timestep IIN) is modified to be equal to the Pervious flow (at IIN) + delayPAQ
				input->PAQi_Routed_m[IIN] += Delay_PAQ_TS_sum;
			}
		}

		//DCIA (directly connected impervious area runoff) computation with Yang and Endreny (2013)
		if (DCIAQ_m > 0.0)
		{
			// for each spot in the vector
			for (int j = 0; j < Time_DCIA_upper_h && IT + j < totalTs; ++j)//
			{
				IIN = IT + j;    // IIN  = timestep + j

				//auto sum = std::accumulate(b, b + stride, 0);
				double first = DCIAQ_m / sqrt(3.14) * sqrt(DCIAQ_RT_B_h) * pow(DCIAQ_it, -1.5) * exp(-DCIAQ_RT_B_h / DCIAQ_it) * exp(-DCIAQ_it / DCIAQ_RT_A_h) * exp(2 * sqrt(DCIAQ_RT_B_h / DCIAQ_RT_A_h));
				if (j == 0)
					DCIAQ_it = 0;
				DCIAQ_it += 1;
				double second = DCIAQ_m / sqrt(3.14) * sqrt(DCIAQ_RT_B_h) * pow(DCIAQ_it, -1.5) * exp(-DCIAQ_RT_B_h / DCIAQ_it) * exp(-DCIAQ_it / DCIAQ_RT_A_h) * exp(2 * sqrt(DCIAQ_RT_B_h / DCIAQ_RT_A_h));
				DCIAQ_it += 1;
				double third = DCIAQ_m / sqrt(3.14) * sqrt(DCIAQ_RT_B_h) * pow(DCIAQ_it, -1.5) * exp(-DCIAQ_RT_B_h / DCIAQ_it) * exp(-DCIAQ_it / DCIAQ_RT_A_h) * exp(2 * sqrt(DCIAQ_RT_B_h / DCIAQ_RT_A_h));
				DCIAQ_it += 1;

				//first second and third defined as RHS of Delay_DCIAQ_TS_sum = *DCIAQ_it + *(DCIAQ_it + 1) + *(DCIAQ_it + 2); 
				double Delay_DCIAQ_TS_sum = first + second + third;

				// Impervious Flow (at timestep IIN) is modified to be equal to the Impervious flow (at IIN) + delayPAQ
				input->DCIAQi_Routed_m[IIN] += Delay_DCIAQ_TS_sum;
			}
		}

		//SubQ_m (baseflow runoff) computation with Yang and Endreny (2013)
		if (SubQ_m > 0.0)
		{
			for (int j = 0; j < Time_Sub_upper_h && IT + j < totalTs; ++j) {

				IIN = IT + j;    // IIN  = timestep + j==

				//if (IIN >= totalTs)  // Parameter24 is the total timesteps
					//break;

				double first = SubQ_m / sqrt(3.14) * sqrt(SSQ_RT_h) * pow(SubQ_it, -1.5) * exp(-SSQ_RT_h / SubQ_it);
				if (j == 0)
					SubQ_it = 0;
				SubQ_it += 1;
				double second = SubQ_m / sqrt(3.14) * sqrt(SSQ_RT_h) * pow(SubQ_it, -1.5) * exp(-SSQ_RT_h / SubQ_it);
				SubQ_it += 1;
				double third = SubQ_m / sqrt(3.14) * sqrt(SSQ_RT_h) * pow(SubQ_it, -1.5) * exp(-SSQ_RT_h / SubQ_it);
				SubQ_it += 1;

				double Delay_SubQ_TS_sum = first + second + third;

				// Base Flow (at timestep IIN) is modified to be equal to the base flow (at IIN) + delayPAQ
				input->SubQi_Routed_m[IIN] += Delay_SubQ_TS_sum;
			}
		}
	}
	//flag = 0, nonrouted flow
	else
	{
		//cte 20200918 developing an alternative routing option to the Yang and Endreny (2013) 2-parameter diffusion hydrograph
		IIN = timeStep;
		//Manning discharge terms
		//PAQ_m is depth (m), W is watershed width (m), n is Manning coefficient, S is surface slope (m/m)
		double W_paq = 484;
		double n_paq = 0.24;
		double S_paq = 0.02;
		double W_dcia = 484;
		double n_dcia = 0.015;
		double S_dcia = 0.02;

		// Manning discharge from Eq 3-3 from SWMM Reference Manual Volume I – Hydrology (Revised) 2016 by L. Rossman & W. Huber | EPA/600/R-15/162A | January 2016 | www2.epa.gov/water-research
		// Q = 1/n * Width * Depth_runoff^(5/3) * Slope^(1/2), where [Width * Depth_runoff] = Area_xs, and [Width * Depth_runoff] * Depth_runoff^(2/3) = Width * Depth_runoff^(5/3) 
		// Need to implement algorithm for creating hydrograph, constrained by available volume per time step, and updating depth
		//input->PAQi_Routed_m[IIN] += (1/n_paq * W_paq * pow(PAQ_m,(5/3)) * pow(S_paq, (1 / 2)));
		//input->DCIAQi_Routed_m[IIN] += (1 / n_dcia * W_dcia * pow(DCIAQ_m, (5 / 3)) * pow(S_dcia, (1 / 2)));
		//input->SubQi_Routed_m[IIN] += SubQ_m;

		//Define catchment runoff depth (m) vectors as time distributed vector values
		//Note: These four input-> vector double variables were created in Inputs.h
		//Pervious area runoff, including runon from impervious area (m)
		input->PAQi_Routed_m[IIN] += PAQ_m;
		//Directly connected impervious area runoff (m)
		input->DCIAQi_Routed_m[IIN] += DCIAQ_m;
		//Subsurface runoff (m)
		input->SubQi_Routed_m[IIN] += SubQ_m;
		//Total runoff (m)
	}

	double Qcheck_m = 0;
	//Remove any anomalously low pervious surface runoff generated by Yang routing algorithm 
	if (input->PAQi_Routed_m[timeStep] > 0 && input->PAQi_Routed_m[timeStep] < 0.0000001) {
		input->PAQi_Routed_m[timeStep] = 0;
	}
	//Set a taper of runoff by dividing it in half
	if (input->PAQi_Routed_m[timeStep] >= 0.0000001 && input->PAQi_Routed_m[timeStep] < 0.00001) {
		Qcheck_m = input->PAQi_Routed_m[timeStep - 1] / 2;
		if (Qcheck_m < 0.0000001) {
			input->PAQi_Routed_m[timeStep] = 0;
		}
		else {
			input->PAQi_Routed_m[timeStep] = Qcheck_m;
		}
	}
	//Remove any anomalously low impervious surface runoff generated by Yang routing algorithm 
	if (input->DCIAQi_Routed_m[timeStep] > 0 && input->DCIAQi_Routed_m[timeStep] < 0.0000001) {
		input->DCIAQi_Routed_m[timeStep] = 0;
	}
	//Set a taper of runoff by dividing it in half
	if (input->DCIAQi_Routed_m[timeStep] > 0.0000001 && input->DCIAQi_Routed_m[timeStep] < 0.00001) {
		Qcheck_m = input->DCIAQi_Routed_m[timeStep - 1] / 2;
		if (Qcheck_m < 0.0000001) {
			input->DCIAQi_Routed_m[timeStep] = 0;
		}
		else {
			input->DCIAQi_Routed_m[timeStep] = Qcheck_m;
		}
	}

	//If timeStep greater than first hour, total runoff is sum of all runoff components
	if (timeStep > 0) {
		//Total runoff is sum of subsurface, impervious surface, pervious surface runoff
		input->Qi_Routed_m[timeStep] = input->SubQi_Routed_m[timeStep] + input->DCIAQi_Routed_m[timeStep] + input->PAQi_Routed_m[timeStep];
	}
	//Else timeStep equals 0 then first timeStep of runoff includes the initial starting baseflow value
	else {
		//Total runoff is sum of initial subsurface runoff, impervious surface, pervious surface runoff
		//Note: Q0_mph variable is from HydroPlusConfig.xml, and with To_m2ph is used to derive AveSMD
		input->Qi_Routed_m[timeStep] = input->RefParams["Q0_mph"] * ratio_hr_to_s * input->SimulationNumericalParams["TimeStep_sec"] +
			input->DCIAQi_Routed_m[timeStep] + input->PAQi_Routed_m[timeStep];

	}

	//Store catchment runoff depth (m) as routed (e.g.,time distributed) vector values
	//total_Runoff_Impervious_R_m (m) DCIAQi_Routed_m
	input->RepoDict["total_Runoff_Impervious_R_m"] += input->DCIAQi_Routed_m[timeStep];
	//total_Runoff_Pervious_R_m (m) PAQi_Routed_m
	input->RepoDict["total_Runoff_Pervious_R_m"] += input->PAQi_Routed_m[timeStep];
	//total_Runoff_Subsurface_R_m (m) SubQi_Routed_m
	input->RepoDict["total_Runoff_Subsurface_R_m"] += input->SubQi_Routed_m[timeStep];
	//total_Runoff_R_m (m) Qi_Routed_m
	input->RepoDict["total_Runoff_R_m"] += input->Qi_Routed_m[timeStep];
}
///////////////////////////////////////////////////////////////////////////////////////////////


/*
////////////////////////////////////////////////////////////////////////////////////////////////
// This is the orinal routing code
// Notes on alpha and beta, from Yang and Endreny (2013):
// When alpha approaches infinity (the kinematic flow celerity is 0), the exponential term goes to 1,
// and the 2-parameter surface flow diffusion hydrograph model (developed by Yang and Endreny) has the same form as
// the 1-parameter subsurface flow diffusion hydrograph model of Criss and Winston [2008b].
// The 2-parameter surface flow diffusion hydrograph model is used for PerviousFlow (PAQ_m) and Impervious Flow (DCIAQ_m).
// The 1-parameter subsurface flow diffusion hydrograph model is used for Baseflow (SubQ_m).
// The time input, j, is used in calculating the discharge from the precip (pulse input) at time j in the future

//This for loop delays runoff from the current timestep, in effect distributing (delaying) some of the current runoff to future timesteps
for (int j=0;j<tupper;j++)
// tupper is set by 10*(b * 7.435 + 4.73).  When b is 1.3, tupper becomes 143, instead of 143.9, due to being static cast as an integer
// tupper is the time it takes for less than X% of Qmax (due to the precip pulse) to be left
// The coefficients in the tupper equation were derived by Yang Yang, examining the relationship of beta and the time to X% of Qmax
{

if(j==0)
{
// Pervious flow uses the two parameter surface flow diffusion hydrograph model
delayPAQ= PAQ_m/sqrt(3.14)*sqrt(beta)*pow(0.1,-1.5)*exp(-beta/0.1)*exp(-0.1/alpha)*exp(2*sqrt(beta/alpha));
// Impervious flow (impervious flow) uses the two parameter surface flow diffusion hydrograph model
delayIAQ= DCIAQ_m/sqrt(3.14)*sqrt(beta)*pow(0.1,-1.5)*exp(-beta/0.1)*exp(-0.1/alpha)*exp(2*sqrt(beta/alpha));
// Base flow uses the one parameter subsurface flow diffusion hydrograph model
delayBaseQ=SubQ_m/sqrt(3.14)*sqrt(b)*pow(0.1,-1.5)*exp(-b/0.1);
	}
else
{
// Pervious flow uses the two parameter surface flow diffusion hydrograph model
delayPAQ= PAQ_m/sqrt(3.14)*sqrt(beta)*pow(j,-1.5)*exp(-beta/j)*exp(-j/alpha)*exp(2*sqrt(beta/alpha));
// Impervious flow (impervious flow) uses the two parameter surface flow diffusion hydrograph model
delayIAQ= DCIAQ_m/sqrt(3.14)*sqrt(beta)*pow(j,-1.5)*exp(-beta/j)*exp(-j/alpha)*exp(2*sqrt(beta/alpha));
// Base flow uses the one parameter subsurface flow diffusion hydrograph model
delayBaseQ=SubQ_m/sqrt(3.14)*sqrt(b)*pow(j,-1.5)*exp(-b/j);
}

IIN = IT + j;    // IIN  = timestep + j
if( IIN >= totalTs)  // Parameter24 is the total timesteps
break;
// If the current timestep + j (the hours into the future we are distributing flow) > the total timesteps, break
// Meaning, there is no point to distributing flow to the hours beyond the model simulation
// Pervious Flow (at timestep IIN) is modified to be equal to the Pervious flow (at IIN) + delayPAQ
repo->runoff.modifyinput->PAQi_Routed_m(repo->runoff.getinput->PAQi_Routed_m(IIN) + delayPAQ, IIN);
// Impervious Flow (at timestep IIN) is modified to be equal to the Impervious flow (at IIN) + delayPAQ
repo->runoff.modifyIAQi(repo->runoff.getIAQi(IIN) + delayIAQ, IIN);
// Base Flow (at timestep IIN) is modified to be equal to the base flow (at IIN) + delayPAQ
repo->runoff.modifyinput->SubQi_Routed_m(repo->runoff.getinput->SubQi_Routed_m(IIN) + delayBaseQ, IIN);

//Total Flow (at timestep IIN) is modified to be equal to the all the flow components (at IIN)
//Old code:
//repo->runoff.modifyinput->Qi_Routed_m(repo->runoff.getinput->Qi_Routed_m(IIN) + (delayPAQ + delayIAQ + delayBaseQ), IIN);
//repo->runoff.setTotalOLQ(repo->runoff.getTotalOLQ() + delayPAQ);
//repo->runoff.setRunTotalBQ(repo->runoff.getRunTotalBQ() + delayBaseQ);
//repo->runoff.setTotalUrImpQ(repo->runoff.getTotalUrImpQ() + delayIAQ);
}
repo->runoff.setTotalQQ(repo->runoff.getRunTotalQQ() + repo->runoff.getinput->Qi_Routed_m(timeStep));
*/
