18544749821

Committed 15 Oct 2025 10:47PM UTC coverage: 89.946% (+45.8%) from 44.166%

Build # 18544749821

Build Type

push

github

Committed by

web-flow

Commit Message

Restructured backend for SolTrace (#63)

* Initial API for refactored SolTrace

* Updates to various APIs; first pass at element/ray source container implementation

* Use existing matrix-vector code as backend for matrix and vector classes; start implementing APIs

* Created new test on power tower example file with 50000 rays. Changed CMake to use cpp 17.

* Added ground results csv to new folder in google-tests directory

* Removed register prefix from variables in mtrand.h. More errors to come.

* Removed lines in CMakeLists causing linux failure, fixed bug in power tower test

* Created parabola.stinput test and a test using powertower optimization on the powertower sample

* Commented out parabola test for now, added nonexistent branch to CI to see how it runs

* Fixed syntax error

* Changed tests based on PR feedback

* Fixed error in CMakeLists that caused failure on windows

* Removed output messages from st_sim_run_test

* Fixed silly mistake in previous commit

* Added comments and removed unused libraries

* Changed naming conventions of tests

* Added high flux solar furnace test changes, changed parabola test file name

* Fixed typo

* Move refactored data to its own directory; add to refactored classes to build

* Add missed files

* Add unit tests for basic linear algebra and container template

* More tests; Start implementation of composite element

* Add more unit tests on element

* Add missed headers

* Remove target from cmake build command in CI

* Initial virtual element implementation; smoke test for virtual element

* Correct aperture spelling; add tests for virtual elements

* Add some unit tests for simulation data

* Move to static library for simulation data (temporary); update cmake files to get correct mac architecture

* Attempt to fix coverage failure

* Attempt to capture inline functions in code coverage; add a few more explicit tests

* Fix build

* Separate storage of CompositeElement and SingleElement; update te... (continued)

Run Details

4357 of 4844 new or added lines in 62 files covered. (89.95%)

4357 of 4844 relevant lines covered (89.95%)

8565975.82 hits per line

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

94.38

/coretrace/simulation_runner/native_runner/process_interaction.cpp


#include "process_interaction.hpp"

// SimulationData headers
#include "simulation_data_export.hpp"

// NativeRunner headers
#include "tracing_errors.hpp"

namespace SolTrace::NativeRunner
{

    void ProcessInteraction(
        // system info
        TSystem *System,
        MTRand &myrng,
        const bool IncludeSunShape,
        const OpticalProperties *optics,
        const bool IncludeErrors,
        // stage info
        const int i,
        // const TStage *Stage,
        const tstage_ptr Stage,
        // const telement_ptr Elem,
        // const int k,
        // ray info
        const uint_fast64_t MultipleHitCount,
        double (&LastDFXYZ)[3],
        // Outputs
        double (&LastCosRaySurfElement)[3],
        int &ErrorFlag,
        double (&CosRayOutElement)[3],
        double (&LastPosRaySurfElement)[3],
        double (&PosRayOutElement)[3],
        int &myrng_counter)
    {
        // Initialize
        double CosIn[3] = {0.0, 0.0, 0.0};
        double CosOut[3] = {0.0, 0.0, 0.0};

        if (!Stage->Virtual)
        {
            // change to account for first hit only in primary stage 8-11-31
            if (IncludeSunShape && i == 0 && MultipleHitCount == 1)
            {
                // Apply sunshape to UNPERTURBED ray at intersection point
                // only apply sunshape error once for primary stage
                CopyVec3(CosIn, LastCosRaySurfElement);
                // sun shape
                Errors(myrng, CosIn, 1, &System->Sun,
                       optics, CosOut, LastDFXYZ);
                CopyVec3(LastCosRaySurfElement, CosOut);
            }

            //{Determine interaction at surface and direction of perturbed ray}
            ErrorFlag = 0;

            // {Apply surface normal errors to surface normal before interaction
            // ray at intersection point - Wendelin 11-23-09}
            if (IncludeErrors)
            {
                CopyVec3(CosIn, CosRayOutElement);
                // surface normal errors
                SurfaceNormalErrors(myrng, LastDFXYZ, optics, CosOut);
                myrng_counter++;
                CopyVec3(LastDFXYZ, CosOut);
            }

            Interaction(myrng, LastPosRaySurfElement, LastCosRaySurfElement,
                        LastDFXYZ, // Stage->ElementList[k]->InteractionType,
                        optics, 630.0, PosRayOutElement, CosRayOutElement,
                        &ErrorFlag);
            myrng_counter++;

            // {Apply specularity optical error to PERTURBED (i.e. after
            // interaction) ray at intersection point}
            if (IncludeErrors)
            {
                // if (optics->error_distribution_type == 'F' ||
                //         optics->error_distribution_type == 'f')
                // {
                //         // Apply diffuse errors relative to surface normal
                //         CopyVec3(CosIn, LastDFXYZ);
                // }
                // else
                // {
                //         // Apply all other errors relative to the specularly-reflected
                //         // direction
                //         CopyVec3(CosIn, CosRayOutElement);
                // }

                // TODO: Not sure what error distribution type 'F' is?
                // Do we need to implement it? For now just use the 'else'
                // clause from the above.
                CopyVec3(CosIn, CosRayOutElement);

                // optical errors
                Errors(myrng, CosIn, 2, &System->Sun,
                       //    Stage->ElementList[k].get(),
                       optics, CosOut, LastDFXYZ);
                myrng_counter++;
                CopyVec3(CosRayOutElement, CosOut);
            }
        }
    }

    inline double sqr(double x) { return (x) * (x); }

    void Interaction(
        MTRand &myrng,
        const double PosXYZ[3],
        const double CosKLM[3],
        const double DFXYZ[3],
        // int InteractionType,
        const OpticalProperties *Opticl,
        double Wavelength,
        double PosOut[3],
        double CosOut[3],
        int *ErrorFlag)
    {
        /* {Purpose: To compute the direction cosines of the ray due to optical interaction
                   at the intersection point of the ray with the surface
             Input - PosXYZ[2] = X,Y,Z coordinates of intersection point.
                     DFXYZ     = direction numbers for the surface normal at the
                                 intersection point (partial derivatives with respect
                                 to X,Y,Z of surface equation)
                     InteractionType = Optical interaction type indicator
                               = 1, refraction
                               = 2, reflection
                               = 3, aperture stop
                               = 4, diffraction, transmission grating
                               = 5, diffraction, reflection grating
                     CosKLM[2] = direction cosines of incident ray
                     Opticl    = record of optical properties
                           .RefractiveIndex[4] = Refractive index of incident and outgoing medium
                                           [0] = real part of incident medium refractive index
                                           [1] = imaginary part of ""
                                           [2] = real part of outgoing medium refractive index
                                           [4] = imaginary part of ""
                           .ApertureStopOrGratingType
                                               for InteractionType = 3, aperture stop
                                                   = 1, slit
                                                   = 2, elliptical
                                               for InteractionType = 4,5 grating
                                                   = 1, planes parallel to Y-Z plane
                                                   = 2, concentric cylinders centered about Z-axis
                           .DiffractionOrder = integral order of diffraction for InteractionTypes=4,5, grating
                           .AB12[4] = coefficients of polynomial specifying grating spacing for InteractionTypes=4,5
                                [0] = lower X limit, ApertureStopOrGratingType = 1
                                      semi-X axis, ApertureStopOrGratingType = 2
                                [1] = lower Y limit, ApertureStopOrGratingType = 1
                                      semi-Y axis, ApertureStopOrGratingType = 2
                                [2] = upper X limit, ApertureStopOrGratingType = 1
                                      unused, ApertureStopOrGratingType = 2
                                [4] = upper Y limit, ApertureStopOrGratingType = 1
                                      unused, ApertureStopOrGratingType = 2
                     Wavelength = wavelength of ray

             Output - PosOut[2] = position of ray after optical interaction
                      CosOut[2] = direction cosines of ray after optical interaction
                      ErrorFlag = Error flag indicating successful interaction
        }*/

        int i = 0;
        double CosUVW[3] = {0.0, 0.0, 0.0};
        int NIter = 0, IType = 0, NMord = 0;
        double Epsilon = 0.0, Refr1 = 0.0, Refr2 = 0.0, RMU = 0.0, RM2 = 0.0;
        double D2 = 0.0, B = 0.0, A = 0.0, A2 = 0.0;
        double Gamn = 0.0, Gamn1 = 0.0;
        double X = 0.0, Y = 0.0, A1 = 0.0, B1 = 0.0, Ellips = 0.0, B2 = 0.0;
        double RK = 0.0, RL = 0.0, RM = 0.0, Denom, U = 0, V = 0, W = 0;
        double Varr = 0, GFactr = 0, Rho2 = 0.0, Rho = 0.0, Term = 0.0, G = 0.0, D = 0.0, XX = 0.0, Ordiff = 0.0, RLamda = 0.0;
        double Rave = 0.0, Rs = 0.0, Rp = 0.0;
        double UnitDFXYZ[3] = {0.0, 0.0, 0.0}, IncidentAngle = 0.0;

        NIter = 10;
        Epsilon = 0.000005;

        *ErrorFlag = 0;
        for (i = 0; i < 3; i++)
            PosOut[i] = PosXYZ[i];

        switch (Opticl->my_type)
        {

            /*{  InteractionType = 1, Refraction
            ===============================================================================}*/
        case InteractionType::REFRACTION:
        {
            // TODO: Check that this grabs the correct/savem values
            // as the commented out bit immediately below.
            // Refr1 = Opticl->RefractiveIndex[0];
            // Refr2 = Opticl->RefractiveIndex[2];
            Refr1 = Opticl->refraction_index_front;
            Refr2 = Opticl->refraction_index_back;
            RMU = Refr1 / Refr2;
            D2 = DOT(DFXYZ, DFXYZ);
            B = (RMU * RMU - 1.0) / D2;
            A = RMU * DOT(CosKLM, DFXYZ) / D2;
            A2 = A * A;
            if (B > A2) // Total internal reflection
            {
                A = DOT(CosKLM, DFXYZ) / DOT(DFXYZ, DFXYZ);
                for (i = 0; i < 3; i++)
                    CosOut[i] = CosKLM[i] - 2.0 * A * DFXYZ[i];
                return;
            }

            // fresnel equations
            UnitDFXYZ[0] = -DFXYZ[0] / sqrt(DOT(DFXYZ, DFXYZ)); // unit surface normals
            UnitDFXYZ[1] = -DFXYZ[1] / sqrt(DOT(DFXYZ, DFXYZ));
            UnitDFXYZ[2] = -DFXYZ[2] / sqrt(DOT(DFXYZ, DFXYZ));
            IncidentAngle = acos(DOT(CosKLM, UnitDFXYZ));
            Rs = sqr(((Refr1 * cos(IncidentAngle) - Refr2 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2)))) /
                     ((Refr1 * cos(IncidentAngle) + Refr2 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2)))));
            Rp = sqr(((Refr1 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2))) - Refr2 * cos(IncidentAngle)) /
                     ((Refr1 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2))) + Refr2 * cos(IncidentAngle)));
            Rave = (Rp + Rs) / 2.0; // average of s and p polarized light; equal parts of both = non-polarized
            if (Rave < myrng())     // transmitted through surface
            {
                Gamn = -B / (2.0 * A);

                // Begin Newton-Raphson loop to converge on correct root.
                bool converged = false;
                for (i = 1; i < NIter; i++)
                {
                    Gamn1 = (Gamn * Gamn - B) / (2.0 * (Gamn + A));
                    if (fabs(Gamn - Gamn1) < Epsilon)
                    {
                        converged = true;
                        break;
                    }

                    Gamn = Gamn1;
                }
                // Failed to converge
                if (converged == false)
                {
                    *ErrorFlag = 12;
                    return;
                }

                // Have converged on Gamma, Compute direction cosines of refracted ray.
                // Label_Converge:
                for (i = 0; i < 3; i++)
                    CosOut[i] = RMU * CosKLM[i] + Gamn1 * DFXYZ[i];
            }
            else // reflected from surface
            {
                A = DOT(CosKLM, DFXYZ) / DOT(DFXYZ, DFXYZ);
                for (i = 0; i < 3; i++)
                    CosOut[i] = CosKLM[i] - 2.0 * A * DFXYZ[i];
            }
            return;
            break;
        }

            /*{  InteractionType = 2, Reflection
            ===============================================================================}*/
        case InteractionType::REFLECTION:
        {
            A = DOT(CosKLM, DFXYZ) / DOT(DFXYZ, DFXYZ);
            // Compute direction cosines for reflected ray
            for (i = 0; i < 3; i++)
                CosOut[i] = CosKLM[i] - 2.0 * A * DFXYZ[i];

            return;
            break;
        }

        //         /*{  InteractionType = 3, Aperture Stop
        //         ===============================================================================}*/
        // case 3:
        // {
        //         X = PosXYZ[0];
        //         Y = PosXYZ[1];
        //         IType = Opticl->ApertureStopOrGratingType;
        //         A1 = Opticl->AB12[0];
        //         B1 = Opticl->AB12[1];

        //         bool ray_missed_aperture = false;
        //         if (IType == 1) // Slit Aperture
        //         {
        //                 A2 = Opticl->AB12[2];
        //                 B2 = Opticl->AB12[3];
        //                 if (X < A1 || X > A2)
        //                 {
        //                         *ErrorFlag = 31;
        //                         ray_missed_aperture = true;
        //                 }
        //                 else
        //                 {
        //                         if (Y >= B1 && Y <= B2)
        //                                 return;

        //                         *ErrorFlag = 31;
        //                         ray_missed_aperture = true;
        //                 }
        //         }

        //         else if (IType == 2) // Elliptical Aperture
        //         {
        //                 Ellips = X * X / (A1 * A1) + Y * Y / (B1 * B1);
        //                 if (Ellips <= 1.0)
        //                         return;
        //                 *ErrorFlag = 32;
        //         }

        //         // RayMissesAperture:
        //         // Ray misses aperture
        //         if (ray_missed_aperture == true)
        //         {
        //                 for (i = 0; i < 3; i++)
        //                         CosOut[i] = 0.0;
        //         }
        //         return;

        //         break;
        // }

        //         /*{  InteractionType = 4,5; Diffraction
        //         ===============================================================================}*/
        // case 4:
        // case 5:
        // {
        //         IType = Opticl->ApertureStopOrGratingType;
        //         NMord = Opticl->DiffractionOrder;
        //         Refr1 = Opticl->RefractiveIndex[0];
        //         Refr2 = Opticl->RefractiveIndex[2];
        //         RMU = Refr1 / Refr2;
        //         D2 = DOT(DFXYZ, DFXYZ);
        //         RK = DFXYZ[0];
        //         RL = DFXYZ[1];
        //         RM = DFXYZ[2];
        //         X = PosXYZ[0];
        //         Y = PosXYZ[1];

        //         if (IType == 1) // Parallel plane grating
        //         {
        //                 Denom = RL * RL + RM * RM;
        //                 U = 1.0 / sqrt(1.0 + RK * RK / Denom);
        //                 V = -RK * RL * U / Denom;
        //                 W = -RK * RM * U / Denom;
        //                 Varr = X;
        //                 GFactr = 1.0 / U;
        //         }

        //         else if (IType == 2) // Concentric Cylinder Grating
        //         {
        //                 Rho2 = X * X + Y * Y;
        //                 Rho = sqrt(Rho2);
        //                 RM2 = RM * RM;
        //                 Term = RL * X - RK * Y;
        //                 G = sqrt(D2 * (RM2 * Rho2 + Term * Term));
        //                 U = (RM2 * X + RL * Term) / G;
        //                 V = (RM2 * Y - RK * Term) / G;
        //                 W = -RM * (RK * X + RL * Y) / G;
        //                 Varr = Rho;
        //                 GFactr = Rho / (X * U + Y * V);
        //         }
        //         // CompDiffInt:         //Compute interaction due to diffraction
        //         CosUVW[0] = U;
        //         CosUVW[1] = V;
        //         CosUVW[2] = W;

        //         D = 0.0;
        //         XX = 1.0;

        //         for (i = 0; i < 4; i++)
        //         {
        //                 D = D + Opticl->AB12[i] * XX;
        //                 XX = XX * Varr;
        //         }

        //         D = D * GFactr;
        //         Ordiff = NMord;
        //         RLamda = Ordiff * Wavelength / (Refr2 * D);
        //         B = (RMU * RMU - 1.0 + RLamda * RLamda - 2.0 * RMU * RLamda * DOT(CosKLM, CosUVW)) / D2;
        //         A = RMU * DOT(CosKLM, DFXYZ) / D2;
        //         A2 = A * A;
        //         if (B > A2) // Total internal reflection
        //         {
        //                 for (i = 0; i < 3; i++)
        //                         CosOut[i] = 0.0;
        //                 *ErrorFlag = 11;
        //                 return;
        //         }

        //         Gamn = -B / (2.0 * A);
        //         if (InteractionType == 5)
        //                 Gamn = -Gamn - 2.0 * A;

        //         // Begin Newton-Raphson loop to converge on correct root.
        //         i = 0;
        //         bool converged = false;
        //         while (i++ < NIter)
        //         {
        //                 Gamn1 = (Gamn * Gamn - B) / (2.0 * (Gamn + A));
        //                 if (fabs(Gamn - Gamn1) < Epsilon)
        //                 {
        //                         converged = true;
        //                         break;
        //                 }
        //                 Gamn = Gamn1;
        //         }
        //         // Failed to converge
        //         if (converged == false)
        //         {
        //                 *ErrorFlag = 12;
        //                 return;
        //         }
        //         // Have converged on Gamn1. Compute direction cosines of diffracted ray.
        //         // CompDCos:
        //         for (i = 0; i < 3; i++)
        //                 CosOut[i] = RMU * CosKLM[i] - RLamda * CosUVW[i] + Gamn1 * DFXYZ[i];

        //         break;
        // }
        default:
            break;
        }
        return;
    }

} // namespace SolTrace::NativeRunner

