Skip to content

Commit

Permalink
make spa threadsafe again (#1232)
Browse files Browse the repository at this point in the history
  • Loading branch information
@dguittet
dguittet authored Oct 30, 2024
1 parent 89554fc commit 77b85d4
Show file tree
Hide file tree
Showing 3 changed files with 151 additions and 115 deletions.
157 changes: 102 additions & 55 deletions shared/lib_irradproc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,10 +49,6 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include "lib_weatherfile.h"

static const int __nday[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};

std::unordered_map<spa_table_key, std::vector<double>> spa_table;
int spa_table_day;

/// Compute the Julian day of year
static int julian(int yr, int month, int day) {
int i = 1, jday = 0, k;
Expand Down Expand Up @@ -803,47 +799,96 @@ double sun_rise_and_set(double *m_rts, double *h_rts, double *delta_prime, doubl
(360.0 * cos(DTOR * (delta_prime[sun])) * cos(DTOR * (latitude)) * sin(DTOR * (h_prime[sun])));
}

void clear_spa_table() {
bool spa_table_key::operator==(const spa_table_key &other) const {
return (jd == other.jd
&& delta_t == other.delta_t
&& pressure == other.pressure
&& temp == other.temp
&& ascension == other.ascension
&& declination == other.declination
);
}

spa_table_key::spa_table_key(double j, double dt, double p, double t, double a, double d):
jd(j), delta_t(dt), ascension(a), declination(d) {
int pressure_bucket = 10;
pressure = ((int)(p + pressure_bucket/2) / pressure_bucket) * pressure_bucket;
pressure = (int)pressure;
int temp_bucket = 5;
temp = ((int)(t + temp_bucket / 2) / temp_bucket) * temp_bucket;
temp = (int)temp;
}

std::size_t std::hash<spa_table_key>::operator()(const spa_table_key& k) const
{
using std::hash;
// Compute individual hash values for first, second, etc and combine them using XOR and bit shifting:
return
((((
((((hash<double>()(k.jd)
^ (hash<double>()(k.delta_t) << 1)) >> 1)
^ (hash<int>()(k.pressure) << 1)) >> 1)
^ (hash<int>()(k.temp) << 1)) >> 1)
^ (hash<double>()(k.ascension) << 1)) >> 1)
^ (hash<double>()(k.declination) << 1)
;
}

solarpos_lookup::solarpos_lookup(){
clear_spa_table();
}

void solarpos_lookup::clear_spa_table() {
spa_table.clear();
spa_table_day = 0;
};

// The algorithm reuses the outputs from the last 3 days or so, so the hash table is emptied every 3 days to reduce size
void roll_spa_table_forward(int day) {
void solarpos_lookup::roll_spa_table_forward(int day) {
if (std::abs(spa_table_day - day) > 3){
spa_table_day = day;
spa_table.clear();
}
};

std::vector<double>* solarpos_lookup::find(spa_table_key spa_key_inputs) {
auto spa_pos = spa_table.end();
spa_pos = spa_table.find(spa_key_inputs);
if (spa_pos != spa_table.end())
return &spa_pos->second;
return nullptr;
}

void solarpos_lookup::insert(spa_table_key spa_key_inputs, std::vector<double> spa_outputs) {
spa_table[spa_key_inputs] = spa_outputs;
}

void
calculate_spa(double jd, double lat, double lng, double alt, double pressure, double temp, double delta_t, double tilt,
double azm_rotation, double ascension_and_declination[2], double needed_values[9]) {
double azm_rotation, double ascension_and_declination[2], double needed_values[9], std::shared_ptr<solarpos_lookup> spa_table) {
//Calculate the Julian and Julian Ephemeris, Day, Century, and Millennium (3.1)
//double jd = julian_day(year, month, day, hour, minute, second, delta_t, tz);
double jc = julian_century(jd); // for 2000 standard epoch
double jde = julian_ephemeris_day(jd,
delta_t); //Adjusted for difference between Earth rotation time and the Terrestrial Time (TT) (derived from observation, reported yearly in Astronomical Almanac)

bool use_table = true;
spa_table_key spa_key_inputs(jd, delta_t, pressure, temp, ascension_and_declination[0], ascension_and_declination[1]);
auto spa_pos = spa_table.end();
if (use_table)
spa_pos = spa_table.find(spa_key_inputs);

if (spa_pos != spa_table.end()){
needed_values[0] = spa_pos->second[0];
needed_values[1] = spa_pos->second[1];
needed_values[2] = spa_pos->second[2];
needed_values[3] = spa_pos->second[3];
needed_values[4] = spa_pos->second[4];
needed_values[5] = spa_pos->second[5];
needed_values[6] = spa_pos->second[6];
needed_values[7] = spa_pos->second[7];
needed_values[8] = spa_pos->second[8];
ascension_and_declination[0] = spa_pos->second[9];
ascension_and_declination[1] = spa_pos->second[10];
return;
spa_table_key spa_key_inputs = {jd, delta_t, pressure, temp, ascension_and_declination[0], ascension_and_declination[1]};
if (spa_table) {
std::vector<double>* results = spa_table->find(spa_key_inputs);
if (results){
needed_values[0] = (*results)[0];
needed_values[1] = (*results)[1];
needed_values[2] = (*results)[2];
needed_values[3] = (*results)[3];
needed_values[4] = (*results)[4];
needed_values[5] = (*results)[5];
needed_values[6] = (*results)[6];
needed_values[7] = (*results)[7];
needed_values[8] = (*results)[8];
ascension_and_declination[0] = (*results)[9];
ascension_and_declination[1] = (*results)[10];
return;
}
}

double jce = julian_ephemeris_century(jde); //for 2000 standard epoch
Expand Down Expand Up @@ -942,11 +987,11 @@ calculate_spa(double jd, double lat, double lng, double alt, double pressure, do
azimuth = M_PI;
}

std::vector<double> spa_outputs = {needed_values[0], needed_values[1], needed_values[2], needed_values[3], needed_values[4], needed_values[5], needed_values[6], needed_values[7], needed_values[8],
ascension_and_declination[0], ascension_and_declination[1]};

if (use_table)
spa_table[spa_key_inputs] = spa_outputs;
if (spa_table) {
std::vector<double> spa_outputs = {needed_values[0], needed_values[1], needed_values[2], needed_values[3], needed_values[4], needed_values[5], needed_values[6], needed_values[7], needed_values[8],
ascension_and_declination[0], ascension_and_declination[1]};
spa_table->insert(spa_key_inputs, spa_outputs);
}

//Calculate the incidence angle for a selected surface (3.16)
//double aoi = surface_incidence_angle(zenith, azimuth_astro, azm_rotation, tilt); //incidence angle for a surface oriented in any direction (degrees)
Expand All @@ -957,7 +1002,7 @@ void
calculate_eot_and_sun_rise_transit_set(double jme, double tz, double alpha, double del_psi, double epsilon, double jd,
int year, int month, int day, double lat, double lng, double alt,
double pressure, double temp, double tilt, double delta_t, double azm_rotation,
double needed_values[4]) {
double needed_values[4], std::shared_ptr<solarpos_lookup> spa_table) {
// Equation of Time (A.1)
double M = sun_mean_longitude(jme); // degrees
double E = eot(M, alpha, del_psi, epsilon); //Equation of Time (degrees)
Expand All @@ -966,7 +1011,7 @@ calculate_eot_and_sun_rise_transit_set(double jme, double tz, double alpha, doub

// Sunrise, Sunset, and Sun Transit (A.2)
int i; //initialize iterator
double sun_declination_and_ascension[2]; // storage for iterating sun declination and ascension results
double sun_declination_and_ascension[2] {0, 0}; // storage for iterating sun declination and ascension results
double needed_values2[13]; //where is this used?
double alpha_array[3]; //storage for alphas in iteration
double delta_array[3]; //storage for deltas in iteration
Expand All @@ -977,15 +1022,15 @@ calculate_eot_and_sun_rise_transit_set(double jme, double tz, double alpha, doub
0); //initialize Julian array with day in question (0 UT (0h0mm0s0tz)
//run calculate_spa to obtain apparent sidereal time at Greenwich at 0 UT (v, degrees)
calculate_spa(jd_array[1], lat, lng, alt, pressure, temp, 67, tilt, azm_rotation, sun_declination_and_ascension,
needed_values_nu);
needed_values_nu, spa_table);
double delta_test = sun_declination_and_ascension[1];
double nu = needed_values_nu[4]; //store apparent sidereal time at Greenwich
jd_array[0] = jd_array[1] - 1; //Julian day for the day prior to the day in question
jd_array[2] = jd_array[1] + 1; //Julian day for the day following the day in question
for (i = 0; i < 3; i++) { // iterate through day behind (0), day in question (1), and day following (2)
//Calculate the spa for each day at 0 TT by setting the delta_T difference between UT and TT to 0
calculate_spa(jd_array[i], lat, lng, alt, pressure, temp, 0, tilt, azm_rotation, sun_declination_and_ascension,
needed_values2);
needed_values2, spa_table);
alpha_array[i] = sun_declination_and_ascension[0]; //store geocentric right ascension for each iteration (degrees)
delta_array[i] = sun_declination_and_ascension[1]; //store geocentric declination for each iteration (degrees)
}
Expand Down Expand Up @@ -1075,7 +1120,8 @@ calculate_eot_and_sun_rise_transit_set(double jme, double tz, double alpha, doub

void
solarpos_spa(int year, int month, int day, int hour, double minute, double second, double lat, double lng, double tz,
double dut1, double alt, double pressure, double temp, double tilt, double azm_rotation, double sunn[9]) {
double dut1, double alt, double pressure, double temp, double tilt, double azm_rotation, double sunn[9],
std::shared_ptr<solarpos_lookup> spa_table) {

int t;
double delta_t;
Expand All @@ -1101,13 +1147,14 @@ solarpos_spa(int year, int month, int day, int hour, double minute, double secon
double needed_values_eot[4]; //preallocate storage for output from calculate_spa
double needed_values_eot_check[4];

roll_spa_table_forward(day);
if (spa_table){
spa_table->roll_spa_table_forward(day);
}
calculate_spa(jd, lat, lng, alt, pressure, temp, delta_t, tilt, azm_rotation, ascension_and_declination,
needed_values_spa); //calculate solar position algorithm values
needed_values_spa, spa_table); //calculate solar position algorithm values
calculate_eot_and_sun_rise_transit_set(needed_values_spa[0], tz, ascension_and_declination[0], needed_values_spa[2],
needed_values_spa[3], jd, year, month, day, lat, lng, alt, pressure, temp,
tilt, delta_t, azm_rotation,
needed_values_eot); //calculate Equation of Time and sunrise/sunset values
tilt, delta_t, azm_rotation, needed_values_eot, spa_table); //calculate Equation of Time and sunrise/sunset values

double __n_days[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};

Expand All @@ -1119,17 +1166,17 @@ solarpos_spa(int year, int month, int day, int hour, double minute, double secon
calculate_eot_and_sun_rise_transit_set(needed_values_spa[0], tz, ascension_and_declination[0],
needed_values_spa[2], needed_values_spa[3], jd, year, month, day + 1,
lat, lng, alt, pressure, temp, tilt, delta_t, azm_rotation,
needed_values_eot_check); //calculate Equation of Time and sunrise/sunset values
needed_values_eot_check, spa_table); //calculate Equation of Time and sunrise/sunset values
else if (month < 12) //on the 1st of the month, need to switch to the last day of previous month
calculate_eot_and_sun_rise_transit_set(needed_values_spa[0], tz, ascension_and_declination[0],
needed_values_spa[2], needed_values_spa[3], jd, year, month + 1, 1,
lat, lng, alt, pressure, temp, tilt, delta_t, azm_rotation,
needed_values_eot_check);
needed_values_eot_check, spa_table);
else //on the first day of the year, need to switch to Dec 31 of last year
calculate_eot_and_sun_rise_transit_set(needed_values_spa[0], tz, ascension_and_declination[0],
needed_values_spa[2], needed_values_spa[3], jd, year + 1, 1, 1, lat,
lng, alt, pressure, temp, tilt, delta_t, azm_rotation,
needed_values_eot_check);
needed_values_eot_check, spa_table);

//on the last day of endless days, sunset is returned as 100 (hour angle too large for calculation), so use today's sunset time as a proxy
needed_values_eot[3] = needed_values_eot_check[3] + 24;
Expand Down Expand Up @@ -1817,7 +1864,7 @@ void irrad::setup() {

irrad::irrad() {
setup();
clear_spa_table();
spa_table = std::make_shared<solarpos_lookup>(solarpos_lookup());
}

irrad::irrad(weather_header hdr,
Expand Down Expand Up @@ -2240,7 +2287,7 @@ int irrad::calc() {
if (!getStoredSolarposOutputs()) {

// calculate sunrise and sunset hours in local standard time for the current day
solarpos_spa(year, month, day, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunAnglesRadians);
solarpos_spa(year, month, day, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunAnglesRadians, spa_table);

double t_sunrise = sunAnglesRadians[4];
double t_sunset = sunAnglesRadians[5];
Expand All @@ -2250,11 +2297,11 @@ int irrad::calc() {
{
double sunanglestemp[9];
if (day > 1) //simply decrement day during month
solarpos_spa(year, month, day - 1, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp);
solarpos_spa(year, month, day - 1, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp, spa_table);
else if (month > 1) //on the 1st of the month, need to switch to the last day of previous month
solarpos_spa(year, month - 1, __nday[month - 2], 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp);
solarpos_spa(year, month - 1, __nday[month - 2], 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp, spa_table);
else //on the first day of the year, need to switch to Dec 31 of last year
solarpos_spa(year - 1, 12, 31, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp);
solarpos_spa(year - 1, 12, 31, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp, spa_table);
//on the last day of endless days, sunset is returned as 100 (hour angle too large for calculation), so use today's sunset time as a proxy
if (sunanglestemp[5] == 100.0)
t_sunset -= 24.0;
Expand All @@ -2268,11 +2315,11 @@ int irrad::calc() {
{
double sunanglestemp[9];
if (day < __nday[month - 1]) //simply increment the day during the month, month is 1-indexed and __nday is 0-indexed
solarpos_spa(year, month, day + 1, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp);
solarpos_spa(year, month, day + 1, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp, spa_table);
else if (month < 12) //on the last day of the month, need to switch to the first day of the next month
solarpos_spa(year, month + 1, 1, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp);
solarpos_spa(year, month + 1, 1, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp, spa_table);
else //on the last day of the year, need to switch to Jan 1 of the next year
solarpos_spa(year + 1, 1, 1, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp);
solarpos_spa(year + 1, 1, 1, 12, 0.0, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunanglestemp, spa_table);
//on the last day of endless days, sunrise would be returned as -100 (hour angle too large for calculations), so use today's sunrise time as a proxy
if (sunanglestemp[4] == -100.0)
t_sunrise += 24.0;
Expand All @@ -2291,7 +2338,7 @@ int irrad::calc() {
timeStepSunPosition[0] = hr_calc;
timeStepSunPosition[1] = (int) min_calc;

solarpos_spa(year, month, day, hr_calc, min_calc, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunAnglesRadians);
solarpos_spa(year, month, day, hr_calc, min_calc, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunAnglesRadians, spa_table);

timeStepSunPosition[2] = 2;
}
Expand All @@ -2304,7 +2351,7 @@ int irrad::calc() {
timeStepSunPosition[0] = hr_calc;
timeStepSunPosition[1] = (int) min_calc;

solarpos_spa(year, month, day, hr_calc, min_calc, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunAnglesRadians);
solarpos_spa(year, month, day, hr_calc, min_calc, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunAnglesRadians, spa_table);

timeStepSunPosition[2] = 3;
}
Expand All @@ -2315,12 +2362,12 @@ int irrad::calc() {
{
timeStepSunPosition[0] = hour;
timeStepSunPosition[1] = (int)minute;
solarpos_spa(year, month, day, hour, minute, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunAnglesRadians);
solarpos_spa(year, month, day, hour, minute, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunAnglesRadians, spa_table);
timeStepSunPosition[2] = 1;
}
else {
// sun is down, assign sundown values
solarpos_spa(year, month, day, hour, minute, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunAnglesRadians);
solarpos_spa(year, month, day, hour, minute, 0.0, latitudeDegrees, longitudeDegrees, timezone, dut1, elevation, pressure, tamb, tiltDegrees, surfaceAzimuthDegrees, sunAnglesRadians, spa_table);
timeStepSunPosition[0] = hour;
timeStepSunPosition[1] = (int) minute;
timeStepSunPosition[2] = 0;
Expand Down
Loading

0 comments on commit 77b85d4

Please sign in to comment.