20348441017

Committed 18 Dec 2025 07:15PM UTC coverage: 89.332% (+0.1%) from 89.184%

Build # 20348441017

Build Type

Pull #94

github

Committed by

web-flow

Commit Message

Merge c719a34ba into a73a760a4

Pull Request Pull Request #94: Round robin validation tests

Run Details

6029 of 6749 relevant lines covered (89.33%)

48291780.77 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,	24,308,261✔
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;	24,308,261✔
38	double Origin[3] = {0.0, 0.0, 0.0},	24,308,261✔
39	Euler[3] = {0.0, 0.0, 0.0};	24,308,261✔
40	double PosIn[3] = {0.0, 0.0, 0.0},	24,308,261✔
41	PosOut[3] = {0.0, 0.0, 0.0};	24,308,261✔
42	DistributionType dist;
43	double delop = 0.0, thetax = 0.0,	24,308,261✔
44	thetay = 0.0, theta2 = 0.0,	24,308,261✔
45	phi = 0.0, theta = 0.0;	24,308,261✔
46	double RRefToLoc[3][3] = {{0.0, 0.0, 0.0},	24,308,261✔
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},	24,308,261✔
50	{0.0, 0.0, 0.0},
51	{0.0, 0.0, 0.0}};
52
53	if (CosIn[2] == 0.0)	24,308,261✔
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]);	24,308,261✔
69	Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));	24,308,261✔
70	}
71
72	Euler[2] = 0.0;	24,308,261✔
73
74	CalculateTransformMatrices(Euler, RRefToLoc, RLocToRef);	24,308,261✔
75
76	// TODO: Add distribution type to optical properties
77	// dist = OptProperties->DistributionType;
78	dist = OptProperties->error_distribution_type;	24,308,261✔
79	// delop = OptProperties->RMSSlopeError / 1000.0;
80	delop = OptProperties->slope_error / 1000.0;	24,308,261✔
81
82	switch (dist)	24,308,261✔
83	{
84	case DistributionType::GAUSSIAN: // case 'g':	24,308,261✔
85	// gaussian distribution
86	thetax = myrng.randNorm(0., delop);	24,308,261✔
87	thetay = myrng.randNorm(0., delop);	24,308,261✔
88
89	theta2 = thetax * thetax + thetay * thetay;	24,308,261✔
90	break;	24,308,261✔
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);	24,308,261✔
109
110	/* {Generate errors in terms of direction cosines in local ray coordinate system} */
111	theta = sqrt(theta2);	24,308,261✔
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	24,308,261✔
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);	24,308,261✔
117	CosOut[1] = sin(theta) * sin(phi);	24,308,261✔
118	CosOut[2] = cos(theta);	24,308,261✔
119
120	for (i = 0; i < 3; i++)	97,233,044✔
121	{
122	PosIn[i] = PosOut[i];	72,924,783✔
123	CosIn[i] = CosOut[i];	72,924,783✔
124	}
125
126	/{Transform perturbed ray back to element system}/
127	TransformToReference(PosIn, CosIn, Origin, RLocToRef, PosOut, CosOut);	24,308,261✔
128	}	24,308,261✔
129
130	void Errors(	47,358,941✔
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};	47,358,941✔
161	double Euler[3] = {0.0, 0.0, 0.0};	47,358,941✔
162	double PosIn[3] = {0.0, 0.0, 0.0};	47,358,941✔
163	double PosOut[3] = {0.0, 0.0, 0.0};	47,358,941✔
164	double delop = 0, thetax = 0, thetay = 0, theta2 = 0, phi = 0, theta = 0, stest = 0;	47,358,941✔
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}};	47,358,941✔
167	double RLocToRef[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};	47,358,941✔
168
169	if (CosIn[2] == 0.0)	47,358,941✔
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]);	47,358,941✔
185	Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));	47,358,941✔
186	}
187
188	Euler[2] = 0.0;	47,358,941✔
189
190	CalculateTransformMatrices(Euler, RRefToLoc, RLocToRef);	47,358,941✔
191
192	unsigned int maxcall = 0;	47,358,941✔
193	// g,p,d
194	if (Source == 1) // sun error	47,358,941✔
195	{
196	delop = Sun->Sigma;	23,050,680✔
197
198	switch (Sun->ShapeIndex)	23,050,680✔
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':	28,836,127✔
208	do
209	{
210	thetax = 2.0 * delop * myrng() - delop;	28,836,127✔
211	thetay = 2.0 * delop * myrng() - delop;	28,836,127✔
212	theta2 = thetax * thetax + thetay * thetay;	28,836,127✔
213	} while (theta2 > (delop * delop));	28,836,127✔
214	//theta = delop * sqrt(myrng()); // Wang et al. 2010 Solar Energy 195 461-474
215	//theta2 = theta * theta;
216	break;	22,640,680✔
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:	159,111✔
245	do
246	{
247	thetax = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	159,111✔
248	thetay = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	159,111✔
249	theta2 = thetax * thetax + thetay * thetay;	159,111✔
250	theta = sqrt(theta2); // wendelin 1-9-12 do the test once on theta NOT individually on thetax and thetay as before	159,111✔
251
252	i = 0;	159,111✔
253	while (i < Sun->SunShapeAngle.size() - 1 && Sun->SunShapeAngle[i] < theta)	2,717,967✔
254	i++;	2,558,856✔
255
256	if (i == 0)	159,111✔
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]) /	159,111✔
260	(Sun->SunShapeAngle[i] - Sun->SunShapeAngle[i - 1]);	159,111✔
261
262	} while ((myrng() > (stest / Sun->MaxIntensity)) \|\| (theta2 > (Sun->MaxAngle * Sun->MaxAngle)));	159,111✔
263	break;	100,000✔
264	default:	×
265	throw std::invalid_argument("Unsupported sun shape in Errors function.");	×
266	}
267	}
268
269	if (Source == 2) // surface error	47,358,941✔
270	{
271	// dist = OptProperties->DistributionType; // errors
272	// // delop = sqrt(4.0*sqr(OptProperties->RMSSlopeError)+sqr(OptProperties->RMSSpecError))/1000.0;
273	delop = OptProperties->specularity_error;	24,308,261✔
274
275	Label_50:	24,308,261✔
276	switch (OptProperties->error_distribution_type)	24,308,261✔
277	{
278	case DistributionType::GAUSSIAN: // case 'g':	24,308,261✔
279	thetax = myrng.randNorm(0., delop);	24,308,261✔
280	thetay = myrng.randNorm(0., delop);	24,308,261✔
281
282	theta2 = thetax * thetax + thetay * thetay;	24,308,261✔
283	break;	24,308,261✔
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);	47,358,941✔
306
307	// {Generate errors in terms of direction cosines in local ray coordinate system}
308	theta = sqrt(theta2);	47,358,941✔
309	theta = theta / 1.e3; // convert from mrad to rad	47,358,941✔
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	47,358,941✔
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);	47,358,941✔
317	CosOut[1] = sin(theta) * sin(phi);	47,358,941✔
318	CosOut[2] = cos(theta);	47,358,941✔
319
320	for (i = 0; i < 3; i++)	189,435,764✔
321	{
322	PosIn[i] = PosOut[i];	142,076,823✔
323	CosIn[i] = CosOut[i];	142,076,823✔
324	}
325
326	//{Transform perturbed ray back to element system}
327	TransformToReference(PosIn, CosIn, Origin, RLocToRef, PosOut, CosOut);	47,358,941✔
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) &&	71,667,202✔
334	(OptProperties->my_type == InteractionType::REFLECTION) &&	24,308,261✔
335	(DOT(CosOut, DFXYZ) < 0) &&	71,667,202✔
336	maxcall++ < 50000)	×
337	{
338	goto Label_50;	×
339	}
340	}	47,358,941✔
341	// End of Procedure--------------------------------------------------------------
342
343	} // namespace SolTrace::NativeRunner

NREL / SolTrace / 20348441017

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