#include "Drainage_StormwaterDevice.h"

//Percolation_Pavement function has two key roles, parameterizing hydraulic conductivity for infiltration and percolation
void Drainage_StormwaterDevice::Percolation_Pavement(Inputs* input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	double PermeabilityReduction_frac = 0.0;
	double Depth_PavementClogged_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Drainage_Limit_m");
	double Time_PavementRegenerationFrequency_hr = beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Drainage_Renewed_d") * Ratio_Hour_to_Day;
	double Total_Drainage_GI_Pavement_m = 0.0;
	double Pavement_HydraulicConductivity_Decayed_mph = 0.0;

	//If Time_PeriodSincePavementRegenerated_hr > Time_PavementRegenerationFrequency_hr then time to regenerate
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Time_PeriodSincePavementRegenerated_hr") > Time_PavementRegenerationFrequency_hr) {

		//Time_PeriodSincePavementRegenerated_hr reset to zero
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Time_PeriodSincePavementRegenerated_hr") = 0;
		//Total_Drainage_GI_Pavement_m3 (m3) reset to zero
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Total_Drainage_GI_Pavement_m3") = 0;
	}

	//Total_Drainage_GI_Pavement_m (m) is Total_Drainage_GI_Pavement_m3 (m3) divided by Area_GI_Unit_Pervious_m2 (m2)
	Total_Drainage_GI_Pavement_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Total_Drainage_GI_Pavement_m3") / beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2");

	//PermeabilityReduction_frac (fraction) is Total_Drainage_GI_Pavement_m divided by Depth_PavementClogged_m
	//Note: Depth_PavementClogged_m is confirmed > 0 in BuildDataOrganizer::Check_HydroPlusConfig_Inputs function
	PermeabilityReduction_frac = Total_Drainage_GI_Pavement_m / Depth_PavementClogged_m;
	//PermeabilityReduction_frac is minimum of PermeabilityReduction_frac and 1
	PermeabilityReduction_frac = MIN(PermeabilityReduction_frac, 1.0);
	//PermeabilityReduction_frac is maximum of PermeabilityReduction_frac and 0
	PermeabilityReduction_frac = MAX(PermeabilityReduction_frac, 0);

	//Pavement_HydraulicConductivity_Decayed_mph is Pavement_HydraulicConductivity_mph reduced by (1-PermeabilityReduction_frac)
	Pavement_HydraulicConductivity_Decayed_mph = beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_HydraulicConductivity_mph") * (1 - PermeabilityReduction_frac);
	//Pavement_HydraulicConductivity_Decayed_mph sent to infiltration functor for GreenAmpt infiltration 
	//Note: InfiltExcessGovernedArea_frac set to 1 in BuildDataOrganizer::Check_HydroPlusConfig_Inputs function to ensure GreenAmpt infiltration
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_HydraulicConductivity_Decayed_mph") = Pavement_HydraulicConductivity_Decayed_mph;

	//Pavement_HydraulicConductivity_Decayed_mpdt is Pavement_HydraulicConductivity_Decayed_mph converted by Ratio_Hour_to_Second * ratio_s_to_timeStep
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_HydraulicConductivity_Decayed_mpdt") = Pavement_HydraulicConductivity_Decayed_mph * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep];

	//If Type equals to PermeablePavement then 
	if (beC->by_key_str(DataDrawer_ID, DataFolder_ID, "Type") == "PermeablePavement") {
		//Percolation_Pavement_m3 (m3) for each time step is product of Pavement_HydraulicConductivity_Decayed_mpdt (m/dt) and Area_GI_Unit_Pervious_m2 (m2)
		//Note: Vault is presumed to have area of Area_GI_Unit_m2 for exfiltration, with water otherwise infiltrating, percolating, and evaporating through Area_GI_Unit_Pervious_m2
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Percolation_Pavement_m3") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_HydraulicConductivity_Decayed_mpdt") *
			beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2");
	}
	//If Type not equal to PermeablePavement then 
	else {
		//Percolation_Pavement_m3 (m3) equals 0
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Percolation_Pavement_m3") = 0.0;
	}

	//Pavement_Drainage_Limit_m (CF): Defines stormwater depth treated before completely clogged
	//CF = Yclog * Pa * (1 + CR) * (1 + VR)  / (VR * (1 - ISF) * T)
	//Note: where Yclog (yr) to completely (100%) clog the pavement, Pa (in) is the annual rainfall amount over the site, 
	//Note: CR is the pavement's capture ratio (area that contributes runoff to the pavement divided by area of the pavement itself), 
	//Note: VR is the system's Void Ratio, ISF is the Impervious Surface Fraction of pavement layer, and 
	//Note: T (in) is the pavement layer Thickness, FClog is fraction clogged. Note phi = e/(e+1) and e=phi/(1-phi)
	//Note: An example, it takes 5 years to completely clog pavement with an annual rainfall is 36 inches/year, a pavement 6 inches thick, 
	//Note: ... a void ratio of 0.2, impervious cover of zero, and no drainage inflow only treating direct rainfall so that CR = 0, then 
	//Note: ... the CF = 5 x 36 x 1 x (1 + 0.2) / 0.2 / 1 / 6 / 1 = 180 in
}

