#ifndef StreamTemperature_H
#define StreamTemperature_H
/*
* HydroPlus is source code and models developed and managed by Theodore Endreny and his students at SUNY ESF, te@esf.edu
* The US Forest Service and Davey Tree have provided strategic software development support through funding and personnel
* Attribution for use of software and code is requested. Please see Executive.cpp for documentation.
*/
#include <iostream>
#include <fstream>
#include <vector>
//_USE_MATH_DEFINES gives access to Pi as M_PI, and must precede #include <cmath> or <math.h>
#define _USE_MATH_DEFINES
#include <cmath>
#include <math.h>
#include <stdio.h>
#include <algorithm>	// min_element, max_element
#include <numeric>		// accumulate
#include <string>
#include <time.h>		// UNSTEADY STATE SITUATION
#include <stdlib.h>
#include "SolarCalculation.h"

//Constant density of water (kg/m^3), approximates value for 20C
#define rho_water  1000
//Constant specific heat capacity of water (J/ kg / C)
#define c_water 4182
//Constant NODATA
#define NODATA -99999

///////////////////////////////////////////////////////////////
//      StreamTemperature class declarations                 //
///////////////////////////////////////////////////////////////



class StreamTemperature
{
public:

	double sectionNr = 0.0;       // Number of the cross section in the reach
	double eastBank = 0.0;        // The height of bank in the eastern side
	double eastTree = 0.0;        // The height of tree in the eastern side
	double eastDistance = 0.0;    // The distance between bank and Canopy in the eastern side
	double westBank = 0.0;        // The height of bank in the western side
	double westTree = 0.0;        // The height of tree in the eastern side
	double westDistance = 0.0;    // The distance between bank and Canopy in the western side

	double dirrectSW = 0.0;        // Direct shortwave radiation
	double diffuseSW = 0.0;        // Diffuse shortwave radiation
	double sectionWidth = 0.0;     // The width of the stream in the desired location
	double downRadLW = 0.0;        // longwave atmospheric radiation to the down
	double upRadLW = 0.0;          // longwave atmospheric radiation to up, we need the net longwave radiation
	double windSpeed = 0.0;        // Wind speed in desired time step based on initial inputs of i-Tree Air
	double sectiondx = 0.0;		 // The length of the river reach
	double waterDepth = 0.0;	     // The depth of the water
	double secArea = 0.0;			 // cross sectional area (m^2)

	double avgVel = 0.0;			 // average velocity of the section


	double VaporPressure_Saturated_kPa = 0.0;
	double sedExchange = 0.0;
	double cl = 0.0;		//Cloudiness
	double humidity = 0.0;

	double netDiffuseE = 0.0;       // Net solar radiation (I'diffuse) in (Eq. #8 in Sun et al., 2015) for Eastern side
	double netDiffuseW = 0.0;       // Net solar radiation (I'diffuse) in (Eq. #8 in Sun et al., 2015) for Western side
	double adjustedNetRadE = 0.0;   // Adjusted direct radiation by topographic shading (I"direct) for Eastern side
	double adjustedNetRadW = 0.0;   // Adjusted direct radiation by topographic shading (I"direct) for Western side

	double effShadowLengthW = 0.0;      // The shadow length as the initial value for shadow (Eq. #6 in Sun et al., 2015) for Western side

	double sedimentHeat = 0.0;        // Heat exchanged between stream and sediment bed (kcal�m-2�s-1)
	double nexThermalEnergyE = 0.0;   // net exchange of thermal energy (Joules�m-1�s-1) for Eastern side
	double nexThermalEnergyW = 0.0;  // net exchange of thermal energy (Joules�m-1�s-1) for Western side

	vector<double> qEast;			// added in May 17
	vector<double> qWest;

	double sedimentTemp = 0.0;       // sediment temperature (C)
	double airTemp_C = 0.0;            // air temperature (C)
	double dewTemp = 0.0;	           // dew point temperature , the input file is in F but should be change to c
	double surfaceTemp = 0.0;         // cross sectionally averaged surface temperature (C)
	double deltaT = 0.0;

	double solarAzimuth = 0.0;
	double solarAltitude = 0.0;
	double streamAzimuth = 0.0;
	double CanopyLAI = 0.0;

