#include "RunoffSummation.h"

#include <math.h>
#include <fstream>

//Note: Runoff_SatExcess_m and Runoff_InfilExcess_m remain at cover type depth
//Note: GI runoff terms arrive as volumes for entire GI area, and are converted to depths for folder area
//Note: Runoff_Subsurface_m (m) is calculated for all folders with pervious area

//RunoffSummation::SumAcrossCoverTypes function sums the runoff leaving the folder, combining cover types to get folder depths
void RunoffSummation::SumAcrossCoverTypes(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	double Runoff_Surface_m = 0.0;
	double Runoff_Subsurface_m = 0.0;
	double Runoff_m = 0.0;
	double Runoff_Vault_m = 0.0;

	double Runoff_SatExcess_m = 0.0;
	double Runoff_InfilExcess_m = 0.0;
	double Runoff_Pervious_m = 0.0;
	double Runoff_Water_m = 0.0;
	double Runoff_Impervious_m = 0.0;

	//Note: Multiplication of variables by fractional area below is to adjust depths from one analysis area to entire folder area
	//Rationale: 1. variables arrive as depths for an analysis area (e.g., pervious area)
	//Rationale: 2. multiplying these depths by their analysis area gives a volume, which is conserved when transfered between different areas
	//Rationale: 3. multiplying these depths by their fractional area is equivalent to multiplying by their analysis area and dividing by the folder area, getting depth for folder area
	//Note: Multiplication is equivalent to multiplying by ratio of analysis area to folder area to convert from /analysis_area to /folder_area

	//ratio_FolderArea_to_CatchmentArea is ratio of folder area to total catchment area
	double ratio_FolderArea_to_CatchmentArea = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2") / input->SimulationNumericalParams["CatchmentArea_m2"];
	//PerviousCover_frac is pervious area fraction
	double PerviousCover_frac = beC->by_key(MapPixel_ID, DataFolder_ID, "PerviousCover_frac");
	//ImperviousCover_frac is impervious area fraction
	double ImperviousCover_frac = beC->by_key(MapPixel_ID, DataFolder_ID, "ImperviousCover_frac");
	//WaterCover_noTreeCanopy_frac is impervious area fraction
	double WaterCover_noTreeCanopy_frac = beC->by_key(MapPixel_ID, DataFolder_ID, "WaterCover_noTreeCanopy_frac");

	//Runoff_SatExcess_m (m) is folder area depth; Runoff_SatExcess_m (m) is cover area depth, converted to folder area depth by multiplication with PerviousCover_frac
	Runoff_SatExcess_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_m") * PerviousCover_frac;
	//Runoff_InfilExcess_m (m) is folder area depth; Runoff_InfilExcess_m (m) is cover area depth, converted to folder area depth by multiplication with PerviousCover_frac
	Runoff_InfilExcess_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_InfilExcess_m") * PerviousCover_frac;
	//Runoff_Pervious_m (m) as folder area depth is sum of Runoff_SatExcess_m (m) and Runoff_InfilExcess_m (m)
	//Note: Runoff_SatExcess_m (m) and Runoff_InfilExcess_m (m) include runon from impervious area; only runoff mechanisms from pervious area
	Runoff_Pervious_m = Runoff_SatExcess_m + Runoff_InfilExcess_m;
	//Runoff_Water_m (m) is folder area depth; Runoff_Water_to_Outlet_m (m) is cover area depth, converted to folder area depth by multiplication with WaterCover_noTreeCanopy_frac
	Runoff_Water_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Water_to_Outlet_m") * WaterCover_noTreeCanopy_frac;

	//Runoff_Impervious_m (m) is folder area depth; Runoff_Impervious_to_Outlet_m (m) converted from cover area depth by multiplication with ImperviousCover_DrainsToOutlet_frac earlier
	//Note: DepressionImpervious::DepressionImpervious function uses the ImperviousCover_DrainsToOutlet_frac conversion to folder area depth
	Runoff_Impervious_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Impervious_to_Outlet_m");

	//If Flag_GI_Simulated is true and Type is not BulkArea, then convert GI runoff volumes (m3) to folder depths (m) for use in AggregateOutput function
	if (input->RepoDict["Flag_GI_Simulated"] == 1 && beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") != "BulkArea" && beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2") > 0) {

		//If Flag_Outflow_to_NextFolder equals 1 then GI drains to next folder, as GI in series
		//Note: Flag_Outflow_to_NextFolder defined in HydroPlusConfig.xml file for GI device that will drain to next GI folder; undefined = 0
		if (beC->has_var("Flag_Outflow_to_NextFolder") && beC->by_key(MapPixel_ID, DataFolder_ID, "Flag_Outflow_to_NextFolder") == 1) {

			//Runoff_Impervious_GItoGI_m3 (m3) is defined as Runoff_Impervious_m3 (m3) coming from upslope GI folder devices
			//Note: Runoff_Impervious_m3 (m3) is volume of water across all GI devices, created in Flux_Manager_StormwaterDevice
			input->RepoDict["Runoff_Impervious_GItoGI_m3"] = input->RepoDict["Runoff_Impervious_m3"];
			//Runoff_Pervious_GItoGI_m3 (m) is defined as Runoff_Pervious_m3 (m3) coming from upslope GI folder devices
			//Note: Runoff_Pervious_m3 (m3) is volume of water across all GI devices, created in Flux_Manager_StormwaterDevice
			input->RepoDict["Runoff_Pervious_GItoGI_m3"] = input->RepoDict["Runoff_Pervious_m3"] / beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2");
			//Runoff_Vault_GItoGI_m3 (m3) is defined as Runoff_Vault_m3 (m3) coming from upslope GI folder devices 
			//Note: Runoff_Vault_m3 (m3) is volume of water across all GI devices, created in Flux_Manager_StormwaterDevice
			input->RepoDict["Runoff_Vault_GItoGI_m3"] = input->RepoDict["Runoff_Vault_m3"];
			//Runoff_Water_GItoGI_m3 (m3) is defined as Runoff_Water_m3 (m3) coming from upslope GI folder devices 
			//Note: Runoff_Water_m3 (m3) is volume of water across all GI devices, created in Flux_Manager_StormwaterDevice
			input->RepoDict["Runoff_Water_GItoGI_m3"] = input->RepoDict["Runoff_Water_m3"];

			//Set all GI to outlet transfer depths to zero given that Flag_Outflow_to_NextFolder = 1
			Runoff_Impervious_m = 0.0;
			Runoff_Pervious_m = 0.0;
			Runoff_SatExcess_m = 0.0;
			Runoff_InfilExcess_m = 0.0;
			Runoff_Vault_m = 0.0;
			Runoff_Water_m = 0.0;
			beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Impervious_to_Outlet_m") = 0;
			beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Water_to_Outlet_m") = 0;
		}
		//Else If Flag_Outflow_to_NextFolder not 1 then GI drains to catchment outlet
		else {
			//Runoff_Impervious_m (m) is defined as Runoff_Impervious_m3 (m3) divided by GI Area_m2 (m2)
			//Note: Runoff_Impervious_m3 (m3) is volume of water across all GI devices, created in Flux_Manager_StormwaterDevice
			Runoff_Impervious_m = input->RepoDict["Runoff_Impervious_m3"] / beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2");
			//Runoff_Pervious_m (m) is defined as Runoff_Pervious_m3 divided by BulkArea Area_m2
			//Note: Runoff_Pervious_m3 (m3) is volume of water across all GI devices, created in Flux_Manager_StormwaterDevice
			Runoff_Pervious_m = input->RepoDict["Runoff_Pervious_m3"] / beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2");
			//Runoff_InfilExcess_m (m) estimated as portion of Runoff_Pervious_m3 controlled by Ratio_RunoffInfilExcess_to_RunoffPervious, computed in Flux_Manager_StormwaterDevice
			Runoff_InfilExcess_m = (input->RepoDict["Runoff_Pervious_m3"] / beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2")) * beC->by_key(MapPixel_ID, DataFolder_ID, "Ratio_RunoffInfilExcess_to_RunoffPervious");
			beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_InfilExcess_m") = Runoff_InfilExcess_m;

			//Note: Consider refactor to allow for Runoff_SatExcess_m based on TI distribution
			//Runoff_InfilExcess_m (m) estimated as portion of Runoff_Pervious_m3 controlled by (1 - Ratio_RunoffInfilExcess_to_RunoffPervious), computed in Flux_Manager_StormwaterDevice
			Runoff_SatExcess_m = (input->RepoDict["Runoff_Pervious_m3"] / beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2")) * (1 - beC->by_key(MapPixel_ID, DataFolder_ID, "Ratio_RunoffInfilExcess_to_RunoffPervious"));
			beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_m") = Runoff_SatExcess_m;
			//Runoff_Water_m (m) is defined as Runoff_Water_m3 (m3) divided by BulkArea Area_m2 (m2)
			//Note: Runoff_Water_m3 (m3) is volume of water across all GI devices, created in Flux_Manager_StormwaterDevice
			Runoff_Water_m = input->RepoDict["Runoff_Water_m3"] / beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2");
			//Runoff_Vault_m (m) is quotient of Runoff_Vault_m3 (m3) and Area_m2; Runoff_Vault_m3 can contain drainage from any cover type, computed in Flux_Manager_StormwaterDevice
			Runoff_Vault_m = input->RepoDict["Runoff_Vault_m3"] / beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2");
			//Runoff_Impervious_m (m) increased by portion of Runoff_Vault_m from ImperviousCover_frac
			Runoff_Impervious_m = Runoff_Impervious_m + Runoff_Vault_m * beC->by_key(MapPixel_ID, DataFolder_ID, "ImperviousCover_frac");
			//Runoff_Pervious_m (m) increased by portion of Runoff_Vault_m from PerviousCover_frac
			Runoff_Pervious_m = Runoff_Pervious_m + Runoff_Vault_m * beC->by_key(MapPixel_ID, DataFolder_ID, "PerviousCover_frac");
			//Runoff_Water_m (m) increased by portion of Runoff_Vault_m from WaterCover_noTreeCanopy_frac
			Runoff_Water_m = Runoff_Water_m + Runoff_Vault_m * beC->by_key(MapPixel_ID, DataFolder_ID, "WaterCover_noTreeCanopy_frac");
			
			//Update beC->by_key variables with GI components:
			//Runoff_Impervious_to_Outlet_m (m) quotient of Runoff_Impervious_m3 (m3) and Area_m2 (m2)
			//Note: Runoff_Impervious_m3 (m3) is volume of water across all GI devices, created in Flux_Manager_StormwaterDevice
			beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Impervious_to_Outlet_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Impervious_m3") / beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2");
			//Runoff_Water_to_Outlet_m (m) quotient of Runoff_Water_m3 (m3) and Area_m2 (m2)
			//Note: Runoff_Water_m3 (m3) is volume of water across all GI devices, created in Flux_Manager_StormwaterDevice
			beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Water_to_Outlet_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Water_m3") / beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2");
			
			//Set all GI to GI transfer volumes to zero given that Flag_Outflow_to_NextFolder = 0
			input->RepoDict["Runoff_Impervious_GItoGI_m3"] = 0;
			input->RepoDict["Runoff_Pervious_GItoGI_m3"] = 0;
			input->RepoDict["Runoff_Vault_GItoGI_m3"] = 0;
			input->RepoDict["Runoff_Water_GItoGI_m3"] = 0;
		}
	}
	//Runoff_Subsurface_m (m) converted to folder area depth, arriving as pervious area depth, Runoff_Subsurface_m (m), from RunoffSubsurface.cpp
	//Note: For 1st timestep, do not use input->InputXml["Discharge_Subsurface_Initial_mph"] for Runoff_Subsurface_m, given this is not in equilibrium with AveSMD_m
	Runoff_Subsurface_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Subsurface_m") * PerviousCover_frac;

	//Runoff_Surface_m (m) as folder area depth total surface runoff, Runoff_Impervious_m + Runoff_Pervious_m + Runoff_Water_m 
	Runoff_Surface_m = Runoff_Impervious_m + Runoff_Pervious_m + Runoff_Water_m;

	//Runoff_m (m) as folder area depth total runoff, Runoff_Surface_m + Runoff_Subsurface_m
	Runoff_m = Runoff_Surface_m + Runoff_Subsurface_m;

	//Runoff as folder area depth (m) sent to AggregateOutput functions and transformed to catchment area depth
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_m") = Runoff_SatExcess_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_InfilExcess_m") = Runoff_InfilExcess_m;
	//Runoff_Impervious_m (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Impervious_m") = Runoff_Impervious_m;
	//Runoff_Pervious_m (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Pervious_m") = Runoff_Pervious_m;
	//Runoff_Water_m (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Water_m") = Runoff_Water_m;
	//Runoff_Subsurface_m (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Subsurface_m") = Runoff_Subsurface_m;
	//Runoff_Surface_m (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Surface_m") = Runoff_Surface_m;
	//Runoff_m (m) = BQ + PAQ + DCIAQ
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_m") = Runoff_m;
}