//Percolation_Soil function receives Storage_GI_Soil_m3 scaled upwards by prior division by Porosity_ThetaSat_ThetaWp_m3pm3
double Drainage_StormwaterDevice::Percolation_Soil(Inputs* input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	//Deficit_Soil_GravitationalWater_m3 (m3/m3) is soil moisture deficit, or capacity to hold more water
	double Deficit_Soil_GravitationalWater_m3 = 0.0;
	double Deficit_Soil_GravitationalWater_m3pm3 = 0.0;
	double Percolation_Soil_m = 0.0;
	double Percolation_Soil_m3 = 0.0;
	double Percolation_Potential_m3 = 0.0;
	double Soil_HydraulicConductivity_m_p_dt = 0.0;

	//If Storage_GI_Soil_m3 (m3) > Storage_Soil_FieldCapacity_m3 (m3) there is Percolation_Soil_m3 (m3) via gravitational drainage
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Soil_m3") > beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Soil_FieldCapacity_m3")) {

		//Deficit_Soil_GravitationalWater_m3 (m3) is Storage_GI_Soil_Max_m3 - Storage_GI_Soil_m3, given Storage_GI_Soil_m3 > Storage_Soil_FieldCapacity_m3(m3)
		Deficit_Soil_GravitationalWater_m3 = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Soil_Max_m3") - beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Soil_m3");

		//Deficit_Soil_GravitationalWater_m (m) is scaled to Deficit_Soil_GravitationalWater_m3pm3 (m3/m3) by division with Volume_GI_SoilLayer_m3 (m3)
		//Note: Deficit_Soil_GravitationalWater_m3pm3 (m3/m3) is volumetric moisture, term used in equation for Percolation_Soil_m
		Deficit_Soil_GravitationalWater_m3pm3 = Deficit_Soil_GravitationalWater_m3 / beC->by_key(DataDrawer_ID, DataFolder_ID, "Volume_GI_SoilLayer_m3");

		//Soil_HydraulicConductivity_m_p_dt (m/dt) is K0_mph (m/hr) adjusted to per time step based on Ratio_Hour_to_Second and SimulationTimeStep_Duration_sec[timeStep]
		Soil_HydraulicConductivity_m_p_dt = beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Ksat_mph") * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep];

		//Percolation_Soil_m (m) is function of with Soil_HydraulicConductivity_m_p_dt, Deficit_Soil_GravitationalWater_m3, Soil_PercolationDecay_Coeff
		//Note: Percolation equation from EPA SWMM code lidproc.c function getSoilPercRate, with behavior described next ..
		//Note: The exponential has negative coefficient to ensure its result ranges as a fraction from 1 to 0, ensuring it will not be greater than 1
		//Note: If Deficit_Soil_GravitationalWater_m3 or Soil_PercolationDecay_Coeff = 0, the exponential result becomes 1 and Percolation_Soil_m3 = Soil_HydraulicConductivity_m_p_dt
		Percolation_Soil_m = Soil_HydraulicConductivity_m_p_dt * exp(-Deficit_Soil_GravitationalWater_m3pm3 * 
			beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_PercolationDecay_Coeff"));
	
		//Percolation_Soil_m3 (m3) is Percolation_Soil_m (m) over pervious cover scaled to folder are by Area_per_GIUnit_m2 and PerviousCover_frac
		//Note: Vault is presumed to have area of Area_GI_Unit_m2 for exfiltration, with water otherwise infiltrating, percolating, and evaporating through Area_GI_Unit_Pervious_m2
		Percolation_Soil_m3 = Percolation_Soil_m * beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2");
		
		//Percolation_Potential_m3 (m3) is available water, Storage_GI_Soil_m3 (m3) minus EvapoTrans_Soil_m3 (m3) that was removed
		Percolation_Potential_m3 = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Soil_m3") - beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Soil_m3");

		//If Percolation_Potential_m3 (m3) > 0 then use it in MIN function
		if (Percolation_Potential_m3 > 0) {
			//Percolation_Soil_m3 (m3) is minimum of Percolation_Soil_m3 (m3) and Percolation_Potential_m3 (m)
			Percolation_Soil_m3 = MIN(Percolation_Soil_m3, Percolation_Potential_m3);
		}
		//Else Percolation_Potential_m3 (m3) < zero and Percolation_Soil_m3 = 0
		else { Percolation_Soil_m3 = 0.0; }
	}
	//Else If Type equal to InfilTrench or RainBarrel then 
	else if (beC->by_key_str(DataDrawer_ID, DataFolder_ID, "Type") == "InfilTrench" || beC->by_key_str(DataDrawer_ID, DataFolder_ID, "Type") == "RainBarrel") {
		Percolation_Soil_m3 = 0.0;
	}
	//Else Storage_GI_Soil_m3(m) <= Storage_Soil_FieldCapacity_m3(m) then Percolation_Soil_m3 (m3) is zero
	else { Percolation_Soil_m3 = 0.0; }

	//Percolation_Soil_m3 (m3) returned as volume for folder area
	return Percolation_Soil_m3;
}

