#include "BufferDynamicCalc.h"
#include "../Outputs/WaterQuality_Output.h"
#include <algorithm>
#include <math.h>

//Author(s) : Li Zhang 
//lzhan126@syr.edu Feb 2021

//#include <boost/math/distributions/normal.hpp>


double BufferDynamicCalc::EmpericalCumulativeDistributionFunction_calc(double eventdata, Inputs* input)
{
    //Number of elements < event data
    double Count_NumberOfEvents=0;
    double EmpericalCumulativeDistributionFunction_calc = 0.00;
    vector<double> HydrologicData_Record = input->DischargeRecord;

    for (auto it = HydrologicData_Record.begin(); it != HydrologicData_Record.end(); ++it) {
        if (*it <= eventdata) { 
            ++Count_NumberOfEvents; 
        }
    }
    //EmpericalCumulativeDistributionFunction_calc range from 0-1
    EmpericalCumulativeDistributionFunction_calc = Count_NumberOfEvents / HydrologicData_Record.size();
    return EmpericalCumulativeDistributionFunction_calc;
}


double BufferDynamicCalc::quantile(const vector<double>& datavector, double Probability_Cummulative_frac, Inputs* input)
{
    //https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/quantile
    int HydrologicData_size = 0;
    int HydrologicVariable_rank_floor = 0;
    double Probability_Adjustment_frac = 0.0;
    double HydrologicVariable_rank_decimal = 0.0;
    double HydrologicVariable_QuantileValue_SourceUnits = 0.0;

    //HydrologicData_size is the sample size
    HydrologicData_size = int(datavector.size());
    //For types 4 through 9, m given below
    //type4: linear interpolation of the empirical cdf
    if (input->BufferTemporalLoadingNumericalParams["QuantileType"] == 4) {
        Probability_Adjustment_frac = 0;
    }
    //type5: piecewise linear function where the knots are the values midway through the steps of the empirical cdf
    if (input->BufferTemporalLoadingNumericalParams["QuantileType"] == 5) {
        Probability_Adjustment_frac = 0.5;
    }
    //type6: used by Minitaband by SPSS.
    if (input->BufferTemporalLoadingNumericalParams["QuantileType"] == 6) {
        Probability_Adjustment_frac = Probability_Cummulative_frac;
    }
    //type7: This is used by S.
    if (input->BufferTemporalLoadingNumericalParams["QuantileType"] == 7) {
        Probability_Adjustment_frac = 1 - Probability_Cummulative_frac;
    }
    //type8: The resulting quantile estimates are approximately median-unbiased regardless of the distribution of x
    if (input->BufferTemporalLoadingNumericalParams["QuantileType"] == 8) {
        Probability_Adjustment_frac = (Probability_Cummulative_frac + 1)/3;
    }
    //type9: The resulting quantile estimates are approximately unbiased for the expected order statistics if x is normally distributed
    if (input->BufferTemporalLoadingNumericalParams["QuantileType"] == 9) {
        Probability_Adjustment_frac = Probability_Cummulative_frac/4 + 3/8;
    }
    
    HydrologicVariable_rank_floor = floor(HydrologicData_size * Probability_Cummulative_frac + Probability_Adjustment_frac);
    HydrologicVariable_rank_decimal = HydrologicData_size * Probability_Cummulative_frac + Probability_Adjustment_frac - HydrologicVariable_rank_floor;
    
    //HydrologicVariable_rank_floor = order number; 
    //insert one more data (0) at the beginning of the vector, so we can use the same order as R code: HydrologicVariable_rank_floor & HydrologicVariable_rank_floor+1
    //if we remove this line, since the vector begin from HydrologicVariable_rank_floor=0, HydrologicVariable_rank_floor & HydrologicVariable_rank_floor+1 in the equation-> HydrologicVariable_rank_floor-1 & HydrologicVariable_rank_floor
    //VectorSorted.insert(VectorSorted.begin(), 0);
    //HydrologicVariable_QuantileValue_SourceUnits = (1 - HydrologicVariable_rank_decimal) * VectorSorted[HydrologicVariable_rank_floor] + HydrologicVariable_rank_decimal * VectorSorted[HydrologicVariable_rank_floor+1];
    if (HydrologicVariable_rank_floor == int(datavector.size())) {
        HydrologicVariable_QuantileValue_SourceUnits = (1 - HydrologicVariable_rank_decimal) * datavector[HydrologicVariable_rank_floor - 1];
    }
    else {
        HydrologicVariable_QuantileValue_SourceUnits = (1 - HydrologicVariable_rank_decimal) * datavector[HydrologicVariable_rank_floor - 1] +
            HydrologicVariable_rank_decimal * datavector[HydrologicVariable_rank_floor];
    }

     return HydrologicVariable_QuantileValue_SourceUnits;
}

