#ifndef ResistanceCal_H
#define ResistanceCal_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 <vector>
#include <string>
#include <math.h>
#include "../Inputs/Inputs.h"
#include "../DataFolder.h"

//Latent heat of vaporization for water is 2257000.0 (J/kg)
#define Lamda_J_p_kg 2257000.0
//Specific heat capacity of air at constant pressure is 1006.5 (J/kg/K) at 310 K 
#define SpecificHeat_DryAir_J_p_kg_K 1006.5
//Density of air is 1.205 (kg/m3)
#define Density_Air_kg_p_m3 1.205
//Constant gravitational acceleration is 9.81 (m/s^2)
#define g_GravitationalAcceleration_m_p_s 9.81
//von Karman constant of 0.41 (model sensitive to hundreths decimal precision of term)
#define k_vonKarman  0.41
//Coeff_Drag Su et al. (2001) define Cd = 0.2 as drag coefficient, Cd, of the foliage elements
//Note: Choudhury and Monteith (1988) Eq 20 use Cd = drag coefficient for individual leaves
#define Coeff_Drag 0.2

//Note: Theory discussed with Dr. Christopher Holst of Kalsruhe Institute of Technology (KIT) Alpine campus in GaPa, developer of PALM-4U model
//Note: Dr. Holst and post-doc ChangXing Lan (UW PhD) provided insights May-June 2023 at KIT that helped parameterize this energy balance model
//Note: They assured me that this simple parameterization is as defensible as other approaches given the tremendous uncertainty in the processes
//Note: Yang et al. (2013) below Eq 33 states no generic formula for the roughness length and zero plane displacement for all landscapes
//Note: Experimental roughness lengths range 0.4-0.7 m for dense low buildings or 0.7-1.5 m for regularly built towns [Wieringa, 1993; Masson, 2000], 
//Note: Value approximates one tenth of the building height, or a weighted average landscape roughness element height 
//Note: Advanced users can specify the values of roughness length and zero plane displacement based on their knowledge of their study area.
//Note: PASATH Eq 33: Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m = log((Zu- Zd_m) / Z0m_m)*log((Zu- Zd_m) / Z0m_m) / (k*k*WindSpeed_m_p_s);
//Note: For PASATH model, it hardcoded Eq 33 with parameters: Zu = 10 m, H = 5 m, Zd_m = 4 = 0.8*H, Z0m_m = 0.5 = 0.1*H
//Note: Height_WindSpeedMeasurement_Standard_m & Height_Landscape_ConstantParameter_m remain constant across all pixels to preserve continuity
//Note: Height_Landscape_ConstantParameter_m has large influence on value of Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m and Tair_K
//Note: Height_Landscape_ConstantParameter_m might be adjusted as parameter useful in representing different scenarios and temperature targets

//Height_WindSpeedMeasurement_Standard_m (m) = 10 m by Yang et al. (2013) within PASATH code
//Note: Height_WindSpeedMeasurement_Standard_m (m) for wind profile above urban canopy to compute aerodynamic resitance to higher mesoscale layer
#define Height_WindSpeedMeasurement_Standard_m 10
//Height_Landscape_ConstantParameter_m (m) = 5 m by Yang et al. (2013) within PASATH code
//Note: Height_Landscape_ConstantParameter_m (m) held constant at 5 m for all simulation areas, all scenarios, as parameterized in the PASATH code
#define Height_Landscape_ConstantParameter_m 5

