#include "SolarCalculator.h"

//References:
//Iqbal, M. (1983). Sun-Earth Astronomical Relationships, Chapter 1. An Introduction to Solar Radiation, Elsevier Science: 408 pages.
//Kumar, L., A. K. Skidmore and E. Knowles (1997). "Modelling topographic variation in solar radiation in a GIS environment." International Journal of Geographical Information Science 11(5): 475-497.
//Reda, I. and A. Andreas (2003). "Solar Position Algorithm for Solar Radiation Applications " National Renewable Energy Laboratory NREL/TP-560-34302. Revised 2008: 55.
//Yang, F., K. Mitchell, Y.-T. Hou, Y. Dai, X. Zeng, Z. Wang and X.-Z. Liang (2008). "Dependence of Land Surface Albedo on Solar Zenith Angle: Observations and Model Parameterization." Journal of Applied Meteorology and Climatology 47(11): 2963-2982.


//SolarCalculator::SolarCalculator defined constructor for use here with input pointer
SolarCalculator::SolarCalculator(Inputs* input) {
    //Initialize variables
    Latitude_NH_pos_dd = 0;
    Longitude_WH_pos_dd = 0;
    DeclinationAngle_Solar_rad = 0;
    Day_Angle_rad = 0;
    Equation_of_Time_HH_point_Mn = 0;
    Standard_Meridian_deg = 0;
    TrueSolarTime_HH_point_Mn = 0;

    HourAngle_Solar_deg = 0;
    HourAngle_Solar_rad = 0;

    Cos_ZenithAngle_Solar = 0;
    ZenithAngle_Solar_rad = 0;

    Sin_AltitudeAngle_Solar = 0;
    AltitudeAngle_Solar_rad = 0;

    AzimuthAngle_Solar_S_0_rad = 0;
    AzimuthAngle_Solar_N_0_rad = 0;
    AzimuthAngle_Solar_N_0_deg = 0;
    Cos_IncidenceAngle_Solar = 0;
}

