solar_position_michalsky.h

Go to the documentation of this file.

#ifndef SOLAR_POSITION_MICHALSKY_H
#define SOLAR_POSITION_MICHALSKY_H
 
#include <cmath>                               // for floor, fmod
#include "../framework/degree_trigonometry.h"  // for atan2_deg, cos_deg, sin_deg, acos_deg
#include "../framework/module.h"
#include "../framework/state_map.h"
 
namespace standardBML
{
class solar_position_michalsky : public direct_module
{
   public:
    solar_position_michalsky(
        state_map const& input_quantities,
        state_map* output_quantities)
        : direct_module{},
 
          // Get references to input quantities
          lat{get_input(input_quantities, "lat")},
          longitude{get_input(input_quantities, "longitude")},
          fractional_doy{get_input(input_quantities, "fractional_doy")},
          time_zone_offset{get_input(input_quantities, "time_zone_offset")},
          year{get_input(input_quantities, "year")},
 
          // Get pointers to output quantities
          cosine_zenith_angle_op{get_op(output_quantities, "cosine_zenith_angle")},
          julian_date_op{get_op(output_quantities, "julian_date")},
          solar_L_op{get_op(output_quantities, "solar_L")},
          solar_g_op{get_op(output_quantities, "solar_g")},
          solar_ell_op{get_op(output_quantities, "solar_ell")},
          solar_ep_op{get_op(output_quantities, "solar_ep")},
          solar_ra_op{get_op(output_quantities, "solar_ra")},
          solar_dec_op{get_op(output_quantities, "solar_dec")},
          gmst_op{get_op(output_quantities, "gmst")},
          lmst_op{get_op(output_quantities, "lmst")},
          lha_op{get_op(output_quantities, "lha")},
          solar_zenith_angle_op{get_op(output_quantities, "solar_zenith_angle")},
          solar_azimuth_angle_op{get_op(output_quantities, "solar_azimuth_angle")}
    {
    }
    static string_vector get_inputs();
    static string_vector get_outputs();
    static std::string get_name() { return "solar_position_michalsky"; }
 
   private:
    // References to input quantities
    double const& lat;
    double const& longitude;
    double const& fractional_doy;
    double const& time_zone_offset;
    double const& year;
 
    // Pointers to output quantities
    double* cosine_zenith_angle_op;
    double* julian_date_op;
    double* solar_L_op;
    double* solar_g_op;
    double* solar_ell_op;
    double* solar_ep_op;
    double* solar_ra_op;
    double* solar_dec_op;
    double* gmst_op;
    double* lmst_op;
    double* lha_op;
    double* solar_zenith_angle_op;
    double* solar_azimuth_angle_op;
 
    // Main operation
    void do_operation() const;
};
 
string_vector solar_position_michalsky::get_inputs()
{
    return {
        "lat",               // degrees (North is positive)
        "longitude",         // degrees (East is positive)
        "fractional_doy",    // days
        "time_zone_offset",  // the offset of the time zone relative to UTC
        "year"               // a year between 1950 and 2050
    };
}
 
string_vector solar_position_michalsky::get_outputs()
{
    return {
        "cosine_zenith_angle",  // dimensionless
        "julian_date",          // days
        "solar_L",              // degrees
        "solar_g",              // degrees
        "solar_ell",            // degrees
        "solar_ep",             // degrees
        "solar_ra",             // degrees
        "solar_dec",            // degrees
        "gmst",                 // hr
        "lmst",                 // hr
        "lha",                  // degrees
        "solar_zenith_angle",   // degrees
        "solar_azimuth_angle"   // degrees
    };
}
 
void solar_position_michalsky::do_operation() const
{
    // Define some constants
    double constexpr deg_to_hr = 1.0 / 15.0;
    double constexpr hr_to_deg = 15.0;
    double constexpr hr_per_day = 24.0;
    double constexpr deg_per_rev = 360.0;
 
    double constexpr jd_ref_1948 = 2432916.5;  // Julian date at midnight on 31 December 1948 (UTC)
    double constexpr jd_ref_2000 = 2451545.0;  // Julian date at noon on 1 January 2000 (UTC)
 
    // Unpack the doy and hour in UTC
    double const fractional_doy_utc = fractional_doy - time_zone_offset / hr_per_day;  // days
    double const doy_utc = std::floor(fractional_doy_utc / hr_per_day);                // days
    double const hour_utc = hr_per_day * (fractional_doy_utc - doy_utc);               // hr
 
    // Calculate the Julian date
    double const delta = year - 1949.0;
    double const leap = floor(0.25 * delta);
    double const jd = jd_ref_1948 + delta * 365.0 + leap + doy_utc + hour_utc / hr_per_day;  // days
 
    // Calculate the "n" date
    double const n = jd - jd_ref_2000;  // days
 
    // Calculate the ecliptic coordinates of the sun
    double const L = fmod(280.460 + 0.9856474 * n, deg_per_rev);                            // degrees
    double const g = fmod(357.528 + 0.9856003 * n, deg_per_rev);                            // degrees
    double const ell = fmod(L + 1.915 * sin_deg(g) + 0.020 * sin_deg(2 * g), deg_per_rev);  // degrees
    double const ep = 23.439 - 0.0000004 * n;                                               // degrees
 
    // Calculate the equatorial celestial coordinates of the sun
    double const ra = atan2_deg(cos_deg(ep) * sin_deg(ell), cos_deg(ell));  // degrees
    double const dec = asin_deg(sin_deg(ep) * sin_deg(ell));                // degrees
 
    // Calculate the sidereal time
    double const gmst = fmod(6.697375 + 0.0657098242 * n + hour_utc, hr_per_day);  // hours
    double const lmst = fmod(gmst + longitude * deg_to_hr, hr_per_day);            // hours
 
    // Convert to local coordinates
    double const lha = fmod(lmst * hr_to_deg - ra, deg_per_rev);  // degrees
 
    double const zen = acos_deg(sin_deg(dec) * sin_deg(lat) +
                                cos_deg(dec) * cos_deg(lat) * cos_deg(lha));  // degrees
 
    double az = asin_deg(-cos_deg(dec) * sin_deg(lha) / cos_deg(90.0 - zen));  // degrees
 
    // Make sure azimuth is in the correct quadrant
    double const el = 90.0 - zen;                                      // degrees
    double const el_critical = asin_deg(sin_deg(dec) / sin_deg(lat));  // degrees
 
    if (el >= el_critical) {
        az = 180.0 - az;  // degrees
    } else if (lha > 0) {
        az = 360.0 + az;  // degrees
    }
 
    // Determine the cosine of the zenith angle
    double const cos_zen = cos_deg(zen);  // dimensionless
 
    // Update the output pointers
    update(cosine_zenith_angle_op, cos_zen);
    update(julian_date_op, jd);
    update(solar_L_op, L);
    update(solar_g_op, g);
    update(solar_ell_op, ell);
    update(solar_ep_op, ep);
    update(solar_ra_op, ra);
    update(solar_dec_op, dec);
    update(gmst_op, gmst);
    update(lmst_op, lmst);
    update(lha_op, lha);
    update(solar_zenith_angle_op, zen);
    update(solar_azimuth_angle_op, az);
}
 
}  // namespace standardBML
#endif