//roughness length (m) for paved impervious, water, or soil by Yang et al. (2013) in PASATH to 0.00137 m, perhaps based on water set to z0 = 0.00137 m Chin (2021) and Thom and Oliver (1977)
//Note: Chin, D. A. (2021). Water Resources Engineering, Fourth Edition. Hoboken, NJ: Pearson Education.
//Note: Thom, A. S., & Oliver, H. R. (1977). On Penman's equation for estimating regional evaporation. Quarterly Journal of the Royal Meteorological Society, 103(436), 345-357. doi:https://doi.org/10.1002/qj.49710343610
#define RoughnessLength_ImpSoil_m 0.00137
//roughness length (m) for classification open, e.g. airport, from Davenport-Wieringa classification, from Table 4.1 in Stull (2000)
#define RoughnessLength_Airport_m 0.03
//RoughnessLength_Water_m = 0.00137 m, based on roughness length (m) for water set to z0 = 0.00137 m Chin (2021) and Thom and Oliver (1977)
//Note: Chin, D. A. (2021). Water Resources Engineering, Fourth Edition. Hoboken, NJ: Pearson Education.
//Note: Thom, A. S., & Oliver, H. R. (1977). On Penman's equation for estimating regional evaporation. Quarterly Journal of the Royal Meteorological Society, 103(436), 345-357. doi:https://doi.org/10.1002/qj.49710343610
#define RoughnessLength_Water_m 0.00137
//RoughnessLength_SVeg_m = 0.03, based on classification open, e.g. grass prarie or farm field, from Davenport-Wieringa classification, from Table 4.1 in Stull (2000)
#define RoughnessLength_SVeg_m 0.03
//RoughnessLength_Tree_m = 0.25 m + 0.75 * Tree_frac, function based on classification rough to closed, e.g. scattered trees to mature forest, from Davenport-Wieringa classification, from Table 4.1 in Stull (2000)
#define RoughnessLength_Tree_m 0.25
//RoughnessLength_NLCD_21_m = 1 m, based on classification closed, e.g. suburban houses, from Davenport-Wieringa classification, from Table 4.1 in Stull (2000)
#define RoughnessLength_NLCD_21_m 1.
//RoughnessLength_NLCD_22_m = 1.8 m, based on 80% between classifications closed and urban center, from Davenport-Wieringa classification, from Table 4.1 in Stull (2000)
#define RoughnessLength_NLCD_22_m 1.8
//RoughnessLength_NLCD_23_m = 2.5 m, based on 0.5 + classification chaotic, e.g. urban center, from Davenport-Wieringa classification, from Table 4.1 in Stull (2000)
#define RoughnessLength_NLCD_23_m 2.5
//RoughnessLength_NLCD_24_m = 4 m, based on 2 + classification chaotic, e.g. urban center, from Davenport-Wieringa classification, from Table 4.1 in Stull (2000)
#define RoughnessLength_NLCD_24_m 4.
//Conversion percent to decimal
#define ratio_percent_to_decimal 100.

//MAX function for finding minimum of x_var and y_var using conditional ternary operator (?:)
//Note: Ternary operator explained: For (x_var<=y_var?x_var:y_var), if x_var<=y_var is true, then x_var is taken, otherwise y_var taken.
#define MAX(x_var,y_var) (((x_var)>=(y_var)) ? (x_var) : (y_var))        

//using namespace std allows for use of standard namespace without the preface 
using namespace std;

//class ResistanceCal
class ResistanceCal
{
private:

	double TreeCover_frac{};
	double ImperviousCover_noTreeCanopy_frac{};
	double PermeablePavementCover_noTreeCanopy_frac{};
	double ShortVegCover_noTreeCanopy_frac{};

	//Aerodynamic resistance above tree
	double Resistance_Aerodynamic_Tree_s_p_m{};
	//Aerodynamic resistance above short vegetation
	double Resistance_Aerodynamic_SVeg_s_p_m{};
	//Aerodynamic resistance above impervious pavement
	double Resistance_Aerodynamic_Impervious_s_p_m{};
	//Aerodynamic resistance above soil
	double Resistance_Aerodynamic_Soil_s_p_m{};
	//Aerodynamic resistance above water
	double Resistance_Aerodynamic_Water_s_p_m{};
	//Aerodynamic resistance above urban canopy, extending to mesoscale air layer
	double Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m{};
	//Surface resistance for impervious 
	double Resistance_Surface_Impervious_s_p_m{};
	//Surface resistance for soil 
	double Resistance_Surface_Soil_s_p_m{};
	//Surface resistance for tree 
	double Resistance_Surface_Tree_s_p_m{};
	//Surface resistance for short vegetation
	double Resistance_Surface_SVeg_s_p_m{};
	//Surface resistance for free water surface 
	double Resistance_Surface_Water_s_p_m{};
	//Aerodynamic conductance for tree
	double Conductance_Aerodynamic_Tree_m_p_s{};
	//Aerodynamic conductance for short vegetation
	double Conductance_Aerodynamic_SVeg_m_p_s{};

