14628243107

Committed 23 Apr 2025 09:13PM UTC coverage: 44.154% (+0.5%) from 43.617%

Build # 14628243107

Build Type

Pull #65

github

Committed by

web-flow

Commit Message

Merge e65736e0a into 96166e291

Pull Request Pull Request #65: Refactor the Trace function (and related functions)

Run Details

363 of 670 new or added lines in 6 files covered. (54.18%)

135 existing lines in 4 files now uncovered.

1941 of 4396 relevant lines covered (44.15%)

10000057.29 hits per line

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

14.29

/coretrace/interaction.cpp


/*******************************************************************************************************
*  Copyright 2018 Alliance for Sustainable Energy, LLC
*
*  NOTICE: This software was developed at least in part by Alliance for Sustainable Energy, LLC
*  ("Alliance") under Contract No. DE-AC36-08GO28308 with the U.S. Department of Energy and the U.S.
*  The Government retains for itself and others acting on its behalf a nonexclusive, paid-up,
*  irrevocable worldwide license in the software to reproduce, prepare derivative works, distribute
*  copies to the public, perform publicly and display publicly, and to permit others to do so.
*
*  Redistribution and use in source and binary forms, with or without modification, are permitted
*  provided that the following conditions are met:
*
*  1. Redistributions of source code must retain the above copyright notice, the above government
*  rights notice, this list of conditions and the following disclaimer.
*
*  2. Redistributions in binary form must reproduce the above copyright notice, the above government
*  rights notice, this list of conditions and the following disclaimer in the documentation and/or
*  other materials provided with the distribution.
*
*  3. The entire corresponding source code of any redistribution, with or without modification, by a
*  research entity, including but not limited to any contracting manager/operator of a United States
*  National Laboratory, any institution of higher learning, and any non-profit organization, must be
*  made publicly available under this license for as long as the redistribution is made available by
*  the research entity.
*
*  4. Redistribution of this software, without modification, must refer to the software by the same
*  designation. Redistribution of a modified version of this software (i) may not refer to the modified
*  version by the same designation, or by any confusingly similar designation, and (ii) must refer to
*  the underlying software originally provided by Alliance as "SolTrace". Except to comply with the 
*  foregoing, the term "SolTrace", or any confusingly similar designation may not be used to refer to 
*  any modified version of this software or any modified version of the underlying software originally 
*  provided by Alliance without the prior written consent of Alliance.
*
*  5. The name of the copyright holder, contributors, the United States Government, the United States
*  Department of Energy, or any of their employees may not be used to endorse or promote products
*  derived from this software without specific prior written permission.
*
*  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
*  IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
*  FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER,
*  CONTRIBUTORS, UNITED STATES GOVERNMENT OR UNITED STATES DEPARTMENT OF ENERGY, NOR ANY OF THEIR
*  EMPLOYEES, BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
*  DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
*  DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
*  IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
*  THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*******************************************************************************************************/


#include <math.h>

#include "types.h"
#include "procs.h"

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

void Interaction(
                        MTRand &myrng,
                        double PosXYZ[3],
                        double CosKLM[3],
                        double DFXYZ[3],
                        int InteractionType,
                        TOpticalProperties *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, 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, 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 (InteractionType)
        {

/*{  InteractionType = 1, Refraction
===============================================================================}*/
        case 1:
        {
                Refr1 = Opticl->RefractiveIndex[0];
                Refr2 = Opticl->RefractiveIndex[2];
                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;
                                        //goto Label_Converge;
                                }

                                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 2:
        {
                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;
                                //goto RayMissesAperture;
                        }
                        else
                        {
                                if (Y >= B1 && Y <= B2) return;

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

                }

                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;
                        //goto CompDiffInt;
                }

                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;
                                //goto CompDCos;
                        }
                        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;
        }
}
//End of Procedure--------------------------------------------------------------

