From 9701d54e38c9acbcfe5bf74451feee51d4fbb31f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jo=C3=A3o=20Marcos=20Mororo=20Costa?= <43460229+mororo250@users.noreply.github.com> Date: Thu, 4 Jun 2020 04:28:36 -0300 Subject: [PATCH] Change code used. --- sandbox/settings/appleseed.studio.xml | 4 + src/appleseed/CMakeLists.txt | 2 + .../environmentedf/hosekenvironmentedf.cpp | 48 +- .../environmentedf/skymodelenvironmentedf.cpp | 27 + .../environmentedf/skymodelenvironmentedf.h | 57 + .../renderer/modeling/light/sunlight.cpp | 25 +- .../renderer/modeling/scene/scene.cpp | 12 +- .../utility/solarpositionalgorithm.cpp | 1291 +++-------------- .../renderer/utility/solarpositionalgorithm.h | 293 ++-- 9 files changed, 383 insertions(+), 1376 deletions(-) create mode 100644 src/appleseed/renderer/modeling/environmentedf/skymodelenvironmentedf.cpp create mode 100644 src/appleseed/renderer/modeling/environmentedf/skymodelenvironmentedf.h diff --git a/sandbox/settings/appleseed.studio.xml b/sandbox/settings/appleseed.studio.xml index cc212381bc..b3754b3b21 100644 --- a/sandbox/settings/appleseed.studio.xml +++ b/sandbox/settings/appleseed.studio.xml @@ -11,6 +11,10 @@ + + + + diff --git a/src/appleseed/CMakeLists.txt b/src/appleseed/CMakeLists.txt index 7ca0c8bbfa..e058eede8f 100644 --- a/src/appleseed/CMakeLists.txt +++ b/src/appleseed/CMakeLists.txt @@ -1829,6 +1829,8 @@ set (renderer_modeling_environmentedf_sources renderer/modeling/environmentedf/oslenvironmentedf.h renderer/modeling/environmentedf/preethamenvironmentedf.cpp renderer/modeling/environmentedf/preethamenvironmentedf.h + renderer/modeling/environmentedf/skymodelenvironmentedf.cpp + renderer/modeling/environmentedf/skymodelenvironmentedf.h renderer/modeling/environmentedf/sphericalcoordinates.h ) list (APPEND appleseed_sources diff --git a/src/appleseed/renderer/modeling/environmentedf/hosekenvironmentedf.cpp b/src/appleseed/renderer/modeling/environmentedf/hosekenvironmentedf.cpp index 93bc44c022..dc6f98d4ab 100644 --- a/src/appleseed/renderer/modeling/environmentedf/hosekenvironmentedf.cpp +++ b/src/appleseed/renderer/modeling/environmentedf/hosekenvironmentedf.cpp @@ -34,7 +34,7 @@ #include "renderer/global/globaltypes.h" #include "renderer/kernel/shading/shadingcontext.h" #include "renderer/modeling/color/colorspace.h" -#include "renderer/modeling/environmentedf/environmentedf.h" +#include "renderer/modeling/environmentedf/skymodelenvironmentedf.h" #include "renderer/modeling/environmentedf/sphericalcoordinates.h" #include "renderer/modeling/input/inputarray.h" #include "renderer/modeling/input/source.h" @@ -93,13 +93,13 @@ namespace const float BaseTurbidity = 2.0f; class HosekEnvironmentEDF - : public EnvironmentEDF + : public SkyModelEnvironmentEDF { public: HosekEnvironmentEDF( const char* name, const ParamArray& params) - : EnvironmentEDF(name, params) + : SkyModelEnvironmentEDF(name, params) { m_inputs.declare("sun_theta", InputFormat::InputFormatFloat); m_inputs.declare("sun_phi", InputFormat::InputFormatFloat); @@ -128,39 +128,18 @@ namespace OnFrameBeginRecorder& recorder, IAbortSwitch* abort_switch) override { - if (!EnvironmentEDF::on_frame_begin(project, parent, recorder, abort_switch)) + if (!SkyModelEnvironmentEDF::on_frame_begin(project, parent, recorder, abort_switch)) return false; // Evaluate uniform values. m_inputs.evaluate_uniforms(&m_uniform_values); - - spa_data test; - test.hour = m_params.get_optional("hour", 12); - test.minute = m_params.get_optional("minute", 0); - test.second = m_params.get_optional("second", 0); - - test.year = m_params.get_optional("year", 2020); - test.month = m_params.get_optional("month", 1); - test.day = m_params.get_optional("day", 1); - - test.timezone = 0; - test.longitude = m_params.get_optional("longitude", 0); - test.latitude = m_params.get_optional("latitude", 0); - - test.azm_rotation = 0.0; - test.elevation = 0.0; - - test.function = SPA_ZA; - - spa_calculate(&test); - - m_uniform_values.m_sun_theta = static_cast(test.zenith); - m_uniform_values.m_sun_phi = static_cast(test.azimuth); // Compute the sun direction. - m_sun_theta = deg_to_rad(m_uniform_values.m_sun_theta); - m_sun_phi = deg_to_rad(m_uniform_values.m_sun_phi); - m_sun_dir = Vector3f::make_unit_vector(m_sun_theta, m_sun_phi); + m_sun_positioner.m_output_value.zenith = deg_to_rad(m_sun_positioner.m_output_value.zenith); + m_sun_positioner.m_output_value.azimuth = deg_to_rad(m_sun_positioner.m_output_value.azimuth); + m_sun_dir = Vector3f::make_unit_vector(static_cast(m_sun_positioner.m_output_value.zenith) + , static_cast(m_sun_positioner.m_output_value.azimuth)); + m_horizon_shift = m_uniform_values.m_horizon_shift; // Precompute the coefficients of the radiance distribution function and // the master luminance value if turbidity is uniform. @@ -174,7 +153,7 @@ namespace compute_coefficients( m_uniform_values.m_turbidity, m_uniform_values.m_ground_albedo, - m_sun_theta, + static_cast(m_sun_positioner.m_output_value.zenith), m_uniform_coeffs, m_uniform_master_Y); } @@ -280,9 +259,6 @@ namespace }; InputValues m_uniform_values; - - float m_sun_theta; // sun zenith angle in radians, 0=zenith - float m_sun_phi; // radians Vector3f m_sun_dir; bool m_uniform_turbidity; @@ -447,7 +423,7 @@ namespace compute_coefficients( turbidity, m_uniform_values.m_ground_albedo, - m_sun_theta, + static_cast(m_sun_positioner.m_output_value.zenith), coeffs, master_Y); @@ -492,7 +468,7 @@ namespace Vector3f shift(Vector3f v) const { - v.y -= m_uniform_values.m_horizon_shift; + v.y -= m_horizon_shift; return normalize(v); } }; diff --git a/src/appleseed/renderer/modeling/environmentedf/skymodelenvironmentedf.cpp b/src/appleseed/renderer/modeling/environmentedf/skymodelenvironmentedf.cpp new file mode 100644 index 0000000000..cc4a3c927b --- /dev/null +++ b/src/appleseed/renderer/modeling/environmentedf/skymodelenvironmentedf.cpp @@ -0,0 +1,27 @@ + +#include "skymodelenvironmentedf.h" + +namespace renderer +{ + SkyModelEnvironmentEDF::SkyModelEnvironmentEDF( + const char* name, + const ParamArray& params) + : EnvironmentEDF(name, params), + m_sun_positioner("Sun Positioner", params) + { + } + + bool SkyModelEnvironmentEDF::on_frame_begin(const Project& project, const BaseGroup* parent, OnFrameBeginRecorder& recorder, foundation::IAbortSwitch* abort_switch) + { + if (!EnvironmentEDF::on_frame_begin(project, parent, recorder, abort_switch)) + return false; + + if (!m_sun_positioner.on_frame_begin(project, parent, recorder, abort_switch)) + return false; + + m_sun_positioner.compute_sun_position(); + + return true; + } + +} // namespace renderer diff --git a/src/appleseed/renderer/modeling/environmentedf/skymodelenvironmentedf.h b/src/appleseed/renderer/modeling/environmentedf/skymodelenvironmentedf.h new file mode 100644 index 0000000000..5c56af385e --- /dev/null +++ b/src/appleseed/renderer/modeling/environmentedf/skymodelenvironmentedf.h @@ -0,0 +1,57 @@ +#pragma once + +// appleseed.renderer headers. +#include "renderer/modeling/environmentedf/environmentedf.h" +#include "renderer/utility/solarpositionalgorithm.h" + +namespace renderer +{ + + // + // An Environment EDF Interface Class For Differents Sky Models implementation. + // + + class APPLESEED_DLLSYMBOL SkyModelEnvironmentEDF + : public EnvironmentEDF + { + public: + + // Constructor. + SkyModelEnvironmentEDF( + const char* name, + const ParamArray& params); + + bool on_frame_begin( + const Project& project, + const BaseGroup* parent, + OnFrameBeginRecorder& recorder, + foundation::IAbortSwitch* abort_switch = nullptr) override; + + float get_sun_theta() const; + + float get_sun_phi() const; + + float get_shift() const; + + protected: + + SunPositioner m_sun_positioner; + float m_horizon_shift; + }; + + inline float SkyModelEnvironmentEDF::get_sun_theta() const + { + return static_cast(m_sun_positioner.m_output_value.zenith); + } + + inline float SkyModelEnvironmentEDF::get_sun_phi() const + { + return static_cast(m_sun_positioner.m_output_value.azimuth); + } + + inline float SkyModelEnvironmentEDF::get_shift() const + { + return m_horizon_shift; + } + +} // namespace renderer diff --git a/src/appleseed/renderer/modeling/light/sunlight.cpp b/src/appleseed/renderer/modeling/light/sunlight.cpp index b2b777af0d..6eac1d3d47 100644 --- a/src/appleseed/renderer/modeling/light/sunlight.cpp +++ b/src/appleseed/renderer/modeling/light/sunlight.cpp @@ -34,7 +34,7 @@ #include "renderer/global/globaltypes.h" #include "renderer/modeling/color/colorspace.h" #include "renderer/modeling/color/wavelengths.h" -#include "renderer/modeling/environmentedf/environmentedf.h" +#include "renderer/modeling/environmentedf/skymodelenvironmentedf.h" #include "renderer/modeling/input/inputarray.h" #include "renderer/modeling/input/source.h" #include "renderer/modeling/light/lighttarget.h" @@ -157,7 +157,7 @@ namespace m_sun_solid_angle = TwoPi() * (1.0f - std::cos(std::atan(SunRadius * m_values.m_size_multiplier / m_values.m_distance))); // If the Sun light is bound to an environment EDF, let it override the Sun's direction and turbidity. - const EnvironmentEDF* env_edf = dynamic_cast(m_inputs.get_entity("environment_edf")); + const SkyModelEnvironmentEDF* env_edf = dynamic_cast(m_inputs.get_entity("environment_edf")); if (env_edf != nullptr) apply_env_edf_overrides(env_edf); @@ -327,24 +327,13 @@ namespace RegularSpectrum31f m_k1; RegularSpectrum31f m_k2; - void apply_env_edf_overrides(const EnvironmentEDF* env_edf) + void apply_env_edf_overrides(const SkyModelEnvironmentEDF* env_edf) { // Use the Sun direction from the EDF if it has one. - const Source* sun_theta_src = env_edf->get_inputs().source("sun_theta"); - const Source* sun_phi_src = env_edf->get_inputs().source("sun_phi"); - const Source* sun_shift_src = env_edf-> get_inputs().source("horizon_shift"); - - if (sun_theta_src != nullptr && - sun_theta_src->is_uniform() && - sun_phi_src != nullptr && - sun_phi_src->is_uniform() && - sun_shift_src != nullptr && - sun_shift_src->is_uniform()) { - float sun_theta, sun_phi, sun_shift; - sun_theta_src->evaluate_uniform(sun_theta); - sun_phi_src->evaluate_uniform(sun_phi); - sun_shift_src->evaluate_uniform(sun_shift); + const float sun_theta = env_edf->get_sun_theta(); + const float sun_phi = env_edf->get_sun_phi(); + const float sun_shift = env_edf->get_shift(); Transformd scratch; const Transformd& env_edf_transform = env_edf->transform_sequence().evaluate(0.0f, scratch); @@ -355,7 +344,7 @@ namespace Matrix4d::make_rotation( Quaterniond::make_rotation( Vector3d(0.0, 0.0, -1.0), // default emission direction of this light - -Vector3d::make_unit_vector(deg_to_rad(sun_theta), deg_to_rad(sun_phi))))) * + -Vector3d::make_unit_vector(sun_theta, sun_phi)))) * env_edf_transform); } diff --git a/src/appleseed/renderer/modeling/scene/scene.cpp b/src/appleseed/renderer/modeling/scene/scene.cpp index 1831f46833..bd1d77ee85 100644 --- a/src/appleseed/renderer/modeling/scene/scene.cpp +++ b/src/appleseed/renderer/modeling/scene/scene.cpp @@ -395,12 +395,6 @@ bool Scene::on_frame_begin( OnFrameBeginRecorder& recorder, IAbortSwitch* abort_switch) { - if (!Entity::on_frame_begin(project, parent, recorder, abort_switch)) - return false; - - if (!BaseGroup::on_frame_begin(project, parent, recorder, abort_switch)) - return false; - bool success = true; success = success && impl->m_default_surface_shader->on_frame_begin(project, this, recorder, abort_switch); success = success && invoke_on_frame_begin(environment_edfs(), project, this, recorder, abort_switch); @@ -411,6 +405,12 @@ bool Scene::on_frame_begin( // Call on_frame_begin() on cameras last because some of them cast rays to sense depth in their autofocus mechanism. success = success && invoke_on_frame_begin(cameras(), project, this, recorder, abort_switch); + if (!Entity::on_frame_begin(project, parent, recorder, abort_switch)) + return false; + + if (!BaseGroup::on_frame_begin(project, parent, recorder, abort_switch)) + return false; + return success; } diff --git a/src/appleseed/renderer/utility/solarpositionalgorithm.cpp b/src/appleseed/renderer/utility/solarpositionalgorithm.cpp index 59eddd9660..a1c39341f0 100644 --- a/src/appleseed/renderer/utility/solarpositionalgorithm.cpp +++ b/src/appleseed/renderer/utility/solarpositionalgorithm.cpp @@ -1,1174 +1,229 @@ -///////////////////////////////////////////// -// Solar Position Algorithm (SPA) // -// for // -// Solar Radiation Application // -// // -// May 12, 2003 // -// // -// Filename: SPA.C // -// // -// Afshin Michael Andreas // -// Afshin.Andreas@NREL.gov (303)384-6383 // -// // -// Metrology Laboratory // -// Solar Radiation Research Laboratory // -// National Renewable Energy Laboratory // -// 15013 Denver W Pkwy, Golden, CO 80401 // -///////////////////////////////////////////// -///////////////////////////////////////////// -// See the SPA.H header file for usage // -// // -// This code is based on the NREL // -// technical report "Solar Position // -// Algorithm for Solar Radiation // -// Application" by I. Reda & A. Andreas // -///////////////////////////////////////////// - -/////////////////////////////////////////////////////////////////////////////////////////////// // -// NOTICE -// Copyright (C) 2008-2011 Alliance for Sustainable Energy, LLC, All Rights Reserved +// This source file is part of appleseed. +// Visit https://appleseedhq.net/ for additional information and resources. // -//The Solar Position Algorithm ("Software") is code in development prepared by employees of the -//Alliance for Sustainable Energy, LLC, (hereinafter the "Contractor"), under Contract No. -//DE-AC36-08GO28308 ("Contract") with the U.S. Department of Energy (the "DOE"). The United -//States Government has been granted for itself and others acting on its behalf a paid-up, non- -//exclusive, irrevocable, worldwide license in the Software to reproduce, prepare derivative -//works, and perform publicly and display publicly. Beginning five (5) years after the date -//permission to assert copyright is obtained from the DOE, and subject to any subsequent five -//(5) year renewals, the United States Government is granted for itself and others acting on -//its behalf a paid-up, non-exclusive, irrevocable, worldwide license in the Software to -//reproduce, prepare derivative works, distribute copies to the public, perform publicly and -//display publicly, and to permit others to do so. If the Contractor ceases to make this -//computer software available, it may be obtained from DOE's Office of Scientific and Technical -//Information's Energy Science and Technology Software Center (ESTSC) at P.O. Box 1020, Oak -//Ridge, TN 37831-1020. THIS SOFTWARE IS PROVIDED BY THE CONTRACTOR "AS IS" AND ANY EXPRESS OR -//IMPLIED WARRANTIES, INCLUDING BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY -//AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE CONTRACTOR OR THE -//U.S. GOVERNMENT BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES -//WHATSOEVER, INCLUDING BUT NOT LIMITED TO CLAIMS ASSOCIATED WITH THE LOSS OF DATA OR PROFITS, -//WHICH MAY RESULT FROM AN ACTION IN CONTRACT, NEGLIGENCE OR OTHER TORTIOUS CLAIM THAT ARISES -//OUT OF OR IN CONNECTION WITH THE ACCESS, USE OR PERFORMANCE OF THIS SOFTWARE. +// This software is released under the MIT license. // -//The Software is being provided for internal, noncommercial purposes only and shall not be -//re-distributed. Please contact the NREL Commercialization and Technology Transfer Office -//for information concerning a commercial license to use the Software, visit: -//http://midcdmz.nrel.gov/spa/ for the contact information. +// Copyright (c) 2020 João Marcos , Jupiter Jazz Limited // -//As a condition of using the Software in an application, the developer of the application -//agrees to reference the use of the Software and make this Notice readily accessible to any -//end-user in a Help|About screen or equivalent manner. +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +// THE SOFTWARE. // -/////////////////////////////////////////////////////////////////////////////////////////////// - -/////////////////////////////////////////////////////////////////////////////////////////////// -// Revised 27-FEB-2004 Andreas -// Added bounds check on inputs and return value for spa_calculate(). -// Revised 10-MAY-2004 Andreas -// Changed temperature bound check minimum from -273.15 to -273 degrees C. -// Revised 17-JUN-2004 Andreas -// Corrected a problem that caused a bogus sunrise/set/transit on the equinox. -// Revised 18-JUN-2004 Andreas -// Added a "function" input variable that allows the selecting of desired outputs. -// Revised 21-JUN-2004 Andreas -// Added 3 new intermediate output values to SPA structure (srha, ssha, & sta). -// Revised 23-JUN-2004 Andreas -// Enumerations for "function" were renamed and 2 were added. -// Prevented bound checks on inputs that are not used (based on function). -// Revised 01-SEP-2004 Andreas -// Changed a local variable from integer to double. -// Revised 12-JUL-2005 Andreas -// Put a limit on the EOT calculation, so that the result is between -20 and 20. -// Revised 26-OCT-2005 Andreas -// Set the atmos. refraction correction to zero, when sun is below horizon. -// Made atmos_refract input a requirement for all "functions". -// Changed atmos_refract bound check from +/- 10 to +/- 5 degrees. -// Revised 07-NOV-2006 Andreas -// Corrected 3 earth periodic terms in the L_TERMS array. -// Corrected 2 earth periodic terms in the R_TERMS array. -// Revised 10-NOV-2006 Andreas -// Corrected a constant used to calculate topocentric sun declination. -// Put a limit on observer hour angle, so result is between 0 and 360. -// Revised 13-NOV-2006 Andreas -// Corrected calculation of topocentric sun declination. -// Converted all floating point inputs in spa structure to doubles. -// Revised 27-FEB-2007 Andreas -// Minor correction made as to when atmos. refraction correction is set to zero. -// Revised 21-JAN-2008 Andreas -// Minor change to two variable declarations. -// Revised 12-JAN-2009 Andreas -// Changed timezone bound check from +/-12 to +/-18 hours. -// Revised 14-JAN-2009 Andreas -// Corrected a constant used to calculate ecliptic mean obliquity. -// Revised 01-APR-2013 Andreas -// Replace floor with new integer function for tech. report consistency, no affect on results. -// Add "utility" function prototypes to header file for use with NREL's SAMPA. -// Rename 4 "utility" function names (remove "sun") for clarity with NREL's SAMPA. -// Added delta_ut1 as required input, which the fractional second difference between UT and UTC. -// Time must be input w/o delta_ut1 adjustment, instead of assuming adjustment was pre-applied. -// Revised 10-JUL-2014 Andreas -// Change second in spa_data structure from an integer to double to allow fractional second -// Revised 08-SEP-2014 Andreas -// Corrected description of azm_rotation in header file -// Limited azimuth180 to range of 0 to 360 deg (instead of -180 to 180) for tech report consistency -// Changed all variables names from azimuth180 to azimuth_astro -// Renamed 2 "utility" function names for consistency -/////////////////////////////////////////////////////////////////////////////////////////////// -#include +// Interface header. #include "solarpositionalgorithm.h" -#define PI 3.1415926535897932384626433832795028841971 -#define SUN_RADIUS 0.26667 - -#define L_COUNT 6 -#define B_COUNT 2 -#define R_COUNT 5 -#define Y_COUNT 63 +// appleseed.foundation headers. +#include "foundation/math/scalar.h" +#include "foundation/utility/job/abortswitch.h" +#include "foundation/utility/uid.h" -#define L_MAX_SUBCOUNT 64 -#define B_MAX_SUBCOUNT 5 -#define R_MAX_SUBCOUNT 40 -enum { TERM_A, TERM_B, TERM_C, TERM_COUNT }; -enum { TERM_X0, TERM_X1, TERM_X2, TERM_X3, TERM_X4, TERM_X_COUNT }; -enum { TERM_PSI_A, TERM_PSI_B, TERM_EPS_C, TERM_EPS_D, TERM_PE_COUNT }; -enum { JD_MINUS, JD_ZERO, JD_PLUS, JD_COUNT }; -enum { SUN_TRANSIT, SUN_RISE, SUN_SET, SUN_COUNT }; +// appleseed.renderer headers. +#include "renderer/modeling/project/project.h" +#include "renderer/modeling/scene/basegroup.h" +#include "renderer/modeling/entity/onframebeginrecorder.h" +#include "renderer/utility/paramarray.h" -#define TERM_Y_COUNT TERM_X_COUNT +// Standard headers. +#include -const int l_subcount[L_COUNT] = { 64,34,20,7,3,1 }; -const int b_subcount[B_COUNT] = { 5,2 }; -const int r_subcount[R_COUNT] = { 40,10,6,2,1 }; +using namespace foundation; -/////////////////////////////////////////////////// -/// Earth Periodic Terms -/////////////////////////////////////////////////// -const double L_TERMS[L_COUNT][L_MAX_SUBCOUNT][TERM_COUNT] = +namespace renderer { - { - {175347046.0,0,0}, - {3341656.0,4.6692568,6283.07585}, - {34894.0,4.6261,12566.1517}, - {3497.0,2.7441,5753.3849}, - {3418.0,2.8289,3.5231}, - {3136.0,3.6277,77713.7715}, - {2676.0,4.4181,7860.4194}, - {2343.0,6.1352,3930.2097}, - {1324.0,0.7425,11506.7698}, - {1273.0,2.0371,529.691}, - {1199.0,1.1096,1577.3435}, - {990,5.233,5884.927}, - {902,2.045,26.298}, - {857,3.508,398.149}, - {780,1.179,5223.694}, - {753,2.533,5507.553}, - {505,4.583,18849.228}, - {492,4.205,775.523}, - {357,2.92,0.067}, - {317,5.849,11790.629}, - {284,1.899,796.298}, - {271,0.315,10977.079}, - {243,0.345,5486.778}, - {206,4.806,2544.314}, - {205,1.869,5573.143}, - {202,2.458,6069.777}, - {156,0.833,213.299}, - {132,3.411,2942.463}, - {126,1.083,20.775}, - {115,0.645,0.98}, - {103,0.636,4694.003}, - {102,0.976,15720.839}, - {102,4.267,7.114}, - {99,6.21,2146.17}, - {98,0.68,155.42}, - {86,5.98,161000.69}, - {85,1.3,6275.96}, - {85,3.67,71430.7}, - {80,1.81,17260.15}, - {79,3.04,12036.46}, - {75,1.76,5088.63}, - {74,3.5,3154.69}, - {74,4.68,801.82}, - {70,0.83,9437.76}, - {62,3.98,8827.39}, - {61,1.82,7084.9}, - {57,2.78,6286.6}, - {56,4.39,14143.5}, - {56,3.47,6279.55}, - {52,0.19,12139.55}, - {52,1.33,1748.02}, - {51,0.28,5856.48}, - {49,0.49,1194.45}, - {41,5.37,8429.24}, - {41,2.4,19651.05}, - {39,6.17,10447.39}, - {37,6.04,10213.29}, - {37,2.57,1059.38}, - {36,1.71,2352.87}, - {36,1.78,6812.77}, - {33,0.59,17789.85}, - {30,0.44,83996.85}, - {30,2.74,1349.87}, - {25,3.16,4690.48} - }, - { - {628331966747.0,0,0}, - {206059.0,2.678235,6283.07585}, - {4303.0,2.6351,12566.1517}, - {425.0,1.59,3.523}, - {119.0,5.796,26.298}, - {109.0,2.966,1577.344}, - {93,2.59,18849.23}, - {72,1.14,529.69}, - {68,1.87,398.15}, - {67,4.41,5507.55}, - {59,2.89,5223.69}, - {56,2.17,155.42}, - {45,0.4,796.3}, - {36,0.47,775.52}, - {29,2.65,7.11}, - {21,5.34,0.98}, - {19,1.85,5486.78}, - {19,4.97,213.3}, - {17,2.99,6275.96}, - {16,0.03,2544.31}, - {16,1.43,2146.17}, - {15,1.21,10977.08}, - {12,2.83,1748.02}, - {12,3.26,5088.63}, - {12,5.27,1194.45}, - {12,2.08,4694}, - {11,0.77,553.57}, - {10,1.3,6286.6}, - {10,4.24,1349.87}, - {9,2.7,242.73}, - {9,5.64,951.72}, - {8,5.3,2352.87}, - {6,2.65,9437.76}, - {6,4.67,4690.48} - }, - { - {52919.0,0,0}, - {8720.0,1.0721,6283.0758}, - {309.0,0.867,12566.152}, - {27,0.05,3.52}, - {16,5.19,26.3}, - {16,3.68,155.42}, - {10,0.76,18849.23}, - {9,2.06,77713.77}, - {7,0.83,775.52}, - {5,4.66,1577.34}, - {4,1.03,7.11}, - {4,3.44,5573.14}, - {3,5.14,796.3}, - {3,6.05,5507.55}, - {3,1.19,242.73}, - {3,6.12,529.69}, - {3,0.31,398.15}, - {3,2.28,553.57}, - {2,4.38,5223.69}, - {2,3.75,0.98} - }, - { - {289.0,5.844,6283.076}, - {35,0,0}, - {17,5.49,12566.15}, - {3,5.2,155.42}, - {1,4.72,3.52}, - {1,5.3,18849.23}, - {1,5.97,242.73} - }, - { - {114.0,3.142,0}, - {8,4.13,6283.08}, - {1,3.84,12566.15} - }, - { - {1,3.14,0} - } -}; -const double B_TERMS[B_COUNT][B_MAX_SUBCOUNT][TERM_COUNT] = -{ + namespace { - {280.0,3.199,84334.662}, - {102.0,5.422,5507.553}, - {80,3.88,5223.69}, - {44,3.7,2352.87}, - {32,4,1577.34} - }, - { - {9,3.9,5507.55}, - {6,1.73,5223.69} + const UniqueID g_class_uid = new_guid(); } -}; -const double R_TERMS[R_COUNT][R_MAX_SUBCOUNT][TERM_COUNT] = -{ - { - {100013989.0,0,0}, - {1670700.0,3.0984635,6283.07585}, - {13956.0,3.05525,12566.1517}, - {3084.0,5.1985,77713.7715}, - {1628.0,1.1739,5753.3849}, - {1576.0,2.8469,7860.4194}, - {925.0,5.453,11506.77}, - {542.0,4.564,3930.21}, - {472.0,3.661,5884.927}, - {346.0,0.964,5507.553}, - {329.0,5.9,5223.694}, - {307.0,0.299,5573.143}, - {243.0,4.273,11790.629}, - {212.0,5.847,1577.344}, - {186.0,5.022,10977.079}, - {175.0,3.012,18849.228}, - {110.0,5.055,5486.778}, - {98,0.89,6069.78}, - {86,5.69,15720.84}, - {86,1.27,161000.69}, - {65,0.27,17260.15}, - {63,0.92,529.69}, - {57,2.01,83996.85}, - {56,5.24,71430.7}, - {49,3.25,2544.31}, - {47,2.58,775.52}, - {45,5.54,9437.76}, - {43,6.01,6275.96}, - {39,5.36,4694}, - {38,2.39,8827.39}, - {37,0.83,19651.05}, - {37,4.9,12139.55}, - {36,1.67,12036.46}, - {35,1.84,2942.46}, - {33,0.24,7084.9}, - {32,0.18,5088.63}, - {32,1.78,398.15}, - {28,1.21,6286.6}, - {28,1.9,6279.55}, - {26,4.59,10447.39} - }, - { - {103019.0,1.10749,6283.07585}, - {1721.0,1.0644,12566.1517}, - {702.0,3.142,0}, - {32,1.02,18849.23}, - {31,2.84,5507.55}, - {25,1.32,5223.69}, - {18,1.42,1577.34}, - {10,5.91,10977.08}, - {9,1.42,6275.96}, - {9,0.27,5486.78} - }, + double get_julian_century(double jd) { - {4359.0,5.7846,6283.0758}, - {124.0,5.579,12566.152}, - {12,3.14,0}, - {9,3.63,77713.77}, - {6,1.87,5573.14}, - {3,5.47,18849.23} - }, - { - {145.0,4.273,6283.076}, - {7,3.92,12566.15} - }, - { - {4,2.56,6283.08} + return (jd - 2451545.0) / 36525.0; } -}; - -//////////////////////////////////////////////////////////////// -/// Periodic Terms for the nutation in longitude and obliquity -//////////////////////////////////////////////////////////////// - -const int Y_TERMS[Y_COUNT][TERM_Y_COUNT] = -{ - {0,0,0,0,1}, - {-2,0,0,2,2}, - {0,0,0,2,2}, - {0,0,0,0,2}, - {0,1,0,0,0}, - {0,0,1,0,0}, - {-2,1,0,2,2}, - {0,0,0,2,1}, - {0,0,1,2,2}, - {-2,-1,0,2,2}, - {-2,0,1,0,0}, - {-2,0,0,2,1}, - {0,0,-1,2,2}, - {2,0,0,0,0}, - {0,0,1,0,1}, - {2,0,-1,2,2}, - {0,0,-1,0,1}, - {0,0,1,2,1}, - {-2,0,2,0,0}, - {0,0,-2,2,1}, - {2,0,0,2,2}, - {0,0,2,2,2}, - {0,0,2,0,0}, - {-2,0,1,2,2}, - {0,0,0,2,0}, - {-2,0,0,2,0}, - {0,0,-1,2,1}, - {0,2,0,0,0}, - {2,0,-1,0,1}, - {-2,2,0,2,2}, - {0,1,0,0,1}, - {-2,0,1,0,1}, - {0,-1,0,0,1}, - {0,0,2,-2,0}, - {2,0,-1,2,1}, - {2,0,1,2,2}, - {0,1,0,2,2}, - {-2,1,1,0,0}, - {0,-1,0,2,2}, - {2,0,0,2,1}, - {2,0,1,0,0}, - {-2,0,2,2,2}, - {-2,0,1,2,1}, - {2,0,-2,0,1}, - {2,0,0,0,1}, - {0,-1,1,0,0}, - {-2,-1,0,2,1}, - {-2,0,0,0,1}, - {0,0,2,2,1}, - {-2,0,2,0,1}, - {-2,1,0,2,1}, - {0,0,1,-2,0}, - {-1,0,1,0,0}, - {-2,1,0,0,0}, - {1,0,0,0,0}, - {0,0,1,2,0}, - {0,0,-2,2,2}, - {-1,-1,1,0,0}, - {0,1,1,0,0}, - {0,-1,1,2,2}, - {2,-1,-1,2,2}, - {0,0,3,2,2}, - {2,-1,0,2,2}, -}; - -const double PE_TERMS[Y_COUNT][TERM_PE_COUNT] = { - {-171996,-174.2,92025,8.9}, - {-13187,-1.6,5736,-3.1}, - {-2274,-0.2,977,-0.5}, - {2062,0.2,-895,0.5}, - {1426,-3.4,54,-0.1}, - {712,0.1,-7,0}, - {-517,1.2,224,-0.6}, - {-386,-0.4,200,0}, - {-301,0,129,-0.1}, - {217,-0.5,-95,0.3}, - {-158,0,0,0}, - {129,0.1,-70,0}, - {123,0,-53,0}, - {63,0,0,0}, - {63,0.1,-33,0}, - {-59,0,26,0}, - {-58,-0.1,32,0}, - {-51,0,27,0}, - {48,0,0,0}, - {46,0,-24,0}, - {-38,0,16,0}, - {-31,0,13,0}, - {29,0,0,0}, - {29,0,-12,0}, - {26,0,0,0}, - {-22,0,0,0}, - {21,0,-10,0}, - {17,-0.1,0,0}, - {16,0,-8,0}, - {-16,0.1,7,0}, - {-15,0,9,0}, - {-13,0,7,0}, - {-12,0,6,0}, - {11,0,0,0}, - {-10,0,5,0}, - {-8,0,3,0}, - {7,0,-3,0}, - {-7,0,0,0}, - {-7,0,3,0}, - {-7,0,3,0}, - {6,0,0,0}, - {6,0,-3,0}, - {6,0,-3,0}, - {-6,0,3,0}, - {-6,0,3,0}, - {5,0,0,0}, - {-5,0,3,0}, - {-5,0,3,0}, - {-5,0,3,0}, - {4,0,0,0}, - {4,0,0,0}, - {4,0,0,0}, - {-4,0,0,0}, - {-4,0,0,0}, - {-4,0,0,0}, - {3,0,0,0}, - {-3,0,0,0}, - {-3,0,0,0}, - {-3,0,0,0}, - {-3,0,0,0}, - {-3,0,0,0}, - {-3,0,0,0}, - {-3,0,0,0}, -}; - -/////////////////////////////////////////////// - -double rad2deg(double radians) -{ - return (180.0 / PI) * radians; -} - -double deg2rad(double degrees) -{ - return (PI / 180.0) * degrees; -} -int integer(double value) -{ - return static_cast(value); -} - -double limit_degrees(double degrees) -{ - double limited; - - degrees /= 360.0; - limited = 360.0 * (degrees - floor(degrees)); - if (limited < 0) limited += 360.0; - - return limited; -} - -double limit_degrees180pm(double degrees) -{ - double limited; - - degrees /= 360.0; - limited = 360.0 * (degrees - floor(degrees)); - if (limited < -180.0) limited += 360.0; - else if (limited > 180.0) limited -= 360.0; - - return limited; -} - -double limit_degrees180(double degrees) -{ - double limited; - - degrees /= 180.0; - limited = 180.0 * (degrees - floor(degrees)); - if (limited < 0) limited += 180.0; + double get_julian_day(int year, int month, int day, double decimal_day) + { + decimal_day += day; - return limited; -} + if (month < 3) + { + month += 12; + year--; + } -double limit_zero2one(double value) -{ - double limited; + double julian_day = std::floor(365.25 * (year + 4716.0)) + + std::floor(30.6001 * (month + 1)) + decimal_day - 1524.5; - limited = value - floor(value); - if (limited < 0) limited += 1.0; + // Test whether to change to gregorian caledar calculation. + if (julian_day > 2299160.0) { + double A = std::floor(year / 100); + julian_day += (2 - A + std::floor(A / 4)); + } - return limited; -} + return julian_day; + } -double limit_minutes(double minutes) -{ - double limited = minutes; + struct SunPositioner::Impl + { + // Input time. + int year; + int month; + int day; + int hour; + int minute; + int second; + int timezone; + + //Iinput location. + double longitude; // Observer longitude (negative west of Greenwich) + double latitude; // Observer latitude (negative south of equator) + double elevation; // Observer elevation [meters] + double azm_rotation; // Surface azimuth rotation. + }; + + // Constructor + SunPositioner::SunPositioner(const char* name, const ParamArray& params) + : Entity(g_class_uid, params), + impl(new Impl()) + { + } - if (limited < -20.0) limited += 1440.0; - else if (limited > 20.0) limited -= 1440.0; + bool SunPositioner::on_frame_begin( + const Project& project, + const BaseGroup* parent, + OnFrameBeginRecorder& recorder, + foundation::IAbortSwitch* abort_switch) + { + if (!Entity::on_frame_begin(project, parent, recorder, abort_switch)) + return false; - return limited; -} + impl->year = m_params.get_optional("year", 2020); + impl->month = m_params.get_optional("month", 1); + impl->day = m_params.get_optional("day", 1); + impl->hour = m_params.get_optional("hour", 12); + impl->minute = m_params.get_optional("minute", 0); + impl->second = m_params.get_optional("second", 0); + impl->timezone = -6; -double dayfrac_to_local_hr(double dayfrac, double timezone) -{ - return 24.0 * limit_zero2one(dayfrac + timezone / 24.0); -} -double third_order_polynomial(double a, double b, double c, double d, double x) -{ - return ((a * x + b) * x + c) * x + d; -} + impl->longitude = deg_to_rad(m_params.get_optional("longitude", 0)); + impl->latitude = deg_to_rad(m_params.get_optional("latitude", 0)); -/////////////////////////////////////////////////////////////////////////////////////////////// -int validate_inputs(spa_data* spa) -{ - if ((spa->year < -2000) || (spa->year > 6000)) return 1; - if ((spa->month < 1) || (spa->month > 12)) return 2; - if ((spa->day < 1) || (spa->day > 31)) return 3; - if ((spa->hour < 0) || (spa->hour > 24)) return 4; - if ((spa->minute < 0) || (spa->minute > 59)) return 5; - if ((spa->second < 0) || (spa->second >= 60)) return 6; - if ((spa->pressure < 0) || (spa->pressure > 5000)) return 12; - if ((spa->temperature <= -273) || (spa->temperature > 6000)) return 13; - if ((spa->delta_ut1 <= -1) || (spa->delta_ut1 >= 1)) return 17; - if ((spa->hour == 24) && (spa->minute > 0)) return 5; - if ((spa->hour == 24) && (spa->second > 0)) return 6; + impl->azm_rotation = 0.0; - if (fabs(spa->delta_t) > 8000) return 7; - if (fabs(spa->timezone) > 18) return 8; - if (fabs(spa->longitude) > 180) return 9; - if (fabs(spa->latitude) > 90) return 10; - if (fabs(spa->atmos_refract) > 5) return 16; - if (spa->elevation < -6500000) return 11; + return true; + } - if ((spa->function == SPA_ZA_INC) || (spa->function == SPA_ALL)) + inline double SunPositioner::mean_obliquity_ecliptic(const double jc) const { - if (fabs(spa->slope) > 360) return 14; - if (fabs(spa->azm_rotation) > 360) return 15; + return 23.0 + (26.0 + (21.448 - jc * (46.815 + jc * (0.00059 - jc * 0.001813))) / 60.0) / 60.0; } - return 0; -} -/////////////////////////////////////////////////////////////////////////////////////////////// -double julian_day(int year, int month, int day, int hour, int minute, double second, double dut1, double tz) -{ - double day_decimal, julian_day, a; - - day_decimal = day + (hour - tz + (minute + (second + dut1) / 60.0) / 60.0) / 24.0; - - if (month < 3) { - month += 12; - year--; + inline double SunPositioner::obliquity_correction(const double jc) const + { + return mean_obliquity_ecliptic(jc) + 0.00256 * std::cos(deg_to_rad(125.04 - 1934.136 * jc)); } - julian_day = integer(365.25 * (year + 4716.0)) + integer(30.6001 * (month + 1)) + day_decimal - 1524.5; - - if (julian_day > 2299160.0) { - a = integer(year / 100); - julian_day += (2 - a + integer(a / 4)); + inline double SunPositioner::mean_sun_anomaly(const double jc) const + { + return 357.52911 + jc * (35999.05029 - 0.0001537 * jc); } - return julian_day; -} - -double julian_century(double jd) -{ - return (jd - 2451545.0) / 36525.0; -} - -double julian_ephemeris_day(double jd, double delta_t) -{ - return jd + delta_t / 86400.0; -} - -double julian_ephemeris_century(double jde) -{ - return (jde - 2451545.0) / 36525.0; -} - -double julian_ephemeris_millennium(double jce) -{ - return (jce / 10.0); -} - -double earth_periodic_term_summation(const double terms[][TERM_COUNT], int count, double jme) -{ - int i; - double sum = 0; - - for (i = 0; i < count; i++) - sum += terms[i][TERM_A] * cos(terms[i][TERM_B] + terms[i][TERM_C] * jme); - - return sum; -} - -double earth_values(double term_sum[], int count, double jme) -{ - int i; - double sum = 0; - - for (i = 0; i < count; i++) - sum += term_sum[i] * pow(jme, i); - - sum /= 1.0e8; - - return sum; -} - -double earth_heliocentric_longitude(double jme) -{ - double sum[L_COUNT]; - int i; - - for (i = 0; i < L_COUNT; i++) - sum[i] = earth_periodic_term_summation(L_TERMS[i], l_subcount[i], jme); - - return limit_degrees(rad2deg(earth_values(sum, L_COUNT, jme))); - -} - -double earth_heliocentric_latitude(double jme) -{ - double sum[B_COUNT]; - int i; - - for (i = 0; i < B_COUNT; i++) - sum[i] = earth_periodic_term_summation(B_TERMS[i], b_subcount[i], jme); - - return rad2deg(earth_values(sum, B_COUNT, jme)); - -} - -double earth_radius_vector(double jme) -{ - double sum[R_COUNT]; - int i; - - for (i = 0; i < R_COUNT; i++) - sum[i] = earth_periodic_term_summation(R_TERMS[i], r_subcount[i], jme); - - return earth_values(sum, R_COUNT, jme); - -} - -double geocentric_longitude(double l) -{ - double theta = l + 180.0; - - if (theta >= 360.0) theta -= 360.0; - - return theta; -} - -double geocentric_latitude(double b) -{ - return -b; -} - -double mean_elongation_moon_sun(double jce) -{ - return third_order_polynomial(1.0 / 189474.0, -0.0019142, 445267.11148, 297.85036, jce); -} - -double mean_anomaly_sun(double jce) -{ - return third_order_polynomial(-1.0 / 300000.0, -0.0001603, 35999.05034, 357.52772, jce); -} - -double mean_anomaly_moon(double jce) -{ - return third_order_polynomial(1.0 / 56250.0, 0.0086972, 477198.867398, 134.96298, jce); -} - -double argument_latitude_moon(double jce) -{ - return third_order_polynomial(1.0 / 327270.0, -0.0036825, 483202.017538, 93.27191, jce); -} - -double ascending_longitude_moon(double jce) -{ - return third_order_polynomial(1.0 / 450000.0, 0.0020708, -1934.136261, 125.04452, jce); -} - -double xy_term_summation(int i, double x[TERM_X_COUNT]) -{ - int j; - double sum = 0; - - for (j = 0; j < TERM_Y_COUNT; j++) - sum += x[j] * Y_TERMS[i][j]; - - return sum; -} - -void nutation_longitude_and_obliquity(double jce, double x[TERM_X_COUNT], double* del_psi, - double* del_epsilon) -{ - int i; - double xy_term_sum, sum_psi = 0, sum_epsilon = 0; - - for (i = 0; i < Y_COUNT; i++) { - xy_term_sum = deg2rad(xy_term_summation(i, x)); - sum_psi += (PE_TERMS[i][TERM_PSI_A] + jce * PE_TERMS[i][TERM_PSI_B]) * sin(xy_term_sum); - sum_epsilon += (PE_TERMS[i][TERM_EPS_C] + jce * PE_TERMS[i][TERM_EPS_D]) * cos(xy_term_sum); + inline double SunPositioner::sun_eqc_of_center(const double jc, const double msa) const + { + return std::sin(deg_to_rad(msa))* (1.914602 - jc * (0.004817 + 0.000014 * jc)) + + std::sin(deg_to_rad(2 * msa)) * (0.019993 - 0.000101 * jc) + + std::sin(deg_to_rad(3 * msa)) * 0.000289; } - *del_psi = sum_psi / 36000000.0; - *del_epsilon = sum_epsilon / 36000000.0; -} - -double ecliptic_mean_obliquity(double jme) -{ - double u = jme / 10.0; - - return 84381.448 + u * (-4680.93 + u * (-1.55 + u * (1999.25 + u * (-51.38 + u * (-249.67 + - u * (-39.05 + u * (7.12 + u * (27.87 + u * (5.79 + u * 2.45))))))))); -} - -double ecliptic_true_obliquity(double delta_epsilon, double epsilon0) -{ - return delta_epsilon + epsilon0 / 3600.0; -} - -double aberration_correction(double r) -{ - return -20.4898 / (3600.0 * r); -} - -double apparent_sun_longitude(double theta, double delta_psi, double delta_tau) -{ - return theta + delta_psi + delta_tau; -} - -double greenwich_mean_sidereal_time(double jd, double jc) -{ - return limit_degrees(280.46061837 + 360.98564736629 * (jd - 2451545.0) + - jc * jc * (0.000387933 - jc / 38710000.0)); -} - -double greenwich_sidereal_time(double nu0, double delta_psi, double epsilon) -{ - return nu0 + delta_psi * cos(deg2rad(epsilon)); -} - -double geocentric_right_ascension(double lamda, double epsilon, double beta) -{ - double lamda_rad = deg2rad(lamda); - double epsilon_rad = deg2rad(epsilon); - - return limit_degrees(rad2deg(atan2(sin(lamda_rad) * cos(epsilon_rad) - - tan(deg2rad(beta)) * sin(epsilon_rad), cos(lamda_rad)))); -} - -double geocentric_declination(double beta, double epsilon, double lamda) -{ - double beta_rad = deg2rad(beta); - double epsilon_rad = deg2rad(epsilon); - - return rad2deg(asin(sin(beta_rad) * cos(epsilon_rad) + - cos(beta_rad) * sin(epsilon_rad) * sin(deg2rad(lamda)))); -} - -double observer_hour_angle(double nu, double longitude, double alpha_deg) -{ - return limit_degrees(nu + longitude - alpha_deg); -} - -double sun_equatorial_horizontal_parallax(double r) -{ - return 8.794 / (3600.0 * r); -} - -void right_ascension_parallax_and_topocentric_dec(double latitude, double elevation, - double xi, double h, double delta, double* delta_alpha, double* delta_prime) -{ - double delta_alpha_rad; - double lat_rad = deg2rad(latitude); - double xi_rad = deg2rad(xi); - double h_rad = deg2rad(h); - double delta_rad = deg2rad(delta); - double u = atan(0.99664719 * tan(lat_rad)); - double y = 0.99664719 * sin(u) + elevation * sin(lat_rad) / 6378140.0; - double x = cos(u) + elevation * cos(lat_rad) / 6378140.0; - - delta_alpha_rad = atan2(-x * sin(xi_rad) * sin(h_rad), - cos(delta_rad) - x * sin(xi_rad) * cos(h_rad)); - - *delta_prime = rad2deg(atan2((sin(delta_rad) - y * sin(xi_rad)) * cos(delta_alpha_rad), - cos(delta_rad) - x * sin(xi_rad) * cos(h_rad))); - - *delta_alpha = rad2deg(delta_alpha_rad); -} - -double topocentric_right_ascension(double alpha_deg, double delta_alpha) -{ - return alpha_deg + delta_alpha; -} - -double topocentric_local_hour_angle(double h, double delta_alpha) -{ - return h - delta_alpha; -} - -double topocentric_elevation_angle(double latitude, double delta_prime, double h_prime) -{ - double lat_rad = deg2rad(latitude); - double delta_prime_rad = deg2rad(delta_prime); - - return rad2deg(asin(sin(lat_rad) * sin(delta_prime_rad) + - cos(lat_rad) * cos(delta_prime_rad) * cos(deg2rad(h_prime)))); -} - -double atmospheric_refraction_correction(double pressure, double temperature, - double atmos_refract, double e0) -{ - double del_e = 0; - - if (e0 >= -1 * (SUN_RADIUS + atmos_refract)) - del_e = (pressure / 1010.0) * (283.0 / (273.0 + temperature)) * - 1.02 / (60.0 * tan(deg2rad(e0 + 10.3 / (e0 + 5.11)))); - - return del_e; -} - -double topocentric_elevation_angle_corrected(double e0, double delta_e) -{ - return e0 + delta_e; -} - -double topocentric_zenith_angle(double e) -{ - return 90.0 - e; -} - -double topocentric_azimuth_angle_astro(double h_prime, double latitude, double delta_prime) -{ - double h_prime_rad = deg2rad(h_prime); - double lat_rad = deg2rad(latitude); - - return limit_degrees(rad2deg(atan2(sin(h_prime_rad), - cos(h_prime_rad) * sin(lat_rad) - tan(deg2rad(delta_prime)) * cos(lat_rad)))); -} - -double topocentric_azimuth_angle(double azimuth_astro) -{ - return limit_degrees(azimuth_astro + 180.0); -} - -double surface_incidence_angle(double zenith, double azimuth_astro, double azm_rotation, - double slope) -{ - double zenith_rad = deg2rad(zenith); - double slope_rad = deg2rad(slope); - - return rad2deg(acos(cos(zenith_rad) * cos(slope_rad) + - sin(slope_rad) * sin(zenith_rad) * cos(deg2rad(azimuth_astro - azm_rotation)))); -} - -double sun_mean_longitude(double jme) -{ - return limit_degrees(280.4664567 + jme * (360007.6982779 + jme * (0.03032028 + - jme * (1 / 49931.0 + jme * (-1 / 15300.0 + jme * (-1 / 2000000.0)))))); -} - -double eot(double m, double alpha, double del_psi, double epsilon) -{ - return limit_minutes(4.0 * (m - 0.0057183 - alpha + del_psi * cos(deg2rad(epsilon)))); -} - -double approx_sun_transit_time(double alpha_zero, double longitude, double nu) -{ - return (alpha_zero - longitude - nu) / 360.0; -} - -double sun_hour_angle_at_rise_set(double latitude, double delta_zero, double h0_prime) -{ - double h0 = -99999; - double latitude_rad = deg2rad(latitude); - double delta_zero_rad = deg2rad(delta_zero); - double argument = (sin(deg2rad(h0_prime)) - sin(latitude_rad) * sin(delta_zero_rad)) / - (cos(latitude_rad) * cos(delta_zero_rad)); - - if (fabs(argument) <= 1) h0 = limit_degrees180(rad2deg(acos(argument))); - - return h0; -} - -void approx_sun_rise_and_set(double* m_rts, double h0) -{ - double h0_dfrac = h0 / 360.0; - - m_rts[SUN_RISE] = limit_zero2one(m_rts[SUN_TRANSIT] - h0_dfrac); - m_rts[SUN_SET] = limit_zero2one(m_rts[SUN_TRANSIT] + h0_dfrac); - m_rts[SUN_TRANSIT] = limit_zero2one(m_rts[SUN_TRANSIT]); -} - -double rts_alpha_delta_prime(double* ad, double n) -{ - double a = ad[JD_ZERO] - ad[JD_MINUS]; - double b = ad[JD_PLUS] - ad[JD_ZERO]; - - if (fabs(a) >= 2.0) a = limit_zero2one(a); - if (fabs(b) >= 2.0) b = limit_zero2one(b); - - return ad[JD_ZERO] + n * (a + b + (b - a) * n) / 2.0; -} - -double rts_sun_altitude(double latitude, double delta_prime, double h_prime) -{ - double latitude_rad = deg2rad(latitude); - double delta_prime_rad = deg2rad(delta_prime); - - return rad2deg(asin(sin(latitude_rad) * sin(delta_prime_rad) + - cos(latitude_rad) * cos(delta_prime_rad) * cos(deg2rad(h_prime)))); -} - -double sun_rise_and_set(double* m_rts, double* h_rts, double* delta_prime, double latitude, - double* h_prime, double h0_prime, int sun) -{ - return m_rts[sun] + (h_rts[sun] - h0_prime) / - (360.0 * cos(deg2rad(delta_prime[sun])) * cos(deg2rad(latitude)) * sin(deg2rad(h_prime[sun]))); -} - -//////////////////////////////////////////////////////////////////////////////////////////////// -// Calculate required SPA parameters to get the right ascension (alpha) and declination (delta) -// Note: JD must be already calculated and in structure -//////////////////////////////////////////////////////////////////////////////////////////////// -void calculate_geocentric_sun_right_ascension_and_declination(spa_data* spa) -{ - double x[TERM_X_COUNT]; - - spa->jc = julian_century(spa->jd); - - spa->jde = julian_ephemeris_day(spa->jd, spa->delta_t); - spa->jce = julian_ephemeris_century(spa->jde); - spa->jme = julian_ephemeris_millennium(spa->jce); - - spa->l = earth_heliocentric_longitude(spa->jme); - spa->b = earth_heliocentric_latitude(spa->jme); - spa->r = earth_radius_vector(spa->jme); - - spa->theta = geocentric_longitude(spa->l); - spa->beta = geocentric_latitude(spa->b); - - x[TERM_X0] = spa->x0 = mean_elongation_moon_sun(spa->jce); - x[TERM_X1] = spa->x1 = mean_anomaly_sun(spa->jce); - x[TERM_X2] = spa->x2 = mean_anomaly_moon(spa->jce); - x[TERM_X3] = spa->x3 = argument_latitude_moon(spa->jce); - x[TERM_X4] = spa->x4 = ascending_longitude_moon(spa->jce); - - nutation_longitude_and_obliquity(spa->jce, x, &(spa->del_psi), &(spa->del_epsilon)); - - spa->epsilon0 = ecliptic_mean_obliquity(spa->jme); - spa->epsilon = ecliptic_true_obliquity(spa->del_epsilon, spa->epsilon0); - - spa->del_tau = aberration_correction(spa->r); - spa->lamda = apparent_sun_longitude(spa->theta, spa->del_psi, spa->del_tau); - spa->nu0 = greenwich_mean_sidereal_time(spa->jd, spa->jc); - spa->nu = greenwich_sidereal_time(spa->nu0, spa->del_psi, spa->epsilon); - - spa->alpha = geocentric_right_ascension(spa->lamda, spa->epsilon, spa->beta); - spa->delta = geocentric_declination(spa->beta, spa->epsilon, spa->lamda); -} - -//////////////////////////////////////////////////////////////////////// -// Calculate Equation of Time (EOT) and Sun Rise, Transit, & Set (RTS) -//////////////////////////////////////////////////////////////////////// - -void calculate_eot_and_sun_rise_transit_set(spa_data* spa) -{ - spa_data sun_rts; - double nu, m, h0, n; - double alpha[JD_COUNT], delta[JD_COUNT]; - double m_rts[SUN_COUNT], nu_rts[SUN_COUNT], h_rts[SUN_COUNT]; - double alpha_prime[SUN_COUNT], delta_prime[SUN_COUNT], h_prime[SUN_COUNT]; - double h0_prime = -1 * (SUN_RADIUS + spa->atmos_refract); - int i; - - sun_rts = *spa; - m = sun_mean_longitude(spa->jme); - spa->eot = eot(m, spa->alpha, spa->del_psi, spa->epsilon); - - sun_rts.second = 0.0; - sun_rts.hour = sun_rts.minute = static_cast(sun_rts.second); - sun_rts.delta_ut1 = sun_rts.timezone = 0.0; - - sun_rts.jd = julian_day(sun_rts.year, sun_rts.month, sun_rts.day, sun_rts.hour, - sun_rts.minute, sun_rts.second, sun_rts.delta_ut1, sun_rts.timezone); - - calculate_geocentric_sun_right_ascension_and_declination(&sun_rts); - nu = sun_rts.nu; - - sun_rts.delta_t = 0; - sun_rts.jd--; - for (i = 0; i < JD_COUNT; i++) { - calculate_geocentric_sun_right_ascension_and_declination(&sun_rts); - alpha[i] = sun_rts.alpha; - delta[i] = sun_rts.delta; - sun_rts.jd++; + inline double SunPositioner::mean_longitude_sun(const double jc) const + { + return std::fmod(280.46646 + jc * (36000.76983 + jc * 0.0003032), 360); } - m_rts[SUN_TRANSIT] = approx_sun_transit_time(alpha[JD_ZERO], spa->longitude, nu); - h0 = sun_hour_angle_at_rise_set(spa->latitude, delta[JD_ZERO], h0_prime); - - if (h0 >= 0) { - - approx_sun_rise_and_set(m_rts, h0); - - for (i = 0; i < SUN_COUNT; i++) { - - nu_rts[i] = nu + 360.985647 * m_rts[i]; - - n = m_rts[i] + spa->delta_t / 86400.0; - alpha_prime[i] = rts_alpha_delta_prime(alpha, n); - delta_prime[i] = rts_alpha_delta_prime(delta, n); - - h_prime[i] = limit_degrees180pm(nu_rts[i] + spa->longitude - alpha_prime[i]); - - h_rts[i] = rts_sun_altitude(spa->latitude, delta_prime[i], h_prime[i]); - } - - spa->srha = h_prime[SUN_RISE]; - spa->ssha = h_prime[SUN_SET]; - spa->sta = h_rts[SUN_TRANSIT]; - - spa->suntransit = dayfrac_to_local_hr(m_rts[SUN_TRANSIT] - h_prime[SUN_TRANSIT] / 360.0, - spa->timezone); - - spa->sunrise = dayfrac_to_local_hr(sun_rise_and_set(m_rts, h_rts, delta_prime, - spa->latitude, h_prime, h0_prime, SUN_RISE), spa->timezone); - - spa->sunset = dayfrac_to_local_hr(sun_rise_and_set(m_rts, h_rts, delta_prime, - spa->latitude, h_prime, h0_prime, SUN_SET), spa->timezone); - + inline double SunPositioner::sun_true_longitude(const double jc, const double mean_long_sun, const double mean_sun_ano) const + { + return mean_long_sun + sun_eqc_of_center(jc, mean_sun_ano); } - else spa->srha = spa->ssha = spa->sta = spa->suntransit = spa->sunrise = spa->sunset = -99999; - -} - -/////////////////////////////////////////////////////////////////////////////////////////// -// Calculate all SPA parameters and put into structure -// Note: All inputs values (listed in header file) must already be in structure -/////////////////////////////////////////////////////////////////////////////////////////// -int spa_calculate(spa_data* spa) -{ - int result; - result = validate_inputs(spa); - - if (result == 0) + inline double SunPositioner::apparent_longitude_of_sun(const double jc, const double mean_long_sun, const double mean_sun_ano) const { - spa->jd = julian_day(spa->year, spa->month, spa->day, spa->hour, - spa->minute, spa->second, spa->delta_ut1, spa->timezone); - - calculate_geocentric_sun_right_ascension_and_declination(spa); - - spa->h = observer_hour_angle(spa->nu, spa->longitude, spa->alpha); - spa->xi = sun_equatorial_horizontal_parallax(spa->r); - - right_ascension_parallax_and_topocentric_dec(spa->latitude, spa->elevation, spa->xi, - spa->h, spa->delta, &(spa->del_alpha), &(spa->delta_prime)); + return sun_true_longitude(jc, mean_long_sun, mean_sun_ano) - 0.00569 - 0.00478 * std::sin(deg_to_rad(125.04 - 1934.136 * jc)); + } - spa->alpha_prime = topocentric_right_ascension(spa->alpha, spa->del_alpha); - spa->h_prime = topocentric_local_hour_angle(spa->h, spa->del_alpha); + inline double SunPositioner::sun_declination(const double obl_cor, const double app_lon_sun) const + { + return std::asin(std::sin(deg_to_rad(obl_cor)) * std::sin(deg_to_rad(app_lon_sun))); + } - spa->e0 = topocentric_elevation_angle(spa->latitude, spa->delta_prime, spa->h_prime); - //spa->del_e = atmospheric_refraction_correction(spa->pressure, spa->temperature, - // spa->atmos_refract, spa->e0); - spa->e = topocentric_elevation_angle_corrected(spa->e0, spa->del_e); + inline double SunPositioner::equation_of_time(const double jc, const double obl_cor, const double mean_long_sun, const double mean_sun_ano) const + { + const double eccent_earth_orbit = 0.016708634 - jc * (0.000042037 + 0.0000001267 * jc); + const double y = std::tan(deg_to_rad(obl_cor / 2.0)) * + std::tan(deg_to_rad(obl_cor / 2.0)); - spa->zenith = topocentric_zenith_angle(spa->e); - spa->azimuth_astro = topocentric_azimuth_angle_astro(spa->h_prime, spa->latitude, - spa->delta_prime); - spa->azimuth = topocentric_azimuth_angle(spa->azimuth_astro); + double eq_of_time = y * std::sin(2 * deg_to_rad(mean_long_sun)); + eq_of_time -= 2 * eccent_earth_orbit * std::sin(deg_to_rad(mean_sun_ano)); + eq_of_time += 4 * eccent_earth_orbit * y * std::sin(deg_to_rad(mean_sun_ano)) * std::cos(2 * deg_to_rad(mean_long_sun)); + eq_of_time -= 0.5 * y * y * std::sin(4 * deg_to_rad(mean_long_sun)); + eq_of_time -= 1.25 * eccent_earth_orbit * eccent_earth_orbit * std::sin(2 * deg_to_rad(mean_sun_ano)); - if ((spa->function == SPA_ZA_INC) || (spa->function == SPA_ALL)) - spa->incidence = surface_incidence_angle(spa->zenith, spa->azimuth_astro, - spa->azm_rotation, spa->slope); + return 4 * rad_to_deg(eq_of_time); + } - if ((spa->function == SPA_ZA_RTS) || (spa->function == SPA_ALL)) - calculate_eot_and_sun_rise_transit_set(spa); + void SunPositioner::compute_sun_position() + { + const double decimal_day = (impl->hour - impl->timezone + (impl->minute + (impl->second) / 60.0) / 60.0) / 24.0; + const double jc = get_julian_century(get_julian_day(impl->year, impl->month, impl->day, decimal_day)); + const double obl_cor = obliquity_correction(jc); + const double mean_long_sun = mean_longitude_sun(jc); + const double mean_sun_ano = mean_sun_anomaly(jc); + const double app_lon_sun = apparent_longitude_of_sun(jc, mean_long_sun, mean_sun_ano); + const double sun_dcli = sun_declination(obl_cor, app_lon_sun); + const double solar_time = + std::fmod((decimal_day * 1440 + equation_of_time(jc, obl_cor, mean_long_sun, mean_sun_ano) + 4 * impl->longitude), 1440); + double hour_angle = solar_time / 4; + hour_angle += (hour_angle < 0) ? 180 : -180; + + // Compute zenith and azimuth. + const double sin_lat = std::sin(impl->latitude); + const double cos_lat = std::cos(impl->latitude); + const double sin_sun_dcli = std::sin(sun_dcli); + const double cos_sun_dcli = std::cos(sun_dcli); + + m_output_value.zenith = std::acos(sin_lat * sin_sun_dcli + cos_lat * cos_sun_dcli * std::cos(deg_to_rad(hour_angle))); + + const double sin_zenith = std::sin(m_output_value.zenith); + const double cos_zenith = std::cos(m_output_value.zenith); + + m_output_value.azimuth = rad_to_deg(std::acos(((sin_lat * cos_zenith) - sin_sun_dcli) / (cos_lat * sin_zenith))); + + if (hour_angle > 0) + m_output_value.azimuth += 180; + else + m_output_value.azimuth = 540 - m_output_value.azimuth; + + m_output_value.zenith = std::fabs(90 - rad_to_deg(m_output_value.zenith)); + m_output_value.azimuth = std::fmod(m_output_value.azimuth, 360); } - return result; -} -/////////////////////////////////////////////////////////////////////////////////////////// +} // namespace renderer diff --git a/src/appleseed/renderer/utility/solarpositionalgorithm.h b/src/appleseed/renderer/utility/solarpositionalgorithm.h index af07837500..b45dbaa69f 100644 --- a/src/appleseed/renderer/utility/solarpositionalgorithm.h +++ b/src/appleseed/renderer/utility/solarpositionalgorithm.h @@ -1,207 +1,104 @@ + + +// +// This source file is part of appleseed. +// Visit https://appleseedhq.net/ for additional information and resources. +// +// This software is released under the MIT license. +// +// Copyright (c) 2020 João Marcos , Jupiter Jazz Limited +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +// THE SOFTWARE. +// + +// This sun positioning algorithm is based on the National Oceanic and Atmospheric +// Administration's (NOAA) Solar Position Calculator which uses the calculations +// presented by Jean Meeus in his book Astronomical Algorithms. + +// The Use of NOAA data and products are in the public domain and may be used +// freely by the public as specify in www.nws.noaa.gov/disclaimer.php + #pragma once +#include "renderer/modeling/entity/entity.h" + +// Forward declarations. +namespace foundation { class IAbortSwitch; } +namespace renderer { class BaseGroup; } +namespace renderer { class OnFrameBeginRecorder; } +namespace renderer { class ParamArray; } +namespace renderer { class Project; } +namespace renderer { class ShadingContext; } -///////////////////////////////////////////// -// HEADER FILE for SPA.C // -// // -// Solar Position Algorithm (SPA) // -// for // -// Solar Radiation Application // -// // -// May 12, 2003 // -// // -// Filename: SPA.H // -// // -// Afshin Michael Andreas // -// afshin_andreas@nrel.gov (303)384-6383 // -// // -// Measurement & Instrumentation Team // -// Solar Radiation Research Laboratory // -// National Renewable Energy Laboratory // -// 1617 Cole Blvd, Golden, CO 80401 // -///////////////////////////////////////////// - -//////////////////////////////////////////////////////////////////////// -// // -// Usage: // -// // -// 1) In calling program, include this header file, // -// by adding this line to the top of file: // -// #include "spa.h" // -// // -// 2) In calling program, declare the SPA structure: // -// spa_data spa; // -// // -// 3) Enter the required input values into SPA structure // -// (input values listed in comments below) // -// // -// 4) Call the SPA calculate function and pass the SPA structure // -// (prototype is declared at the end of this header file): // -// spa_calculate(&spa); // -// // -// Selected output values (listed in comments below) will be // -// computed and returned in the passed SPA structure. Output // -// will based on function code selected from enumeration below. // -// // -// Note: A non-zero return code from spa_calculate() indicates that // -// one of the input values did not pass simple bounds tests. // -// The valid input ranges and return error codes are also // -// listed below. // -// // -//////////////////////////////////////////////////////////////////////// - -#ifndef __solar_position_algorithm_header -#define __solar_position_algorithm_header - - -//enumeration for function codes to select desired final outputs from SPA -enum { - SPA_ZA, //calculate zenith and azimuth - SPA_ZA_INC, //calculate zenith, azimuth, and incidence - SPA_ZA_RTS, //calculate zenith, azimuth, and sun rise/transit/set values - SPA_ALL, //calculate all SPA output values -}; - -typedef struct +namespace renderer { - //----------------------INPUT VALUES------------------------ - int year; // 4-digit year, valid range: -2000 to 6000, error code: 1 - int month; // 2-digit month, valid range: 1 to 12, error code: 2 - int day; // 2-digit day, valid range: 1 to 31, error code: 3 - int hour; // Observer local hour, valid range: 0 to 24, error code: 4 - int minute; // Observer local minute, valid range: 0 to 59, error code: 5 - double second; // Observer local second, valid range: 0 to <60, error code: 6 + class APPLESEED_DLLSYMBOL SunPositioner + :public Entity + { + public: + + struct output_values + { + double zenith; // Topocentric zenith angle [degrees] + double azimuth; // Topocentric azimuth angle (eastward from north) [for navigators and solar radiation] + + double solar_noon; + double sunrise; // Local sunrise time (+/- 30 seconds). + double sunset; // Local sunset time (+/- 30 seconds). + } m_output_value; + + // Constructor + SunPositioner( + const char* name, + const ParamArray& params); + + void compute_sun_position(); + + bool on_frame_begin( + const Project& project, + const BaseGroup* parent, + OnFrameBeginRecorder& recorder, + foundation::IAbortSwitch* abort_switch) override; + + void SunPositioner::release() + { + delete this; + } + + private: + struct Impl; + Impl* impl; + + double mean_obliquity_ecliptic(const double jc) const; + double obliquity_correction(const double jc) const; + + double mean_sun_anomaly(const double jc) const; + double mean_longitude_sun(const double jc) const; + double sun_eqc_of_center(const double jc, const double msa) const; + double sun_true_longitude(const double jc, const double mean_long_sun, const double mean_sun_ano) const; + double apparent_longitude_of_sun(const double jc, const double mean_long_sun, const double mean_sun_ano) const; - double delta_ut1; // Fractional second difference between UTC and UT which is used - // to adjust UTC for earth's irregular rotation rate and is derived - // from observation only and is reported in this bulletin: - // http://maia.usno.navy.mil/ser7/ser7.dat, - // where delta_ut1 = DUT1 - // valid range: -1 to 1 second (exclusive), error code 17 + double sun_declination(const double obj_cor, const double app_lon_sun) const; - double delta_t; // Difference between earth rotation time and terrestrial time - // It is derived from observation only and is reported in this - // bulletin: http://maia.usno.navy.mil/ser7/ser7.dat, - // where delta_t = 32.184 + (TAI-UTC) - DUT1 - // valid range: -8000 to 8000 seconds, error code: 7 + double equation_of_time(const double jc, const double obl_cor, const double mean_long_sun, const double mean_sun_ano) const; - double timezone; // Observer time zone (negative west of Greenwich) - // valid range: -18 to 18 hours, error code: 8 - - double longitude; // Observer longitude (negative west of Greenwich) - // valid range: -180 to 180 degrees, error code: 9 - - double latitude; // Observer latitude (negative south of equator) - // valid range: -90 to 90 degrees, error code: 10 - - double elevation; // Observer elevation [meters] - // valid range: -6500000 or higher meters, error code: 11 - - double pressure; // Annual average local pressure [millibars] - // valid range: 0 to 5000 millibars, error code: 12 - - double temperature; // Annual average local temperature [degrees Celsius] - // valid range: -273 to 6000 degrees Celsius, error code; 13 - - double slope; // Surface slope (measured from the horizontal plane) - // valid range: -360 to 360 degrees, error code: 14 - - double azm_rotation; // Surface azimuth rotation (measured from south to projection of - // surface normal on horizontal plane, negative east) - // valid range: -360 to 360 degrees, error code: 15 - - double atmos_refract;// Atmospheric refraction at sunrise and sunset (0.5667 deg is typical) - // valid range: -5 to 5 degrees, error code: 16 - - int function; // Switch to choose functions for desired output (from enumeration) - - //-----------------Intermediate OUTPUT VALUES-------------------- - - double jd; //Julian day - double jc; //Julian century - - double jde; //Julian ephemeris day - double jce; //Julian ephemeris century - double jme; //Julian ephemeris millennium - - double l; //earth heliocentric longitude [degrees] - double b; //earth heliocentric latitude [degrees] - double r; //earth radius vector [Astronomical Units, AU] - - double theta; //geocentric longitude [degrees] - double beta; //geocentric latitude [degrees] - - double x0; //mean elongation (moon-sun) [degrees] - double x1; //mean anomaly (sun) [degrees] - double x2; //mean anomaly (moon) [degrees] - double x3; //argument latitude (moon) [degrees] - double x4; //ascending longitude (moon) [degrees] - - double del_psi; //nutation longitude [degrees] - double del_epsilon; //nutation obliquity [degrees] - double epsilon0; //ecliptic mean obliquity [arc seconds] - double epsilon; //ecliptic true obliquity [degrees] - - double del_tau; //aberration correction [degrees] - double lamda; //apparent sun longitude [degrees] - double nu0; //Greenwich mean sidereal time [degrees] - double nu; //Greenwich sidereal time [degrees] - - double alpha; //geocentric sun right ascension [degrees] - double delta; //geocentric sun declination [degrees] - - double h; //observer hour angle [degrees] - double xi; //sun equatorial horizontal parallax [degrees] - double del_alpha; //sun right ascension parallax [degrees] - double delta_prime; //topocentric sun declination [degrees] - double alpha_prime; //topocentric sun right ascension [degrees] - double h_prime; //topocentric local hour angle [degrees] - - double e0; //topocentric elevation angle (uncorrected) [degrees] - double del_e; //atmospheric refraction correction [degrees] - double e; //topocentric elevation angle (corrected) [degrees] - - double eot; //equation of time [minutes] - double srha; //sunrise hour angle [degrees] - double ssha; //sunset hour angle [degrees] - double sta; //sun transit altitude [degrees] - - //---------------------Final OUTPUT VALUES------------------------ - - double zenith; //topocentric zenith angle [degrees] - double azimuth_astro;//topocentric azimuth angle (westward from south) [for astronomers] - double azimuth; //topocentric azimuth angle (eastward from north) [for navigators and solar radiation] - double incidence; //surface incidence angle [degrees] - - double suntransit; //local sun transit time (or solar noon) [fractional hour] - double sunrise; //local sunrise time (+/- 30 seconds) [fractional hour] - double sunset; //local sunset time (+/- 30 seconds) [fractional hour] - -} spa_data; + }; -//-------------- Utility functions for other applications (such as NREL's SAMPA) -------------- -double deg2rad(double degrees); -double rad2deg(double radians); -double limit_degrees(double degrees); -double third_order_polynomial(double a, double b, double c, double d, double x); -double geocentric_right_ascension(double lamda, double epsilon, double beta); -double geocentric_declination(double beta, double epsilon, double lamda); -double observer_hour_angle(double nu, double longitude, double alpha_deg); -void right_ascension_parallax_and_topocentric_dec(double latitude, double elevation, - double xi, double h, double delta, double* delta_alpha, double* delta_prime); -double topocentric_right_ascension(double alpha_deg, double delta_alpha); -double topocentric_local_hour_angle(double h, double delta_alpha); -double topocentric_elevation_angle(double latitude, double delta_prime, double h_prime); -double atmospheric_refraction_correction(double pressure, double temperature, - double atmos_refract, double e0); -double topocentric_elevation_angle_corrected(double e0, double delta_e); -double topocentric_zenith_angle(double e); -double topocentric_azimuth_angle_astro(double h_prime, double latitude, double delta_prime); -double topocentric_azimuth_angle(double azimuth_astro); - - -//Calculate SPA output values (in structure) based on input values passed in structure -int spa_calculate(spa_data* spa); - -#endif +} // namespace renderer