#include "Removal_StormwaterPollutants.h"
#include "../Buffer/BufferSpatialCalc.h"
#include <iomanip>
#include <string>
#include <fstream>

//Author(s) : Li Zhang, modified by Professor Endreny
//lzhan126@syr.edu March 2021, 2022

using namespace std;
map<string, double> variables;

//Removal_StormwaterPollutants::Pollutant_Removal_via_GI function computes pollutant load (g) for use in TimeSeries_Budget_Pollutants_GI.csv file
//Note:	Name_Release_GI_Total_Pollutant_Mean_g can become negative at any time step, which is evidence of pollutant release and is recorded as negative removal
//Note: Pollutant_Name and Pollutant_ConcentrationMean_mg_p_L vectors used below are from Pollutants.csv, read in Inputs.cpp into vector 
//Note: beC->by_key(MapPixel_ID, DataFolder_ID, Name_Removal_Rate_Pollutant_percent) contains removal efficiencies from HydroPlusConfig.xml elements 
//Note: HydroPlusConfig.xml element names and Pollutant.dat variable names must be identical strings; e.g., TSS, BOD, COD, TP, SolP, TKN, NO2_3, Cu, Pb, Zn
void Removal_StormwaterPollutants::Pollutant_Removal_via_GI(Inputs* input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	//create variable names 3 types of flow and 10 types of pollutants 
	for (int pollutant_ID = 0; pollutant_ID < input->Pollutant_Name.size(); pollutant_ID++) {
		//Create string to represent variable associated with Pollutant_Name
		string Name_Load_GI_Inflow_Pollutant_Mean_g = "Load_GI_Inflow_" + input->Pollutant_Name[pollutant_ID] + "_Mean_g";
		string Name_Load_GI_Runoff_Surface_Pollutant_Mean_g = "Load_GI_Runoff_Surface_" + input->Pollutant_Name[pollutant_ID] + "_Mean_g";
		string Name_Load_GI_Runoff_Vault_Pollutant_Mean_g = "Load_GI_Runoff_Vault_" + input->Pollutant_Name[pollutant_ID] + "_Mean_g";
		string Name_Load_GI_Drainage_VadoseZone_Pollutant_Mean_g = "Load_GI_Drainage_VadoseZone_" + input->Pollutant_Name[pollutant_ID] + "_Mean_g";
		string Name_Release_GI_Total_Pollutant_Mean_g = "Release_GI_Total_" + input->Pollutant_Name[pollutant_ID] + "_Mean_g";
		string Name_Removal_Rate_Pollutant_percent = input->Pollutant_Name[pollutant_ID] + "_RemovalEfficiency_percent";

		//Pollutant inflow mass (g), product (GI_Precipitation_m3, GI_Inflow_from_BulkArea_m3, Ratio_L_to_m3, Pollutant_ConcentrationMean_mg_p_L (mg/L) * Ratio_g_to_mg)
		//Note: Actual variable name for beC->by_key is string defined above
		beC->by_key(DataDrawer_ID, DataFolder_ID, Name_Load_GI_Inflow_Pollutant_Mean_g) = ((beC->by_key(DataDrawer_ID, DataFolder_ID, "GI_Precipitation_m3") + beC->by_key(DataDrawer_ID, DataFolder_ID, "GI_Inflow_from_BulkArea_m3")) * Ratio_L_to_m3 * input->Pollutant_ConcentrationMean_mg_p_L[pollutant_ID] * Ratio_g_to_mg);

		//Pollutant outflow mass (g), product (GI_Runoff_Surface_m (m), Ratio_m3_to_L, Pollutant_ConcentrationMean_mg_p_L (mg/L) * Ratio_g_to_mg)
		//Note: Retained pollutant is computed as that not removed, using last term [1 - (Removal_Efficiency_percent * Ratio_Decimal_to_Percent)]
		//Note: Removal_Efficiency_percent from HydroPlusConfig.xml is in the variable beC->by_key(MapPixel_ID, DataFolder_ID, Name_Removal_Rate_Pollutant_percent)
		//Note: Actual variable name for beC->by_key is string defined above
		//Li 12252024 pollutant concentration won't be removed for surface outflow
		beC->by_key(DataDrawer_ID, DataFolder_ID, Name_Load_GI_Runoff_Surface_Pollutant_Mean_g) = beC->by_key(DataDrawer_ID, DataFolder_ID, "GI_Runoff_Surface_m3") * Ratio_L_to_m3 *
			input->Pollutant_ConcentrationMean_mg_p_L[pollutant_ID] * Ratio_g_to_mg * 1;

		//Name_Load_GI_Runoff_Vault_Pollutant_Mean_g (g) is based on GI_Runoff_Vault_m3 (m3) and Pollutant_ConcentrationMean_mg_p_L and removal fraction
		beC->by_key(DataDrawer_ID, DataFolder_ID, Name_Load_GI_Runoff_Vault_Pollutant_Mean_g) = (beC->by_key(DataDrawer_ID, DataFolder_ID, "GI_Runoff_Vault_m3") * Ratio_L_to_m3 *
			input->Pollutant_ConcentrationMean_mg_p_L[pollutant_ID] * Ratio_g_to_mg *
			//(1 - beC->by_key(MapPixel_ID, DataFolder_ID, Name_Removal_Rate_Pollutant_percent) is (1 - removal fraction), or the pollutant release as a load
			(1 - beC->by_key(DataDrawer_ID, DataFolder_ID, Name_Removal_Rate_Pollutant_percent) * Ratio_Decimal_to_Percent));

		//Name_Load_GI_Drainage_VadoseZone_Pollutant_Mean_g (g) is based on GI_Drainage_VadoseZone_m3 (m3) and Pollutant_ConcentrationMean_mg_p_L and removal fraction
		//Note: GI_Drainage_VadoseZone_m3 (m3) is Exfiltration_m3, Percolation_Soil_m3, or Percolation_Pavement_m3, based on GI device layers in  Drainage_StormwaterDevice::Drainage_VadoseZone
		beC->by_key(DataDrawer_ID, DataFolder_ID, Name_Load_GI_Drainage_VadoseZone_Pollutant_Mean_g) = (beC->by_key(DataDrawer_ID, DataFolder_ID, "GI_Drainage_VadoseZone_m3") * Ratio_L_to_m3 *
			input->Pollutant_ConcentrationMean_mg_p_L[pollutant_ID] * Ratio_g_to_mg *
			//(1 - beC->by_key(MapPixel_ID, DataFolder_ID, Name_Removal_Rate_Pollutant_percent) is (1 - removal fraction), or the pollutant release as a load
			(1 - beC->by_key(DataDrawer_ID, DataFolder_ID, Name_Removal_Rate_Pollutant_percent) * Ratio_Decimal_to_Percent));

		//Name_Release_GI_Total_Pollutant_Mean_g (g) is pollutant removal mass (g) as difference between Inflow Name_Load_GI_Inflow_Pollutant_Mean_g (g) and ...
		//Note: ...Outflows: Name_Load_GI_Runoff_Surface_Pollutant_Mean_g + Name_Load_GI_Runoff_Vault_Pollutant_Mean_g + Name_Load_GI_Drainage_VadoseZone_Pollutant_Mean_g
		//Note: Name_Release_GI_Total_Pollutant_Mean_g can become negative at any time step, which is evidence of pollutant release and is recorded as negative removal
		beC->by_key(DataDrawer_ID, DataFolder_ID, Name_Release_GI_Total_Pollutant_Mean_g) = beC->by_key(DataDrawer_ID, DataFolder_ID, Name_Load_GI_Runoff_Surface_Pollutant_Mean_g) +
			beC->by_key(DataDrawer_ID, DataFolder_ID, Name_Load_GI_Runoff_Vault_Pollutant_Mean_g) + beC->by_key(DataDrawer_ID, DataFolder_ID, Name_Load_GI_Drainage_VadoseZone_Pollutant_Mean_g);

	}
}

