18984697163


#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};
        // char dist = 'g';
        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 / 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 * 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,	601,899✔
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;	601,899✔
38	double Origin[3] = {0.0, 0.0, 0.0},	601,899✔
39	Euler[3] = {0.0, 0.0, 0.0};	601,899✔
40	double PosIn[3] = {0.0, 0.0, 0.0},	601,899✔
41	PosOut[3] = {0.0, 0.0, 0.0};	601,899✔
42	DistributionType dist;
43	double delop = 0.0, thetax = 0.0,	601,899✔
44	thetay = 0.0, theta2 = 0.0,	601,899✔
45	phi = 0.0, theta = 0.0;	601,899✔
46	double RRefToLoc[3][3] = {{0.0, 0.0, 0.0},	601,899✔
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},	601,899✔
50	{0.0, 0.0, 0.0},
51	{0.0, 0.0, 0.0}};
52
53	if (CosIn[2] == 0.0)	601,899✔
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;	×
NEW 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]);	601,899✔
69	Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));	601,899✔
70	}
71
72	Euler[2] = 0.0;	601,899✔
73
74	CalculateTransformMatrices(Euler, RRefToLoc, RLocToRef);	601,899✔
75
76	// TODO: Add distribution type to optical properties
77	// dist = OptProperties->DistributionType;
78	dist = OptProperties->error_distribution_type;	601,899✔
79	// delop = OptProperties->RMSSlopeError / 1000.0;
80	delop = OptProperties->slope_error / 1000.0;	601,899✔
81
82	switch (dist)	601,899✔
83	{
84	case DistributionType::GAUSSIAN: // case 'g':	601,899✔
85	// gaussian distribution
86	thetax = myrng.randNorm(0., delop);	601,899✔
87	thetay = myrng.randNorm(0., delop);	601,899✔
88
89	theta2 = thetax * thetax + thetay * thetay;	601,899✔
90	break;	601,899✔
NEW 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));	×
UNCOV 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);	601,899✔
109
110	/* {Generate errors in terms of direction cosines in local ray coordinate system} */
111	theta = sqrt(theta2);	601,899✔
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	601,899✔
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);	601,899✔
117	CosOut[1] = sin(theta) * sin(phi);	601,899✔
118	CosOut[2] = cos(theta);	601,899✔
119
120	for (i = 0; i < 3; i++)	2,407,596✔
121	{
122	PosIn[i] = PosOut[i];	1,805,697✔
123	CosIn[i] = CosOut[i];	1,805,697✔
124	}
125
126	/{Transform perturbed ray back to element system}/
127	TransformToReference(PosIn, CosIn, Origin, RLocToRef, PosOut, CosOut);	601,899✔
128	}	601,899✔
129
130	void Errors(	1,161,087✔
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,161,087✔
161	double Euler[3] = {0.0, 0.0, 0.0};	1,161,087✔
162	double PosIn[3] = {0.0, 0.0, 0.0};	1,161,087✔
163	double PosOut[3] = {0.0, 0.0, 0.0};	1,161,087✔
164	// char dist = 'g';
165	double delop = 0, thetax = 0, thetay = 0, theta2 = 0, phi = 0, theta = 0, stest = 0;	1,161,087✔
166	uint_fast64_t i;
167	double RRefToLoc[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};	1,161,087✔
168	double RLocToRef[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};	1,161,087✔
169
170	if (CosIn[2] == 0.0)	1,161,087✔
171	{
172	if (CosIn[0] == 0.0)	×
173	{
174	Euler[0] = 0.0;	×
175	Euler[1] = PI / 2.0;	×
176	}
177	else
178	{
179	Euler[0] = PI / 2.0;	×
NEW 180	Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));	×
181	}
182	}
183	else
184	{
185	Euler[0] = atan2(CosIn[0], CosIn[2]);	1,161,087✔
186	Euler[1] = atan2(CosIn[1], sqrt(CosIn[0] * CosIn[0] + CosIn[2] * CosIn[2]));	1,161,087✔
187	}
188
189	Euler[2] = 0.0;	1,161,087✔
190
191	CalculateTransformMatrices(Euler, RRefToLoc, RLocToRef);	1,161,087✔
192
193	unsigned int maxcall = 0;	1,161,087✔
194	// g,p,d
195	if (Source == 1) // sun error	1,161,087✔
196	{
197	delop = Sun->Sigma / 1000.0;	559,188✔
198
199	switch (Sun->ShapeIndex)	559,188✔
200	{
201	case SunShape::GAUSSIAN: // case 'g':	119,974✔
202	thetax = myrng.randNorm(0., delop);	119,974✔
203	thetay = myrng.randNorm(0., delop);	119,974✔
204
205	theta2 = thetax * thetax + thetay * thetay;	119,974✔
206	break;	119,974✔
207
208	case SunShape::PILLBOX: // case 'p':	177,215✔
209	do
210	{
211	thetax = 2.0 * delop * myrng() - delop;	177,215✔
212	thetay = 2.0 * delop * myrng() - delop;	177,215✔
213	theta2 = thetax * thetax + thetay * thetay;	177,215✔
214	} while (theta2 > (delop * delop));	177,215✔
215	break;	139,214✔
216
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
227	theta2 = theta2 / 1.e6; // convert from mrad^2 to rad^2	100,000✔
228	break;	100,000✔
229
230	case SunShape::BUIE_CSR:	11,850,588✔
231	// This sun model has long tails so this might take more iterations
232	// TODO: add an option to set the max angle (thereby reducing the tail)
233	do
234	{
235	thetax = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	11,850,588✔
236	thetay = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	11,850,588✔
237	theta2 = thetax * thetax + thetay * thetay;	11,850,588✔
238	theta = sqrt(theta2);	11,850,588✔
239
240	if (std::abs(theta) <= 4.65) // within solar disc	11,850,588✔
241	stest = cos(0.326 * theta) / cos(0.308 * theta);	105,365✔
242	else // within circumsolar region
243	stest = std::exp(Sun->buie_kappa) * std::pow(std::abs(theta), Sun->buie_gamma);	11,745,223✔
244
245	} while ((myrng() > (stest / Sun->MaxIntensity)) \|\| (theta2 > (Sun->MaxAngle * Sun->MaxAngle)));	11,850,588✔
246
247	theta2 = theta2 / 1.e6; // convert from mrad^2 to rad^2	100,000✔
248	break;	100,000✔
249
250	case SunShape::USER_DEFINED:	159,111✔
251	do
252	{
253	thetax = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	159,111✔
254	thetay = 2.0 * Sun->MaxAngle * myrng() - Sun->MaxAngle;	159,111✔
255	theta2 = thetax * thetax + thetay * thetay;	159,111✔
256	theta = sqrt(theta2); // wendelin 1-9-12 do the test once on theta NOT individually on thetax and thetay as before	159,111✔
257
258	i = 0;	159,111✔
259	while (i < Sun->SunShapeAngle.size() - 1 && Sun->SunShapeAngle[i] < theta)	2,717,967✔
260	i++;	2,558,856✔
261
262	if (i == 0)	159,111✔
NEW 263	stest = Sun->SunShapeIntensity[0];	×
264	else // linear interpolation (switched from average) 12-20-11 wendelin
265	stest = Sun->SunShapeIntensity[i - 1] + (Sun->SunShapeIntensity[i] - Sun->SunShapeIntensity[i - 1]) * (theta - Sun->SunShapeAngle[i - 1]) /	159,111✔
266	(Sun->SunShapeAngle[i] - Sun->SunShapeAngle[i - 1]);	159,111✔
267
268	} while ((myrng() > (stest / Sun->MaxIntensity)) \|\| (theta2 > (Sun->MaxAngle * Sun->MaxAngle)));	159,111✔
269
270	theta2 = theta2 / 1.e6; // convert from mrad^2 to rad^2	100,000✔
271	break;	100,000✔
272
NEW 273	default:	×
274	// TODO: Add error message here.
275	//throw std::exception("Unsupported sun shape in Errors function.");
NEW 276	break;	×
277	}
278	}
279
280	if (Source == 2) // surface error	1,161,087✔
281	{
282	// dist = OptProperties->DistributionType; // errors
283	// // delop = sqrt(4.0*sqr(OptProperties->RMSSlopeError)+sqr(OptProperties->RMSSpecError))/1000.0;
284	delop = OptProperties->specularity_error / 1000.0;	601,899✔
285
286	Label_50:	607,743✔
287	switch (OptProperties->error_distribution_type)	607,743✔
288	{
289	case DistributionType::GAUSSIAN: // case 'g':	607,743✔
290	thetax = myrng.randNorm(0., delop);	607,743✔
291	thetay = myrng.randNorm(0., delop);	607,743✔
292
293	theta2 = thetax * thetax + thetay * thetay;	607,743✔
294	break;	607,743✔
295
NEW 296	case DistributionType::PILLBOX: // case 'p':	×
297	do
298	{
NEW 299	thetax = 2.0 * delop * myrng() - delop;	×
NEW 300	thetay = 2.0 * delop * myrng() - delop;	×
NEW 301	theta2 = thetax * thetax + thetay * thetay;	×
NEW 302	} while (theta2 > (delop * delop));	×
NEW 303	break;	×
304
NEW 305	case DistributionType::DIFFUSE:	×
NEW 306	theta2 = pow(asin(sqrt(myrng())), 2);	×
NEW 307	break;	×
308
NEW 309	default:	×
310	// TODO: Add error message here.
NEW 311	break;	×
312	}
313	}
314
315	// {Transform to local coordinate system of ray to set up rotation matrices for coordinate and inverse transforms}
316	TransformToLocal(PosIn, CosIn, Origin, RRefToLoc, PosOut, CosOut);	1,166,931✔
317
318	// {Generate errors in terms of direction cosines in local ray coordinate system}
319	theta = sqrt(theta2);	1,166,931✔
320
321	// phi = atan2(thetay, thetax); //This function appears to present irregularities that bias results incorrectly for small values of thetay or thetax
322	phi = myrng() * 2.0 * PI; // Therefore have chosen to randomize phi rather than calculate from randomized theta components	1,166,931✔
323	// obtained from the distribution. The two approaches are equivalent save for this issue with arctan2. wendelin 01-12-11
324
325	CosOut[0] = sin(theta) * cos(phi);	1,166,931✔
326	CosOut[1] = sin(theta) * sin(phi);	1,166,931✔
327	CosOut[2] = cos(theta);	1,166,931✔
328
329	for (i = 0; i < 3; i++)	4,667,724✔
330	{
331	PosIn[i] = PosOut[i];	3,500,793✔
332	CosIn[i] = CosOut[i];	3,500,793✔
333	}
334
335	//{Transform perturbed ray back to element system}
336	TransformToReference(PosIn, CosIn, Origin, RLocToRef, PosOut, CosOut);	1,166,931✔
337
338	// TODO: Remove goto, should we always do dot product check? // We could move this out of the function and into the caller.
339
340	/*{If reflection error application and new ray direction (after errors) physically goes through opaque surface,
341	then go back and get new perturbation 06-12-07}*/
342	if ((Source == 2) &&	1,774,674✔
343	(OptProperties->my_type == InteractionType::REFLECTION) &&	607,743✔
344	(DOT(CosOut, DFXYZ) < 0) &&	1,780,518✔
345	maxcall++ < 50000)	5,844✔
346	{
347	goto Label_50;	5,844✔
348	}
349	}	1,161,087✔
350	// End of Procedure--------------------------------------------------------------
351
352	} // namespace SolTrace::NativeRunner

NREL / SolTrace / 18984697163

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