1
2	#include "process_interaction.hpp"
3
4	// SimulationData headers
5	#include "simulation_data_export.hpp"
6
7	// NativeRunner headers
8	#include "tracing_errors.hpp"
9
10	namespace SolTrace::NativeRunner
11	{
12
13	void ProcessInteraction(	634,824✔
14	// system info
15	TSystem *System,
16	MTRand &myrng,
17	const bool IncludeSunShape,
18	const OpticalProperties *optics,
19	const bool IncludeErrors,
20	// stage info
21	const int i,
22	// const TStage *Stage,
23	const tstage_ptr Stage,
24	// const telement_ptr Elem,
25	// const int k,
26	// ray info
27	const uint_fast64_t MultipleHitCount,
28	double (&LastDFXYZ)[3],
29	// Outputs
30	double (&LastCosRaySurfElement)[3],
31	int &ErrorFlag,
32	double (&CosRayOutElement)[3],
33	double (&LastPosRaySurfElement)[3],
34	double (&PosRayOutElement)[3],
35	int &myrng_counter)
36	{
37	// Initialize
38	double CosIn[3] = {0.0, 0.0, 0.0};	634,824✔
39	double CosOut[3] = {0.0, 0.0, 0.0};	634,824✔
40
41	if (!Stage->Virtual)	634,824✔
42	{
43	// change to account for first hit only in primary stage 8-11-31
44	if (IncludeSunShape && i == 0 && MultipleHitCount == 1)	634,824✔
45	{
46	// Apply sunshape to UNPERTURBED ray at intersection point
47	// only apply sunshape error once for primary stage
48	CopyVec3(CosIn, LastCosRaySurfElement);	159,188✔
49	// sun shape
50	Errors(myrng, CosIn, 1, &System->Sun,	159,188✔
51	optics, CosOut, LastDFXYZ);
52	CopyVec3(LastCosRaySurfElement, CosOut);	159,188✔
53	}
54
55	//{Determine interaction at surface and direction of perturbed ray}
56	ErrorFlag = 0;	634,824✔
57
58	// {Apply surface normal errors to surface normal before interaction
59	// ray at intersection point - Wendelin 11-23-09}
60	if (IncludeErrors)	634,824✔
61	{
62	CopyVec3(CosIn, CosRayOutElement);	201,899✔
63	// surface normal errors
64	SurfaceNormalErrors(myrng, LastDFXYZ, optics, CosOut);	201,899✔
65	myrng_counter++;	201,899✔
66	CopyVec3(LastDFXYZ, CosOut);	201,899✔
67	}
68
69	Interaction(myrng, LastPosRaySurfElement, LastCosRaySurfElement,	634,824✔
70	LastDFXYZ, // Stage->ElementList[k]->InteractionType,
71	optics, 630.0, PosRayOutElement, CosRayOutElement,
72	&ErrorFlag);
73	myrng_counter++;	634,824✔
74
75	// {Apply specularity optical error to PERTURBED (i.e. after
76	// interaction) ray at intersection point}
77	if (IncludeErrors)	634,824✔
78	{
79	// if (optics->error_distribution_type == 'F' \|\|
80	// optics->error_distribution_type == 'f')
81	// {
82	// // Apply diffuse errors relative to surface normal
83	// CopyVec3(CosIn, LastDFXYZ);
84	// }
85	// else
86	// {
87	// // Apply all other errors relative to the specularly-reflected
88	// // direction
89	// CopyVec3(CosIn, CosRayOutElement);
90	// }
91
92	// TODO: Not sure what error distribution type 'F' is?
93	// Do we need to implement it? For now just use the 'else'
94	// clause from the above.
95	CopyVec3(CosIn, CosRayOutElement);	201,899✔
96
97	// optical errors
98	Errors(myrng, CosIn, 2, &System->Sun,	201,899✔
99	// Stage->ElementList[k].get(),
100	optics, CosOut, LastDFXYZ);
101	myrng_counter++;	201,899✔
102	CopyVec3(CosRayOutElement, CosOut);	201,899✔
103	}
104	}
105	}	634,824✔
106
107	inline double sqr(double x) { return (x) * (x); }	335,436✔
108
109	void Interaction(	634,824✔
110	MTRand &myrng,
111	const double PosXYZ[3],
112	const double CosKLM[3],
113	const double DFXYZ[3],
114	// int InteractionType,
115	const OpticalProperties *Opticl,
116	double Wavelength,
117	double PosOut[3],
118	double CosOut[3],
119	int *ErrorFlag)
120	{
121	/* {Purpose: To compute the direction cosines of the ray due to optical interaction
122	at the intersection point of the ray with the surface
123	Input - PosXYZ[2] = X,Y,Z coordinates of intersection point.
124	DFXYZ = direction numbers for the surface normal at the
125	intersection point (partial derivatives with respect
126	to X,Y,Z of surface equation)
127	InteractionType = Optical interaction type indicator
128	= 1, refraction
129	= 2, reflection
130	= 3, aperture stop
131	= 4, diffraction, transmission grating
132	= 5, diffraction, reflection grating
133	CosKLM[2] = direction cosines of incident ray
134	Opticl = record of optical properties
135	.RefractiveIndex[4] = Refractive index of incident and outgoing medium
136	[0] = real part of incident medium refractive index
137	[1] = imaginary part of ""
138	[2] = real part of outgoing medium refractive index
139	[4] = imaginary part of ""
140	.ApertureStopOrGratingType
141	for InteractionType = 3, aperture stop
142	= 1, slit
143	= 2, elliptical
144	for InteractionType = 4,5 grating
145	= 1, planes parallel to Y-Z plane
146	= 2, concentric cylinders centered about Z-axis
147	.DiffractionOrder = integral order of diffraction for InteractionTypes=4,5, grating
148	.AB12[4] = coefficients of polynomial specifying grating spacing for InteractionTypes=4,5
149	[0] = lower X limit, ApertureStopOrGratingType = 1
150	semi-X axis, ApertureStopOrGratingType = 2
151	[1] = lower Y limit, ApertureStopOrGratingType = 1
152	semi-Y axis, ApertureStopOrGratingType = 2
153	[2] = upper X limit, ApertureStopOrGratingType = 1
154	unused, ApertureStopOrGratingType = 2
155	[4] = upper Y limit, ApertureStopOrGratingType = 1
156	unused, ApertureStopOrGratingType = 2
157	Wavelength = wavelength of ray
158
159	Output - PosOut[2] = position of ray after optical interaction
160	CosOut[2] = direction cosines of ray after optical interaction
161	ErrorFlag = Error flag indicating successful interaction
162	}*/
163
164	int i = 0;	634,824✔
165	double CosUVW[3] = {0.0, 0.0, 0.0};	634,824✔
166	int NIter = 0, IType = 0, NMord = 0;	634,824✔
167	double Epsilon = 0.0, Refr1 = 0.0, Refr2 = 0.0, RMU = 0.0, RM2 = 0.0;	634,824✔
168	double D2 = 0.0, B = 0.0, A = 0.0, A2 = 0.0;	634,824✔
169	double Gamn = 0.0, Gamn1 = 0.0;	634,824✔
170	double X = 0.0, Y = 0.0, A1 = 0.0, B1 = 0.0, Ellips = 0.0, B2 = 0.0;	634,824✔
171	double RK = 0.0, RL = 0.0, RM = 0.0, Denom, U = 0, V = 0, W = 0;	634,824✔
172	double Varr = 0, GFactr = 0, Rho2 = 0.0, Rho = 0.0, Term = 0.0, G = 0.0, D = 0.0, XX = 0.0, Ordiff = 0.0, RLamda = 0.0;	634,824✔
173	double Rave = 0.0, Rs = 0.0, Rp = 0.0;	634,824✔
174	double UnitDFXYZ[3] = {0.0, 0.0, 0.0}, IncidentAngle = 0.0;	634,824✔
175
176	NIter = 10;	634,824✔
177	Epsilon = 0.000005;	634,824✔
178
179	*ErrorFlag = 0;	634,824✔
180	for (i = 0; i < 3; i++)	2,539,296✔
181	PosOut[i] = PosXYZ[i];	1,904,472✔
182
183	switch (Opticl->my_type)	634,824✔
184	{
185
186	/*{ InteractionType = 1, Refraction
187	===============================================================================}*/
188	case InteractionType::REFRACTION:	61,685✔
189	{
190	// TODO: Check that this grabs the correct/savem values
191	// as the commented out bit immediately below.
192	// Refr1 = Opticl->RefractiveIndex[0];
193	// Refr2 = Opticl->RefractiveIndex[2];
194	Refr1 = Opticl->refraction_index_front;	61,685✔
195	Refr2 = Opticl->refraction_index_back;	61,685✔
196	RMU = Refr1 / Refr2;	61,685✔
197	D2 = DOT(DFXYZ, DFXYZ);	61,685✔
198	B = (RMU * RMU - 1.0) / D2;	61,685✔
199	A = RMU * DOT(CosKLM, DFXYZ) / D2;	61,685✔
200	A2 = A * A;	61,685✔
201	if (B > A2) // Total internal reflection	61,685✔
202	{
203	A = DOT(CosKLM, DFXYZ) / DOT(DFXYZ, DFXYZ);	5,779✔
204	for (i = 0; i < 3; i++)	23,116✔
205	CosOut[i] = CosKLM[i] - 2.0 * A * DFXYZ[i];	17,337✔
206	return;	5,779✔
207	}
208
209	// fresnel equations
210	UnitDFXYZ[0] = -DFXYZ[0] / sqrt(DOT(DFXYZ, DFXYZ)); // unit surface normals	55,906✔
211	UnitDFXYZ[1] = -DFXYZ[1] / sqrt(DOT(DFXYZ, DFXYZ));	55,906✔
212	UnitDFXYZ[2] = -DFXYZ[2] / sqrt(DOT(DFXYZ, DFXYZ));	55,906✔
213	IncidentAngle = acos(DOT(CosKLM, UnitDFXYZ));	55,906✔
214	Rs = sqr(((Refr1 * cos(IncidentAngle) - Refr2 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2)))) /	55,906✔
215	((Refr1 * cos(IncidentAngle) + Refr2 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2)))));	55,906✔
216	Rp = sqr(((Refr1 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2))) - Refr2 * cos(IncidentAngle)) /	55,906✔
217	((Refr1 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2))) + Refr2 * cos(IncidentAngle)));	55,906✔
218	Rave = (Rp + Rs) / 2.0; // average of s and p polarized light; equal parts of both = non-polarized	55,906✔
219	if (Rave < myrng()) // transmitted through surface	55,906✔
220	{
221	Gamn = -B / (2.0 * A);	52,014✔
222
223	// Begin Newton-Raphson loop to converge on correct root.
224	bool converged = false;	52,014✔
225	for (i = 1; i < NIter; i++)	197,754✔
226	{
227	Gamn1 = (Gamn * Gamn - B) / (2.0 * (Gamn + A));	197,754✔
228	if (fabs(Gamn - Gamn1) < Epsilon)	197,754✔
229	{
230	converged = true;	52,014✔
231	break;	52,014✔
232	}
233
234	Gamn = Gamn1;	145,740✔
235	}
236	// Failed to converge
237	if (converged == false)	52,014✔
238	{
NEW 239	*ErrorFlag = 12;	×
NEW 240	return;	×
241	}
242
243	// Have converged on Gamma, Compute direction cosines of refracted ray.
244	// Label_Converge:
245	for (i = 0; i < 3; i++)	208,056✔
246	CosOut[i] = RMU * CosKLM[i] + Gamn1 * DFXYZ[i];	156,042✔
247	}
248	else // reflected from surface
249	{
250	A = DOT(CosKLM, DFXYZ) / DOT(DFXYZ, DFXYZ);	3,892✔
251	for (i = 0; i < 3; i++)	15,568✔
252	CosOut[i] = CosKLM[i] - 2.0 * A * DFXYZ[i];	11,676✔
253	}
254	return;	55,906✔
255	break;
256	}
257
258	/*{ InteractionType = 2, Reflection
259	===============================================================================}*/
260	case InteractionType::REFLECTION:	573,139✔
261	{
262	A = DOT(CosKLM, DFXYZ) / DOT(DFXYZ, DFXYZ);	573,139✔
263	// Compute direction cosines for reflected ray
264	for (i = 0; i < 3; i++)	2,292,556✔
265	CosOut[i] = CosKLM[i] - 2.0 * A * DFXYZ[i];	1,719,417✔
266
267	return;	573,139✔
268	break;
269	}
270
271	// /*{ InteractionType = 3, Aperture Stop
272	// ===============================================================================}*/
273	// case 3:
274	// {
275	// X = PosXYZ[0];
276	// Y = PosXYZ[1];
277	// IType = Opticl->ApertureStopOrGratingType;
278	// A1 = Opticl->AB12[0];
279	// B1 = Opticl->AB12[1];
280
281	// bool ray_missed_aperture = false;
282	// if (IType == 1) // Slit Aperture
283	// {
284	// A2 = Opticl->AB12[2];
285	// B2 = Opticl->AB12[3];
286	// if (X < A1 \|\| X > A2)
287	// {
288	// *ErrorFlag = 31;
289	// ray_missed_aperture = true;
290	// }
291	// else
292	// {
293	// if (Y >= B1 && Y <= B2)
294	// return;
295
296	// *ErrorFlag = 31;
297	// ray_missed_aperture = true;
298	// }
299	// }
300
301	// else if (IType == 2) // Elliptical Aperture
302	// {
303	// Ellips = X * X / (A1 * A1) + Y * Y / (B1 * B1);
304	// if (Ellips <= 1.0)
305	// return;
306	// *ErrorFlag = 32;
307	// }
308
309	// // RayMissesAperture:
310	// // Ray misses aperture
311	// if (ray_missed_aperture == true)
312	// {
313	// for (i = 0; i < 3; i++)
314	// CosOut[i] = 0.0;
315	// }
316	// return;
317
318	// break;
319	// }
320
321	// /*{ InteractionType = 4,5; Diffraction
322	// ===============================================================================}*/
323	// case 4:
324	// case 5:
325	// {
326	// IType = Opticl->ApertureStopOrGratingType;
327	// NMord = Opticl->DiffractionOrder;
328	// Refr1 = Opticl->RefractiveIndex[0];
329	// Refr2 = Opticl->RefractiveIndex[2];
330	// RMU = Refr1 / Refr2;
331	// D2 = DOT(DFXYZ, DFXYZ);
332	// RK = DFXYZ[0];
333	// RL = DFXYZ[1];
334	// RM = DFXYZ[2];
335	// X = PosXYZ[0];
336	// Y = PosXYZ[1];
337
338	// if (IType == 1) // Parallel plane grating
339	// {
340	// Denom = RL * RL + RM * RM;
341	// U = 1.0 / sqrt(1.0 + RK * RK / Denom);
342	// V = -RK * RL * U / Denom;
343	// W = -RK * RM * U / Denom;
344	// Varr = X;
345	// GFactr = 1.0 / U;
346	// }
347
348	// else if (IType == 2) // Concentric Cylinder Grating
349	// {
350	// Rho2 = X * X + Y * Y;
351	// Rho = sqrt(Rho2);
352	// RM2 = RM * RM;
353	// Term = RL * X - RK * Y;
354	// G = sqrt(D2 * (RM2 * Rho2 + Term * Term));
355	// U = (RM2 * X + RL * Term) / G;
356	// V = (RM2 * Y - RK * Term) / G;
357	// W = -RM * (RK * X + RL * Y) / G;
358	// Varr = Rho;
359	// GFactr = Rho / (X * U + Y * V);
360	// }
361	// // CompDiffInt: //Compute interaction due to diffraction
362	// CosUVW[0] = U;
363	// CosUVW[1] = V;
364	// CosUVW[2] = W;
365
366	// D = 0.0;
367	// XX = 1.0;
368
369	// for (i = 0; i < 4; i++)
370	// {
371	// D = D + Opticl->AB12[i] * XX;
372	// XX = XX * Varr;
373	// }
374
375	// D = D * GFactr;
376	// Ordiff = NMord;
377	// RLamda = Ordiff * Wavelength / (Refr2 * D);
378	// B = (RMU * RMU - 1.0 + RLamda * RLamda - 2.0 * RMU * RLamda * DOT(CosKLM, CosUVW)) / D2;
379	// A = RMU * DOT(CosKLM, DFXYZ) / D2;
380	// A2 = A * A;
381	// if (B > A2) // Total internal reflection
382	// {
383	// for (i = 0; i < 3; i++)
384	// CosOut[i] = 0.0;
385	// *ErrorFlag = 11;
386	// return;
387	// }
388
389	// Gamn = -B / (2.0 * A);
390	// if (InteractionType == 5)
391	// Gamn = -Gamn - 2.0 * A;
392
393	// // Begin Newton-Raphson loop to converge on correct root.
394	// i = 0;
395	// bool converged = false;
396	// while (i++ < NIter)
397	// {
398	// Gamn1 = (Gamn * Gamn - B) / (2.0 * (Gamn + A));
399	// if (fabs(Gamn - Gamn1) < Epsilon)
400	// {
401	// converged = true;
402	// break;
403	// }
404	// Gamn = Gamn1;
405	// }
406	// // Failed to converge
407	// if (converged == false)
408	// {
409	// *ErrorFlag = 12;
410	// return;
411	// }
412	// // Have converged on Gamn1. Compute direction cosines of diffracted ray.
413	// // CompDCos:
414	// for (i = 0; i < 3; i++)
415	// CosOut[i] = RMU * CosKLM[i] - RLamda * CosUVW[i] + Gamn1 * DFXYZ[i];
416
417	// break;
418	// }
NEW 419	default:	×
NEW 420	break;	×
421	}
NEW 422	return;	×
423	}
424
425	} // namespace SolTrace::NativeRunner

NREL / SolTrace / 18544749821

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