void Removal_StormwaterPollutants::reduce_EC_EMC_after_GI(Inputs* input) {
	
	int nRows = Inputs::nRows;
	int nCols = Inputs::nCols;
	double NODATA_code = Inputs::NODATA_code;
	//GI processed catchment totals for all time steps
	string Name_Load_GI_Inflow_Pollutant_Mean_Sum_kg;
	string Name_Release_GI_Total_Pollutant_Mean_Sum_kg;

	

	//create variable names 3 types of flow and 10 types of pollutants 
	for (int pollutant_ID = 0; pollutant_ID < input->Pollutant_Name.size(); pollutant_ID++) {
		//Create string to represent variable associated with Pollutant_Name
		//string Name_Removal_Rate_Pollutant_percent = input->Pollutant_Name[pollutant_ID] + "_RemovalEfficiency_percent";
		//Create string of variable name within beC->vec_by_key, comprised of flow type label in quotes, Pollutant_Name vector, and "_Mean_Sum_kg"
		Name_Load_GI_Inflow_Pollutant_Mean_Sum_kg = "Load_GI_Inflow_" + input->Pollutant_Name[pollutant_ID] + "_Mean_Sum_kg";
		Name_Release_GI_Total_Pollutant_Mean_Sum_kg = "Release_GI_Total_" + input->Pollutant_Name[pollutant_ID] + "_Mean_Sum_kg";

		double removalrate_frac = (input->RepoDict[Name_Load_GI_Inflow_Pollutant_Mean_Sum_kg] - input->RepoDict[Name_Release_GI_Total_Pollutant_Mean_Sum_kg])
			/ input->RepoDict[Name_Load_GI_Inflow_Pollutant_Mean_Sum_kg];

		//for 3 types of pollutant, TP, TSS and TKN (assume TKN = TN)
		if (input->Pollutant_Name[pollutant_ID] == "TP") {
			for (int i = 0; i < nRows * nCols; i++) {
				//ContributingAreaPixels_map will have contributing pixel > 0 , others are nodata
				if (BufferSpatialCalc::ContributingAreaPixels_map_writeout[i] > 0) {
					//only apply the removal rate to non-nodata pixels EC
					if (BufferSpatialCalc::total_ec_tp_wtd[i] != NODATA_code) {
						BufferSpatialCalc::total_ec_tp_wtd[i] = BufferSpatialCalc::total_ec_tp_wtd[i] * (1 - removalrate_frac);
					}
					//only apply the removal rate to non-nodata pixels EMC
					if (BufferSpatialCalc::total_emc_tp_wtd[i] != NODATA_code) {
						BufferSpatialCalc::total_emc_tp_wtd[i] = BufferSpatialCalc::total_emc_tp_wtd[i] * (1 - removalrate_frac);
					}
				}
			}
		}
		if (input->Pollutant_Name[pollutant_ID] == "TSS") {
			for (int i = 0; i < nRows * nCols; i++) {
				//ContributingAreaPixels_map will have contributing pixel > 0 , others are nodata
				if (BufferSpatialCalc::ContributingAreaPixels_map_writeout[i] > 0) {
					//only apply the removal rate to non-nodata pixels EC
					if (BufferSpatialCalc::total_ec_tss_wtd[i] != NODATA_code) {
						BufferSpatialCalc::total_ec_tss_wtd[i] = BufferSpatialCalc::total_ec_tss_wtd[i] * (1 - removalrate_frac);
					}
					//only apply the removal rate to non-nodata pixels EMC
					if (BufferSpatialCalc::total_emc_tss_wtd[i] != NODATA_code) {
						BufferSpatialCalc::total_emc_tss_wtd[i] = BufferSpatialCalc::total_emc_tss_wtd[i] * (1 - removalrate_frac);
					}
				}
			}
		}

		if (input->Pollutant_Name[pollutant_ID] == "TKN") {
			for (int i = 0; i < nRows * nCols; i++) {
				//ContributingAreaPixels_map will have contributing pixel > 0 , others are nodata
				if (BufferSpatialCalc::ContributingAreaPixels_map_writeout[i] > 0) {
					//only apply the removal rate to non-nodata pixels EC
					if (BufferSpatialCalc::total_ec_tn_wtd[i] != NODATA_code) {
						BufferSpatialCalc::total_ec_tn_wtd[i] = BufferSpatialCalc::total_ec_tn_wtd[i] * (1 - removalrate_frac);
					}
					//only apply the removal rate to non-nodata pixels EMC
					if (BufferSpatialCalc::total_emc_tn_wtd[i] != NODATA_code) {
						BufferSpatialCalc::total_emc_tn_wtd[i] = BufferSpatialCalc::total_emc_tn_wtd[i] * (1 - removalrate_frac);
					}
				}
			}
		}
	}
}