21223558770

Committed 21 Jan 2026 07:50PM UTC coverage: 87.837% (+0.02%) from 87.815%

Build # 21223558770

Build Type

Pull #94

github

Committed by

web-flow

Commit Message

Merge cc11fa7f4 into e78d2bfb0

Pull Request Pull Request #94: Round robin validation tests

Run Details

6153 of 7005 relevant lines covered (87.84%)

7296207.18 hits per line

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

77.63

/coretrace/simulation_runner/native_runner/tracing_errors.cpp


#include "tracing_errors.hpp"

#include "simulation_data_export.hpp"

namespace SolTrace::NativeRunner {

    // TODO: These two functions are basically the same. Refactor to avoid code duplication.
        // NOTES:
        // SurfaceNormalErrors()                                                Errors()
    //     - uses slope error                                                        - uses specularity error
    //     - applies to surface normal (input)                    - applies to ray direction (input)
    //     - no diffuse option                                                  - has diffuse option
        //                                                                                                        - has an dot product check (goto) for surface reflection (requires DFXYZ)
    //                                                                                                        - handles sun shape errors
        //      
    // Plan: Break up surface and sun shape error handling into separate functions
        //     - remove the special case of diffuse surfaces before Errors()
    //     - Reduce the random number calls. I.e., sample theta directly rather than thetax and thetay -> this will break tests because the number of RNG calls will change.

void SurfaceNormalErrors(MTRand &myrng,
                                                 double CosIn[3],
                                                 const OpticalProperties *OptProperties,
                                                 double CosOut[3]) noexcept(false) // throw(nanexcept)
{

        /*{Purpose:  To add error terms to the surface normal vector at the surface in question

                           Input - Seed    = Seed for RNG
                                           CosIn   = Direction cosine vector of surface normal to which errors will be applied.
                                           Element = Element data record
                                           DFXYZ   = surface normal vector at interaction point

                           Output - CosOut  = Output direction cosine vector of surface normal after error terms have been included
                                           }*/

        int i = 0;
        double Origin[3] = {0.0, 0.0, 0.0},
                   Euler[3] = {0.0, 0.0, 0.0};
        double PosIn[3] = {0.0, 0.0, 0.0},
                   PosOut[3] = {0.0, 0.0, 0.0};
        DistributionType dist;
        double delop = 0.0, thetax = 0.0,
                   thetay = 0.0, theta2 = 0.0,
                   phi = 0.0, theta = 0.0;
        double RRefToLoc[3][3] = {{0.0, 0.0, 0.0},
                                                          {0.0, 0.0, 0.0},
                                                          {0.0, 0.0, 0.0}};
        double RLocToRef[3][3] = {{0.0, 0.0, 0.0},
                                                          {0.0, 0.0, 0.0},
                                                          {0.0, 0.0, 0.0}};

        if (CosIn[2] == 0.0)
        {
                if (CosIn[0] == 0.0)
                {
                        Euler[0] = 0.0;
                        Euler[1] = PI / 2.0;
                }
                else
                {
                        Euler[0] = PI / 2.0;
                        Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));
                }
        }
        else
        {
                Euler[0] = atan2(CosIn[0], CosIn[2]);
                Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));
        }

        Euler[2] = 0.0;

        CalculateTransformMatrices(Euler, RRefToLoc, RLocToRef);

        // TODO: Add distribution type to optical properties
        // dist = OptProperties->DistributionType;
        dist = OptProperties->error_distribution_type;
        // delop = OptProperties->RMSSlopeError / 1000.0;
        delop = OptProperties->slope_error / 1000.0;

        switch (dist)
        {
        case DistributionType::GAUSSIAN:                // case 'g':
                // gaussian distribution
                thetax = myrng.randNorm(0., delop);
                thetay = myrng.randNorm(0., delop);

                theta2 = thetax * thetax + thetay * thetay;
                break;
        case DistributionType::PILLBOX:                        // case 'p':
                // pillbox distribution
                do
                {
                        thetax = 2.0 * delop * myrng() - delop;
                        thetay = 2.0 * delop * myrng() - delop;
                        theta2 = thetax * thetax + thetay * thetay;
                } while (theta2 > (delop * delop));
                break;
        default:
                // TODO: Need an error here.
                break;
        }

        /* {Transform to local coordinate system of ray to set up rotation matrices for coord and inverse
           transforms} */

        TransformToLocal(PosIn, CosIn, Origin, RRefToLoc, PosOut, CosOut);

        /* {Generate errors in terms of direction cosines in local ray coordinate system} */
        theta = sqrt(theta2);
        // phi = atan2(thetay, thetax); //This function appears to  present irregularities that bias results incorrectly for small values of thetay or thetax
        phi = myrng() * 2.0 * PI; // Therefore have chosen to randomize phi rather than calculate from randomized theta components
                                  //  obtained from the distribution. The two approaches are equivalent save for this issue with arctan2.      wendelin 01-12-11

        CosOut[0] = sin(theta) * cos(phi);
        CosOut[1] = sin(theta) * sin(phi);
        CosOut[2] = cos(theta);

        for (i = 0; i < 3; i++)
        {
                PosIn[i] = PosOut[i];
                CosIn[i] = CosOut[i];
        }

        /*{Transform perturbed ray back to element system}*/
        TransformToReference(PosIn, CosIn, Origin, RLocToRef, PosOut, CosOut);
}