vector< vector<double> > BufferDynamicCalc::pollutant_landcover_list;
void BufferDynamicCalc::Interpolation(Inputs* input) {
    double Q1st, Q2nd, Qlast, Qlast2nd,slope_Q_begin, slope_Q_end,EC_begin, EC_end;
    vector<double> HydrologicData_Record = input->DischargeRecord;
    vector<double> DataRecordSorted(HydrologicData_Record.size());
    //sort the discharge record
    partial_sort_copy(HydrologicData_Record.begin(), HydrologicData_Record.end() , DataRecordSorted.begin(), DataRecordSorted.end());
    //use the quantile function to get the HydrologicVariable_QuantileValue_SourceUnits at the tails
    Q1st= quantile(DataRecordSorted, 0, input);//0th
    Qlast = quantile(DataRecordSorted, 1, input);//100th
    Q2nd = quantile(DataRecordSorted, 0.05, input);//5th
    Qlast2nd = quantile(DataRecordSorted, 0.95, input);//95th
    //find the slope at the tails
    if (Q2nd == 0) {
        slope_Q_begin = 0;
    }
    else slope_Q_begin = Q1st / Q2nd;

    if (Qlast2nd == 0) {
        slope_Q_end = 0;
    }
    else slope_Q_end = Qlast / Qlast2nd;
    //creating list of vectors
    //Total Phosphorus
    //Urban
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TP"]);//21
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TP"]);//22
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TP"]);//23
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TP"]);//24
    //Forest
    pollutant_landcover_list.push_back(input->EC_vec_map["Forest_TP"]);//41
    pollutant_landcover_list.push_back(input->EC_vec_map["Forest_TP"]);//42
    pollutant_landcover_list.push_back(input->EC_vec_map["Forest_TP"]);//43
    //Rangeland
    pollutant_landcover_list.push_back(input->EC_vec_map["Rangeland_TP"]);//71
    //Pasture_hay
    pollutant_landcover_list.push_back(input->EC_vec_map["Pasture_hay_TP"]);//81
    //Cropland
    pollutant_landcover_list.push_back(input->EC_vec_map["Cropland_TP"]);//82

    //Total Nitrogen 
    //Urban
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TN"]);//21
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TN"]);//22
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TN"]);//23
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TN"]);//24
    //Forest
    pollutant_landcover_list.push_back(input->EC_vec_map["Forest_TN"]);//41
    pollutant_landcover_list.push_back(input->EC_vec_map["Forest_TN"]);//42
    pollutant_landcover_list.push_back(input->EC_vec_map["Forest_TN"]);//43
    //Rangeland
    pollutant_landcover_list.push_back(input->EC_vec_map["Rangeland_TN"]);//71
    //Pasture_hay
    pollutant_landcover_list.push_back(input->EC_vec_map["Pasture_hay_TN"]);//81
    //Cropland
    pollutant_landcover_list.push_back(input->EC_vec_map["Cropland_TN"]);//82

    //total suspended solids
    //Urban
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TSS"]);//21
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TSS"]);//22
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TSS"]);//23
    pollutant_landcover_list.push_back(input->EC_vec_map["Urban_TSS"]);//24
    //Forest
    pollutant_landcover_list.push_back(input->EC_vec_map["Forest_TSS"]);//41
    pollutant_landcover_list.push_back(input->EC_vec_map["Forest_TSS"]);//42
    pollutant_landcover_list.push_back(input->EC_vec_map["Forest_TSS"]);//43
    //Rangeland
    pollutant_landcover_list.push_back(input->EC_vec_map["Rangeland_TSS"]);//71
    //Pasture_hay
    pollutant_landcover_list.push_back(input->EC_vec_map["Pasture_hay_TSS"]);//81
    //Cropland
    pollutant_landcover_list.push_back(input->EC_vec_map["Cropland_TSS"]);//82


    for (auto list_it = pollutant_landcover_list.begin(); list_it != pollutant_landcover_list.end(); ++list_it) {
        //vector<double> &V = *list_it;
        int num_insert = 0;
        EC_begin = (*list_it)[0] * slope_Q_begin;
        EC_end = (*list_it).back() * slope_Q_end;
        (*list_it).insert((*list_it).begin(), EC_begin);
        (*list_it).insert((*list_it).end(), EC_end);
        //now we have 0th, 5th, 10th, 25th, 50th, 75th, 90th, 95th, 100th (9 in total)
        //insert 12 new values using linear interpolation, to have 5th interval (21 in total)
        int interval = 5;
        //filling the 1st gap
        double x_increment1 = 15;//25-10
        double y_increment1 = (*list_it)[3] - (*list_it)[2];
        double slope1 = y_increment1 / x_increment1;
        for (auto i = 0; i < x_increment1 / interval-1; ++i) {//insert 2 times
            //inserts EC at (3+i)-th location, for this loop, we add 2 values, add at begin() + 3 and begin() + 4
            auto it = (*list_it).insert((*list_it).begin() + 3+i, (*list_it)[2]+slope1* interval*(i+1));
        }
        num_insert += x_increment1 / interval - 1;
        
        //filling the 2nd gap
        double x_increment2 = 25;//50-25
        double y_increment2 = (*list_it)[3 + num_insert+1] - (*list_it)[3+num_insert ];
        double slope2 = y_increment2 / x_increment2;
        for (auto i = 0; i < x_increment2 / interval - 1; ++i) {//insert 4 times
            auto it = (*list_it).insert((*list_it).begin() + 4+ num_insert + i, (*list_it)[3+ num_insert]+ slope2 * interval * (i + 1));
        }
        num_insert += x_increment2 / interval - 1;

        //filling the 3rd gap
        double x_increment3 = 25;//75-50
        double y_increment3 = (*list_it)[4+ num_insert+1] - (*list_it)[4 + num_insert];
        double slope3 = y_increment3 / x_increment3;
        for (auto i = 0; i < x_increment3 / interval - 1; ++i) {//insert 4 times
            auto it = (*list_it).insert((*list_it).begin() + 5 + num_insert + i, (*list_it)[4+ num_insert] + slope3 * interval * (i + 1));
        }
        num_insert += x_increment3 / interval - 1;

        //filling the 4rd gap
        double x_increment4 = 15;//90-75
        double y_increment4 = (*list_it)[5 + num_insert + 1] - (*list_it)[5 + num_insert];
        double slope4 = y_increment4 / x_increment4;
        for (auto i = 0; i < x_increment4 / interval - 1; ++i) {//insert 2 times
            auto it = (*list_it).insert((*list_it).begin() + 6 + num_insert + i, (*list_it)[5 + num_insert] + slope4 * interval * (i + 1));
        }
        num_insert += x_increment4 / interval - 1;
        
        /*cout << Qlast << "," << Qlast2nd << endl;
        vector<double> temp = *list_it;
        for (auto w = temp.begin(); w !=temp.end(); w++) {
            
            cout << *w << ", ";
        }
        cout << endl;*/
    }
}

