SolarFluxVariableCustomSinusoidal.h

Go to the documentation of this file.

#ifndef SOLAR_FLUX_CUSTOM_SINUSOIDAL_H
#define SOLAR_FLUX_CUSTOM_SINUSOIDAL_H
 
#define VERSION "6"
 
#include <cmath>
#include <iostream>
#include <memory>
#include <numbers>
#include <string>
#include <valarray>
 
#include "../../src/GenericSubmodule.h"
#include "../../src/ModuleFactory.h"
#include "../../src/PhysicalConstants.h"
 
template <typename T>
class SolarFluxVariableCustomSinusoidal : public GenericSubmodule<T>
{
    private:
        static bool m_registered;
        std::string m_ini_filepath{"solar_flux/SolarFluxVariableCustomSinusoidal.ini"};
        T m_period{};                 // Period of the day-night-cycle in s
        T m_heliocentric_distance{};  // Distance between body and sun in AU
        T m_declination{};            // Declination of body in radians
        T m_latitude{};               // Latitude of the point on body in radians
        T m_slope_angle{};            // Angle of the plane against the surface normal in radians
        T m_azimuth_topo_grad{};  // Azimuth of the topographic gradient, where the zero meridian is
                                  // South.
        T m_solar_const{};        // Solar constant
 
        T calculateIncidenceAngle();
 
    public:
        SolarFluxVariableCustomSinusoidal(SimulationClassBase<T>* sim);
        bool setup(std::vector<std::shared_ptr<GenericSubmodule<T>>> all_submodules)
            override;  // loads module into necessary module chains (auto load for calculation
                       // chains)
        bool exec(std::string_view param) override;  // main functionality of the module
 
        bool init() override;
        bool preTimeStep() override;
        bool postTimeStep() override { return true; };
        bool output() override { return true; };
 
        void calculateSolarFlux();
 
        void setFieldPtr(std::shared_ptr<std::valarray<T>> field_ptr) override
        {
            this->module_field = field_ptr;
        }
 
        static std::shared_ptr<GenericSubmodule<T>> createMethode(SimulationClassBase<T>* sim)
        {
            return std::make_shared<SolarFluxVariableCustomSinusoidal<T>>(sim);
        }
        static std::string getName() { return "SolarFluxVariableCustomSinusoidal"; }
        std::string_view getNameLocal() const override
        {
            return "SolarFluxVariableCustomSinusoidal";
        };
        std::vector<std::string> getDependencies() const override
        {
            return this->ini_file_data.getStringVectorParameters("dependencies");
        };
};
 
// ================= Implementation =================
 
template <typename T>
SolarFluxVariableCustomSinusoidal<T>::SolarFluxVariableCustomSinusoidal(SimulationClassBase<T>* sim)
    : GenericSubmodule<T>(sim)
{
    try
    {
        std::string ini_folder_path{
            this->sim->m_simulation_config.getStringParameters("ini_folder_path")};
        this->ini_file_data.loadUserInput(ini_folder_path + m_ini_filepath);
        this->m_generic_submodules = this->ini_file_data.getStringVectorParameters("submodules");
    }
    catch (const BadInput& e)
    {
        std::cerr << e.what() << '\n';
    }
}
 
template <typename T>
bool SolarFluxVariableCustomSinusoidal<T>::setup(
    std::vector<std::shared_ptr<GenericSubmodule<T>>> all_submodules)
{
    m_period = this->ini_file_data.getDoubleParameters("period");
    m_declination = this->ini_file_data.getDoubleParameters("declination") * std::numbers::pi / 180;
    m_latitude = this->ini_file_data.getDoubleParameters("latitude") * std::numbers::pi / 180;
    m_slope_angle = this->ini_file_data.getDoubleParameters("slope_angle") * std::numbers::pi / 180;
    m_azimuth_topo_grad =
        this->ini_file_data.getDoubleParameters("azimuth_topo_grad") * std::numbers::pi / 180;
    m_solar_const = PhysicalConstants::solar_const;
    return true;
}
 
template <typename T>
bool SolarFluxVariableCustomSinusoidal<T>::exec(std::string_view param)
{
    if (param == "InitChain")
    {
        return init();
    }
    if (param == "PreTimeStepChain")
    {
        return preTimeStep();
    }
    if (param == "PostTimeStepChain")
    {
        return postTimeStep();
    }
    if (param == "OutputChain")
    {
        return output();
    }
    return false;
}
 
template <typename T>
T SolarFluxVariableCustomSinusoidal<T>::calculateIncidenceAngle()
{
    T hour_angle{2 * std::numbers::pi * std::fmod(this->sim->elapsed_time / m_period, 1.0)};
    // beta is the angle of 90 minus incidence angle for horizontal surfaces. Equation 2.1 in
    // Schörghofer N. (2025). Planetary Code Collection User Guide.
    T sinbeta{std::cos(m_latitude) * std::cos(m_declination) * std::cos(hour_angle)
              + std::sin(m_latitude) * std::sin(m_declination)};
    T cosbeta{std::sqrt(1 - std::pow(sinbeta, 2))};
    T cos_az_sun{(std::sin(m_declination) - std::sin(m_latitude) * sinbeta)
                 / (std::cos(m_latitude) * cosbeta)};
    // roundoff in both directions
    if (cos_az_sun > 1.0)
    {
        cos_az_sun = 1.0;
    }
    else if (cos_az_sun < -1.0)
    {
        cos_az_sun = -1.0;
    }
    T azSun{std::acos(cos_az_sun)};
    if (std::sin(hour_angle) >= 0)
    {
        azSun = 2 * std::numbers::pi - azSun;
    }
    // Equation 2.2 in Schörghofer N. (2025). Planetary Code Collection User Guide. theta = 90 minus
    // incidence angle for sloped surface
    T sintheta = std::cos(m_slope_angle) * sinbeta
                 - std::sin(m_slope_angle) * cosbeta * std::cos(azSun - m_azimuth_topo_grad);
    if (cosbeta == 0)
    {
        sintheta = std::cos(m_slope_angle) * sinbeta;
    }
    // sintheta < 0 = sun below local horizon, sinbeta < 0 = sun below horizon at infinity
    if (sintheta < 0.0 || sinbeta < 0.0)
    {
        sintheta = 0.0;
    }
    return sintheta;
}
 
template <typename T>
void SolarFluxVariableCustomSinusoidal<T>::calculateSolarFlux()
{
    m_heliocentric_distance = this->sim->getField("HeliocentricDistance")[0];
    T sintheta{calculateIncidenceAngle()};
    const std::valarray<T>& albedo{this->sim->getField("Albedo")};
    for (int i = 0; i < albedo.size(); ++i)
    {
        (*this->module_field)[i] =
            (1 - albedo[i]) * sintheta * m_solar_const / std::pow(m_heliocentric_distance, 2);
    }
}
 
template <typename T>
bool SolarFluxVariableCustomSinusoidal<T>::init()
{
    calculateSolarFlux();
    return true;
}
 
template <typename T>
bool SolarFluxVariableCustomSinusoidal<T>::preTimeStep()
{
    calculateSolarFlux();
    return true;
}
 
template <typename T>
bool SolarFluxVariableCustomSinusoidal<T>::m_registered =
    ModuleFactory<T>::registerModule(getName(), createMethode);
 
#endif // SOLAR_FLUX_CUSTOM_SINUSOIDAL_H