18757590279

Committed 23 Oct 2025 06:04PM UTC coverage: 88.943% (-1.0%) from 89.946%

Build # 18757590279

Build Type

Pull #75

github

Committed by

web-flow

Commit Message

Merge 1c15e794a into f6f121007

Pull Request Pull Request #75: New sun models

Run Details

72 of 133 new or added lines in 8 files covered. (54.14%)

5 existing lines in 2 files now uncovered.

4384 of 4929 relevant lines covered (88.94%)

8585596.28 hits per line

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

66.45

/coretrace/simulation_runner/native_runner/tracing_errors.cpp


#include "tracing_errors.hpp"

#include "simulation_data_export.hpp"

namespace SolTrace::NativeRunner {

// TODO: These two functions have a lot of overlap
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, delop3 = 0.0, thetax = 0.0,
                   thetay = 0.0, ttheta = 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;

        int nninner = 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 * 3.1415926535897932385; // 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};
        // char dist = 'g';
        double delop = 0, delop3 = 0, thetax = 0, thetay = 0, ttheta = 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 / 1000.0;

                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));
                        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)));

                        theta2 = theta2 / 1.e6; // convert from mrad^2 to rad^2
                        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)));

                        theta2 = theta2 / 1.e6; // convert from mrad^2 to rad^2
                        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)));

            theta2 = theta2 / 1.e6;        // convert from mrad^2 to rad^2
                        break;

                default:
                        // TODO: Add error message here.
            //throw std::exception("Unsupported sun shape in Errors function.");
                        break;
                }
        }

        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 / 1000.0;

        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);

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

NREL / SolTrace / 18757590279

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