double BufferDynamicCalc::LCEventload(const vector<double>& Landuse_EC_vectors, double Eventdata, double cell_num, Inputs* input) {
    double LandCover_ExportCoefficient_kgphapyr = 0.0;
    double Event_EmpericalCumulativeDistributionFunction = 0.0;
    double Event_Total_ExportCoefficient_kgphapyr = 0.0;
    double LC_area_fraction = 0.0;
    vector<double> HydrologicData_Record = input->DischargeRecord;
    //get the EmpericalCumulativeDistributionFunction_calc for this event, function returns cumulative distribution probability between 0 and 1, it is (1- exceedance probability)
    //DataRecord is a vector that reading the daily discharge, Eventdata is a scaler of that day's discharge value, cfs
    Event_EmpericalCumulativeDistributionFunction = EmpericalCumulativeDistributionFunction_calc(Eventdata, input);
    //LandCover_ExportCoefficient_kgphapyr use the quantile function to find the EC of this event in this landcover
    LandCover_ExportCoefficient_kgphapyr = quantile(Landuse_EC_vectors, Event_EmpericalCumulativeDistributionFunction, input);
    //EC_L,t * (AL/A), the contribution of pollutant from each landcover type
    LC_area_fraction = cell_num / input->cell_count_Nowater;
    Event_Total_ExportCoefficient_kgphapyr = LandCover_ExportCoefficient_kgphapyr * LC_area_fraction;
    return Event_Total_ExportCoefficient_kgphapyr; 
}