	double solarZenith = 0.0;
	double SWrefl = 0.0;
	double depthMeasBed = 0.0;
	int TotalTime = 0;
	int TotalDist = 0;
	int SolvingMethod = 0;				// The methodology for solving the SW radiation flux, 1 is based on sky view factor and 2 is based on fixed values like 0.25 and 0.75 as the SF and sky view
	int SedimentMethod = 0;				// The methodology selection for sediment flux for explicit finite difference method
	vector<vector<double>> solarSW_Dir_LI2D, solarSW_Diff_LI2D, aitTemp_LI2D, rela_humidity_LI2D, windSp_LI2D, cloudiness_LI2D;
	vector<double> hourlyTime, SWreflVec, hourlyTime_unSteady,hourlyTime_min;
	vector<double> adjustedDirSW_E, adjustedDirSW_W, adjustedDiffSW_E, adjustedDiffSW_W, atmLW, lcLW, backLW, latent, sensible, bed;
	vector<double> BankH_E_LI, TreeH_E_LI, BankH_W_LI, TreeH_W_LI, BuildH_E_LI, BuildH_W_LI, CanDist_E_LI,
		CanDist_W_LI, BanDist_E_LI, BanDist_W_LI, BuildDist_E_LI, BuildDist_W_LI, Buff_E_LI, Buff_W_LI, elev_LI,
		StrAzimuth_LI;
	vector<double> colNum_Dist, rowNum_Dist, LongValue_Dist, latValue_Dist;
	vector<double> SolAzimuth_LI, SolAltitude_LI;
	vector<double> SF_LI;
	vector<double> SkyToView_LI;
	vector<double> secWidth_LI, secDepth_LI, secSlope_LI, secArea_LI, Z_LI;
	vector<vector<double>> heat_flux_2D, SWDir_2D, SWDiff_2D, SW_2D, LW_2D, atmLW_2D, LCLW_2D, backLW_2D, Latent_2D, Sensible_2D, Bed_2D;
	vector<vector<double>> T;  // T: temperature array


	vector<double> fluxxx, tNew;
	vector<double> SkyV_E, SkyV_W, VSAEast, TSAEast, BSAEast, VSAWest, TSAWest, BSAWest;

	// ANALYSING THE ENTERIOR RESULTS AND VECTOR - REZA
	vector<vector<double>> a_c, b_c, c_c, o, p, q, o_c, p_c, q_c, d, g, k, m, a, b, c;
	vector<vector<double>> depth_total, width_total, area_total, wp_total, ksed1;
	vector<vector<double>> t_minBased;
	//added vectors for minutely outputs
	vector<vector<double>> heat_flux_2D_min, SWDir_2D_min, SWDiff_2D_min, SW_2D_min, LW_2D_min, atmLW_2D_min, LCLW_2D_min, backLW_2D_min, Latent_2D_min, Sensible_2D_min, Bed_2D_min;
	vector<vector<double>> delta_2D_min, depthAvgALL_2D_min, AreaALL_2D_min, VelocityALL_2D_min, WidthALL_2D_min, RetentionT_2D_min, Pw_ALL_2D_min, WSElev_ALL_2D_min, Q_Str_2D_min, HyporheicExchange_2D_min, Dispersion_2D_min,
		FluxConduction_2D_min, depth_Uns_2D_min, Pw_unsteady_2D_min, width_UnsTotal_2D_min, depth_AvgUns_2D_min, Hyporheic_Exchange_2D_min, Rh_unsteady_2D_min;

	vector<vector<double>> Q_half_min1, volume1, discharge_time1;


	vector<vector<double>> A;
	vector<double> x, b_copy;