//Exfiltration_Vault function computes potential Exfiltration_m3 (m3) value of drainage from vault zone
double Drainage_StormwaterDevice::Exfiltration_Vault(Inputs* input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	double Exfiltration_Max_m3 = Constant_1E_positive10;
	double Exfiltration_m = 0.0;
	double Exfiltration_m3 = 0.0;
	double Exfiltration_Rate_mpdt = 0.0;
	double Total_GI_Drainage_VadoseZone_m = 0.0;
	double Exfiltration_Reduction_frac = 0.0;
	double Exfiltration_Limit_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Exfiltration_Limit_m");
	//HydraulicConductivity_Catchment_mph (m/h) is defined with HydroPlusConfig.xml element Soil_Ksat_mph (m/h) for BulkArea reference folder
	double HydraulicConductivity_Catchment_mph = input->InputXml["Soil_Ksat_mph"];
	//ExfiltrationRate_Vault_mpdt (m/h) is defined with HydroPlusConfig.xml element Vault_ExfiltrationRate_mph (m/h) for GI folder
	double ExfiltrationRate_Vault_mph = beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_ExfiltrationRate_mph");

	//Exfiltration_Rate_mpdt (m/dt) is product of ExfiltrationRate_Vault_mph (m/h) and Ratio_Hour_to_Second and SimulationTimeStep_Duration_sec[timeStep]
	//Note: Default is to define Exfiltration_Rate_mpdt using ExfiltrationRate_Vault_mph, and if this is zero, then use HydraulicConductivity_Catchment_mph
	Exfiltration_Rate_mpdt = ExfiltrationRate_Vault_mph * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep];

	//If ExfiltrationRate_Vault_mph equals zero, and Flag_Exfiltration_to_Catchment = 1, then use HydraulicConductivity_Catchment_mph for exfiltration
	if (Exfiltration_Rate_mpdt <= 0) {
		//Exfiltration_Rate_mpdt (m/dt) is product of HydraulicConductivity_Catchment_mph (m/h) and Ratio_Hour_to_Second and SimulationTimeStep_Duration_sec[timeStep]
		Exfiltration_Rate_mpdt = HydraulicConductivity_Catchment_mph * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep];
	}


	//If Exfiltration_Rate_mpdt <= zero, then Exfiltration_m3 (m3) equals zero
	if (Exfiltration_Rate_mpdt <= 0.0) {
		Exfiltration_m3 = 0.0;
		return Exfiltration_m3;
	}

	//Total_GI_Drainage_VadoseZone_m (m) is Total_GI_Drainage_VadoseZone_m3 (m3) divided by Area_GI_Unit_Pervious_m2 (m2)
	//Note: Vault is presumed to have area of Area_GI_Unit_m2 for exfiltration, with water otherwise infiltrating, percolating, and evaporating through Area_GI_Unit_Pervious_m2
	Total_GI_Drainage_VadoseZone_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Total_GI_Drainage_VadoseZone_m3") / beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_m2");

	//Exfiltration_Reduction_frac  (fraction) is Total_GI_Drainage_VadoseZone_m divided by Exfiltration_Limit_m
	//Note: Exfiltration_Limit_m is confirmed > 0 in BuildDataOrganizer::Check_HydroPlusConfig_Inputs function
	Exfiltration_Reduction_frac = Total_GI_Drainage_VadoseZone_m / Exfiltration_Limit_m;
	//Exfiltration_Reduction_frac is minimum of Exfiltration_Reduction_frac and 1
	Exfiltration_Reduction_frac = MIN(Exfiltration_Reduction_frac, 1.0);
	//Exfiltration_Reduction_frac is maximum of Exfiltration_Reduction_frac and 0
	Exfiltration_Reduction_frac = MAX(Exfiltration_Reduction_frac, 0);
	
	//Exfiltration_m (m) is product of (1-Exfiltration_Reduction_frac), Exfiltration_Rate_mpdt (m/dt)
	Exfiltration_m = (1.0 - Exfiltration_Reduction_frac) * Exfiltration_Rate_mpdt;

	//Exfiltration_m3 (m3) is Exfiltration_m (m) scaled to folder volume by product with Area_per_GIUnit_m2 and PerviousCover_frac
	//Note: Vault is presumed to have area of Area_GI_Unit_m2 for exfiltration, with water otherwise infiltrating, percolating, and evaporating through Area_GI_Unit_Pervious_m2
	Exfiltration_m3 = Exfiltration_m * beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_m2");

	//Exfiltration_m3 (m3) is MAX of Exfiltration_m3 and zero
	Exfiltration_m3 = MAX(Exfiltration_m3, 0);
	
	//If Flag_Exfiltration_to_Catchment not equal to 1, then GI Type has zero exfiltration to catchment
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Flag_Exfiltration_to_Catchment")!= 1) {
		//Exfiltration_m3 (m3) equals 0
		Exfiltration_m3 = 0.0;
	}

	//return Exfiltration_m3 (m3) for removal from GI and transfer to BulkArea
	return Exfiltration_m3;
}

