<commit_before>/*
  The incident shortwave radiation model is an algebraic model with dependencies.

  Generated via evaluator_generator with:
Aspect modified shortwave radiation is determined by a factor which
is multiplied by the 'incoming radiation incident on a flat surface'
to determine the 'incoming radiation incident on a sloping surface of
a given aspect' as a function of latitude, slope, aspect, and Julian
day of the year, and time of day.

Note that some careful checking and experimentation has found that, in
general, the daily average incident radiation times the 12-noon aspect
modifier correlates reasonably well with the daily average of the
product of the hourly incident radiation and the hourly aspect
modifier.  It is notably better than the daily average radiation times
the daily average aspect modifier.



  Authors: <NAME> (<EMAIL>)
*/
#include <cmath>

#include "Teuchos_ParameterList.hpp"
#include "dbc.hh"
#include "exceptions.hh"
#include "errors.hh"
#include "incident_shortwave_radiation_model.hh"

namespace Amanzi {
namespace SurfaceBalance {
namespace Relations {

// Constructor from ParameterList
IncidentShortwaveRadiationModel::IncidentShortwaveRadiationModel(Teuchos::ParameterList& plist)
{
  InitializeFromPlist_(plist);
}


// Initialize parameters
void
IncidentShortwaveRadiationModel::InitializeFromPlist_(Teuchos::ParameterList& plist)
{
  daily_avg_ = plist.get<bool>("daily averaged", true);

  lat_ = plist.get<double>("latitude [degrees]");
  if (lat_ < -90 || lat_ > 90) {
    Errors::Message msg("IncidentShortwaveRadiationModel: \"domain-averaged latitude [degrees]\" not in valid range [-90,90]");
    Exceptions::amanzi_throw(msg);
  }

  doy0_ = plist.get<int>("day of year at time 0 [Julian days]", 0);
  if (doy0_ < 0 || doy0_ > 364) {
    Errors::Message msg("IncidentShortwaveRadiationModel: \"day of year at time 0 [Julian days]\" not in valid range [0,364]");
    Exceptions::amanzi_throw(msg);
  }
}


// main method
double
IncidentShortwaveRadiationModel::IncidentShortwaveRadiation(double slope, double aspect, double qSWin, double time) const
{
  double time_days = time / 86400.0;
  double doy = std::fmod((double)doy0_ + time_days, (double)365);

  int doy_i = std::lround(doy);
  if (doy_i == 365) {
    // can round up!
    doy_i = 0;
    doy = doy - 365.0;
  }

  double rad = 0.0;
  if (daily_avg_) {
    double hour = 12;
    // to keep this function smooth, we interpolate between neighboring days
    if (doy_i < doy) {
      double rad_i = Impl::Radiation(slope, aspect, doy_i, hour, lat_, qSWin);
      int doy_ii = doy_i + 1;
      if (doy_ii > 364) doy_ii = 0;
      double rad_ii = Impl::Radiation(slope, aspect, doy_ii, hour, lat_, qSWin);
      rad = rad_i + (doy - doy_i) * rad_ii;
    } else {
      double rad_i = Impl::Radiation(slope, aspect, doy_i, hour, lat_, qSWin);
      int doy_ii = doy_i - 1;
      if (doy_ii < 0) doy_ii = 364;
      double rad_ii = Impl::Radiation(slope, aspect, doy_ii, hour, lat_, qSWin);
      rad = rad_i + (doy_i - doy) * rad_ii;
    }
  } else {
    double hour = 12.0 + 24 * (doy - doy_i);
    rad = Impl::Radiation(slope, aspect, doy_i, hour, lat_, qSWin);
  }

  return rad;
}

double
IncidentShortwaveRadiationModel::DIncidentShortwaveRadiationDSlope(double slope, double aspect, double qSWin, double time) const
{
  AMANZI_ASSERT(false);
  return 0;
}

double
IncidentShortwaveRadiationModel::DIncidentShortwaveRadiationDAspect(double slope, double aspect, double qSWin, double time) const
{
  AMANZI_ASSERT(false);
  return 0;
}

double
IncidentShortwaveRadiationModel::DIncidentShortwaveRadiationDIncomingShortwaveRadiation(double slope, double aspect, double qSWin, double time) const
{
  return IncidentShortwaveRadiation(slope, aspect, qSWin, time) / qSWin;
}

namespace Impl {

    /*Declination angle

    Parameters
    ----------
    doy : int
      Julian day of the year

    Returns
    -------
    delta : double
      The angle, in radians
    */
double DeclinationAngle(double doy) {
  return 23.45 * M_PI / 180.0 * std::cos(M_2_PI / 365 * (172 - doy));
}