	// HEC-RAS method
	double dx_Interval_m, dtime_Interval_min;
	vector<double> HyporheicVelocity_HR, DeltaT_HR, calcT_Temp;
	vector<double> Q_HR_LI, MinChElev_HR_LI, WSElev_HR_LI, VelChnl_HR_LI, FlowArea_HR_LI, TopWidth_HR_LI, WSSlope_HR_LI, WP_HR_LI;
	vector<vector<double>> Q_Str_HR, WidthALL_HR, DepthALL_HR, DepthAvgALL_HR, VelocityALL_HR,
		AreaALL_HR, HeatFluxRetentionTime_p_dt, Pw_ALL_HR, WSElev_ALL_HR, HyporheicExchange_HR, Dispersion_HR, FluxConduction_HR, bedTemp_HR;
	vector<vector<double>> Q_Inflow_HR, T_Inflow_HR;
	vector<vector<double>> Q_In_HR, T_In_HR, DeltaTotal_HR;
	vector<vector<double>> Flux_Evaporation_HR, Flux_Conduction_HR, Flux_Convection_HR, FLux_AtmLW_HR, Flux_LCLW_HR,
		Flux_BackLW_HR, Flux_LW_HR, Flux_DirSW_HR, Flux_DiffSW_HR, Flux_SW_HR, Flux_Total_HRdist, Temp_HRdist;
	vector<vector<double>> Flux_Evaporation_Nodes, Flux_Conduction_Nodes, Flux_Convection_Nodes, FLux_AtmLW_Nodes, Flux_LCLW_Nodes,
		Flux_BackLW_Nodes, Flux_LW_Nodes, Flux_DirSW_Nodes, Flux_DiffSW_Nodes, Flux_SW_Nodes, FluxTotal_HR, HeatFluxADR_Temperature_C;
	vector<double> Hyporheic_Exchange_HR, Pw_HR_1D, Hor_Con_HR, Hor_Con_FD, ho_HR, ho_FD, hl_HR, hl_FD;
	float average_PW_HR, averageHor_Con_HR, averageHor_ho_HR, averageHor_hl_HR;
	int time_HR, Node_HR;
	double Flux_HEF;
	int theNode_HR, theDistance_HR;
	int calculation;  // 1 = steady; 2 = unsteady; 3 = HEC-RAS

	// UNSTEADY STATE VARIABLES =============================================== (Aug. 17)
	int TotalTime_min, TotalDist_stable, TotalTime_Uns_FD, TotalTime_Uns, TotalTime_Uns_CN;
	int time_Uns, Node_Uns;
	int flag_BC;
	vector<double> delta_t, S1, S2;
	vector<double> SkyToView_UnsE, SkyToView_UnsW, SF_LI_Uns,
		SWreflVec_Uns, SolAzimuth_LI_Uns, SolAltitude_LI_Uns;
	vector<int> Flag_SkipNode;
	vector<vector<double>> Q, Temp_unsteady, depth_unsteady, AreaX, velocity_unsteady, inflow_Rate, inflow_Temp_Uns;
	vector<double> Q_minBased, inflow_Q_2_min, inflow_Q_3_min, Tobs_minBased, inflow_temp_2_temp, inflow_temp_3_temp,
		width_unsteady, Pw_unsteady, Rh_unsteady, Z_unsteady, n_unsteady, s_unsteady, width_UnsTotal, discharge_Unsteady, obs_xo_LI;
	vector<double> Q_GW, Temp_GW, depth_AvgUns;
	vector<vector<double>> Q_In, T_In, T_trib_shift;
	vector<double> conductivity_Uns, particleSize_Uns, embeddedness_Uns, Hyporheic_Exchange, BedTn_LI_Uns, hyp_Slope_Uns;
	vector<double> dispersion, T_unsTemporary;
	double Froude_Densiometric;
	double evap_Rate = 0;
	double DirSW_UnsE, DirSW_UnsW, DiffSW_UnsE, DiffSW_UnsW, LW_Atm_Uns, LW_LC_Uns, LW_Back_Uns, LW_Uns,
		latent_Uns, DirSW_Uns_temp, DiffSW_Uns_temp, SW_Uns, Sensible_Uns, bed_Uns;
	vector<vector<double>> Flux_Total, Flux_Conduction, Flux_Evaporation, Flux_Convection, FLux_AtmLW, Flux_LCLW, Flux_BackLW, Flux_LW, Flux_DirSW,
		Flux_DiffSW, Flux_SW, T_Uns, depth_Uns, Porosity, Pw_unsteady1, width_UnsTotal1, depth_AvgUns1, Hyporheic_Exchange1,Rh_unsteady1;

	vector<vector<double>> hour_Flux_Total, hour_Flux_Conduction, hour_Flux_Evaporation, hour_Flux_Convection, hour_FLux_AtmLW,
		hour_Flux_LCLW, hour_Flux_BackLW, hour_Flux_LW, hour_Flux_DirSW,
		hour_Flux_DiffSW, hour_Flux_SW, hour_T_Uns, hour_depth_Uns, hour_Pw_unsteady, hour_width_UnsTotal, hour_depth_AvgUns,
		hour_Hyporheic_Exchange, hour_Rh_unsteady;