//Drainage_VadoseZone function determines if percolation replaces exfiltration as drainage to the native 	
//Note: Groundwater_Mounding function modifies Drainage_VadoseZone_m3 with changes in saturated mounding beneath the GI device 
void Drainage_StormwaterDevice::Drainage_VadoseZone(Inputs* input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	//If Vault_Thickness_m > Constant_1E_negative10 And Flag_Exfiltration_to_Catchment equals 1 or True then 
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Thickness_m") > Constant_1E_negative10 && beC->by_key(DataDrawer_ID, DataFolder_ID, "Flag_Exfiltration_to_Catchment")== 1) {
		//Drainage_VadoseZone_m3 (m3) is any Exfiltration_m3 (m3)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Drainage_VadoseZone_m3") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Exfiltration_m3");
	}
	//Else If Vault_Thickness_m <= Constant_1E_negative10 And Soil_Thickness_m > Constant_1E_negative10 And 
	//Note: Flag_Exfiltration_to_Catchment equals 1 or True then 
	else if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Thickness_m") <= Constant_1E_negative10 && beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m") > Constant_1E_negative10 &&
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Flag_Exfiltration_to_Catchment")== 1) {
		//Drainage_VadoseZone_m3 (m3) is any Percolation_Soil_m3 (m3)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Drainage_VadoseZone_m3") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Percolation_Soil_m3");
	}
	//Else If Vault_Thickness_m And Soil_Thickness_m <= Constant_1E_negative10 And Pavement_Thickness_m > Constant_1E_negative10 And 
	//Note: Flag_Exfiltration_to_Catchment equals 1 or True then 
	else if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Thickness_m") <= Constant_1E_negative10 && beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m") <= Constant_1E_negative10 &&
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") > Constant_1E_negative10 && beC->by_key(DataDrawer_ID, DataFolder_ID, "Flag_Exfiltration_to_Catchment")== 1) {
		//Drainage_VadoseZone_m3 (m3) is any Percolation_Pavement_m3 (m3)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Drainage_VadoseZone_m3") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Percolation_Pavement_m3");
	}
	//Else Flag_Exfiltration_to_Catchment equals 0 or False 
	else {
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Drainage_VadoseZone_m3 = 0.0");
	}
}