	//Coefficient_alpha Choudhury and Monteith (1988) Eq 23 alpha is attenuation coefficient for momentum; Defined as alpha=2 in Table 2; Yang defined ak_attenuation = 2.25
	double Coefficient_alpha{};
	//K_H_soil from Eq 24 of Choudhury and Monteith (1988) as function relating turbulent transfer coefficient, friction velocity, and height; used in Eq 25
	double K_H_soil{};

	//Coefficient_a Choudhury and Monteith (1988) Eq 26 and 29; Defined as a=0.01 s/m below Eq 26; Yang in PASATH defined a_attenuation 0.00662
	double Coefficient_a_spm{};
	//Coefficient_alphaPrime Choudhury and Monteith (1988) Eq 27 and 29 alphaPrime is attenuation coefficient for wind speed; Defined as alphaPrime=3 in Table 2
	//Note: alphaPrime term (denoted alpha) as function of LAI in https://rdrr.io/github/lhmet-forks/bigleaf/man/Gb.Choudhury.html: αlpha = 4.39 - 3.97*exp(-0.258*LAI)
	double Coefficient_alphaPrime{};


	double Soil_SaturationPoint_m3pm3{};
	double Soil_WiltingPoint_m3pm3{};
	double Soil_FieldCapacity_m3pm3{};
	double Coeff_SurfaceResistance_Tree_s_p_m{};
	double Coeff_SurfaceResistance_SVeg_s_p_m{};
	//Richardson number proposed but not calculated, see Liang et al. (1994) on calculation
	double Ri{};
	//Atmopshere stability correction number proposed but not used
	double F{};
	double Height_WindMeasurement_m{};
	double Leaf_Width_Tree_m{};
	double Leaf_Width_SVeg_m{};

	double Height_Water_m{};
	double Height_Soil_m{};
	double Height_SVeg_m{};
	double Height_Tree_m{};
	double Height_Urban21_m{};
	double Height_Urban22_m{};
	double Height_Urban23_m{};
	double Height_Urban24_m{};

public:

	//ResistanceCal(Inputs* input) constructor must match .h ResistanceCal class name
	//Note: Variables defined in the constructor are shared across all void functions, otherwise limited to void function
	ResistanceCal(Inputs* input);
	//~ResistanceCal() {} destructor that ensures resources used by the object can be released 
	~ResistanceCal() {};

	//friend class HeatFluxCal gives it access to the private and protected members of ResistanceCal
	friend class HeatFluxCal;

	//get_F function returns F; could be made constant with const
	double get_F() { return F; }

	//ResistanceCal::AerodynamicResistance_CanopyToMesoLayer function computes aerodynamic resistances (s/m) for between canopy layer and mesoscale layer
	void AerodynamicResistance_CanopyToMesoLayer(Inputs* input, CompactRagged* beC, DataOrganizer* organizer, WeatherProcessor* WeatherPro, int MapPixel_ID, int DataFolder_ID, double WindSpeed_reference_m_p_s, double Height_LandCover_Average_m, double Radiation_Shortwave_Direct_W_p_m2, double Radiation_Shortwave_Diffuse_W_p_m2, int timeStep);
	//ResistanceCal::AerodynamicResistance_OpenArea function computes aerodynamic resistances (s/m) for open area without urban canyons
	void AerodynamicResistance_OpenArea(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, double WindSpeed_reference_m_p_s, double LAI_Tree_m2pm2, double LAI_SVeg_m2pm2,	double Height_LandCover_Average_m, int timeStep);
	//ResistanceCal::AerodynamicResistance_UrbanCanyon function computes aerodynamic resistances (s/m) for urban canyon area, using aspect ratio
	void AerodynamicResistance_UrbanCanyon(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, double WindSpeed_reference_m_p_s, double LAI_Tree_m2pm2, double LAI_SVeg_m2pm2,	double AspectRatio_Canyon_frac, double Height_LandCover_Average_m, double SpecificHeat_HumidAir_JpkgK);
	//ResistanceCal::SurfaceResistance_LandCover function computes surface resistances (s/m) for soil moisture based land cover
	void SurfaceResistance_LandCover(CompactRagged* beC,int MapPixel_ID, int DataFolder_ID, double Saturation_SoilStorage_frac, double LAI_Tree_m2pm2, double LAI_SVeg_m2pm2);
};

#endif