	vector<vector<double>> inflow_Discharge;
	vector<vector<double>> inflow_Temp;
	vector<string> BedSedType;
	vector<double> convertedCode;
	vector<double> DistBedMeas, BedDepthMeas, BedT0, BedTn, K_sedType_LI;

	double slope= 0.0;
	double aspect= 0.0;

	double SenHeat= 0.0;
	double SenHeatRes= 0.0;
	double skyViewGeneral= 0.0;
	double skyViewFacW= 0.0;
	double skyViewFacE= 0.0;

	double treeP= 0.0;
	double soilP= 0.0;
	double shortVegP= 0.0;
	double imperP= 0.0;
	double waterP= 0.0;

	vector<double> totAirTemp;  // A vector for accessing the generated air temperature for the targeted cell in Cell.h and TemperaturePrime class

	Inputs* input;

	SolarCalculation* solarCal;

	// For calculating the average deltaT
	//StreamTemperature(Inputs* input, SolarCalculation* solarCal) to obtain pointers
	StreamTemperature(Inputs* input, SolarCalculation* solarCal);
	StreamTemperature() {}
	//~StreamTemperature() {}
	void clear();

	// For calling the input files from the shading class and the required parameters which is already exist.
	//void initialVariable(int ts, int id);	
	// Calculating the sky view factor adjusted by riparian vegetation buffer considering Eastern info
	// and then calculating the "net solar radiation" using Eastern data (Eq. #8 in Sun et al., 2015)
	void netSolarDiffRadE(const vector<double>&, int);

	// Calculating the sky view factor adjusted by riparian vegetation buffer considering Western info
	// and then calculating the "net solar radiation" using Western data (Eq. #8 in Sun et al., 2015)
	void netSolarDiffRadW(const vector<double>&, int);

	/* Here using the equation in Sun et al., (2015), I am going to calculate the "effective shading width (Xeff)"
	   This function covers equations 1, 3, 4, 5 and 6 and also table 1 of the mentioned paper.
	   Like the previous set of functions, I am going to define two similar functions for east and west side.*/
	double effectiveShadingE(const vector<double>&, int, int);
	double effectiveShadingW(const vector<double>&, int, int);


	/* Here, using two following function, I am going to calculate the adjusted direct radiation by topographic shading.
	   This term (I"direct) is represented in equation #07 in Sun et al., (2015) */
	   // modified by Reza for testing the Dr. Lautz's work: HFlux
	void adjustedDirRadE(const vector<double>&, int);
	void adjustedDirRadW(const vector<double>&, int);

	void Thomas_Numerical_Solver(const vector<double>& a, vector<double> b, const vector<double>& c, const vector<double>& rhs, vector<double>& x);
	/* For the stream segments of the heat exchange between stream and sediment bed (qsed)
	Chapra et al. (2008) - QUAL2K model  */
	//double streamBedTransfer(int);

	/* Calculating the surface temperature. Surface temperature is a function of the surface sensible heat and resistance, as well as
	the air temperature, specific heat, and density. So we can calculate the Surface temperature based on the mentioned terms.
	According to the Yang et al., (2013) --> Eqn. 3, 5 and 7 we can argue that:
	Ts = (H*rH + airDensity * Specific Heat of air) / airDensity * Specific Heat of air */
	//double surTemperature(int, double, double);

	/* Sky view factor:
	 calculating the sky view factor above the riparian vegetation
	 for Eastern bank
	 Reference: DHSVM model, CanopyShading.c CPP Line #197 */
	 //double CalcCanopySkyViewE();

	 // Sky view factor: for Western bank
	 //double CalcCanopySkyViewW();

	 /* The net transfer of thermal energy (qw) due to shortwave solar radiation, longwave atmospheric radiation,
	 blackbody radiation from the water surface, latent and sensible heat conduction and sediment heat flux
	 (e.g., Chen et al, 1998; Chapra et al, 2008; Yearsley, 2009)*/
	 //void netExchangeE(int, double, double);    //East side
	 //void netExchangeW(int);    // West side

