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