1.3.77
 
Loading...
Searching...
No Matches
selfTest.cpp
1#include "SolarPosition.h"
2
3#define DOCTEST_CONFIG_IMPLEMENT
4#include <doctest.h>
5#include "doctest_utils.h"
6
7using namespace helios;
8
9TEST_CASE("SolarPosition sun position Boulder") {
10 Context context_s;
11
12 DOCTEST_CHECK_NOTHROW(context_s.setDate(make_Date(1, 1, 2000)));
13 DOCTEST_CHECK_NOTHROW(context_s.setTime(make_Time(10, 30, 0)));
14
15 SolarPosition sp(7, 40.1250f, 105.2369f, &context_s);
16 float theta_s = sp.getSunElevation() * 180.f / M_PI;
17 float phi_s = sp.getSunAzimuth() * 180.f / M_PI;
18
19 DOCTEST_CHECK(std::fabs(theta_s - 29.49f) <= 10.0f);
20 DOCTEST_CHECK(std::fabs(phi_s - 154.18f) <= 5.0f);
21}
22
23TEST_CASE("SolarPosition fractional UTC offset") {
24 // A fractional UTC offset (e.g. +05:30 zones like India, Helios convention -5.5) must be
25 // honored rather than truncated to a whole hour. The local standard time meridian scales
26 // linearly with the offset, so the half-hour offset should fall between the two neighboring
27 // integer offsets and be very close to their average.
28 Context context_s;
29 DOCTEST_CHECK_NOTHROW(context_s.setDate(make_Date(80, 2023))); // near equinox
30 DOCTEST_CHECK_NOTHROW(context_s.setTime(make_Time(12, 0, 0)));
31
32 const float latitude = 28.6139f; // New Delhi
33 const float longitude = -77.2090f; // Helios convention: +West, so New Delhi is negative
34
35 SolarPosition sp_lo(5.f, latitude, longitude, &context_s);
36 SolarPosition sp_mid(5.5f, latitude, longitude, &context_s);
37 SolarPosition sp_hi(6.f, latitude, longitude, &context_s);
38
39 float az_lo = sp_lo.getSunAzimuth();
40 float az_mid = sp_mid.getSunAzimuth();
41 float az_hi = sp_hi.getSunAzimuth();
42
43 // The half-hour offset must differ from both whole-hour offsets (i.e. it was not truncated to
44 // a whole hour) and must lie strictly between them, since the local standard time meridian — and
45 // hence the computed sun azimuth — varies monotonically with the UTC offset over this 1-hour span.
46 DOCTEST_CHECK(std::fabs(az_mid - az_lo) > 1e-4f);
47 DOCTEST_CHECK(std::fabs(az_mid - az_hi) > 1e-4f);
48 DOCTEST_CHECK(((az_lo < az_mid && az_mid < az_hi) || (az_hi < az_mid && az_mid < az_lo)));
49}
50
51TEST_CASE("SolarPosition ambient longwave model") {
52 Context context_s;
53 DOCTEST_CHECK_NOTHROW(context_s.setDate(make_Date(5, 5, 2003)));
54 DOCTEST_CHECK_NOTHROW(context_s.setTime(make_Time(9, 10, 0)));
55
56 SolarPosition sp(6, 36.5289f, 97.4439f, &context_s);
57
58 // Set atmospheric conditions
59 sp.setAtmosphericConditions(101325.f, 290.f, 0.5f, 0.02f);
60
61 float LW;
62 DOCTEST_CHECK_NOTHROW(LW = sp.getAmbientLongwaveFlux());
63
64 DOCTEST_CHECK(doctest::Approx(310.03192f).epsilon(1e-6f) == LW);
65}
66
67TEST_CASE("SolarPosition sunrise and sunset") {
68 Context context_s;
69 DOCTEST_CHECK_NOTHROW(context_s.setDate(make_Date(1, 1, 2023)));
70 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
71
72 Time sunrise;
73 DOCTEST_CHECK_NOTHROW(sunrise = sp.getSunriseTime());
74 Time sunset;
75 DOCTEST_CHECK_NOTHROW(sunset = sp.getSunsetTime());
76
77 DOCTEST_CHECK(!(sunrise.hour == 0 && sunrise.minute == 0));
78 DOCTEST_CHECK(!(sunset.hour == 0 && sunset.minute == 0));
79}
80
81TEST_CASE("SolarPosition sun direction vector") {
82 Context context_s;
83 DOCTEST_CHECK_NOTHROW(context_s.setDate(make_Date(1, 1, 2023)));
84 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
85
86 vec3 dir;
87 DOCTEST_CHECK_NOTHROW(dir = sp.getSunDirectionVector());
88 DOCTEST_CHECK(dir.x != 0.f);
89 DOCTEST_CHECK(dir.y != 0.f);
90 DOCTEST_CHECK(dir.z != 0.f);
91}
92
93TEST_CASE("SolarPosition sun direction spherical") {
94 Context context_s;
95 DOCTEST_CHECK_NOTHROW(context_s.setDate(make_Date(1, 1, 2023)));
96 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
97
98 SphericalCoord dir;
99 DOCTEST_CHECK_NOTHROW(dir = sp.getSunDirectionSpherical());
100 DOCTEST_CHECK(dir.elevation > 0.f);
101 DOCTEST_CHECK(dir.azimuth > 0.f);
102}
103
104TEST_CASE("SolarPosition flux and fractions") {
105 Context context_s;
106 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
107
108 // Set atmospheric conditions
109 sp.setAtmosphericConditions(101325.f, 300.f, 0.5f, 0.02f);
110
111 float flux;
112 DOCTEST_CHECK_NOTHROW(flux = sp.getSolarFlux());
113 DOCTEST_CHECK(flux > 0.f);
114
115 float diffuse_fraction;
116 DOCTEST_CHECK_NOTHROW(diffuse_fraction = sp.getDiffuseFraction());
117 DOCTEST_CHECK(diffuse_fraction >= 0.f);
118 DOCTEST_CHECK(diffuse_fraction <= 1.f);
119
120 float flux_par;
121 DOCTEST_CHECK_NOTHROW(flux_par = sp.getSolarFluxPAR());
122 DOCTEST_CHECK(flux_par > 0.f);
123
124 float flux_nir;
125 DOCTEST_CHECK_NOTHROW(flux_nir = sp.getSolarFluxNIR());
126 DOCTEST_CHECK(flux_nir > 0.f);
127}
128
129TEST_CASE("SolarPosition elevation, zenith, azimuth") {
130 Context context_s;
131 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
132
133 float elevation;
134 DOCTEST_CHECK_NOTHROW(elevation = sp.getSunElevation());
135 DOCTEST_CHECK(elevation >= 0.f);
136 DOCTEST_CHECK(elevation <= M_PI / 2.f);
137
138 float zenith;
139 DOCTEST_CHECK_NOTHROW(zenith = sp.getSunZenith());
140 DOCTEST_CHECK(zenith >= 0.f);
141 DOCTEST_CHECK(zenith <= M_PI);
142
143 float azimuth;
144 DOCTEST_CHECK_NOTHROW(azimuth = sp.getSunAzimuth());
145 DOCTEST_CHECK(azimuth >= 0.f);
146 DOCTEST_CHECK(azimuth <= 2.f * M_PI);
147}
148
149TEST_CASE("SolarPosition turbidity calibration") {
150 Context context_s;
151 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
152 std::string label = "test_flux_timeseries";
153
154 if (!context_s.doesTimeseriesVariableExist(label.c_str())) {
155 return; // skip test if data does not exist
156 }
157
158 float turbidity;
159 std::string captured_warnings;
160 {
161 capture_cerr cerr_buffer;
162 DOCTEST_CHECK_NOTHROW(turbidity = sp.calibrateTurbidityFromTimeseries(label));
163 captured_warnings = cerr_buffer.get_captured_output();
164 } // Capture goes out of scope before assertions
165
166 DOCTEST_CHECK(turbidity > 0.f);
167
168 // Turbidity calibration may produce fzero warnings due to the nature of the optimization problem
169 // This is expected behavior and should not be considered a test failure
170 // Just verify the function completes and returns a valid result
171}
172
173TEST_CASE("SolarPosition invalid lat/long") {
174 Context context_s;
175
176 bool had_output_1, had_output_2;
177 {
178 capture_cerr cerr_buffer;
179 SolarPosition sp_1(7, -100.f, 105.2369f, &context_s);
180 had_output_1 = cerr_buffer.has_output();
181
182 cerr_buffer.clear();
183 SolarPosition sp_2(7, 40.125f, -200.f, &context_s);
184 had_output_2 = cerr_buffer.has_output();
185 } // capture goes out of scope before assertions
186 DOCTEST_CHECK(had_output_1);
187 DOCTEST_CHECK(had_output_2);
188}
189
190TEST_CASE("SolarPosition invalid solar angle") {
191 Context context_s;
192 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
193
194 // Set atmospheric conditions
195 sp.setAtmosphericConditions(101325.f, 300.f, 0.5f, 0.02f);
196
197 DOCTEST_CHECK_NOTHROW(sp.setSunDirection(make_SphericalCoord(0.75 * M_PI, M_PI / 2.f)));
198
199 float flux;
200 DOCTEST_CHECK_NOTHROW(flux = sp.getSolarFlux());
201 DOCTEST_CHECK(flux == 0.f);
202}
203
204
205TEST_CASE("SolarPosition solor position overridden") {
206 Context context_s;
207 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
208
209 DOCTEST_CHECK_NOTHROW(sp.setSunDirection(make_SphericalCoord(M_PI / 4.f, M_PI / 2.f)));
210
211 float elevation;
212 DOCTEST_CHECK_NOTHROW(elevation = sp.getSunElevation());
213 DOCTEST_CHECK(elevation >= 0.f);
214 DOCTEST_CHECK(elevation <= M_PI / 2.f);
215
216 float zenith;
217 DOCTEST_CHECK_NOTHROW(zenith = sp.getSunZenith());
218 DOCTEST_CHECK(zenith >= 0.f);
219 DOCTEST_CHECK(zenith <= M_PI);
220
221 float azimuth;
222 DOCTEST_CHECK_NOTHROW(azimuth = sp.getSunAzimuth());
223 DOCTEST_CHECK(azimuth >= 0.f);
224 DOCTEST_CHECK(azimuth <= 2.f * M_PI);
225
226 vec3 sun_vector;
227 DOCTEST_CHECK_NOTHROW(sun_vector = sp.getSunDirectionVector());
228
229 SphericalCoord sun_spherical;
230 DOCTEST_CHECK_NOTHROW(sun_spherical = sp.getSunDirectionSpherical());
231}
232
233TEST_CASE("SolarPosition cloud calibration") {
234 Context context_s;
235 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
236
237 DOCTEST_CHECK_NOTHROW(context_s.loadTabularTimeseriesData("lib/testdata/cimis.csv", {"CIMIS"}, ","));
238
239 context_s.setDate(make_Date(13, 7, 2023));
240 context_s.setTime(make_Time(12, 0, 0));
241
242 // Set atmospheric conditions
243 sp.setAtmosphericConditions(101325.f, 300.f, 0.5f, 0.02f);
244
245 DOCTEST_CHECK_NOTHROW(sp.enableCloudCalibration("net_radiation"));
246
247 float flux;
248 DOCTEST_CHECK_NOTHROW(flux = sp.getSolarFlux());
249 DOCTEST_CHECK(flux > 0.f);
250
251 float diffuse_fraction;
252 DOCTEST_CHECK_NOTHROW(diffuse_fraction = sp.getDiffuseFraction());
253 DOCTEST_CHECK(diffuse_fraction >= 0.f);
254 DOCTEST_CHECK(diffuse_fraction <= 1.f);
255
256 float flux_par;
257 DOCTEST_CHECK_NOTHROW(flux_par = sp.getSolarFluxPAR());
258 DOCTEST_CHECK(flux_par > 0.f);
259
260 float flux_nir;
261 DOCTEST_CHECK_NOTHROW(flux_nir = sp.getSolarFluxNIR());
262 DOCTEST_CHECK(flux_nir > 0.f);
263
264 DOCTEST_CHECK_NOTHROW(sp.disableCloudCalibration());
265
266 {
267 capture_cerr cerr_buffer;
268 DOCTEST_CHECK_THROWS_AS(sp.enableCloudCalibration("non_existent_timeseries"), std::runtime_error);
269 }
270}
271
272TEST_CASE("SolarPosition turbidity calculation") {
273 Context context_s;
274 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
275
276 DOCTEST_CHECK_NOTHROW(context_s.loadTabularTimeseriesData("lib/testdata/cimis.csv", {"CIMIS"}, ","));
277
278 float turbidity;
279 {
280 // Turbidity calibration may produce fzero warnings - capture them
281 capture_cerr cerr_buffer;
282 DOCTEST_CHECK_NOTHROW(turbidity = sp.calibrateTurbidityFromTimeseries("net_radiation"));
283 } // capture goes out of scope before assertion
284 DOCTEST_CHECK(turbidity > 0.f);
285
286 {
287 capture_cerr cerr_buffer;
288 DOCTEST_CHECK_THROWS_AS(turbidity = sp.calibrateTurbidityFromTimeseries("non_existent_timeseries"), std::runtime_error);
289 }
290}
291
292TEST_CASE("SolarPosition setAtmosphericConditions valid inputs") {
293 Context context_s;
294 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
295
296 // Test setting valid atmospheric conditions
297 DOCTEST_CHECK_NOTHROW(sp.setAtmosphericConditions(101325.f, 300.f, 0.5f, 0.02f));
298
299 // Verify global data was set correctly
300 float pressure, temperature, humidity, turbidity;
301 DOCTEST_CHECK_NOTHROW(context_s.getGlobalData("atmosphere_pressure_Pa", pressure));
302 DOCTEST_CHECK_NOTHROW(context_s.getGlobalData("atmosphere_temperature_K", temperature));
303 DOCTEST_CHECK_NOTHROW(context_s.getGlobalData("atmosphere_humidity_rel", humidity));
304 DOCTEST_CHECK_NOTHROW(context_s.getGlobalData("atmosphere_turbidity", turbidity));
305
306 DOCTEST_CHECK(doctest::Approx(101325.f).epsilon(1e-6f) == pressure);
307 DOCTEST_CHECK(doctest::Approx(300.f).epsilon(1e-6f) == temperature);
308 DOCTEST_CHECK(doctest::Approx(0.5f).epsilon(1e-6f) == humidity);
309 DOCTEST_CHECK(doctest::Approx(0.02f).epsilon(1e-6f) == turbidity);
310}
311
312TEST_CASE("SolarPosition setAtmosphericConditions validation") {
313 Context context_s;
314 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
315
316 // Test invalid pressure (negative)
317 DOCTEST_CHECK_THROWS_AS(sp.setAtmosphericConditions(-1000.f, 300.f, 0.5f, 0.02f), std::runtime_error);
318
319 // Test invalid pressure (zero)
320 DOCTEST_CHECK_THROWS_AS(sp.setAtmosphericConditions(0.f, 300.f, 0.5f, 0.02f), std::runtime_error);
321
322 // Test invalid temperature (negative)
323 DOCTEST_CHECK_THROWS_AS(sp.setAtmosphericConditions(101325.f, -10.f, 0.5f, 0.02f), std::runtime_error);
324
325 // Test invalid temperature (zero)
326 DOCTEST_CHECK_THROWS_AS(sp.setAtmosphericConditions(101325.f, 0.f, 0.5f, 0.02f), std::runtime_error);
327
328 // Test invalid humidity (negative)
329 DOCTEST_CHECK_THROWS_AS(sp.setAtmosphericConditions(101325.f, 300.f, -0.1f, 0.02f), std::runtime_error);
330
331 // Test invalid humidity (> 1)
332 DOCTEST_CHECK_THROWS_AS(sp.setAtmosphericConditions(101325.f, 300.f, 1.5f, 0.02f), std::runtime_error);
333
334 // Test invalid turbidity (negative)
335 DOCTEST_CHECK_THROWS_AS(sp.setAtmosphericConditions(101325.f, 300.f, 0.5f, -0.01f), std::runtime_error);
336
337 // Test boundary values (should succeed)
338 DOCTEST_CHECK_NOTHROW(sp.setAtmosphericConditions(0.001f, 0.001f, 0.f, 0.f));
339 DOCTEST_CHECK_NOTHROW(sp.setAtmosphericConditions(200000.f, 400.f, 1.f, 1.f));
340}
341
342TEST_CASE("SolarPosition getAtmosphericConditions retrieval") {
343 Context context_s;
344 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
345
346 // Set atmospheric conditions
347 sp.setAtmosphericConditions(98000.f, 295.f, 0.65f, 0.05f);
348
349 // Retrieve atmospheric conditions
350 float pressure, temperature, humidity, turbidity;
351 DOCTEST_CHECK_NOTHROW(sp.getAtmosphericConditions(pressure, temperature, humidity, turbidity));
352
353 // Verify retrieved values match what was set
354 DOCTEST_CHECK(doctest::Approx(98000.f).epsilon(1e-6f) == pressure);
355 DOCTEST_CHECK(doctest::Approx(295.f).epsilon(1e-6f) == temperature);
356 DOCTEST_CHECK(doctest::Approx(0.65f).epsilon(1e-6f) == humidity);
357 DOCTEST_CHECK(doctest::Approx(0.05f).epsilon(1e-6f) == turbidity);
358}
359
360TEST_CASE("SolarPosition getAtmosphericConditions defaults") {
361 Context context_s;
362 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
363
364 // Get atmospheric conditions without setting them first
365 float pressure, temperature, humidity, turbidity;
366 bool had_warning;
367 {
368 capture_cerr cerr_buffer;
369 sp.getAtmosphericConditions(pressure, temperature, humidity, turbidity);
370 had_warning = cerr_buffer.has_output();
371 } // capture goes out of scope before assertions
372
373 // Verify warning was issued
374 DOCTEST_CHECK(had_warning);
375
376 // Verify default values are used
377 DOCTEST_CHECK(doctest::Approx(101325.f).epsilon(1e-6f) == pressure);
378 DOCTEST_CHECK(doctest::Approx(300.f).epsilon(1e-6f) == temperature);
379 DOCTEST_CHECK(doctest::Approx(0.5f).epsilon(1e-6f) == humidity);
380 DOCTEST_CHECK(doctest::Approx(0.02f).epsilon(1e-6f) == turbidity);
381}
382
383TEST_CASE("SolarPosition parameter-free flux methods") {
384 Context context_s;
385 context_s.setDate(make_Date(1, 6, 2023));
386 context_s.setTime(make_Time(12, 0, 0));
387
388 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
389
390 // Set atmospheric conditions
391 sp.setAtmosphericConditions(101325.f, 300.f, 0.5f, 0.02f);
392
393 // Test parameter-free getSolarFlux()
394 float flux;
395 DOCTEST_CHECK_NOTHROW(flux = sp.getSolarFlux());
396 DOCTEST_CHECK(flux > 0.f);
397
398 // Test parameter-free getSolarFluxPAR()
399 float flux_par;
400 DOCTEST_CHECK_NOTHROW(flux_par = sp.getSolarFluxPAR());
401 DOCTEST_CHECK(flux_par > 0.f);
402
403 // Test parameter-free getSolarFluxNIR()
404 float flux_nir;
405 DOCTEST_CHECK_NOTHROW(flux_nir = sp.getSolarFluxNIR());
406 DOCTEST_CHECK(flux_nir > 0.f);
407
408 // Test parameter-free getDiffuseFraction()
409 float diffuse_fraction;
410 DOCTEST_CHECK_NOTHROW(diffuse_fraction = sp.getDiffuseFraction());
411 DOCTEST_CHECK(diffuse_fraction >= 0.f);
412 DOCTEST_CHECK(diffuse_fraction <= 1.f);
413
414 // Test parameter-free getAmbientLongwaveFlux()
415 float lw_flux;
416 DOCTEST_CHECK_NOTHROW(lw_flux = sp.getAmbientLongwaveFlux());
417 DOCTEST_CHECK(lw_flux > 0.f);
418}
419
420TEST_CASE("SolarPosition parameter-free methods with defaults") {
421 Context context_s;
422 context_s.setDate(make_Date(1, 6, 2023));
423 context_s.setTime(make_Time(12, 0, 0));
424
425 SolarPosition sp(7, 40.125f, 105.2369f, &context_s);
426
427 // Call parameter-free methods without setting atmospheric conditions
428 // Should use default values
429 float flux;
430 DOCTEST_CHECK_NOTHROW(flux = sp.getSolarFlux());
431 DOCTEST_CHECK(flux > 0.f);
432
433 // Verify defaults are being used
434 float pressure, temperature, humidity, turbidity;
435 {
436 capture_cerr cerr_buffer;
437 sp.getAtmosphericConditions(pressure, temperature, humidity, turbidity);
438 } // capture goes out of scope before assertions
439 DOCTEST_CHECK(doctest::Approx(101325.f).epsilon(1e-6f) == pressure);
440 DOCTEST_CHECK(doctest::Approx(300.f).epsilon(1e-6f) == temperature);
441 DOCTEST_CHECK(doctest::Approx(0.5f).epsilon(1e-6f) == humidity);
442 DOCTEST_CHECK(doctest::Approx(0.02f).epsilon(1e-6f) == turbidity);
443}
444
445TEST_CASE("SSolar-GOA spectral irradiance") {
446 Context context_s;
447
448 DOCTEST_CHECK_NOTHROW(context_s.setDate(make_Date(16, 7, 2023)));
449 DOCTEST_CHECK_NOTHROW(context_s.setTime(make_Time(12, 0, 0)));
450
451 SolarPosition sp(0, 36.93f, 3.33f, &context_s);
452 sp.setAtmosphericConditions(87700.f, 298.f, 0.5f, 0.026f);
453
454 DOCTEST_CHECK_NOTHROW(sp.calculateGlobalSolarSpectrum("global_test"));
455 DOCTEST_CHECK_NOTHROW(sp.calculateDirectSolarSpectrum("direct_test"));
456 DOCTEST_CHECK_NOTHROW(sp.calculateDiffuseSolarSpectrum("diffuse_test"));
457
458 std::vector<vec2> global_spectrum, direct_spectrum, diffuse_spectrum;
459 DOCTEST_CHECK_NOTHROW(context_s.getGlobalData("global_test", global_spectrum));
460 DOCTEST_CHECK_NOTHROW(context_s.getGlobalData("direct_test", direct_spectrum));
461 DOCTEST_CHECK_NOTHROW(context_s.getGlobalData("diffuse_test", diffuse_spectrum));
462
463 DOCTEST_CHECK(global_spectrum.size() == 2301);
464 DOCTEST_CHECK(direct_spectrum.size() == 2301);
465 DOCTEST_CHECK(diffuse_spectrum.size() == 2301);
466
467 DOCTEST_CHECK(doctest::Approx(300.f).epsilon(0.1f) == global_spectrum.front().x);
468 DOCTEST_CHECK(doctest::Approx(2600.f).epsilon(0.1f) == global_spectrum.back().x);
469
470 for (const auto &point: global_spectrum) {
471 DOCTEST_CHECK(point.y >= 0.f);
472 DOCTEST_CHECK(std::isfinite(point.y));
473 }
474
475 // Verify integrated PAR flux is reasonable for clear sky at solar noon
476 float par_flux = 0.f;
477 for (size_t i = 1; i < global_spectrum.size(); ++i) {
478 float wl = global_spectrum[i].x;
479 if (wl >= 400.f && wl <= 700.f) {
480 float dw = global_spectrum[i].x - global_spectrum[i - 1].x;
481 float avg_irr = 0.5f * (global_spectrum[i].y + global_spectrum[i - 1].y);
482 par_flux += avg_irr * dw;
483 }
484 }
485 DOCTEST_CHECK(par_flux > 400.f);
486 DOCTEST_CHECK(par_flux < 500.f);
487}
488
489TEST_CASE("SSolar-GOA spectral resolution") {
490 Context context_s;
491
492 DOCTEST_CHECK_NOTHROW(context_s.setDate(make_Date(16, 7, 2023)));
493 DOCTEST_CHECK_NOTHROW(context_s.setTime(make_Time(12, 0, 0)));
494
495 SolarPosition sp(0, 36.93f, 3.33f, &context_s);
496 sp.setAtmosphericConditions(87700.f, 298.f, 0.5f, 0.026f);
497
498 // Test default 1 nm resolution
499 DOCTEST_CHECK_NOTHROW(sp.calculateGlobalSolarSpectrum("res_1nm", 1.0f));
500 std::vector<vec2> spectrum_1nm;
501 DOCTEST_CHECK_NOTHROW(context_s.getGlobalData("res_1nm", spectrum_1nm));
502 DOCTEST_CHECK(spectrum_1nm.size() == 2301);
503
504 // Test 10 nm resolution
505 DOCTEST_CHECK_NOTHROW(sp.calculateGlobalSolarSpectrum("res_10nm", 10.0f));
506 std::vector<vec2> spectrum_10nm;
507 DOCTEST_CHECK_NOTHROW(context_s.getGlobalData("res_10nm", spectrum_10nm));
508 DOCTEST_CHECK(spectrum_10nm.size() == 231); // (2600-300)/10 + 1 = 231
509
510 // Test 50 nm resolution
511 DOCTEST_CHECK_NOTHROW(sp.calculateGlobalSolarSpectrum("res_50nm", 50.0f));
512 std::vector<vec2> spectrum_50nm;
513 DOCTEST_CHECK_NOTHROW(context_s.getGlobalData("res_50nm", spectrum_50nm));
514 DOCTEST_CHECK(spectrum_50nm.size() == 47); // (2600-300)/50 + 1 = 47
515
516 // Verify wavelength spacing
517 DOCTEST_CHECK(doctest::Approx(300.f).epsilon(0.5f) == spectrum_10nm.front().x);
518 DOCTEST_CHECK(doctest::Approx(310.f).epsilon(0.5f) == spectrum_10nm[1].x);
519
520 // Verify irradiance values are still reasonable
521 for (const auto &point: spectrum_10nm) {
522 DOCTEST_CHECK(point.y >= 0.f);
523 DOCTEST_CHECK(std::isfinite(point.y));
524 }
525}
526
527TEST_CASE("SSolar-GOA validation against Python reference") {
528 Context context_s;
529
530 DOCTEST_CHECK_NOTHROW(context_s.setDate(make_Date(16, 7, 2023)));
531 DOCTEST_CHECK_NOTHROW(context_s.setTime(make_Time(12, 0, 0)));
532
533 SolarPosition sp(0, 36.93f, 3.33f, &context_s);
534 sp.setAtmosphericConditions(87700.f, 298.f, 0.5f, 0.026f);
535
536 DOCTEST_CHECK_NOTHROW(sp.calculateGlobalSolarSpectrum("validation"));
537
538 std::vector<vec2> global_spectrum;
539 DOCTEST_CHECK_NOTHROW(context_s.getGlobalData("validation", global_spectrum));
540
541 // Compare against Python reference if available
542 std::ifstream ref_file("plugins/solarposition/tests/validate_reference_global.txt");
543 if (ref_file.is_open()) {
544 std::string header_line;
545 std::getline(ref_file, header_line); // Skip header
546
547 std::vector<float> ref_wavelengths, ref_irradiances;
548 std::string line;
549 while (std::getline(ref_file, line)) {
550 if (line.empty() || line[0] == '#')
551 continue;
552
553 std::istringstream iss(line);
554 float wl, irr;
555 if (iss >> wl >> irr) {
556 ref_wavelengths.push_back(wl);
557 ref_irradiances.push_back(irr);
558 }
559 }
560 ref_file.close();
561
562 DOCTEST_CHECK(ref_wavelengths.size() == 2301);
563
564 float max_rel_error = 0.f;
565 float sum_sq_error = 0.f;
566 size_t n_compared = 0;
567
568 for (size_t i = 0; i < std::min(global_spectrum.size(), ref_wavelengths.size()); ++i) {
569 DOCTEST_CHECK(doctest::Approx(ref_wavelengths[i]).epsilon(1e-6f) == global_spectrum[i].x);
570
571 float cpp_irr = global_spectrum[i].y;
572 float ref_irr = ref_irradiances[i];
573 float abs_error = std::fabs(cpp_irr - ref_irr);
574 float rel_error = abs_error / (ref_irr + 1e-10f);
575
576 max_rel_error = std::max(max_rel_error, rel_error);
577 sum_sq_error += abs_error * abs_error;
578 n_compared++;
579 }
580
581 float rms_error = std::sqrt(sum_sq_error / n_compared);
582
583 DOCTEST_CHECK(max_rel_error < 0.01f);
584 DOCTEST_CHECK(rms_error < 0.01f);
585
586 } else {
587 DOCTEST_WARN("Python reference file not found - run validate_detailed.py to enable detailed validation");
588 }
589}
590
591int SolarPosition::selfTest(int argc, char **argv) {
592 return helios::runDoctestWithValidation(argc, argv);
593}
594
595// ===== Prague Sky Model Tests =====
596
597TEST_CASE("SolarPosition - Prague model initialization") {
598 Context context;
599 SolarPosition solar(&context);
600
601 DOCTEST_CHECK(!solar.isPragueSkyModelEnabled());
602 DOCTEST_CHECK_NOTHROW(solar.enablePragueSkyModel());
603 DOCTEST_CHECK(solar.isPragueSkyModelEnabled());
604}
605
606TEST_CASE("SolarPosition - Prague angular parameter fitting - clear sky") {
607 Context context;
608 SolarPosition solar(&context);
609 solar.enablePragueSkyModel();
610 solar.setSunDirection(make_SphericalCoord(0.5236f, 0.0f)); // 60° elevation
611 solar.setAtmosphericConditions(101325.f, 300.f, 0.5f, 0.05f); // Clear sky turbidity
612 solar.updatePragueSkyModel();
613
614 std::vector<float> params;
615 DOCTEST_CHECK_NOTHROW(context.getGlobalData("prague_sky_spectral_params", params));
616
617 DOCTEST_CHECK(params.size() == 225 * 6);
618
619 // Check 550nm parameters (index 38: (550-360)/5 = 38)
620 int idx = 38 * 6;
621 float wavelength = params[idx + 0];
622 float L_zenith = params[idx + 1];
623 float circ_str = params[idx + 2];
624 float circ_width = params[idx + 3];
625 float horiz_bright = params[idx + 4];
626 float norm = params[idx + 5];
627
628 DOCTEST_CHECK(wavelength == doctest::Approx(550.0f).epsilon(1.0f));
629 DOCTEST_CHECK(L_zenith > 0.0f);
630 DOCTEST_CHECK(L_zenith < 0.5f);
631 DOCTEST_CHECK(circ_str >= 0.0f);
632 DOCTEST_CHECK(circ_str <= 20.0f); // Max clamp value
633 DOCTEST_CHECK(circ_width >= 5.0f);
634 DOCTEST_CHECK(circ_width <= 60.0f);
635 DOCTEST_CHECK(horiz_bright >= 1.0f);
636 DOCTEST_CHECK(horiz_bright < 5.0f);
637 DOCTEST_CHECK(norm > 0.0f);
638 DOCTEST_CHECK(norm < 2.0f);
639
640 // Verify global data validity flag
641 int valid = 0;
642 DOCTEST_CHECK_NOTHROW(context.getGlobalData("prague_sky_valid", valid));
643 DOCTEST_CHECK(valid == 1);
644}
645
646TEST_CASE("SolarPosition - Prague lazy evaluation") {
647 Context context;
648 SolarPosition solar(&context);
649 solar.enablePragueSkyModel();
650 solar.setSunDirection(make_SphericalCoord(0.5236f, 0.0f)); // 60° elevation
651 solar.setAtmosphericConditions(101325.f, 300.f, 0.5f, 0.1f);
652 solar.updatePragueSkyModel();
653
654 // Small turbidity change - should not need update
655 solar.setAtmosphericConditions(101325.f, 300.f, 0.5f, 0.101f);
656 DOCTEST_CHECK(!solar.pragueSkyModelNeedsUpdate(0.33f));
657
658 // Large turbidity change - should need update
659 solar.setAtmosphericConditions(101325.f, 300.f, 0.5f, 0.15f);
660 DOCTEST_CHECK(solar.pragueSkyModelNeedsUpdate(0.33f));
661
662 // Sun direction change
663 solar.setSunDirection(make_SphericalCoord(0.6109f, 0.0f)); // 55° elevation (change > 5°)
664 DOCTEST_CHECK(solar.pragueSkyModelNeedsUpdate(0.33f));
665
666 // Albedo change
667 solar.setSunDirection(make_SphericalCoord(0.5236f, 0.0f)); // Back to 60°
668 solar.setAtmosphericConditions(101325.f, 300.f, 0.5f, 0.1f); // Reset turbidity
669 solar.updatePragueSkyModel(); // Update with new sun direction
670 DOCTEST_CHECK(solar.pragueSkyModelNeedsUpdate(0.5f)); // Albedo 0.33 -> 0.5
671}
672
673TEST_CASE("SolarPosition - Prague performance benchmark") {
674 Context context;
675 SolarPosition solar(&context);
676 solar.enablePragueSkyModel();
677 solar.setSunDirection(make_SphericalCoord(0.5236f, 0.0f)); // 60° elevation
678 solar.setAtmosphericConditions(101325.f, 300.f, 0.5f, 0.1f);
679
680 auto start = std::chrono::high_resolution_clock::now();
681 solar.updatePragueSkyModel();
682 auto end = std::chrono::high_resolution_clock::now();
683
684 auto duration_ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
685
686 DOCTEST_CHECK(duration_ms.count() < 10000); // <10 seconds
687}
688
689TEST_CASE("SolarPosition - Prague error handling") {
690 Context context;
691 SolarPosition solar(&context);
692
693 // Try to update without enabling
694 DOCTEST_CHECK_THROWS_WITH_AS(solar.updatePragueSkyModel(), "ERROR (SolarPosition::updatePragueSkyModel): Prague model not enabled. Call enablePragueSkyModel() first.", std::runtime_error);
695
696 // Check needs update returns false when not enabled
697 DOCTEST_CHECK(!solar.pragueSkyModelNeedsUpdate(0.33f));
698}