	 // The last equation for calculating the cross-sectionally averaged water temperature difference (dt) (degree C)
	void avgTemperature();
	void solarReflectivity();
	void longWaveRadiation(const vector<double>&, const vector<double>&, const vector<double>&, const vector<double>&, int);
	void latentHeat(const vector<double>&, const vector<double>&, const vector<double>&, const vector<double>&, int);
	void sensibleHeat(const vector<double>&, const vector<double>&, const vector<double>&, const vector<double>&, int);
	void bedSediment(const vector<vector<double>>&, const vector<vector<double>>&, const vector<double>&, const vector<double>&,
		const vector<double>&, const vector<double>&, int);
	void bedSedimentHS(const vector<double>&, int, int);;
	void minToHour(int, int, vector<int>&);
	double interpolate(vector<double>&, vector<double>&, double, bool);
	vector<double> minute(vector<vector<double>>&, int);
	vector<double> gauss(vector<vector<double>>&);
	double maxValue(double, double, double);


	vector<double> interp1(vector<double>&, vector<double>&, vector<double>&);	//Nearest interpolation

	// UNSTEADY STATE FUNCTIONS =============================================== (Aug. 17)
	void hydroStability();
	void initialConditions(int);
	void hydraulics(int, int);
	void McCormick_Temp1(const vector<double>&, int);
	void McCormick_Temp2(const vector<double>&, int);
	void mixing(int, int, int);
	void adjustedDirRadE_Uns(const vector<double>&, int, int);
	void adjustedDirRadW_Uns(const vector<double>&, int, int);
	double effectiveShadingE_Uns(const vector<double>&, int, int);
	double effectiveShadingW_Uns(const vector<double>&, int, int);
	void netSolarDiffRadE_Uns(const vector<double>&, int, int);
	void netSolarDiffRadW_Uns(const vector<double>&, int, int);
	void longWaveRadiation_Uns(const vector<double>&, const vector<double>&, const vector<double>&, int, int);
	void latentHeat_Uns(const vector<double>&, const vector<double>&, const vector<double>&, int, int);
	void sensibleHeat_Uns(const vector<double>&, const vector<double>&, const vector<double>&, int, int);
	void bedSediment_Uns(const vector<double>&, const vector<double>&, const vector<double>&, int, int);

	void skyViewE();
	void skyViewW();

	void ConductionFlux_HR(vector<vector<double>>&, const vector<double>&, int, int, int);

	int getTotalTime() { return TotalTime; }
	double getTotalDist() { return TotalDist; }
	vector<vector<double>> getModT() { return T; }
	vector<vector<double>> getSWDir_2D() { return SWDir_2D; }
	vector<vector<double>> getSWDiff_2D() { return SWDiff_2D; }
	vector<vector<double>> getAtmLW() { return atmLW_2D; }
	vector<vector<double>> getLcLW() { return LCLW_2D; }
	vector<vector<double>> getBackLW() { return backLW_2D; }
	vector<vector<double>> getLatent() { return Latent_2D; }
	vector<vector<double>> getSensible() { return Sensible_2D; }
	vector<vector<double>> getBed() { return Bed_2D; }
	vector<vector<double>> getSW() { return SW_2D; }
	vector<vector<double>> getLW() { return LW_2D; }
	vector<vector<double>> getHeatFlux() { return heat_flux_2D; }
	vector<vector<double>> getT_MinBased() { return t_minBased; }

	vector<vector<double>> get_TotalFlux() { return Flux_Total; }
	vector<vector<double>> get_Conduction() { return Flux_Conduction; }
	vector<vector<double>> get_Evaporation() { return Flux_Evaporation; }
	vector<vector<double>> get_Convection() { return Flux_Convection; }
	vector<vector<double>> get_AtmLW() { return FLux_AtmLW; }
	vector<vector<double>> get_LCLW() { return Flux_LCLW; }
	vector<vector<double>> get_BackLW() { return Flux_BackLW; }
	vector<vector<double>> get_Flux_LW() { return Flux_LW; }
	vector<vector<double>> get_Flux_DirSW() { return Flux_DirSW; }
	vector<vector<double>> get_Flux_DiffSW() { return Flux_DiffSW; }
	vector<vector<double>> get_Flux_SW() { return Flux_SW; }
	vector<vector<double>> get_T_Uns() { return T_Uns; }
	vector<vector<double>> get_depth_Uns() { return depth_Uns; }