1
2	/*******************************************************************************************************
3	* Copyright 2018 Alliance for Sustainable Energy, LLC
4	*
5	* NOTICE: This software was developed at least in part by Alliance for Sustainable Energy, LLC
6	* ("Alliance") under Contract No. DE-AC36-08GO28308 with the U.S. Department of Energy and the U.S.
7	* The Government retains for itself and others acting on its behalf a nonexclusive, paid-up,
8	* irrevocable worldwide license in the software to reproduce, prepare derivative works, distribute
9	* copies to the public, perform publicly and display publicly, and to permit others to do so.
10	*
11	* Redistribution and use in source and binary forms, with or without modification, are permitted
12	* provided that the following conditions are met:
13	*
14	* 1. Redistributions of source code must retain the above copyright notice, the above government
15	* rights notice, this list of conditions and the following disclaimer.
16	*
17	* 2. Redistributions in binary form must reproduce the above copyright notice, the above government
18	* rights notice, this list of conditions and the following disclaimer in the documentation and/or
19	* other materials provided with the distribution.
20	*
21	* 3. The entire corresponding source code of any redistribution, with or without modification, by a
22	* research entity, including but not limited to any contracting manager/operator of a United States
23	* National Laboratory, any institution of higher learning, and any non-profit organization, must be
24	* made publicly available under this license for as long as the redistribution is made available by
25	* the research entity.
26	*
27	* 4. Redistribution of this software, without modification, must refer to the software by the same
28	* designation. Redistribution of a modified version of this software (i) may not refer to the modified
29	* version by the same designation, or by any confusingly similar designation, and (ii) must refer to
30	* the underlying software originally provided by Alliance as "SolTrace". Except to comply with the
31	* foregoing, the term "SolTrace", or any confusingly similar designation may not be used to refer to
32	* any modified version of this software or any modified version of the underlying software originally
33	* provided by Alliance without the prior written consent of Alliance.
34	*
35	* 5. The name of the copyright holder, contributors, the United States Government, the United States
36	* Department of Energy, or any of their employees may not be used to endorse or promote products
37	* derived from this software without specific prior written permission.
38	*
39	* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
40	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
41	* FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER,
42	* CONTRIBUTORS, UNITED STATES GOVERNMENT OR UNITED STATES DEPARTMENT OF ENERGY, NOR ANY OF THEIR
43	* EMPLOYEES, BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
44	* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
45	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
46	* IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
47	* THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
48	*******************************************************************************************************/
49
50
51	#include <math.h>
52
53	#include "types.h"
54	#include "procs.h"
55
56	inline double sqr(double x) { return (x)*(x); }	×
57
58	void Interaction(	114,857✔
59	MTRand &myrng,
60	double PosXYZ[3],
61	double CosKLM[3],
62	double DFXYZ[3],
63	int InteractionType,
64	TOpticalProperties *Opticl,
65	double Wavelength,
66	double PosOut[3],
67	double CosOut[3],
68	int *ErrorFlag )
69	{
70	/* {Purpose: To compute the direction cosines of the ray due to optical interaction
71	at the intersection point of the ray with the surface
72	Input - PosXYZ[2] = X,Y,Z coordinates of intersection point.
73	DFXYZ = direction numbers for the surface normal at the
74	intersection point (partial derivatives with respect
75	to X,Y,Z of surface equation)
76	InteractionType = Optical interaction type indicator
77	= 1, refraction
78	= 2, reflection
79	= 3, aperture stop
80	= 4, diffraction, transmission grating
81	= 5, diffraction, reflection grating
82	CosKLM[2] = direction cosines of incident ray
83	Opticl = record of optical properties
84	.RefractiveIndex[4] = Refractive index of incident and outgoing medium
85	[0] = real part of incident medium refractive index
86	[1] = imaginary part of ""
87	[2] = real part of outgoing medium refractive index
88	[4] = imaginary part of ""
89	.ApertureStopOrGratingType
90	for InteractionType = 3, aperture stop
91	= 1, slit
92	= 2, elliptical
93	for InteractionType = 4,5 grating
94	= 1, planes parallel to Y-Z plane
95	= 2, concentric cylinders centered about Z-axis
96	.DiffractionOrder = integral order of diffraction for InteractionTypes=4,5, grating
97	.AB12[4] = coefficients of polynomial specifying grating spacing for InteractionTypes=4,5
98	[0] = lower X limit, ApertureStopOrGratingType = 1
99	semi-X axis, ApertureStopOrGratingType = 2
100	[1] = lower Y limit, ApertureStopOrGratingType = 1
101	semi-Y axis, ApertureStopOrGratingType = 2
102	[2] = upper X limit, ApertureStopOrGratingType = 1
103	unused, ApertureStopOrGratingType = 2
104	[4] = upper Y limit, ApertureStopOrGratingType = 1
105	unused, ApertureStopOrGratingType = 2
106	Wavelength = wavelength of ray
107
108	Output - PosOut[2] = position of ray after optical interaction
109	CosOut[2] = direction cosines of ray after optical interaction
110	ErrorFlag = Error flag indicating successful interaction
111	}*/
112
113	int i = 0;	114,857✔
114	double CosUVW[3] = {0.0, 0.0, 0.0};	114,857✔
115	int NIter = 0, IType = 0, NMord = 0;	114,857✔
116	double Epsilon = 0.0, Refr1 = 0.0, Refr2 = 0.0, RMU = 0.0, RM2 = 0.0, D2 = 0.0, B = 0.0, A = 0.0, A2 = 0.0;	114,857✔
117	double Gamn = 0.0, Gamn1 = 0.0;	114,857✔
118	double X = 0.0,Y = 0.0,A1 = 0.0,B1 = 0.0,Ellips = 0.0,B2 = 0.0;	114,857✔
119	double RK = 0.0,RL = 0.0,RM = 0.0,Denom,U=0,V=0,W=0;	114,857✔
120	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;	114,857✔
121	double Rave = 0.0, Rs = 0.0, Rp = 0.0, UnitDFXYZ[3] = {0.0,0.0,0.0}, IncidentAngle = 0.0;	114,857✔
122
123	NIter = 10;	114,857✔
124	Epsilon = 0.000005;	114,857✔
125
126	*ErrorFlag = 0;	114,857✔
127	for (i=0;i<3;i++)	459,428✔
128	PosOut[i] = PosXYZ[i];	344,571✔
129
130	switch (InteractionType)	114,857✔
131	{
132
133	/*{ InteractionType = 1, Refraction
134	===============================================================================}*/
135	case 1:	×
136	{
NEW 137	Refr1 = Opticl->RefractiveIndex[0];	×
NEW 138	Refr2 = Opticl->RefractiveIndex[2];	×
NEW 139	RMU = Refr1 / Refr2;	×
NEW 140	D2 = DOT(DFXYZ, DFXYZ);	×
NEW 141	B = (RMU * RMU - 1.0) / D2;	×
NEW 142	A = RMU * DOT(CosKLM, DFXYZ) / D2;	×
NEW 143	A2 = A * A;	×
NEW 144	if (B > A2) //Total internal reflection	×
145	{
NEW 146	A = DOT(CosKLM, DFXYZ) / DOT(DFXYZ, DFXYZ);	×
NEW 147	for (i = 0; i < 3; i++)	×
NEW 148	CosOut[i] = CosKLM[i] - 2.0 * A * DFXYZ[i];	×
149	return;	114,857✔
150	}
151
152	//fresnel equations
NEW 153	UnitDFXYZ[0] = -DFXYZ[0] / sqrt(DOT(DFXYZ, DFXYZ)); //unit surface normals	×
NEW 154	UnitDFXYZ[1] = -DFXYZ[1] / sqrt(DOT(DFXYZ, DFXYZ));	×
NEW 155	UnitDFXYZ[2] = -DFXYZ[2] / sqrt(DOT(DFXYZ, DFXYZ));	×
NEW 156	IncidentAngle = acos(DOT(CosKLM, UnitDFXYZ));	×
NEW 157	Rs = sqr(((Refr1 * cos(IncidentAngle) - Refr2 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2)))) /	×
NEW 158	((Refr1 * cos(IncidentAngle) + Refr2 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2)))));	×
NEW 159	Rp = sqr(((Refr1 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2))) - Refr2 * cos(IncidentAngle)) /	×
NEW 160	((Refr1 * sqrt(1 - sqr(Refr1 * sin(IncidentAngle) / Refr2))) + Refr2 * cos(IncidentAngle)));	×
NEW 161	Rave = (Rp + Rs) / 2.0; //average of s and p polarized light; equal parts of both = non-polarized	×
NEW 162	if (Rave < myrng()) //transmitted through surface	×
163	{
NEW 164	Gamn = -B / (2.0 * A);	×
165
166	//Begin Newton-Raphson loop to converge on correct root.
NEW 167	bool converged = false;	×
NEW 168	for (i = 1; i < NIter; i++)	×
169	{
NEW 170	Gamn1 = (Gamn * Gamn - B) / (2.0 * (Gamn + A));	×
NEW 171	if (fabs(Gamn - Gamn1) < Epsilon)	×
172	{
NEW 173	converged = true;	×
NEW 174	break;	×
175	//goto Label_Converge;
176	}
177
NEW 178	Gamn = Gamn1;	×
179	}
180	//Failed to converge
NEW 181	if (converged == false)	×
182	{
NEW 183	*ErrorFlag = 12;	×
NEW 184	return;	×
185	}
186
187	//Have converged on Gamma, Compute direction cosines of refracted ray.
188	//Label_Converge:
NEW 189	for (i = 0; i < 3; i++)	×
NEW 190	CosOut[i] = RMU * CosKLM[i] + Gamn1 * DFXYZ[i];	×
191	}
192	else //reflected from surface
193	{
NEW 194	A = DOT(CosKLM, DFXYZ) / DOT(DFXYZ, DFXYZ);	×
NEW 195	for (i = 0; i < 3; i++)	×
NEW 196	CosOut[i] = CosKLM[i] - 2.0 * A * DFXYZ[i];	×
197	}
NEW 198	return;	×
199	break;
200	}
201
202	/*{ InteractionType = 2, Reflection
203	===============================================================================}*/
204	case 2:	114,857✔
205	{
206	A = DOT(CosKLM, DFXYZ) / DOT(DFXYZ, DFXYZ);	114,857✔
207	//Compute direction cosines for reflected ray
208	for (i = 0; i < 3; i++)	459,428✔
209	CosOut[i] = CosKLM[i] - 2.0 * A * DFXYZ[i];	344,571✔
210
211	return;	114,857✔
212	break;
213	}
214
215	/*{ InteractionType = 3, Aperture Stop
216	===============================================================================}*/
217	case 3:	×
218	{
NEW 219	X = PosXYZ[0];	×
NEW 220	Y = PosXYZ[1];	×
NEW 221	IType = Opticl->ApertureStopOrGratingType;	×
NEW 222	A1 = Opticl->AB12[0];	×
NEW 223	B1 = Opticl->AB12[1];	×
224
NEW 225	bool ray_missed_aperture = false;	×
NEW 226	if (IType == 1) //Slit Aperture	×
227	{
NEW 228	A2 = Opticl->AB12[2];	×
NEW 229	B2 = Opticl->AB12[3];	×
NEW 230	if (X < A1 \|\| X > A2)	×
231	{
232	*ErrorFlag = 31;	×
NEW 233	ray_missed_aperture = true;	×
234	//goto RayMissesAperture;
235	}
236	else
237	{
NEW 238	if (Y >= B1 && Y <= B2) return;	×
239
NEW 240	*ErrorFlag = 31;	×
NEW 241	ray_missed_aperture = true;	×
242	//goto RayMissesAperture;
243	}
244
245	}
246
NEW 247	else if (IType == 2) //Elliptical Aperture	×
248	{
NEW 249	Ellips = X * X / (A1 * A1) + Y * Y / (B1 * B1);	×
NEW 250	if (Ellips <= 1.0) return;	×
NEW 251	*ErrorFlag = 32;	×
252	}
253
254	//RayMissesAperture:
255	//Ray misses aperture
NEW 256	if (ray_missed_aperture == true)	×
257	{
NEW 258	for (i = 0; i < 3; i++)	×
UNCOV 259	CosOut[i] = 0.0;	×
260	}
NEW 261	return;	×
262
263	break;
264	}
265
266	/*{ InteractionType = 4,5; Diffraction
267	===============================================================================}*/
268	case 4:	×
269	case 5:
270	{
NEW 271	IType = Opticl->ApertureStopOrGratingType;	×
NEW 272	NMord = Opticl->DiffractionOrder;	×
NEW 273	Refr1 = Opticl->RefractiveIndex[0];	×
NEW 274	Refr2 = Opticl->RefractiveIndex[2];	×
NEW 275	RMU = Refr1 / Refr2;	×
NEW 276	D2 = DOT(DFXYZ, DFXYZ);	×
NEW 277	RK = DFXYZ[0];	×
NEW 278	RL = DFXYZ[1];	×
NEW 279	RM = DFXYZ[2];	×
NEW 280	X = PosXYZ[0];	×
NEW 281	Y = PosXYZ[1];	×
282
NEW 283	if (IType == 1) //Parallel plane grating	×
284	{
NEW 285	Denom = RL * RL + RM * RM;	×
NEW 286	U = 1.0 / sqrt(1.0 + RK * RK / Denom);	×
NEW 287	V = -RK * RL * U / Denom;	×
NEW 288	W = -RK * RM * U / Denom;	×
NEW 289	Varr = X;	×
NEW 290	GFactr = 1.0 / U;	×
291	//goto CompDiffInt;
292	}
293
NEW 294	else if (IType == 2) //Concentric Cylinder Grating	×
295	{
NEW 296	Rho2 = X * X + Y * Y;	×
NEW 297	Rho = sqrt(Rho2);	×
NEW 298	RM2 = RM * RM;	×
NEW 299	Term = RL * X - RK * Y;	×
NEW 300	G = sqrt(D2 * (RM2 * Rho2 + Term * Term));	×
NEW 301	U = (RM2 * X + RL * Term) / G;	×
NEW 302	V = (RM2 * Y - RK * Term) / G;	×
NEW 303	W = -RM * (RK * X + RL * Y) / G;	×
NEW 304	Varr = Rho;	×
NEW 305	GFactr = Rho / (X * U + Y * V);	×
306	}
307	//CompDiffInt: //Compute interaction due to diffraction
NEW 308	CosUVW[0] = U;	×
NEW 309	CosUVW[1] = V;	×
NEW 310	CosUVW[2] = W;	×
311
NEW 312	D = 0.0;	×
NEW 313	XX = 1.0;	×
314
NEW 315	for (i = 0; i < 4; i++)	×
316	{
NEW 317	D = D + Opticl->AB12[i] * XX;	×
NEW 318	XX = XX * Varr;	×
319	}
320
NEW 321	D = D * GFactr;	×
NEW 322	Ordiff = NMord;	×
NEW 323	RLamda = Ordiff * Wavelength / (Refr2 * D);	×
NEW 324	B = (RMU * RMU - 1.0 + RLamda * RLamda - 2.0 * RMU * RLamda * DOT(CosKLM, CosUVW)) / D2;	×
NEW 325	A = RMU * DOT(CosKLM, DFXYZ) / D2;	×
NEW 326	A2 = A * A;	×
NEW 327	if (B > A2) //Total internal reflection	×
328	{
NEW 329	for (i = 0; i < 3; i++)	×
NEW 330	CosOut[i] = 0.0;	×
NEW 331	*ErrorFlag = 11;	×
NEW 332	return;	×
333	}
334
NEW 335	Gamn = -B / (2.0 * A);	×
NEW 336	if (InteractionType == 5)	×
NEW 337	Gamn = -Gamn - 2.0 * A;	×
338
339	//Begin Newton-Raphson loop to converge on correct root.
NEW 340	i = 0;	×
NEW 341	bool converged = false;	×
NEW 342	while (i++ < NIter)	×
343	{
NEW 344	Gamn1 = (Gamn * Gamn - B) / (2.0 * (Gamn + A));	×
NEW 345	if (fabs(Gamn - Gamn1) < Epsilon)	×
346	{
NEW 347	converged = true;	×
NEW 348	break;	×
349	//goto CompDCos;
350	}
NEW 351	Gamn = Gamn1;	×
352	}
353	//Failed to converge
NEW 354	if (converged == false)	×
355	{
356	*ErrorFlag = 12;	×
357	return;	×
358	}
359	//Have converged on Gamn1. Compute direction cosines of diffracted ray.
360	//CompDCos:
NEW 361	for (i = 0; i < 3; i++)	×
NEW 362	CosOut[i] = RMU * CosKLM[i] - RLamda * CosUVW[i] + Gamn1 * DFXYZ[i];	×
363
364	break;	×
365	}
366	default:	×
367	break;	×
368	}
369	}
370	//End of Procedure--------------------------------------------------------------

NREL / SolTrace / 14628243107

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