void Errors(
        MTRand &myrng,
        double CosIn[3],
        int Source,
        TSun *Sun,
        // telement_ptr Element,
        const OpticalProperties *OptProperties,
        // TElement *Element,
        // TOpticalProperties *OptProperties,
        double CosOut[3],
        double DFXYZ[3])
{
        /*{Purpose:  To add error terms to the perturbed ray at the surface in question

                           Input - Seed    = Seed for RNG
                                           CosIn   = Direction cosine vector of ray to which errors will be applied.
                                                                  If Source below is 1 (i.e. sunshape) then this ray vector is before interaction with element surface
                                                                  If Source below is 2 (i.e. surface error) then this ray vector is after interaction with element surface
                                                                        (i.e. reflected ray or transmitted ray)

                                           Source  = Source indicator flag
                                                           = 1 for Sunshape error (Can be gaussian, pillbox or profile data distribution)
                                                           = 2 for surface errors (Can be gaussian or pillbox distribution)
                                           Sun     = Sun data record
                                           Element = Element data record
                                           DFXYZ   = surface normal vector at interaction point

                           Output - CosOut  = Output direction cosine vector of ray after error terms have been included
                                           }*/

        double Origin[3] = {0.0, 0.0, 0.0};
        double Euler[3] = {0.0, 0.0, 0.0};
        double PosIn[3] = {0.0, 0.0, 0.0};
        double PosOut[3] = {0.0, 0.0, 0.0};
        double delop = 0, thetax = 0, thetay = 0, theta2 = 0, phi = 0, theta = 0, stest = 0;
        uint_fast64_t i;
        double RRefToLoc[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
        double RLocToRef[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};

        if (CosIn[2] == 0.0)
        {
                if (CosIn[0] == 0.0)
                {
                        Euler[0] = 0.0;
                        Euler[1] = PI / 2.0;
                }
                else
                {
                        Euler[0] = PI / 2.0;
                        Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));
                }
        }
        else
        {
                Euler[0] = atan2(CosIn[0], CosIn[2]);
                Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));
        }

        Euler[2] = 0.0;

        CalculateTransformMatrices(Euler, RRefToLoc, RLocToRef);

        unsigned int maxcall = 0;
        // g,p,d
        if (Source == 1)  // sun error
        {
                delop = Sun->Sigma;

                switch (Sun->ShapeIndex)
                {
                case SunShape::GAUSSIAN:                        // case 'g':
                        thetax = myrng.randNorm(0., delop);
                        thetay = myrng.randNorm(0., delop);

                        theta2 = thetax * thetax + thetay * thetay;
                        break;

                case SunShape::PILLBOX:                                // case 'p':
                        do
                        {
                                thetax = 2.0 * delop * myrng() - delop;
                                thetay = 2.0 * delop * myrng() - delop;
                                theta2 = thetax * thetax + thetay * thetay;
                        } while (theta2 > (delop * delop));
                        //theta = delop * sqrt(myrng()); // Wang et al. 2010 Solar Energy 195 461-474
                        //theta2 = theta * theta;
                        break;
                case SunShape::LIMBDARKENED:
                        do {
                                thetax = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;
                                thetay = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;
                                theta2 = thetax * thetax + thetay * thetay;
                                theta = sqrt(theta2);

                                stest = 1.0 - 0.5138 * std::pow((theta / Sun->MaxAngle), 4);
                        } while ((myrng() > (stest / Sun->MaxIntensity)) || (theta2 > (Sun->MaxAngle * Sun->MaxAngle)));
                        break;
                case SunShape::BUIE_CSR:
                        // This sun model has long tails so this might take more iterations
                        // TODO: add an option to set the max angle (thereby reducing the tail)
                        do 
                        {
                                thetax = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;
                                thetay = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;
                                theta2 = thetax * thetax + thetay * thetay;
                                theta = sqrt(theta2);

                                if (std::abs(theta) <= 4.65) // within solar disc
                                        stest = cos(0.326 * theta) / cos(0.308 * theta);
                                else // within circumsolar region
                                        stest = std::exp(Sun->buie_kappa) * std::pow(std::abs(theta), Sun->buie_gamma);

                        } while ((myrng() > (stest / Sun->MaxIntensity)) || (theta2 > (Sun->MaxAngle * Sun->MaxAngle)));
                        break;
                case SunShape::USER_DEFINED:
                        do
                        {
                                thetax = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;
                                thetay = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;
                                theta2 = thetax * thetax + thetay * thetay;
                                theta = sqrt(theta2); // wendelin 1-9-12  do the test once on theta NOT individually on thetax and thetay as before

                                i = 0;
                                while (i < Sun->SunShapeAngle.size() - 1 && Sun->SunShapeAngle[i] < theta)
                                        i++;

                                if (i == 0)
                                        stest = Sun->SunShapeIntensity[0];
                                else // linear interpolation (switched from average) 12-20-11 wendelin
                                        stest = Sun->SunShapeIntensity[i - 1] + (Sun->SunShapeIntensity[i] - Sun->SunShapeIntensity[i - 1]) * (theta - Sun->SunShapeAngle[i - 1]) /
                                        (Sun->SunShapeAngle[i] - Sun->SunShapeAngle[i - 1]);

                        } while ((myrng() > (stest / Sun->MaxIntensity)) || (theta2 > (Sun->MaxAngle * Sun->MaxAngle)));
                        break;
                default:
                        throw std::invalid_argument("Unsupported sun shape in Errors function.");
                }
        }

        if (Source == 2)        // surface error
        {
                // dist = OptProperties->DistributionType; // errors
                // // delop = sqrt(4.0*sqr(OptProperties->RMSSlopeError)+sqr(OptProperties->RMSSpecError))/1000.0;
                delop = OptProperties->specularity_error;

        Label_50:
                switch (OptProperties->error_distribution_type)
                {
                case DistributionType::GAUSSIAN:                        // case 'g':
                        thetax = myrng.randNorm(0., delop);
                        thetay = myrng.randNorm(0., delop);

                        theta2 = thetax * thetax + thetay * thetay;
                        break;

                case DistributionType::PILLBOX:                                // case 'p':
                        do
                        {
                                thetax = 2.0 * delop * myrng() - delop;
                                thetay = 2.0 * delop * myrng() - delop;
                                theta2 = thetax * thetax + thetay * thetay;
                        } while (theta2 > (delop * delop));
                        break;

                case DistributionType::DIFFUSE:
                        theta2 = pow(asin(sqrt(myrng())), 2);
                        break;

                default:
                        // TODO: Add error message here.
                        break;
                }
        }

        // {Transform to local coordinate system of ray to set up rotation matrices for coordinate and inverse transforms}
        TransformToLocal(PosIn, CosIn, Origin, RRefToLoc, PosOut, CosOut);

        // {Generate errors in terms of direction cosines in local ray coordinate system}
        theta = sqrt(theta2);
        theta = theta / 1.e3; // convert from mrad to rad


        // phi = atan2(thetay, thetax); //This function appears to  present irregularities that bias results incorrectly for small values of thetay or thetax
        phi = myrng() * 2.0 * PI; // Therefore have chosen to randomize phi rather than calculate from randomized theta components
                                                          //  obtained from the distribution. The two approaches are equivalent save for this issue with arctan2.      wendelin 01-12-11

        CosOut[0] = sin(theta) * cos(phi);
        CosOut[1] = sin(theta) * sin(phi);
        CosOut[2] = cos(theta);

        for (i = 0; i < 3; i++)
        {
                PosIn[i] = PosOut[i];
                CosIn[i] = CosOut[i];
        }

        //{Transform perturbed ray back to element system}
        TransformToReference(PosIn, CosIn, Origin, RLocToRef, PosOut, CosOut);

        // TODO: Remove goto, should we always do dot product check? // We could move this out of the function and into the caller.

        /*{If reflection error application and new ray direction (after errors) physically goes through opaque surface,
    then go back and get new perturbation 06-12-07}*/                
        if ((Source == 2) &&
                (OptProperties->my_type == InteractionType::REFLECTION) &&
                (DOT(CosOut, DFXYZ) < 0) &&
                maxcall++ < 50000)
        {
                goto Label_50;
        }
}
// End of Procedure--------------------------------------------------------------

} // namespace SolTrace::NativeRunner