//SolarCalculator_Inputs function computes solar incidence angle using Reda et al. (2003) NREL SPA Algorithm
//Note: NREL SPA is National Renewable Energy Laboratory Solar Position Algorithm and SolPos Calculator 
//Note: Results of this algorithm agree with NREL SPA calculator and SolPos calculator, https://midcdmz.nrel.gov/solpos/spa.html 
//Note: Incidence angle equations are equivalent between Reda et al., (2003) Eq 47 and Iqbal (1983) Eq 1.6.5b
//Note: Transition to Reda et al. (2003) did not change results obtained from method used by Yang et al. (2013)
//Note: Bracketed RHS term in Eq 25 of Yang et al. (2013) is cosine solar incidence angle, referencing Kumar et al. (1997)
//Note: Cosine incidence angle of Eq 16 in Kumar et al. (1997) uses formulation different than Yang et al. (2013)
void SolarCalculator::SolarCalculator_Inputs(Inputs* input, CompactRagged* beC, WeatherProcessor* WeatherPro, string StationID_string, int MapPixel_ID, int DataFolder_ID, int timeStep, int Flag_simulateReferenceStation) {

    //YYYYMMDD is input->SimulationDate_YYYYMMDD[timeStep]
    int Date_YYYYMMDD = input->SimulationDate_YYYYMMDD[timeStep];
    //Date_JulianDay = input->Date_YYYYMMDD_to_JulianDay(YYYYMMDD)
    int Date_JulianDay = input->Date_YYYYMMDD_to_JulianDay(Date_YYYYMMDD);

    //Day_Angle_rad (radians)
    //Note: Consider comment to provide Hirabayashi and Endreny equation number 
    Day_Angle_rad = (360 * (Date_JulianDay - 1) / 365.0) * Ratio_Radian_to_Degree;
    //Equation_of_Time_HH_point_Mn (decimal minutes), equation of time
    //Note: Consider comment to provide Hirabayashi and Endreny equation number 
    Equation_of_Time_HH_point_Mn = (0.000075 + 0.001868 * cos(Day_Angle_rad) - 0.032077 * sin(Day_Angle_rad)
        - 0.014615 * cos(2 * Day_Angle_rad) - 0.040849 * sin(2 * Day_Angle_rad)) * 229.18;
    //DeclinationAngle_Solar_rad (radians), solar declination angle
    DeclinationAngle_Solar_rad = (0.006918 - 0.399912 * cos(Day_Angle_rad) + 0.070257 * sin(Day_Angle_rad)
        - 0.006758 * cos(2 * Day_Angle_rad) + 0.000907 * sin(2 * Day_Angle_rad)
        - 0.002697 * cos(3 * Day_Angle_rad) + 0.00148 * sin(3 * Day_Angle_rad));

    //If Flag_LatitudeLongitude_Computed is true then enter for obtaining Latitude_Longitude_vectorPair_dd values
    //Note: Flag set in Inputs.cpp if HydroPlusConfig.xml is Albers, or (>= 32601 and <= 32660)
    if (input->Flag_LatitudeLongitude_Computed == true) {
        //If MapPixel_ID < Inputs::nCols * Inputs::nRows then in map
        if (MapPixel_ID < Inputs::nCols * Inputs::nRows) {
            //Note: Conversion between MapPixel_ID and row & column pair, with MapPixel_ID starting at 0, row & col starting at 0
            //Note: Eqs: row=MapPixel_ID / Inputs::nCols; col=MapPixel_ID % Inputs::nCols; MapPixel_ID = (row * nCols) + col
            int row = MapPixel_ID / Inputs::nCols;
            //col = MapPixel_ID % Inputs::nCols, derived from MapPixel_ID, first col indexed to 0
            int col = MapPixel_ID % Inputs::nCols;

            //Latitude_NH_pos_dd (dd) = Latitude_Longitude_vectorPair_dd[row][col].first
            //Note: Latitude_Longitude_vectorPair_dd indexed to 0, row and col indexed to 0
            Latitude_NH_pos_dd = input->Latitude_Longitude_vectorPair_dd[row][col].first;
            //Longitude_WH_pos_dd (dd) = Latitude_Longitude_vectorPair_dd[row][col].second * -1; 
            //Note: Negative sign makes Western Hemisphere +
            //Note: Latitude_Longitude_vectorPair_dd indexed to 0, row and col indexed to 0
            Longitude_WH_pos_dd = input->Latitude_Longitude_vectorPair_dd[row][col].second * -1;
        }
        //Else multiple stations is active and data appended to vector
        else {
            //Latitude_NH_pos_dd = input->Latitude_Longitude_vectorPair_dd[MapPixel_ID][0].first
            //Note: Reference station appended at end: use MapPixel_ID directly
            Latitude_NH_pos_dd = input->Latitude_Longitude_vectorPair_dd[MapPixel_ID][0].first;
            //Longitude_WH_pos_dd = input->Latitude_Longitude_vectorPair_dd[MapPixel_ID][0].second * -1
            //Note: Negative sign makes Western Hemisphere +
            Longitude_WH_pos_dd = input->Latitude_Longitude_vectorPair_dd[MapPixel_ID][0].second * -1;
        }
    }
    //Else set use Radiation.csv Latitude_rad and HourAngle_rad values below and do not use Longitude_WH_pos_dd
    else {
        //Latitude_NH_pos_dd is Latitude_rad * Ratio_Degree_to_Radian
        Latitude_NH_pos_dd = input->Latitude_rad * Ratio_Degree_to_Radian;
        //Longitude_WH_pos_dd set to prime meridian but not used
        Longitude_WH_pos_dd = 0;
    }

    //To handle both Northern/Southern hemisphere, accept negative latitude (for Southern)
    //Latitude_NH_pos_rad (radians), latitude
    double Latitude_NH_pos_rad = Latitude_NH_pos_dd * Ratio_Radian_to_Degree;
    double Longitude_WH_neg_dd = -1 * Longitude_WH_pos_dd;
    //GMT_Offset_int = Longitude_WH_neg_dd / 15, negative values in western hemisphere
    //Note: Earth divided into 24 time zones, each separated into 15 degrees, where 24 * 15 deg = 360 deg
    int GMT_Offset_int = int(Longitude_WH_neg_dd / 15.0);

    //Standard_Meridian_deg (degrees) takes 15 deg for each GMT offset hour
    Standard_Meridian_deg = 15.0 * GMT_Offset_int;

    //HH_MM_SS is input->SimulationTime_HMS[timeStep]
    string HH_MM_SS = input->SimulationTime_HMS[timeStep];
    //Date_Hour_HH extracted from HH_MM_SS; values must be 1 to 24
    int Date_Hour_HH = stoi(HH_MM_SS.substr(0, 2)) + 1;

    //TrueSolarTime_HH_point_Mn (hour with decimal minutes), true solar time
    TrueSolarTime_HH_point_Mn = Date_Hour_HH + (4 * (abs(Standard_Meridian_deg) - abs(Longitude_WH_pos_dd)) + Equation_of_Time_HH_point_Mn) / 60.0;
    //hourAngle (degrees)
    HourAngle_Solar_deg = 15.0 * TrueSolarTime_HH_point_Mn - 180;
    //hourAngle (radians)
    HourAngle_Solar_rad = HourAngle_Solar_deg * Ratio_Radian_to_Degree;

    //If Flag_LatitudeLongitude_Computed is not true then enter to define HourAngle_Solar_rad with Radiation.csv input
    //Note: Flag set in Inputs.cpp if HydroPlusConfig.xml Albers, 3034, or (>= 32601 and <= 32660)
    if (input->Flag_LatitudeLongitude_Computed != true) {
        HourAngle_Solar_rad = input->HourAngle_Solar_rad[timeStep];
    }

    //If Flag_MultipleStations and Flag_simulateReferenceStation true then read in WeatherPro-> station variables for flux calculations
    if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1 && Flag_simulateReferenceStation == 1) {
        //StationID_string used to define strings to define WeatherPro-> variables
        string keynamepart1 = StationID_string;
        string keyname_DeclinationAngle_Solar_rad = keynamepart1 + "_" + "DeclinationAngle_Solar_rad";
        string keyname_Latitude_rad = keynamepart1 + "_" + "Latitude_rad";
        string keyname_HourAngle_Solar_rad = keynamepart1 + "_" + "HourAngle_Solar_rad";
        //DeclinationAngle_Solar_rad (radian) is WeatherPro->RadiationMap[keyname_DeclinationAngle_Solar_rad][timeStep]
        DeclinationAngle_Solar_rad = WeatherPro->RadiationMap[keyname_DeclinationAngle_Solar_rad][timeStep];
        //Latitude_NH_pos_rad (radian) is WeatherPro->LatitudeMap[keyname_Latitude_rad]
        Latitude_NH_pos_rad = WeatherPro->LatitudeMap[keyname_Latitude_rad];
        //HourAngle_Solar_rad (radian) is WeatherPro->RadiationMap[keyname_HourAngle_Solar_rad][timeStep]
        HourAngle_Solar_rad = WeatherPro->RadiationMap[keyname_HourAngle_Solar_rad][timeStep];
    }

    //ElevationAngle_Solar_rad (radians) is asin(sin(Latitude_NH_pos_rad) * sin(DeclinationAngle_Solar_rad) + cos(DeclinationAngle_Solar_rad) * cos(Latitude_NH_pos_rad) * cos(HourAngle_Solar_rad))
    //Note: Based on Eq 41 of Reda et al. (2003) and the NREL SPA
    //Note: Calculation is complement of calculation of Cos_ZenithAngle_Solar from Eq 9-12 Vanek et al (2012)
    double ElevationAngle_Solar_rad = asin(sin(Latitude_NH_pos_rad) * sin(DeclinationAngle_Solar_rad) + cos(DeclinationAngle_Solar_rad) * cos(Latitude_NH_pos_rad) * cos(HourAngle_Solar_rad));

    //ElevationAngle_Solar_deg (degree) is ElevationAngle_Solar_rad * Ratio_Degree_to_Radian
    double ElevationAngle_Solar_deg = ElevationAngle_Solar_rad * Ratio_Degree_to_Radian;

    //Correction_Atmospheric_Refraction_deg initialized to zero
    double Correction_Atmospheric_Refraction_deg = 0.0;
    //AtmosphericPressure_hPa is input->AtmPres_kPa[timeStep] * 10.0
    double AtmosphericPressure_hPa = input->AtmPres_kPa[timeStep] * 10.0;
    //Tair_C (C) is converted from Tair_C (C)
    double Tair_C = input->Tair_C[timeStep];

    //Correction_Atmospheric_Refraction_deg is (AtmosphericPressure_hPa / 1010.0) * (283.0 / (273 + Tair_C)) * 1.02 / (60.0 * tan(ElevationAngle_Solar_rad + 10.3 / (ElevationAngle_Solar_deg + 5.11)))
    //Note: Based on Eq 42 of Reda et al. (2003) and the NREL SPA to account for air mass refraction 
    //Note: tan function directly converts ElevationAngle_Solar_rad, and indirectly converts ElevationAngle_Solar_deg within quotient
    //Note: Eq 42 more generalized than correction in WeatherPrep::CalcSolarZenithAngle and by Liljegren et al. (2008) WBGT code
    Correction_Atmospheric_Refraction_deg = (AtmosphericPressure_hPa / 1010.0) * (283.0 / (273 + Tair_C)) * 1.02 / (60.0 * tan(ElevationAngle_Solar_rad + 10.3 / (ElevationAngle_Solar_deg + 5.11)));

    //ZenithAngle_Solar_rad is (90 - ElevationAngle_Solar_deg) * Ratio_Radian_to_Degree
    //Note: Based on Eq 44 of Reda et al. (2003), which calculated these terms in degrees, conversion to radians in Eq 12 
    ZenithAngle_Solar_rad = (90.0 - ElevationAngle_Solar_deg) * Ratio_Radian_to_Degree;

    //Cos_ZenithAngle_Solar is cos(ZenithAngle_Solar_rad), which takes cosine of ZenithAngle_Solar_rad
    //Note: Cos_ZenithAngle_Solar is used throughout code for adjustment of incoming solar radiation on surface
    Cos_ZenithAngle_Solar = cos(ZenithAngle_Solar_rad);

    //If Cos_ZenithAngle_Solar > 1 then set to 1
    if (Cos_ZenithAngle_Solar > 1) {
        Cos_ZenithAngle_Solar = 1.0;
    }
    //Else If Cos_ZenithAngle_Solar < -1 then set to -1
    else if (Cos_ZenithAngle_Solar < -1) {
        Cos_ZenithAngle_Solar = -1.0;
    }

    //Sin_AltitudeAngle_Solar from Eq 2 of Kumar et al. (1997)
    Sin_AltitudeAngle_Solar = sin(DeclinationAngle_Solar_rad) * sin(Latitude_NH_pos_rad) + cos(DeclinationAngle_Solar_rad) * cos(Latitude_NH_pos_rad) * cos(HourAngle_Solar_rad);
    //AltitudeAngle_Solar_rad is asine of Sin_AltitudeAngle_Solar
    AltitudeAngle_Solar_rad = asin(Sin_AltitudeAngle_Solar);

    //AzimuthAngle_Solar_S_0_rad is atan2(sin(HourAngle_Solar_rad), cos(HourAngle_Solar_rad) * sin(Latitude_NH_pos_rad) - tan(DeclinationAngle_Solar_rad) * cos(Latitude_NH_pos_rad))
    //Note: Based on Eq 45 of Reda et al. (2003) NREL SPA, with conversion to N=0, keeping 360 degree rotation clockwise
    //Note: Eq 45 computes AzimuthAngle_Solar_S_0_rad as South pointing = 0, with 360 degree rotation clockwise
    //Note: This Eq 45 is preferred over Eq 9-13 of Vanek et al (2012) and Eq 3 of Kumar et al. (1997) of PASATH approach
    AzimuthAngle_Solar_S_0_rad = atan2(sin(HourAngle_Solar_rad), cos(HourAngle_Solar_rad) * sin(Latitude_NH_pos_rad) - tan(DeclinationAngle_Solar_rad) * cos(Latitude_NH_pos_rad));

    //AzimuthAngle_Solar_N_0_rad is computed with Eq 46 of Reda et al. (2003), with Pole or N=0, 360 degree rotation clockwise
    AzimuthAngle_Solar_N_0_rad = AzimuthAngle_Solar_S_0_rad + M_PI;
    //AzimuthAngle_Solar_N_0_deg is converted from radians as AzimuthAngle_Solar_N_0_rad * Ratio_Degree_to_Radian 
    AzimuthAngle_Solar_N_0_deg = AzimuthAngle_Solar_N_0_rad * Ratio_Degree_to_Radian;
    //AzimuthAngle_Solar_N_0_deg is AzimuthAngle_Solar_N_0_deg / 360.0, constrained to 360 deg 
    //Note: Based on algorithm Step 3.2.6 of Reda et al. (2003), with Pole or N=0, 360 degree rotation clockwise
    AzimuthAngle_Solar_N_0_deg = AzimuthAngle_Solar_N_0_deg / 360.0;
    //AzimuthAngle_Solar_N_0_deg is 360.0 * (AzimuthAngle_Solar_N_0_deg - floor(AzimuthAngle_Solar_N_0_deg))
    AzimuthAngle_Solar_N_0_deg = 360.0 * (AzimuthAngle_Solar_N_0_deg - floor(AzimuthAngle_Solar_N_0_deg));
    //If AzimuthAngle_Solar_N_0_deg < 0.0 then add 360.0
    if (AzimuthAngle_Solar_N_0_deg < 0.0) {
        AzimuthAngle_Solar_N_0_deg = AzimuthAngle_Solar_N_0_deg + 360.0;
    }

    //End of Algorithm Step 3.2.6 of Reda et al. (2003)
    //AzimuthAngle_Solar_N_0_rad is AzimuthAngle_Solar_N_0_deg * Ratio_Radian_to_Degree
    //Note: AzimuthAngle_Solar_N_0_rad uses N=0, 360 degree clockwise rotation
    AzimuthAngle_Solar_N_0_rad = AzimuthAngle_Solar_N_0_deg * Ratio_Radian_to_Degree;

    //Cos_IncidenceAngle_Solar is cos(ZenithAngle_Solar_rad) * cos(input->SlopeGround_rad[MapPixel_ID]) + sin(ZenithAngle_Solar_rad) * sin(input->SlopeGround_rad[MapPixel_ID]) * cos(AzimuthAngle_Solar_N_0_rad - input->AspectGround_N_0_rad[MapPixel_ID])
    //Note: Based on Eq 47 of Reda et al. (2003), w/ AzimuthAngle_Solar_N_0_rad and AspectGround_N_0_rad, N=0, 360 degree rotation clockwise
    //Note: Eq 47 is equivalent to Iqbal (1983) Eq 1.6.5b, where Iqbal (1983) cited for calculations in WeatherPrep
    Cos_IncidenceAngle_Solar = cos(ZenithAngle_Solar_rad) * cos(input->SlopeGround_rad[MapPixel_ID]) + sin(ZenithAngle_Solar_rad) * sin(input->SlopeGround_rad[MapPixel_ID]) * cos(AzimuthAngle_Solar_N_0_rad - input->AspectGround_N_0_rad[MapPixel_ID]);

    //Write to folder variables for TemperatureOutputWriter 
    //Note: Confirming Equivalence between solar angles of Yang et al (2013) and Reda et al (2003) 
    beC->by_key(MapPixel_ID, DataFolder_ID, "Cos_IncidenceAngle_Solar") = Cos_IncidenceAngle_Solar;
    beC->by_key(MapPixel_ID, DataFolder_ID, "ZenithAngle_Solar_rad") = ZenithAngle_Solar_rad;
    beC->by_key(MapPixel_ID, DataFolder_ID, "ElevationAngle_Solar_rad") = ElevationAngle_Solar_rad;
    beC->by_key(MapPixel_ID, DataFolder_ID, "AzimuthAngle_Solar_N_0_rad") = AzimuthAngle_Solar_N_0_rad;
}