	vector<vector<double>> get_hour_TotalFlux() { return hour_Flux_Total; }
	vector<vector<double>> get_hour_Conduction() { return hour_Flux_Conduction; }
	vector<vector<double>> get_hour_Evaporation() { return hour_Flux_Evaporation; }
	vector<vector<double>> get_hour_Convection() { return hour_Flux_Convection; }
	vector<vector<double>> get_hour_AtmLW() { return hour_FLux_AtmLW; }
	vector<vector<double>> get_hour_LCLW() { return hour_Flux_LCLW; }
	vector<vector<double>> get_hour_BackLW() { return hour_Flux_BackLW; }
	vector<vector<double>> get_hour_Flux_LW() { return hour_Flux_LW; }
	vector<vector<double>> get_hour_Flux_DirSW() { return hour_Flux_DirSW; }
	vector<vector<double>> get_hour_Flux_DiffSW() { return hour_Flux_DiffSW; }
	vector<vector<double>> get_hour_Flux_SW() { return hour_Flux_SW; }
	vector<vector<double>> get_hour_T_Uns() { return hour_T_Uns; }
	vector<vector<double>> get_hour_depth_Uns() { return hour_depth_Uns; }

	vector<double> getHypEx() { return Hyporheic_Exchange; }
	vector<double> getQ_GW() { return Q_GW; }

	double getSkyViewFactorW(int i) { return SkyV_W[i]; }
	double getSkyViewFactorE(int i) { return SkyV_E[i]; }


	double getSenHeat(int) { return SenHeat; }
	double getSenHeatRes(int) { return SenHeatRes; }
	double getAirTemp(int) { return airTemp_C; }
	double getSurfaceTemp(int) { return surfaceTemp; }
	double getSedimentTransprt(int) { return sedimentHeat; }
	double getSkyView(int i) { return SkyV_E[i]; }

	double getNexThermalEnergyE(int) { return nexThermalEnergyE; }
	double getNexThermalEnergyW(int) { return nexThermalEnergyW; }


	double getDeltaTemp(int) { return deltaT; }


	// get functions for writing the HEC-RAS related 2d vectors  
	vector<vector<double>> get_EvapHRdist() { return Flux_Evaporation_HR; }
	vector<vector<double>> get_CondHRdist() { return Flux_Conduction_HR; }
	vector<vector<double>> get_ConvHRdist() { return Flux_Convection_HR; }
	vector<vector<double>> get_AtmLWHRdist() { return FLux_AtmLW_HR; }
	vector<vector<double>> get_LCLWHRdist() { return Flux_LCLW_HR; }
	vector<vector<double>> get_BackLWHRdist() { return Flux_BackLW_HR; }
	vector<vector<double>> get_LWHRdist() { return Flux_LW_HR; }
	vector<vector<double>> get_DirSWHRdist() { return Flux_DirSW_HR; }
	vector<vector<double>> get_DiffSWHRdist() { return Flux_DiffSW_HR; }
	vector<vector<double>> get_SWHRdist() { return Flux_SW_HR; }
	vector<vector<double>> get_TotalFluxHRdist() { return Flux_Total_HRdist; }
	vector<vector<double>> get_TempHRdist() { return Temp_HRdist; }

	vector<vector<double>> get_EvapHRnodes() { return Flux_Evaporation_Nodes; }
	vector<vector<double>> get_CondHRnodes() { return Flux_Conduction_Nodes; }
	vector<vector<double>> get_ConvHRnodes() { return Flux_Convection_Nodes; }
	vector<vector<double>> get_get_AtmLWHRnodes() { return FLux_AtmLW_Nodes; }
	vector<vector<double>> get_LCLWHRnodes() { return Flux_LCLW_Nodes; }
	vector<vector<double>> get_BackLWHRnodes() { return Flux_BackLW_Nodes; }
	vector<vector<double>> get_LWHRnodes() { return Flux_LW_Nodes; }
	vector<vector<double>> get_DirSWHRnodes() { return Flux_DirSW_Nodes; }
	vector<vector<double>> get_DiffSWHRnodes() { return Flux_DiffSW_Nodes; }
	vector<vector<double>> get_SWHRnodes() { return Flux_SW_Nodes; }
	vector<vector<double>> get_TotalFluxHRnodes() { return FluxTotal_HR; }
	vector<vector<double>> get_TempHRnodes() { return HeatFluxADR_Temperature_C; }
};

#endif