1
2	#include "tracing_errors.hpp"
3
4	#include "simulation_data_export.hpp"
5
6	namespace SolTrace::NativeRunner {
7
8	// TODO: These two functions are basically the same. Refactor to avoid code duplication.
9	// NOTES:
10	// SurfaceNormalErrors() Errors()
11	// - uses slope error - uses specularity error
12	// - applies to surface normal (input) - applies to ray direction (input)
13	// - no diffuse option - has diffuse option
14	// - has an dot product check (goto) for surface reflection (requires DFXYZ)
15	// - handles sun shape errors
16	//
17	// Plan: Break up surface and sun shape error handling into separate functions
18	// - remove the special case of diffuse surfaces before Errors()
19	// - Reduce the random number calls. I.e., sample theta directly rather than thetax and thetay -> this will break tests because the number of RNG calls will change.
20
21	void SurfaceNormalErrors(MTRand &myrng,	739,668✔
22	double CosIn[3],
23	const OpticalProperties *OptProperties,
24	double CosOut[3]) noexcept(false) // throw(nanexcept)
25	{
26
27	/*{Purpose: To add error terms to the surface normal vector at the surface in question
28
29	Input - Seed = Seed for RNG
30	CosIn = Direction cosine vector of surface normal to which errors will be applied.
31	Element = Element data record
32	DFXYZ = surface normal vector at interaction point
33
34	Output - CosOut = Output direction cosine vector of surface normal after error terms have been included
35	}*/
36
37	int i = 0;	739,668✔
38	double Origin[3] = {0.0, 0.0, 0.0},	739,668✔
39	Euler[3] = {0.0, 0.0, 0.0};	739,668✔
40	double PosIn[3] = {0.0, 0.0, 0.0},	739,668✔
41	PosOut[3] = {0.0, 0.0, 0.0};	739,668✔
42	DistributionType dist;
43	double delop = 0.0, thetax = 0.0,	739,668✔
44	thetay = 0.0, theta2 = 0.0,	739,668✔
45	phi = 0.0, theta = 0.0;	739,668✔
46	double RRefToLoc[3][3] = {{0.0, 0.0, 0.0},	739,668✔
47	{0.0, 0.0, 0.0},
48	{0.0, 0.0, 0.0}};
49	double RLocToRef[3][3] = {{0.0, 0.0, 0.0},	739,668✔
50	{0.0, 0.0, 0.0},
51	{0.0, 0.0, 0.0}};
52
53	if (CosIn[2] == 0.0)	739,668✔
54	{
55	if (CosIn[0] == 0.0)	×
56	{
57	Euler[0] = 0.0;	×
58	Euler[1] = PI / 2.0;	×
59	}
60	else
61	{
62	Euler[0] = PI / 2.0;	×
63	Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));	×
64	}
65	}
66	else
67	{
68	Euler[0] = atan2(CosIn[0], CosIn[2]);	739,668✔
69	Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));	739,668✔
70	}
71
72	Euler[2] = 0.0;	739,668✔
73
74	CalculateTransformMatrices(Euler, RRefToLoc, RLocToRef);	739,668✔
75
76	// TODO: Add distribution type to optical properties
77	// dist = OptProperties->DistributionType;
78	dist = OptProperties->error_distribution_type;	739,668✔
79	// delop = OptProperties->RMSSlopeError / 1000.0;
80	delop = OptProperties->slope_error / 1000.0;	739,668✔
81
82	switch (dist)	739,668✔
83	{
84	case DistributionType::GAUSSIAN: // case 'g':	739,668✔
85	// gaussian distribution
86	thetax = myrng.randNorm(0., delop);	739,668✔
87	thetay = myrng.randNorm(0., delop);	739,668✔
88
89	theta2 = thetax * thetax + thetay * thetay;	739,668✔
90	break;	739,668✔
91	case DistributionType::PILLBOX: // case 'p':	×
92	// pillbox distribution
93	do
94	{
95	thetax = 2.0 * delop * myrng() - delop;	×
96	thetay = 2.0 * delop * myrng() - delop;	×
97	theta2 = thetax * thetax + thetay * thetay;	×
98	} while (theta2 > (delop * delop));	×
99	break;	×
100	default:	×
101	// TODO: Need an error here.
102	break;	×
103	}
104
105	/* {Transform to local coordinate system of ray to set up rotation matrices for coord and inverse
106	transforms} */
107
108	TransformToLocal(PosIn, CosIn, Origin, RRefToLoc, PosOut, CosOut);	739,668✔
109
110	/* {Generate errors in terms of direction cosines in local ray coordinate system} */
111	theta = sqrt(theta2);	739,668✔
112	// phi = atan2(thetay, thetax); //This function appears to present irregularities that bias results incorrectly for small values of thetay or thetax
113	phi = myrng() * 2.0 * PI; // Therefore have chosen to randomize phi rather than calculate from randomized theta components	739,668✔
114	// obtained from the distribution. The two approaches are equivalent save for this issue with arctan2. wendelin 01-12-11
115
116	CosOut[0] = sin(theta) * cos(phi);	739,668✔
117	CosOut[1] = sin(theta) * sin(phi);	739,668✔
118	CosOut[2] = cos(theta);	739,668✔
119
120	for (i = 0; i < 3; i++)	2,958,672✔
121	{
122	PosIn[i] = PosOut[i];	2,219,004✔
123	CosIn[i] = CosOut[i];	2,219,004✔
124	}
125
126	/{Transform perturbed ray back to element system}/
127	TransformToReference(PosIn, CosIn, Origin, RLocToRef, PosOut, CosOut);	739,668✔
128	}	739,668✔
129
130	void Errors(	1,297,574✔
131	MTRand &myrng,
132	double CosIn[3],
133	int Source,
134	TSun *Sun,
135	// telement_ptr Element,
136	const OpticalProperties *OptProperties,
137	// TElement *Element,
138	// TOpticalProperties *OptProperties,
139	double CosOut[3],
140	double DFXYZ[3])
141	{
142	/*{Purpose: To add error terms to the perturbed ray at the surface in question
143
144	Input - Seed = Seed for RNG
145	CosIn = Direction cosine vector of ray to which errors will be applied.
146	If Source below is 1 (i.e. sunshape) then this ray vector is before interaction with element surface
147	If Source below is 2 (i.e. surface error) then this ray vector is after interaction with element surface
148	(i.e. reflected ray or transmitted ray)
149
150	Source = Source indicator flag
151	= 1 for Sunshape error (Can be gaussian, pillbox or profile data distribution)
152	= 2 for surface errors (Can be gaussian or pillbox distribution)
153	Sun = Sun data record
154	Element = Element data record
155	DFXYZ = surface normal vector at interaction point
156
157	Output - CosOut = Output direction cosine vector of ray after error terms have been included
158	}*/
159
160	double Origin[3] = {0.0, 0.0, 0.0};	1,297,574✔
161	double Euler[3] = {0.0, 0.0, 0.0};	1,297,574✔
162	double PosIn[3] = {0.0, 0.0, 0.0};	1,297,574✔
163	double PosOut[3] = {0.0, 0.0, 0.0};	1,297,574✔
164	double delop = 0, thetax = 0, thetay = 0, theta2 = 0, phi = 0, theta = 0, stest = 0;	1,297,574✔
165	uint_fast64_t i;
166	double RRefToLoc[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};	1,297,574✔
167	double RLocToRef[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};	1,297,574✔
168
169	if (CosIn[2] == 0.0)	1,297,574✔
170	{
171	if (CosIn[0] == 0.0)	×
172	{
173	Euler[0] = 0.0;	×
174	Euler[1] = PI / 2.0;	×
175	}
176	else
177	{
178	Euler[0] = PI / 2.0;	×
179	Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));	×
180	}
181	}
182	else
183	{
184	Euler[0] = atan2(CosIn[0], CosIn[2]);	1,297,574✔
185	Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));	1,297,574✔
186	}
187
188	Euler[2] = 0.0;	1,297,574✔
189
190	CalculateTransformMatrices(Euler, RRefToLoc, RLocToRef);	1,297,574✔
191
192	unsigned int maxcall = 0;	1,297,574✔
193	// g,p,d
194	if (Source == 1) // sun error	1,297,574✔
195	{
196	delop = Sun->Sigma;	557,906✔
197
198	switch (Sun->ShapeIndex)	557,906✔
199	{
200	case SunShape::GAUSSIAN: // case 'g':	110,000✔
201	thetax = myrng.randNorm(0., delop);	110,000✔
202	thetay = myrng.randNorm(0., delop);	110,000✔
203
204	theta2 = thetax * thetax + thetay * thetay;	110,000✔
205	break;	110,000✔
206
207	case SunShape::PILLBOX: // case 'p':	177,255✔
208	do
209	{
210	thetax = 2.0 * delop * myrng() - delop;	177,255✔
211	thetay = 2.0 * delop * myrng() - delop;	177,255✔
212	theta2 = thetax * thetax + thetay * thetay;	177,255✔
213	} while (theta2 > (delop * delop));	177,255✔
214	//theta = delop * sqrt(myrng()); // Wang et al. 2010 Solar Energy 195 461-474
215	//theta2 = theta * theta;
216	break;	139,209✔
217	case SunShape::LIMBDARKENED:	154,013✔
218	do {
219	thetax = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	154,013✔
220	thetay = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	154,013✔
221	theta2 = thetax * thetax + thetay * thetay;	154,013✔
222	theta = sqrt(theta2);	154,013✔
223
224	stest = 1.0 - 0.5138 * std::pow((theta / Sun->MaxAngle), 4);	154,013✔
225	} while ((myrng() > (stest / Sun->MaxIntensity)) \|\| (theta2 > (Sun->MaxAngle * Sun->MaxAngle)));	154,013✔
226	break;	100,000✔
227	case SunShape::BUIE_CSR:	11,850,588✔
228	// This sun model has long tails so this might take more iterations
229	// TODO: add an option to set the max angle (thereby reducing the tail)
230	do
231	{
232	thetax = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	11,850,588✔
233	thetay = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	11,850,588✔
234	theta2 = thetax * thetax + thetay * thetay;	11,850,588✔
235	theta = sqrt(theta2);	11,850,588✔
236
237	if (std::abs(theta) <= 4.65) // within solar disc	11,850,588✔
238	stest = cos(0.326 * theta) / cos(0.308 * theta);	105,365✔
239	else // within circumsolar region
240	stest = std::exp(Sun->buie_kappa) * std::pow(std::abs(theta), Sun->buie_gamma);	11,745,223✔
241
242	} while ((myrng() > (stest / Sun->MaxIntensity)) \|\| (theta2 > (Sun->MaxAngle * Sun->MaxAngle)));	11,850,588✔
243	break;	100,000✔
244	case SunShape::USER_DEFINED:	172,662✔
245	do
246	{
247	thetax = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	172,662✔
248	thetay = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	172,662✔
249	theta2 = thetax * thetax + thetay * thetay;	172,662✔
250	theta = sqrt(theta2); // wendelin 1-9-12 do the test once on theta NOT individually on thetax and thetay as before	172,662✔
251
252	i = 0;	172,662✔
253	while (i < Sun->SunShapeAngle.size() - 1 && Sun->SunShapeAngle[i] < theta)	3,324,729✔
254	i++;	3,152,067✔
255
256	if (i == 0)	172,662✔
257	stest = Sun->SunShapeIntensity[0];	×
258	else // linear interpolation (switched from average) 12-20-11 wendelin
259	stest = Sun->SunShapeIntensity[i - 1] + (Sun->SunShapeIntensity[i] - Sun->SunShapeIntensity[i - 1]) * (theta - Sun->SunShapeAngle[i - 1]) /	172,662✔
260	(Sun->SunShapeAngle[i] - Sun->SunShapeAngle[i - 1]);	172,662✔
261
262	} while ((myrng() > (stest / Sun->MaxIntensity)) \|\| (theta2 > (Sun->MaxAngle * Sun->MaxAngle)));	172,662✔
263	break;	108,697✔
264	default:	×
265	throw std::invalid_argument("Unsupported sun shape in Errors function.");	×
266	}
267	}
268
269	if (Source == 2) // surface error	1,297,574✔
270	{
271	// dist = OptProperties->DistributionType; // errors
272	// // delop = sqrt(4.0*sqr(OptProperties->RMSSlopeError)+sqr(OptProperties->RMSSpecError))/1000.0;
273	delop = OptProperties->specularity_error;	739,668✔
274
275	Label_50:	739,668✔
276	switch (OptProperties->error_distribution_type)	739,668✔
277	{
278	case DistributionType::GAUSSIAN: // case 'g':	739,668✔
279	thetax = myrng.randNorm(0., delop);	739,668✔
280	thetay = myrng.randNorm(0., delop);	739,668✔
281
282	theta2 = thetax * thetax + thetay * thetay;	739,668✔
283	break;	739,668✔
284
285	case DistributionType::PILLBOX: // case 'p':	×
286	do
287	{
288	thetax = 2.0 * delop * myrng() - delop;	×
289	thetay = 2.0 * delop * myrng() - delop;	×
290	theta2 = thetax * thetax + thetay * thetay;	×
291	} while (theta2 > (delop * delop));	×
292	break;	×
293
294	case DistributionType::DIFFUSE:	×
295	theta2 = pow(asin(sqrt(myrng())), 2);	×
296	break;	×
297
298	default:	×
299	// TODO: Add error message here.
300	break;	×
301	}
302	}
303
304	// {Transform to local coordinate system of ray to set up rotation matrices for coordinate and inverse transforms}
305	TransformToLocal(PosIn, CosIn, Origin, RRefToLoc, PosOut, CosOut);	1,297,574✔
306
307	// {Generate errors in terms of direction cosines in local ray coordinate system}
308	theta = sqrt(theta2);	1,297,574✔
309	theta = theta / 1.e3; // convert from mrad to rad	1,297,574✔
310
311
312	// phi = atan2(thetay, thetax); //This function appears to present irregularities that bias results incorrectly for small values of thetay or thetax
313	phi = myrng() * 2.0 * PI; // Therefore have chosen to randomize phi rather than calculate from randomized theta components	1,297,574✔
314	// obtained from the distribution. The two approaches are equivalent save for this issue with arctan2. wendelin 01-12-11
315
316	CosOut[0] = sin(theta) * cos(phi);	1,297,574✔
317	CosOut[1] = sin(theta) * sin(phi);	1,297,574✔
318	CosOut[2] = cos(theta);	1,297,574✔
319
320	for (i = 0; i < 3; i++)	5,190,296✔
321	{
322	PosIn[i] = PosOut[i];	3,892,722✔
323	CosIn[i] = CosOut[i];	3,892,722✔
324	}
325
326	//{Transform perturbed ray back to element system}
327	TransformToReference(PosIn, CosIn, Origin, RLocToRef, PosOut, CosOut);	1,297,574✔
328
329	// TODO: Remove goto, should we always do dot product check? // We could move this out of the function and into the caller.
330
331	/*{If reflection error application and new ray direction (after errors) physically goes through opaque surface,
332	then go back and get new perturbation 06-12-07}*/
333	if ((Source == 2) &&	2,037,242✔
334	(OptProperties->my_type == InteractionType::REFLECTION) &&	739,668✔
335	(DOT(CosOut, DFXYZ) < 0) &&	2,037,242✔
336	maxcall++ < 50000)	×
337	{
338	goto Label_50;	×
339	}
340	}	1,297,574✔
341	// End of Procedure--------------------------------------------------------------
342
343	} // namespace SolTrace::NativeRunner

NREL / SolTrace / 21223558770

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