//Ratio_AlbedoDirectCZA_to_AlbedoDirectCZA60 function implements Eq 15 of Yang et al. (2008)
//Note: Yang et al. (2008) explains albedo dependencies were tested for many cover types, but more testing is needed in forests
//Note: Albedo for direct shortwave radiation shown to vary with cosine of solar zenith angle for many cover types
//Note: Ratio_AlbedoDirCZA_to_AlbedoDirCZA60 is for horizontal surfaces; for walls ratio_vertical = 1/Ratio_AlbedoDirCZA_to_AlbedoDirCZA60
//Note: ... ratio_vertical = 1/Ratio_AlbedoDirCZA_to_AlbedoDirCZA60 is equivalent of inverting Eq 15 of Yang et al. (2008)
double SolarCalculator::Ratio_AlbedoDirectCZA_to_AlbedoDirectCZA60(double Cos_ZenithAngle_Solar) {
    //Coeff_Yang_2008_Eq15_a defined as 0.775 in Eq 15 of Yang et al. (2008)
    double Coeff_Yang_2008_Eq15_a = 0.775;
    //Coeff_Yang_2008_Eq15_a defined as 0.775 in Eq 15 of Yang et al. (2008)
    double Coeff_Yang_2008_Eq15_b = 1.55;
    //Ratio_AlbedoDirCZA_to_AlbedoDirCZA60 is (1 + Coeff_Yang_2008_Eq15_a) / (1 + Coeff_Yang_2008_Eq15_b * Cos_ZenithAngle_Solar)
    //Note: Ratio_AlbedoDirCZA_to_AlbedoDirCZA60 is the left hand side of Eq 15 of Yang et al. (2008)
    double Ratio_AlbedoDirCZA_to_AlbedoDirCZA60 = (1 + Coeff_Yang_2008_Eq15_a) / (1 + Coeff_Yang_2008_Eq15_b * Cos_ZenithAngle_Solar);

    //return Ratio_AlbedoDirCZA_to_AlbedoDirCZA60 
    return Ratio_AlbedoDirCZA_to_AlbedoDirCZA60;
}