//YearAverage function will loop through all years in simulation
void BufferDynamicCalc::YearAverage(Inputs* input) {

    vector<double> HydrologicData_Record = input->DischargeRecord;
    vector<int> YearRecord = input->DischargeYear;
    vector<int> MonthRecord = input->DischargeMonth;
    vector<int> DayRecord = input->DischargeDay;

    double Eventdata = 0.000;
    vector<double> uniqueyear;
    vector<double>AverageYearLoad_TP;
    vector<double>AverageYearLoad_TN;
    vector<double>AverageYearLoad_TSS;

    double Area_Watershed_ha = 0.0;
    double cell_size_m = 0.0;
    double cell_count = 0.0;
    
    //vectors saving data for each time step, to write extened output
    vector<double> TS_Load_TP_entered, TS_Load_TN_entered, TS_Load_TSS_entered,
        TS_Load_TP_delivered, TS_Load_TN_delivered, TS_Load_TSS_delivered,
        TS_Storage_TP_enter, TS_Storage_TN_enter, TS_Storage_TSS_enter,
        TS_Storage_TP_NA, TS_Storage_TN_NA, TS_Storage_TSS_NA,
        TS_Storage_TP_release, TS_Storage_TN_release, TS_Storage_TSS_release,
        TS_Storage_TP_kg, TS_Storage_TN_kg, TS_Storage_TSS_kg,
        TS_daily_TP_load, TS_daily_TN_load, TS_daily_TSS_load,
        TS_Delivery_ratio, TS_ActiveRiverLength_m, TS_ECDF,
        TS_date_year, TS_date_month, TS_date_day, TS_datarecord;


    double river_length_m = input->BufferTemporalLoadingNumericalParams["River_length_km"] * 1000;
    double ActiveRiverLength_Coefficient_C = input->BufferTemporalLoadingNumericalParams["ActiveRiverLength_Coefficient_C"];
    double DischargeStorageRelease_frac = input->BufferTemporalLoadingNumericalParams["DischargeStorageRelease_frac"];
    double Base_delivery_ratio = input->BufferTemporalLoadingNumericalParams["Base_delivery_ratio_frac"];//sediment delivery ratio 
    double Lambda = input->BufferTemporalLoadingNumericalParams["Lambda"];
    double Active_Storage_days = input->BufferTemporalLoadingNumericalParams["Active_Storage_days"];
    //ActiveStorage_initial_TP_kg_p_day is the initial estimate of active storage rate
    double ActiveStorage_initial_TP_kg_p_day = input->BufferTemporalLoadingNumericalParams["ActiveStorage_initial_TP_kg_p_day"];
    double ActiveStorage_initial_TN_kg_p_day = input->BufferTemporalLoadingNumericalParams["ActiveStorage_initial_TN_kg_p_day"];
    double ActiveStorage_initial_TSS_kg_p_day = input->BufferTemporalLoadingNumericalParams["ActiveStorage_initial_TSS_kg_p_day"];
    //initialize
    double ActiveStorage_TP_kg = ActiveStorage_initial_TP_kg_p_day * Active_Storage_days;
    double ActiveStorage_TN_kg = ActiveStorage_initial_TN_kg_p_day * Active_Storage_days;
    double ActiveStorage_TSS_kg = ActiveStorage_initial_TSS_kg_p_day * Active_Storage_days;
    //Storage at Active_Storage_days
    double ActiveStorage_final_TP_p_day, ActiveStorage_final_TN_p_day, ActiveStorage_final_TSS_p_day;
    
    //Linear Interpolation
    Interpolation(input);

    //watershed Area_Watershed_ha is input->cell_count_Nowater * (Inputs::Length_Pixel_Side_m * Inputs::Length_Pixel_Side_m) * Ratio_ha_to_m2
    Area_Watershed_ha = input->cell_count_Nowater * (Inputs::Length_Pixel_Side_m * Inputs::Length_Pixel_Side_m) * Ratio_ha_to_m2;
    //DateParts_Gregorian initialized as vector integer
    int DateParts_Gregorian[4];
    //input->ExtractGregorianDateParts(input->SimulationNumericalParams["StartDate_YYYYMMDD"], DateParts_Gregorian); gets Year_Start
    input->ExtractGregorianDateParts(input->SimulationNumericalParams["StartDate_YYYYMMDD"], DateParts_Gregorian);
    //Year_Start_YYYY = DateParts_Gregorian[0]
    int Year_Start_YYYY = DateParts_Gregorian[0];
    //input->ExtractGregorianDateParts(input->SimulationNumericalParams["StopDate_YYYYMMDD"], DateParts_Gregorian); gets Year_Stop
    input->ExtractGregorianDateParts(input->SimulationNumericalParams["StopDate_YYYYMMDD"], DateParts_Gregorian);
    //Year_Stop_YYYY = DateParts_Gregorian[0]
    int Year_Stop_YYYY = DateParts_Gregorian[0];


    //for each Year_YYYY extracted from StartDate_YYYYMMDD to StopDate_YYYYMMDD
    for (int Year_YYYY = Year_Start_YYYY; Year_YYYY < Year_Stop_YYYY + 1; ++Year_YYYY) {
        //
        vector<double>datarecord_year;
        //set an iterator at the beginning of the discharge data vector, record_it function will iterate with year_it      
        auto year_it = YearRecord.begin();
        auto month_it = MonthRecord.begin();
        auto day_it = DayRecord.begin();
        for (auto record_it = HydrologicData_Record.begin(); record_it != HydrologicData_Record.end(); ++record_it) {
            if (*year_it == Year_YYYY) {
                datarecord_year.push_back(*record_it);
                //saving values to vectors, to write out the extened output
                TS_date_year.push_back(*year_it);
                TS_date_month.push_back(*month_it);
                TS_date_day.push_back(*day_it);
                TS_datarecord.push_back(*record_it);
            }
            ++year_it;
            ++month_it;
            ++day_it;
        }
        
        //sum of the load for each Year_YYYY
        double yearly_TP_load_kg = 0.0, yearly_TN_load_kg = 0.0, yearly_TSS_load_kg = 0.0;
        //for each event in this Year_YYYY
        //if use all data record to build the probability distribution, rather than using event > th, change DataRecord_th to DataRecord for the follwing line within the loop 
        for (auto record_day = datarecord_year.begin(); record_day != datarecord_year.end(); ++record_day) {
            double daily_TP_load_kg_d, daily_TN_load_kg_d, daily_TSS_load_kg_d;
            Eventdata = *record_day;
            //all of these parameters are temp, will be initialized for each day 
            double ECTP21, ECTP22, ECTP23, ECTP24, ECTP41, ECTP42, ECTP43, ECTP71, ECTP81, ECTP82;
            double ECTN21, ECTN22, ECTN23 , ECTN24, ECTN41, ECTN42, ECTN43, ECTN71, ECTN81, ECTN82;
            double ECTSS21, ECTSS22 , ECTSS23, ECTSS24, ECTSS41, ECTSS42, ECTSS43, ECTSS71, ECTSS81, ECTSS82;
            double ActiveRiverLength_m, Delivery_ratio, Length_ratio;
            double Load_TP_delivered_kg_d, Load_TN_delivered_kg_d, Load_TSS_delivered_kg_d;
            double Load_TP_entered_kg_d, Load_TN_entered_kg_d, Load_TSS_entered_kg_d;
            double Storage_TP_enter_kg_d, Storage_TN_enter_kg_d, Storage_TSS_enter_kg_d;
            double Storage_TP_NA_kg_d, Storage_TN_NA_kg_d, Storage_TSS_NA_kg_d;
            double Storage_TP_release_kg_d, Storage_TN_release_kg_d, Storage_TSS_release_kg_d;
            
            ECTP21 = LCEventload(pollutant_landcover_list[0], Eventdata, input->Cell_number[0], input);
            ECTP22 = LCEventload(pollutant_landcover_list[1], Eventdata, input->Cell_number[1], input);
            ECTP23 = LCEventload(pollutant_landcover_list[2], Eventdata, input->Cell_number[2], input);
            ECTP24 = LCEventload(pollutant_landcover_list[3], Eventdata, input->Cell_number[3], input);
            ECTP41 = LCEventload(pollutant_landcover_list[4], Eventdata, input->Cell_number[4], input);
            ECTP42 = LCEventload(pollutant_landcover_list[5], Eventdata, input->Cell_number[5], input);
            ECTP43 = LCEventload(pollutant_landcover_list[6], Eventdata, input->Cell_number[6], input);
            ECTP71 = LCEventload(pollutant_landcover_list[7], Eventdata, input->Cell_number[7], input);
            ECTP81 = LCEventload(pollutant_landcover_list[8], Eventdata, input->Cell_number[8], input);
            ECTP82 = LCEventload(pollutant_landcover_list[9], Eventdata, input->Cell_number[9], input);

            ECTN21 = LCEventload(pollutant_landcover_list[10], Eventdata, input->Cell_number[0], input);
            ECTN22 = LCEventload(pollutant_landcover_list[11], Eventdata, input->Cell_number[1], input);
            ECTN23 = LCEventload(pollutant_landcover_list[12], Eventdata, input->Cell_number[2], input);
            ECTN24 = LCEventload(pollutant_landcover_list[13], Eventdata, input->Cell_number[3], input);
            ECTN41 = LCEventload(pollutant_landcover_list[14], Eventdata, input->Cell_number[4], input);
            ECTN42 = LCEventload(pollutant_landcover_list[15], Eventdata, input->Cell_number[5], input);
            ECTN43 = LCEventload(pollutant_landcover_list[16], Eventdata, input->Cell_number[6], input);
            ECTN71 = LCEventload(pollutant_landcover_list[17], Eventdata, input->Cell_number[7], input);
            ECTN81 = LCEventload(pollutant_landcover_list[18], Eventdata, input->Cell_number[8], input);
            ECTN82 = LCEventload(pollutant_landcover_list[19], Eventdata, input->Cell_number[9], input);

            ECTSS21 = LCEventload(pollutant_landcover_list[20], Eventdata, input->Cell_number[0], input);
            ECTSS22 = LCEventload(pollutant_landcover_list[21], Eventdata, input->Cell_number[1], input);
            ECTSS23 = LCEventload(pollutant_landcover_list[22], Eventdata, input->Cell_number[2], input);
            ECTSS24 = LCEventload(pollutant_landcover_list[23], Eventdata, input->Cell_number[3], input);
            ECTSS41 = LCEventload(pollutant_landcover_list[14], Eventdata, input->Cell_number[4], input);
            ECTSS42 = LCEventload(pollutant_landcover_list[25], Eventdata, input->Cell_number[5], input);
            ECTSS43 = LCEventload(pollutant_landcover_list[26], Eventdata, input->Cell_number[6], input);
            ECTSS71 = LCEventload(pollutant_landcover_list[27], Eventdata, input->Cell_number[7], input);
            ECTSS81 = LCEventload(pollutant_landcover_list[28], Eventdata, input->Cell_number[8], input);
            ECTSS82 = LCEventload(pollutant_landcover_list[29], Eventdata, input->Cell_number[9], input);
           
            //Load entering the drainage network, convert from kg/yr/ha to kg/day
            Load_TP_entered_kg_d = (ECTP21 + ECTP22 + ECTP23 + ECTP24 + ECTP41 + ECTP42 + ECTP43 + ECTP71 + ECTP81 + ECTP82) / datarecord_year.size() * Area_Watershed_ha;
            Load_TN_entered_kg_d = (ECTN21 + ECTN22 + ECTP23 + ECTN24 + ECTN41 + ECTN42 + ECTN43 + ECTN71 + ECTN81 + ECTN82) / datarecord_year.size() * Area_Watershed_ha;
            Load_TSS_entered_kg_d = (ECTSS21 + ECTSS22 + ECTSS23 + ECTSS24 + ECTSS41 + ECTSS42 + ECTSS43 + ECTSS71 + ECTSS81 + ECTSS82) / datarecord_year.size() * Area_Watershed_ha;

            //Dynamic ActiveRiverLength_m Eq
            if (DischargeStorageRelease_frac == 1){
                ActiveRiverLength_m = 0;
            }
            else {
                ActiveRiverLength_m = ActiveRiverLength_Coefficient_C * river_length_m / (1 - DischargeStorageRelease_frac) * exp(-Lambda * (1 + (1 - EmpericalCumulativeDistributionFunction_calc(Eventdata, input)) / (1 - DischargeStorageRelease_frac)));
            }
            
            //Dynamic Delivery_ratio Eq
            if (ActiveRiverLength_m > river_length_m) {
                Length_ratio = 1;
                Delivery_ratio = Length_ratio;
            }
            else {
                Length_ratio = ActiveRiverLength_m / river_length_m;
                if (DischargeStorageRelease_frac == 1) {
                    Delivery_ratio = Base_delivery_ratio;
                }
                else{
                    Delivery_ratio = Base_delivery_ratio + Lambda * exp(-Lambda * (1 + (1 - EmpericalCumulativeDistributionFunction_calc(Eventdata, input)) / (1 - DischargeStorageRelease_frac)) + 1);
                }
            }
            //Clamp to ensure variable is between 0 and 1 
            Delivery_ratio = clamp(Delivery_ratio, 0.0, 1.0);

            //Load released to the basin outlet
            Load_TP_delivered_kg_d = Load_TP_entered_kg_d * Delivery_ratio;
            Load_TN_delivered_kg_d = Load_TN_entered_kg_d * Delivery_ratio;
            Load_TSS_delivered_kg_d = Load_TSS_entered_kg_d * Delivery_ratio;
             
            //if run storage
            if (input->BufferTemporalLoadingNumericalParams["Flag_RunStorage"] == 1) {
                //pollutant load entering storage within the drainage network 
                Storage_TP_enter_kg_d = Load_TP_entered_kg_d - Load_TP_delivered_kg_d;
                Storage_TN_enter_kg_d = Load_TN_entered_kg_d - Load_TN_delivered_kg_d;
                Storage_TSS_enter_kg_d = Load_TSS_entered_kg_d - Load_TSS_delivered_kg_d;
            
                //The value of the unavailable storage for each time step 
                Storage_TP_NA_kg_d = ActiveStorage_TP_kg / Active_Storage_days;
                Storage_TN_NA_kg_d = ActiveStorage_TN_kg / Active_Storage_days;
                Storage_TSS_NA_kg_d = ActiveStorage_TSS_kg / Active_Storage_days;

                //1st time update the storage
                ActiveStorage_TP_kg = ActiveStorage_TP_kg + Storage_TP_enter_kg_d - Storage_TP_NA_kg_d;
                ActiveStorage_TN_kg = ActiveStorage_TN_kg + Storage_TN_enter_kg_d - Storage_TN_NA_kg_d;
                ActiveStorage_TSS_kg = ActiveStorage_TSS_kg + Storage_TSS_enter_kg_d - Storage_TSS_NA_kg_d;
                
                // if p>threshold, release some pollutant from the storage
                if (EmpericalCumulativeDistributionFunction_calc(Eventdata, input) > DischargeStorageRelease_frac) {
                    Storage_TP_release_kg_d = ActiveStorage_TP_kg * Delivery_ratio * Length_ratio;
                    Storage_TN_release_kg_d = ActiveStorage_TN_kg * Delivery_ratio * Length_ratio;
                    Storage_TSS_release_kg_d = ActiveStorage_TSS_kg * Delivery_ratio * Length_ratio;
                
                    if (Storage_TP_release_kg_d > ActiveStorage_TP_kg) {
                        Storage_TP_release_kg_d = ActiveStorage_TP_kg;
                    }
                    if (Storage_TN_release_kg_d > ActiveStorage_TN_kg) {
                        Storage_TN_release_kg_d = ActiveStorage_TN_kg;
                    }
                    if (Storage_TSS_release_kg_d > ActiveStorage_TSS_kg) {
                        Storage_TSS_release_kg_d = ActiveStorage_TSS_kg;
                    }
                    //2nd time update the storage
                    ActiveStorage_TP_kg = ActiveStorage_TP_kg - Storage_TP_release_kg_d;
                    ActiveStorage_TN_kg = ActiveStorage_TN_kg - Storage_TN_release_kg_d;
                    ActiveStorage_TSS_kg = ActiveStorage_TSS_kg - Storage_TSS_release_kg_d;
                }
                else {
                    Storage_TP_release_kg_d = 0;
                    Storage_TN_release_kg_d = 0;
                    Storage_TSS_release_kg_d = 0;
                }
                //load for each day
                daily_TP_load_kg_d = Load_TP_delivered_kg_d + Storage_TP_release_kg_d;
                daily_TN_load_kg_d = Load_TN_delivered_kg_d + Storage_TN_release_kg_d;
                daily_TSS_load_kg_d = Load_TSS_delivered_kg_d + Storage_TSS_release_kg_d;
            }
            else {// if not run storage, daily load will be load delivered
                daily_TP_load_kg_d = Load_TP_delivered_kg_d;
                daily_TN_load_kg_d = Load_TN_delivered_kg_d;
                daily_TSS_load_kg_d = Load_TSS_delivered_kg_d;
            }

            
            //accumulating load for a Year_YYYY
            yearly_TP_load_kg += daily_TP_load_kg_d;
            yearly_TN_load_kg += daily_TN_load_kg_d;
            yearly_TSS_load_kg += daily_TSS_load_kg_d;
            
            //saving discharge event data
            TS_ECDF.push_back(EmpericalCumulativeDistributionFunction_calc(Eventdata, input));
            TS_Delivery_ratio.push_back(Delivery_ratio);
            TS_ActiveRiverLength_m.push_back(ActiveRiverLength_m);

            //saving data into vectors for each time step
            TS_Load_TP_entered.push_back(Load_TP_entered_kg_d);
            TS_Load_TN_entered.push_back(Load_TN_entered_kg_d);
            TS_Load_TSS_entered.push_back(Load_TSS_entered_kg_d);
            
            TS_Load_TP_delivered.push_back(Load_TP_delivered_kg_d);
            TS_Load_TN_delivered.push_back(Load_TN_delivered_kg_d);
            TS_Load_TSS_delivered.push_back(Load_TSS_delivered_kg_d);

            //if run storage
            if (input->BufferTemporalLoadingNumericalParams["Flag_RunStorage"] == 1) {
                TS_Storage_TP_enter.push_back(Storage_TP_enter_kg_d);
                TS_Storage_TN_enter.push_back(Storage_TN_enter_kg_d);
                TS_Storage_TSS_enter.push_back(Storage_TSS_enter_kg_d);

                TS_Storage_TP_NA.push_back(Storage_TP_NA_kg_d);
                TS_Storage_TN_NA.push_back(Storage_TN_NA_kg_d);
                TS_Storage_TSS_NA.push_back(Storage_TSS_NA_kg_d);

                TS_Storage_TP_release.push_back(Storage_TP_release_kg_d);
                TS_Storage_TN_release.push_back(Storage_TN_release_kg_d);
                TS_Storage_TSS_release.push_back(Storage_TSS_release_kg_d);

                TS_Storage_TP_kg.push_back(ActiveStorage_TP_kg);
                TS_Storage_TN_kg.push_back(ActiveStorage_TN_kg);
                TS_Storage_TSS_kg.push_back(ActiveStorage_TSS_kg);
            }
            

            TS_daily_TP_load.push_back(daily_TP_load_kg_d);
            TS_daily_TN_load.push_back(daily_TN_load_kg_d);
            TS_daily_TSS_load.push_back(daily_TSS_load_kg_d);
        }
        //saving yearly load to vector, kg/yr
        AverageYearLoad_TP.push_back(yearly_TP_load_kg);
        AverageYearLoad_TN.push_back(yearly_TN_load_kg);
        AverageYearLoad_TSS.push_back(yearly_TSS_load_kg);

    }


    vector <pair<string, vector<double>>> result_TP_map, result_TN_map, result_TSS_map;

    //If BufferTemporalLoadingNumericalParams["Flag_RunStorage"] == 1 then report storage results
    if (input->BufferTemporalLoadingNumericalParams["Flag_RunStorage"] == 1) {

        // Helper lambda function to print pollutant storage summary
        auto report_storage_gap = [](const string& name_pollutant, double initial_storage_kg, double final_storage_kg, double gap_frac) {
            cout << "Pollutant " << name_pollutant << " had a final average storage rate of " << final_storage_kg << " kg/day" << endl;
            cout << "... and an initial storage rate set to " << initial_storage_kg << " kg/day" << endl;
            //If gap_frac > 0.05 then suggest update of initial storage
            if (abs(gap_frac) > 0.05) {
                cout << "... The " << gap_frac * 100 << " % difference suggest the HydroPlusConfig.xml initial storage is updated with the final value." << endl;
            }
        };

        //for user manully spin up to improve the initial guess of active storage
        double gap_P_frac, gap_N_frac, gap_SS_frac;
        // Compute final values
        ActiveStorage_final_TP_p_day = TS_Storage_TP_kg[Active_Storage_days] / Active_Storage_days;
        ActiveStorage_final_TN_p_day = TS_Storage_TN_kg[Active_Storage_days] / Active_Storage_days;
        ActiveStorage_final_TSS_p_day = TS_Storage_TSS_kg[Active_Storage_days] / Active_Storage_days;

        // Compute percent gaps
        gap_P_frac = (ActiveStorage_initial_TP_kg_p_day - ActiveStorage_final_TP_p_day) / ActiveStorage_final_TP_p_day;
        gap_N_frac = (ActiveStorage_initial_TN_kg_p_day - ActiveStorage_final_TN_p_day) / ActiveStorage_final_TN_p_day;
        gap_SS_frac = (ActiveStorage_initial_TSS_kg_p_day - ActiveStorage_final_TSS_p_day) / ActiveStorage_final_TSS_p_day;

        // Print reports
        report_storage_gap("TP", ActiveStorage_initial_TP_kg_p_day, ActiveStorage_final_TP_p_day, gap_P_frac);
        report_storage_gap("TN", ActiveStorage_initial_TN_kg_p_day, ActiveStorage_final_TN_p_day, gap_N_frac);
        report_storage_gap("TSS", ActiveStorage_initial_TSS_kg_p_day, ActiveStorage_final_TSS_p_day, gap_SS_frac);
       
        //result_TP_map prepared
        result_TP_map = {
            {"Year", TS_date_year},
            {"Month", TS_date_month},
            {"Day", TS_date_day},
            {"Discharge", TS_datarecord},
            {"p_cumulative", TS_ECDF},
            {"Delivery_ratio", TS_Delivery_ratio},
            {"ActiveRiverLength_m", TS_ActiveRiverLength_m},
            {"Load_TP_entered_kg_d", TS_Load_TP_entered},
            {"Load_TP_delivered_kg_d", TS_Load_TP_delivered},
            {"Storage_TP_enter_kg_d", TS_Storage_TP_enter}, 
            {"Storage_TP_NA_kg_d", TS_Storage_TP_NA},
            {"Storage_TP_release_kg_d", TS_Storage_TP_release},
            {"ActiveStorage_TP_kg", TS_Storage_TP_kg},
            {"Daily_TP_load_kg_d", TS_daily_TP_load}
        };
        //result_TN_map prepared
       result_TN_map = {
            {"Year", TS_date_year},
            {"Month", TS_date_month},
            {"Day", TS_date_day},
            {"Discharge", TS_datarecord},
            {"p_cumulative", TS_ECDF},
            {"Delivery_ratio", TS_Delivery_ratio},
            {"ActiveRiverLength_m", TS_ActiveRiverLength_m},
            {"Load_TN_entered_kg_d", TS_Load_TN_entered},
            {"Load_TN_delivered_kg_d", TS_Load_TN_delivered},
            {"Storage_TN_enter_kg_d", TS_Storage_TN_enter},
            {"Storage_TN_NA_kg_d", TS_Storage_TN_NA},
            {"Storage_TN_release_kg_d", TS_Storage_TN_release},
            {"ActiveStorage_TN_kg", TS_Storage_TN_kg},
            {"Daily_TN_load_kg_d", TS_daily_TN_load},
            
       };
       //result_TSS_map prepared
       result_TSS_map = {
            {"Year", TS_date_year},
            {"Month", TS_date_month},
            {"Day", TS_date_day},
            {"Discharge", TS_datarecord},
            {"p_cumulative", TS_ECDF},
            {"Delivery_ratio", TS_Delivery_ratio},
            {"ActiveRiverLength_m", TS_ActiveRiverLength_m},
            {"Load_TSS_entered_kg_d", TS_Load_TSS_entered},
            {"Load_TSS_delivered_kg_d", TS_Load_TSS_delivered},
            {"Storage_TSS_enter_kg_d", TS_Storage_TSS_enter},
            {"Storage_TSS_NA_kg_d", TS_Storage_TSS_NA},
            {"Storage_TSS_release_kg_d", TS_Storage_TSS_release},
            {"ActiveStorage_TSS_kg", TS_Storage_TSS_kg},
            {"Daily_TSS_load_kg_d", TS_daily_TSS_load}
       };
    }
    //Else if not running storage
    else {
        result_TP_map = {
            {"Year", TS_date_year},
            {"Month", TS_date_month},
            {"Day", TS_date_day},
            {"Discharge", TS_datarecord},
            {"p_cumulative", TS_ECDF},
            {"Delivery_ratio", TS_Delivery_ratio},
            {"ActiveRiverLength_m", TS_ActiveRiverLength_m},

            {"Load_TP_entered_kg_d", TS_Load_TP_entered},
            {"Load_TP_delivered_kg_d", TS_Load_TP_delivered},

            {"Daily_TP_load_kg_d", TS_daily_TP_load}
        };
        result_TN_map = {
            {"Year", TS_date_year},
            {"Month", TS_date_month},
            {"Day", TS_date_day},
            {"Discharge", TS_datarecord},
            {"p_cumulative", TS_ECDF},
            {"Delivery_ratio", TS_Delivery_ratio},
            {"ActiveRiverLength_m", TS_ActiveRiverLength_m},

            {"Load_TN_entered_kg_d", TS_Load_TN_entered},
            {"Load_TN_delivered_kg_d", TS_Load_TN_delivered},

            {"Daily_TN_load_kg_d", TS_daily_TN_load},

        };
        result_TSS_map = {
            {"Year", TS_date_year},
            {"Month", TS_date_month},
            {"Day", TS_date_day},
            {"Discharge", TS_datarecord},
            {"p_cumulative", TS_ECDF},
            {"Delivery_ratio", TS_Delivery_ratio},
            {"ActiveRiverLength_m", TS_ActiveRiverLength_m},

            {"Load_TSS_entered_kg_d", TS_Load_TSS_entered},
            {"Load_TSS_delivered_kg_d", TS_Load_TSS_delivered},

            {"Daily_TSS_load_kg_d", TS_daily_TSS_load}
        };
    }
    
    uniqueyear = TS_date_year;
    sort(uniqueyear.begin(), uniqueyear.end());
    uniqueyear.erase(unique(uniqueyear.begin(), uniqueyear.end()), uniqueyear.end());

    //WaterQuality_Output::writeLoad_BufferDynamic function will process total phosphorus TP values
    WaterQuality_Output::writeLoad_BufferDynamic(AverageYearLoad_TP, uniqueyear, input, "TP");
    //WaterQuality_Output::writeLoad_BufferDynamic function will process total nitrogen TN values
    WaterQuality_Output::writeLoad_BufferDynamic(AverageYearLoad_TN, uniqueyear, input, "TN");
    //WaterQuality_Output::writeLoad_BufferDynamic function will process total suspended solids TSS values
    WaterQuality_Output::writeLoad_BufferDynamic(AverageYearLoad_TSS, uniqueyear, input, "TSS");

    if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
        //writing extened output
        WaterQuality_Output::writeLoad_BufferDynamic_extended(result_TP_map, input, "TP");
        WaterQuality_Output::writeLoad_BufferDynamic_extended(result_TN_map, input, "TN");
        WaterQuality_Output::writeLoad_BufferDynamic_extended(result_TSS_map, input, "TSS");
    }
    cout << "Writing output to: " << input->SimulationStringParams["OutputFolder_Path"] << endl;
    cout << "===" << endl;
}