    /*Angle of the sun as a function of the time of day

    Parameters
    ----------
    hour : double
       Time, in a 24 hour clock in [0,24)

    Returns
    -------
    tau : angle in radians
    */
double HourAngle(double hour)
{
  return (hour + 12) * M_PI / 12.0;
}
<commit_msg>bugfix: incorrect value in solar declination angle calculation<commit_after>/*
  The incident shortwave radiation model is an algebraic model with dependencies.

  Generated via evaluator_generator with:
Aspect modified shortwave radiation is determined by a factor which
is multiplied by the 'incoming radiation incident on a flat surface'
to determine the 'incoming radiation incident on a sloping surface of
a given aspect' as a function of latitude, slope, aspect, and Julian
day of the year, and time of day.

Note that some careful checking and experimentation has found that, in
general, the daily average incident radiation times the 12-noon aspect
modifier correlates reasonably well with the daily average of the
product of the hourly incident radiation and the hourly aspect
modifier.  It is notably better than the daily average radiation times
the daily average aspect modifier.



  Authors: <NAME> (<EMAIL>)
*/
#include <cmath>

#include "Teuchos_ParameterList.hpp"
#include "dbc.hh"
#include "exceptions.hh"
#include "errors.hh"
#include "incident_shortwave_radiation_model.hh"

namespace Amanzi {
namespace SurfaceBalance {
namespace Relations {

// Constructor from ParameterList
IncidentShortwaveRadiationModel::IncidentShortwaveRadiationModel(Teuchos::ParameterList& plist)
{
  InitializeFromPlist_(plist);
}


// Initialize parameters
void
IncidentShortwaveRadiationModel::InitializeFromPlist_(Teuchos::ParameterList& plist)
{
  daily_avg_ = plist.get<bool>("daily averaged", true);

  lat_ = plist.get<double>("latitude [degrees]");
  if (lat_ < -90 || lat_ > 90) {
    Errors::Message msg("IncidentShortwaveRadiationModel: \"domain-averaged latitude [degrees]\" not in valid range [-90,90]");
    Exceptions::amanzi_throw(msg);
  }

  doy0_ = plist.get<int>("day of year at time 0 [Julian days]", 0);
  if (doy0_ < 0 || doy0_ > 364) {
    Errors::Message msg("IncidentShortwaveRadiationModel: \"day of year at time 0 [Julian days]\" not in valid range [0,364]");
    Exceptions::amanzi_throw(msg);
  }
}


// main method
double
IncidentShortwaveRadiationModel::IncidentShortwaveRadiation(double slope, double aspect, double qSWin, double time) const
{
  double time_days = time / 86400.0;
  double doy = std::fmod((double)doy0_ + time_days, (double)365);

  int doy_i = std::lround(doy);
  if (doy_i == 365) {
    // can round up!
    doy_i = 0;
    doy = doy - 365.0;
  }

  double rad = 0.0;
  if (daily_avg_) {
    double hour = 12;
    // to keep this function smooth, we interpolate between neighboring days
    if (doy_i < doy) {
      double rad_i = Impl::Radiation(slope, aspect, doy_i, hour, lat_, qSWin);
      int doy_ii = doy_i + 1;
      if (doy_ii > 364) doy_ii = 0;
      double rad_ii = Impl::Radiation(slope, aspect, doy_ii, hour, lat_, qSWin);
      rad = rad_i + (doy - doy_i) * rad_ii;
    } else {
      double rad_i = Impl::Radiation(slope, aspect, doy_i, hour, lat_, qSWin);
      int doy_ii = doy_i - 1;
      if (doy_ii < 0) doy_ii = 364;
      double rad_ii = Impl::Radiation(slope, aspect, doy_ii, hour, lat_, qSWin);
      rad = rad_i + (doy_i - doy) * rad_ii;
    }
  } else {
    double hour = 12.0 + 24 * (doy - doy_i);
    rad = Impl::Radiation(slope, aspect, doy_i, hour, lat_, qSWin);
  }

  return rad;
}

double
IncidentShortwaveRadiationModel::DIncidentShortwaveRadiationDSlope(double slope, double aspect, double qSWin, double time) const
{
  AMANZI_ASSERT(false);
  return 0;
}

double
IncidentShortwaveRadiationModel::DIncidentShortwaveRadiationDAspect(double slope, double aspect, double qSWin, double time) const
{
  AMANZI_ASSERT(false);
  return 0;
}

double
IncidentShortwaveRadiationModel::DIncidentShortwaveRadiationDIncomingShortwaveRadiation(double slope, double aspect, double qSWin, double time) const
{
  return IncidentShortwaveRadiation(slope, aspect, qSWin, time) / qSWin;
}

namespace Impl {

    /*Declination angle

    Parameters
    ----------
    doy : int
      Julian day of the year

    Returns
    -------
    delta : double
      The angle, in radians
    */
double DeclinationAngle(double doy) {
  return 23.45 * M_PI / 180.0 * std::cos(2 * M_PI / 365 * (172 - doy));
}

    /*Angle of the sun as a function of the time of day

    Parameters
    ----------
    hour : double
       Time, in a 24 hour clock in [0,24)

    Returns
    -------
    tau : angle in radians
    */
double HourAngle(double hour)
{
  return (hour + 12) * M_PI / 12.0;
}

{'fraction_non_alphanumeric': 0.10538243626062323, 'fraction_numerical': 0.019943342776203965, 'mean_word_length': 3.7248394004282654, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1571086', 'n_tokens_mistral': 3111, 'n_tokens_neox': 2886, 'n_words': 912}