19380488106

Committed 14 Nov 2025 11:30PM UTC coverage: 88.752% (-1.3%) from 90.052%

Build # 19380488106

Build Type

push

github

Committed by

web-flow

Commit Message

Sun position calcs (#83)

* adding solar position calculators

* adding initial sun position tests

* Remove unused solar calculator utility modules

Deleted lib_pv_incidence_modifier, lib_util, and lib_weatherfile modules from solar_position_calculators, and updated CMakeLists and includes accordingly. These files are no longer required for current solar position calculations (SPA).

* Add SolarPositionCalculator class and tests

Introduces SolarPositionCalculator with multiple calculation methods (LEGACY, DUFFIE, SOLPOS, SPA_ORIGINAL, SPA), input validation, and output retrieval. Updates CMakeLists and headers for integration, refactors sun position logic, and adds comprehensive unit tests for validation and cross-method consistency.

* fixing some compiler errors and warnings

* Guard M_PI definition with preprocessor check

Added a preprocessor guard to define M_PI only if it is not already defined, preventing potential redefinition errors. DTOR and RTOD constants are now always defined without constexpr.

* Update google-tests/unit-tests/simulation_data/sun_position_test.cpp

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

* addressing PR comments

* Update sun_position_test.cpp

updating test value

* Update expected sun vector values in test

Adjusted the expected_sun_x value and relaxed the precision for sun_x comparison in the CrossValidationTest to reflect updated calculation results.

* minor changes

* Initialize elevation in SunPositionCalculator

* Update solar_position_calculator.cpp

---------

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Taylor Brown <60201147+taylorbrown75@users.noreply.github.com>

Coverage Stats

996 of 1195 new or added lines in 4 files covered. (83.35%)

5492 of 6188 relevant lines covered (88.75%)

8572678.39 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

94.44

/coretrace/simulation_data/solar_position_calculators/basic_sun_position.cpp


#include "basic_sun_position.hpp"
#include "constants.hpp"

namespace SolTrace::Data {

int st_sun_position(double lat, double day, double hour,
                    double *x, double *y, double *z)
{
    /*
    computes the sun vector xyz given arguments
    lat : [deg] latitude
    day : [] day of the year
    hour : [hour] solar time. 12.00 corresponds to sun at maximum elevation and does not necessarily match local time

    xyz coordinate system:
        x: +west
        y: +zenith
        z: +north
    */
    double Elevation, Azimuth, Zenith;

    legacy_sun_position(lat, day, hour, &Azimuth, &Zenith);
    Elevation = 90.0 - Zenith;

    // TODO: Update coordinate system
    *x = -sin(Azimuth * D2R) * cos(Elevation * D2R);
    *y = sin(Elevation * D2R);
    *z = cos(Azimuth * D2R) * cos(Elevation * D2R);

    return 1;
}

void legacy_sun_position(double lat, double day, double hour, double* azimuth, double* zenith)
{
    /*
    Computes the sun Azimuth and Zenith angles using the legacy SolTrace method. This method is similar to that presented in
    Duffie, J.A. and Beckman, W.A., "Solar Engineering of Thermal Processes," 4th Edition, Wiley, 2013. but not identical.
        lat : [deg] latitude
        day : [] day of the year
        hour : [hour] solar time. 12.00 corresponds to sun at maximum elevation and does not necessarily match local time

    Returns:
        azimuth : [deg] azimuth angle, measured east from north
        zenith : [deg] zenith angle

    */

    double Declination, HourAngle, Elevation, Azimuth;

    Declination = R2D * asin(0.39795 * cos(0.98563 * D2R * (day - 173)));
    HourAngle = 15 * (hour - 12);
    Elevation = R2D * asin(sin(Declination * D2R) * sin(lat * D2R) + cos(Declination * D2R) * cos(HourAngle * D2R) * cos(lat * D2R));
    Azimuth = R2D * acos((sin(D2R * Declination) * cos(D2R * lat) - cos(D2R * Declination) * sin(D2R * lat) * cos(D2R * HourAngle)) / cos(D2R * Elevation) + 0.0000000001);
    if (sin(HourAngle * D2R) > 0.0)
        Azimuth = 360 - Azimuth;

    *azimuth = Azimuth;
    *zenith = 90.0 - Elevation;
}

void duffie_sun_position(double lat, double lng, double tz, double day, double hour, double* azimuth, double* zenith)
{
    /*
    Computes the sun Azimuth and Zenith angles using the method presented in
    Duffie, J.A. and Beckman, W.A., "Solar Engineering of Thermal Processes," 4th Edition, Wiley, 2013.

    lat : [deg] latitude
    day : [] day of the year
    hour : [hour] solar time. 12.00 corresponds to sun at maximum elevation and does not necessarily match local time

    Returns:
    azimuth : [deg] azimuth angle, measured east from north
    zenith : [deg] zenith angle
    */
    // TODO: we could improve this by 
    //  - calculate a fractional day of year
    //  - Use More accurate equation for declination


    double B = (day - 1) * 360. / 365.; // [degrees] Equation 1.4.2
    double E = 229.2 * (0.000075 + 0.001868 * cos(D2R * B) - 0.032077 * sin(D2R * B) - 0.014615 * cos(D2R * 2 * B) - 0.04089 * sin(D2R * 2 * B)); // [minutes] Equation 1.5.3 
    double solar_time = hour + (lng / 15.0 - tz) + E / 60.0;  // [hour] solar time, Equation 1.5.2

    double HourAngle = 15 * (solar_time - 12);
    //double Declination = 23.45 * sin(D2R * (360.0 * (284 + day) / 365.0));
    double Declination = R2D * (0.006918 - 0.399912 * cos(D2R * B) + 0.070257 * sin(D2R * B) - 0.006758 * cos(D2R * 2 * B) + 0.000907 * sin(D2R * 2 * B) - 0.002697 * cos(D2R * 3 * B) + 0.00148 * sin(D2R * 3 * B));
    double Zenith = R2D * acos(cos(D2R * lat) * cos(D2R * Declination) * cos(D2R * HourAngle) + sin(D2R * lat) * sin(D2R * Declination));
    double ratio = (cos(D2R * Zenith) * sin(D2R * lat) - sin(D2R * Declination)) / (sin(D2R * Zenith) * cos(D2R * lat));

    double Azimuth;
    // adjusting numerical error
    if (ratio > 1.0 && (ratio - 1.0) < 1.e-6)
        Azimuth = 0.0;
    else if (ratio < -1.0 && (-1.0 - ratio) < 1.e-6)
        Azimuth = 180.0;
    else 
        Azimuth = R2D * fabs(acos(ratio));

    if (HourAngle < 0)
        Azimuth *= -1.0;
    Azimuth += 180.0;    // Shift to typical North reference (east positive angle)

    *azimuth = Azimuth;
    *zenith = Zenith;
}



} // namespace SolTrace::Data

1
2	#include "basic_sun_position.hpp"
3	#include "constants.hpp"
4
5	namespace SolTrace::Data {
6
7	int st_sun_position(double lat, double day, double hour,	2✔
8	double x, double y, double *z)
9	{
10	/*
11	computes the sun vector xyz given arguments
12	lat : [deg] latitude
13	day : [] day of the year
14	hour : [hour] solar time. 12.00 corresponds to sun at maximum elevation and does not necessarily match local time
15
16	xyz coordinate system:
17	x: +west
18	y: +zenith
19	z: +north
20	*/
21	double Elevation, Azimuth, Zenith;
22
23	legacy_sun_position(lat, day, hour, &Azimuth, &Zenith);	2✔
24	Elevation = 90.0 - Zenith;	2✔
25
26	// TODO: Update coordinate system
27	x = -sin(Azimuth D2R) * cos(Elevation * D2R);	2✔
28	y = sin(Elevation D2R);	2✔
29	z = cos(Azimuth D2R) * cos(Elevation * D2R);	2✔
30
31	return 1;	2✔
32	}
33
34	void legacy_sun_position(double lat, double day, double hour, double* azimuth, double* zenith)	4✔
35	{
36	/*
37	Computes the sun Azimuth and Zenith angles using the legacy SolTrace method. This method is similar to that presented in
38	Duffie, J.A. and Beckman, W.A., "Solar Engineering of Thermal Processes," 4th Edition, Wiley, 2013. but not identical.
39	lat : [deg] latitude
40	day : [] day of the year
41	hour : [hour] solar time. 12.00 corresponds to sun at maximum elevation and does not necessarily match local time
42
43	Returns:
44	azimuth : [deg] azimuth angle, measured east from north
45	zenith : [deg] zenith angle
46
47	*/
48
49	double Declination, HourAngle, Elevation, Azimuth;
50
51	Declination = R2D * asin(0.39795 * cos(0.98563 * D2R * (day - 173)));	4✔
52	HourAngle = 15 * (hour - 12);	4✔
53	Elevation = R2D * asin(sin(Declination * D2R) * sin(lat * D2R) + cos(Declination * D2R) * cos(HourAngle * D2R) * cos(lat * D2R));	4✔
54	Azimuth = R2D * acos((sin(D2R * Declination) * cos(D2R * lat) - cos(D2R * Declination) * sin(D2R * lat) * cos(D2R * HourAngle)) / cos(D2R * Elevation) + 0.0000000001);	4✔
55	if (sin(HourAngle * D2R) > 0.0)	4✔
56	Azimuth = 360 - Azimuth;	1✔
57
58	*azimuth = Azimuth;	4✔
59	*zenith = 90.0 - Elevation;	4✔
60	}	4✔
61
62	void duffie_sun_position(double lat, double lng, double tz, double day, double hour, double* azimuth, double* zenith)	541✔
63	{
64	/*
65	Computes the sun Azimuth and Zenith angles using the method presented in
66	Duffie, J.A. and Beckman, W.A., "Solar Engineering of Thermal Processes," 4th Edition, Wiley, 2013.
67
68	lat : [deg] latitude
69	day : [] day of the year
70	hour : [hour] solar time. 12.00 corresponds to sun at maximum elevation and does not necessarily match local time
71
72	Returns:
73	azimuth : [deg] azimuth angle, measured east from north
74	zenith : [deg] zenith angle
75	*/
76	// TODO: we could improve this by
77	// - calculate a fractional day of year
78	// - Use More accurate equation for declination
79
80
81	double B = (day - 1) * 360. / 365.; // [degrees] Equation 1.4.2	541✔
82	double E = 229.2 * (0.000075 + 0.001868 * cos(D2R * B) - 0.032077 * sin(D2R * B) - 0.014615 * cos(D2R * 2 * B) - 0.04089 * sin(D2R * 2 * B)); // [minutes] Equation 1.5.3	541✔
83	double solar_time = hour + (lng / 15.0 - tz) + E / 60.0; // [hour] solar time, Equation 1.5.2	541✔
84
85	double HourAngle = 15 * (solar_time - 12);	541✔
86	//double Declination = 23.45 * sin(D2R * (360.0 * (284 + day) / 365.0));
87	double Declination = R2D * (0.006918 - 0.399912 * cos(D2R * B) + 0.070257 * sin(D2R * B) - 0.006758 * cos(D2R * 2 * B) + 0.000907 * sin(D2R * 2 * B) - 0.002697 * cos(D2R * 3 * B) + 0.00148 * sin(D2R * 3 * B));	541✔
88	double Zenith = R2D * acos(cos(D2R * lat) * cos(D2R * Declination) * cos(D2R * HourAngle) + sin(D2R * lat) * sin(D2R * Declination));	541✔
89	double ratio = (cos(D2R * Zenith) * sin(D2R * lat) - sin(D2R * Declination)) / (sin(D2R * Zenith) * cos(D2R * lat));	541✔
90
91	double Azimuth;
92	// adjusting numerical error
93	if (ratio > 1.0 && (ratio - 1.0) < 1.e-6)	541✔
NEW 94	Azimuth = 0.0;	×
95	else if (ratio < -1.0 && (-1.0 - ratio) < 1.e-6)	541✔
NEW 96	Azimuth = 180.0;	×
97	else
98	Azimuth = R2D * fabs(acos(ratio));	541✔
99
100	if (HourAngle < 0)	541✔
101	Azimuth *= -1.0;	239✔
102	Azimuth += 180.0; // Shift to typical North reference (east positive angle)	541✔
103
104	*azimuth = Azimuth;	541✔
105	*zenith = Zenith;	541✔
106	}	541✔
107
108
109
110	} // namespace SolTrace::Data

NREL / SolTrace / 19380488106

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous