Designing Optical Systems

Manipulate[Deploy[LensLayout[{wavelength1,wavelength2,sourceAngle,sourceXY1,sourceXY2,lensAngle,thicknessAndAperture1,thicknessAndAperture2,curvature1,curvature2,lensLocation1,lensLocation2}]],
{{wavelength1,0.532,"(Subscript[λ, 1])"},.2,1,ControlType->Manipulator},{{wavelength2,0.532,"(Subscript[λ, 2])"},.2,1,ControlType->Manipulator},{{sourceAngle,{0.,0.}},{-Pi,-Pi},{Pi,Pi}},{{sourceXY1,{0.,0.}},{0,-80},{200,80}},{{sourceXY2,{0.,10.}},{0,-80},{200,80}},{{lensAngle,{0.,0.}},{-Pi,-Pi},{Pi,Pi}},{{thicknessAndAperture1,{10.,50.}},{0,1},{100,100}},{{thicknessAndAperture2,{10.,50.}},{0,1},{100,100}},{{curvature1,{0.01,-0.01}},{-0.03,-0.03},{0.03,0.03}},{{curvature2,{0.01,-0.01}},{-0.03,-0.03},{0.03,0.03}},{{lensLocation1,{50.,0.}},{0,-80},{200,80}},{{lensLocation2,{75.,0.}},{0,-80},{200,80}},ControlType->None,
AutorunSequencing->{1,2,9,10},
SynchronousInitialization->False,
Initialization:>(
micronsToColor[m_]:=RGBColor[Abs[-0.3183098861837907*ArcTan[63.75 - 0.25*Abs[-500. + 2.*Mod[100. + 500.*m, 500.]]] + 
    0.3183098861837907*ArcTan[43. - 0.1*Abs[-500. + 2.*Mod[100. + 500.*m, 500.]]]], 
  Abs[0.3183098861837907*ArcTan[60. - 0.2*Abs[-500. + 2.*Mod[100. + 500.*m, 500.]]] - 
    0.3183098861837907*ArcTan[18. - 0.1*Abs[-500. + 2.*Mod[100. + 500.*m, 500.]]]], 
  Abs[-0.3183098861837907*ArcTan[10. - 0.2*Abs[-500. + 2.*Mod[100. + 500.*m, 500.]]] + 
    0.3183098861837907*ArcTan[22.5 - 0.125*Abs[-500. + 2.*Mod[100. + 500.*m, 500.]]]]];
(*rayfunction-wedge*)
rayfunction = Compile[{as1, wl1, xs1, ys1, aa, ab, apa, apb, ca1, ca2, cb1, cb2, ta, tb, xa1, xb1, ya1, yb1}, 
 Module[{activeregion, allcreatedamplitudes, allcreatedrays, 
   allowzeroraylength, alternateindex, amplitude, axismagnitude, 
   axismatrix, clonedata, clonegenerationnumber, cloneinput, 
   clonenumber, cloneposition, clonereportedrays, clonerotation, 
   clonescale, clonescalevector, clonesorthistogram, clonetraceinput, 
   clonetraceoutput, closestdistances, closestposition, closestpositions, 
   closestraylength, completedray, componentnumber, cosangle, 
   crossmatrix, deflecteddummyrayspresent, deflectedlocalrays, deflectedrays, 
   deflectedrayspresent, dummymatrix, dummymatrixlist, dummyvector, 
   frontid, generationlimit, generationnumber, holeintersections, 
   holeQ, holeQlist, hor, index, index1, 
   index2, inrays, intensity, intensityscale, 
   intersectionnumber, intersections, intersectionslist, 
   intersectionspresent, localdummyoutrayspresent, localoutrays, 
   localoutrayspresent, magnitude, maxclonegenerationnumber, maximumposition, 
   minimumraylength, minresonateraylength, n12, newcomponentnumber, 
   newdeflectedrays, newdistance, newsurfaceid, newsurfacenumber, 
   nonintersectingrays, normalizedpolarization, oldgenerationnumber, 
   opticallength, originalray, originalraytilt, originalrotation, 
   originalsurfnorm, outrays, parallelamplitude, PAxis, 
   perpendicularamplitude, phaseshift, pol1, pol2, 
   pol3, previoustilt, propagate, propagatedrays, r1, 
   r2, r3, r4, r5, r6, r7, 
   r8, r9, ray1, ray2, ray3, RayL, 
   rayslot, raysreported, raytiltx, raytiltz, 
   reporteddummyrayspresent, reportedlocalrays, reportedrays, 
   reportedrayspresent, resonateintersections, resonatenumber, 
   resonateraylength, rex, rey, rez, rotationaxis, 
   rsx, RSx, rsy, RSy, rsz, RSz, RTx, RTy, RTz, rx, 
   rxhlp, ry, ryhlp, rz, rzhlp, SCx, SCy, 
   sid, sinangle, sno1, sno2, sno3, 
   sno4, surfaceid, surfacenumber, surfacerotationmatrix, 
   surfacetranslationvector, surfnorm, symr1, symr2, 
   symr3, symr4, symr5, symr6, symr7, 
   symr8, symr9, tempclonedata, terminatedeflection, 
   thresholdintensity, tiltmagnitude, transmatrix, TV1, TV10, TV2, TV3, TV4, TV5, TV6, TV7, TV8, 
   TV9, u, unfilteredintersections, unit, uu, ver, 
   wavelength, XAxis, xunit, YAxis, yunit, 
   ZAxis, zunit}, activeregion = 1.; 
   allcreatedamplitudes = {{0. + 0.*I, 0. + 0.*I, 0. + 0.*I}}; allcreatedrays = {{1.}}; 
   allowzeroraylength = True; alternateindex = 1.; amplitude = {1. + 0.*I, 0. + 0.*I, 0. + 0.*I}; 
   axismagnitude = 1.; axismatrix = {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}; clonedata = {{1.}}; 
   clonegenerationnumber = 1.; cloneinput = {1.}; clonenumber = 1; 
   cloneposition = {0., 0., 0.}; clonereportedrays = {{1.}}; clonerotation = {{1.}}; 
   clonescale = 1.; clonescalevector = {1.}; clonesorthistogram = {1}; 
   clonetraceinput = {{1.}}; clonetraceoutput = {{1.}}; closestdistances = {1.}; 
   closestposition = 0; closestpositions = {0}; closestraylength = 1.; 
   completedray = {1.}; componentnumber = -1.; cosangle = 1.; 
   crossmatrix = {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}; deflecteddummyrayspresent = False; 
   deflectedlocalrays = {{1.}}; deflectedrays = {{1.}}; deflectedrayspresent = False; 
   dummymatrix = {{1.}}; dummymatrixlist = {{{1.}}}; dummyvector = {1.}; 
   frontid = 1.; generationlimit = 200.; generationnumber = 1.; 
   holeintersections = {{-1., 0., 0.}}; holeQ = False; holeQlist = {0}; 
   hor = 1. + 0.*I; index = 1.; index1 = 1.; index2 = 1.; 
   inrays = {{1.}}; intensity = 1.; intensityscale = 1.; 
   intersectionnumber = 1.; intersections = {{-1., 0., 0.}}; intersectionslist = {{{-1., 0., 0.}}}; 
   intersectionspresent = False; localdummyoutrayspresent = False; localoutrays = {{1.}}; 
   localoutrayspresent = False; magnitude = 1.; maxclonegenerationnumber = 1.; 
   maximumposition = 1; minimumraylength = 0.0001; minresonateraylength = 1.; 
   n12 = 1.; newcomponentnumber = -1.; newdeflectedrays = {{1.}}; 
   newdistance = 1.; newsurfaceid = -1.; newsurfacenumber = -1.; 
   nonintersectingrays = {{1.}}; normalizedpolarization = {1., 1.}; oldgenerationnumber = 1.; 
   opticallength = 1.; originalray = {1.}; originalraytilt = {1.}; 
   originalrotation = {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}; originalsurfnorm = {{1.}}; 
   outrays = {{1.}}; parallelamplitude = 1. + 0.*I; PAxis = {1., 0., 0.}; 
   perpendicularamplitude = 1. + 0.*I; phaseshift = 0.; pol1 = 1.; pol2 = 1.; 
   pol3 = 1.; previoustilt = {1., 0., 0.}; propagate = True; 
   propagatedrays = {{1.}}; r1 = 1.; r2 = 1.; r3 = 1.; 
   r4 = 1.; r5 = 1.; r6 = 1.; r7 = 1.; r8 = 1.; 
   r9 = 1.; ray1 = 1.; ray2 = 1.; ray3 = 1.; RayL = 1.; 
   rayslot = 1.; raysreported = False; raytiltx = 1.; raytiltz = 1.; 
   reporteddummyrayspresent = False; reportedlocalrays = {{1.}}; reportedrays = {{1.}}; 
   reportedrayspresent = False; resonateintersections = {{-1., 0., 0.}}; resonatenumber = 0; 
   resonateraylength = 1.; rex = 1.; rey = 1.; rez = 1.; 
   rsx = 1.; RSx = 1.; rsy = 1.; RSy = 1.; rsz = 1.; RSz = 1.; RTx = 1.; RTy = 1.; RTz = 1.; 
   rx = 1.; rxhlp = 1.; ry = 1.; ryhlp = 1.; rz = 1.; 
   rzhlp = 1.; SCx = 1.; SCy = 1.; sid = 1.; sinangle = 1.; sno1 = 1.; 
   sno2 = 1.; sno3 = 1.; sno4 = 1.; surfaceid = -1.; 
   surfacenumber = -1.; surfacerotationmatrix = {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}; 
   surfacetranslationvector = {0., 0., 0.}; surfnorm = {{1.}}; symr1 = 1.; 
   symr2 = 1.; symr3 = 1.; symr4 = 1.; symr5 = 1.; 
   symr6 = 1.; symr7 = 1.; symr8 = 1.; symr9 = 1.; 
   tempclonedata = {{1.}}; terminatedeflection = False; thresholdintensity = 1.; 
   tiltmagnitude = 1.; transmatrix = {{{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}}; TV1 = 1.; TV10 = 1.; TV2 = 1.; 
   TV3 = 1.; TV4 = 1.; TV5 = 1.; TV6 = 1.; TV7 = 1.; TV8 = 1.; TV9 = 1.; u = {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}; 
   unfilteredintersections = {{-1., 0., 0.}}; unit = {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}; 
   uu = {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}; ver = 0. + 0.*I; wavelength = 1.; 
   XAxis = {1., 0., 0.}; xunit = 1.; YAxis = {1., 0., 0.}; yunit = 1.; 
   ZAxis = {1., 0., 0.}; zunit = 1.; propagatedrays = 
    Join[{symr1, symr2, symr3, symr4, symr5, 
        symr6, symr7, symr8, symr9} = 
       {Cos[as1*Degree], Sin[as1*Degree], 0, -Sin[as1*Degree], Cos[as1*Degree], 0, 0, 0, 1}; 
      (({rsx, rsy, rsz, rex, rey, 
           rez, r1, r2, r3, r4, 
           r5, r6, r7, r8, r9, 
           wavelength, intensity, intensityscale, ray1, 
           index, opticallength, phaseshift, componentnumber, 
           surfacenumber, surfaceid, ray2, ray3, sid, 
           frontid, sno1, sno2, sno3, sno4, 
           rayslot, oldgenerationnumber, intersectionnumber, pol1, 
           pol2, pol3} = #1; wavelength = wl1; 
         {TV1, TV2, TV3} = {rsx*symr1 + rsy*symr4 + 
            rsz*symr7 + xs1, rsx*symr2 + 
            rsy*symr5 + rsz*symr8 + ys1, 
           rsx*symr3 + rsy*symr6 + 
            rsz*symr9}; {rsx, rsy, rsz} = 
          {TV1, TV2, TV3}; {TV1, TV2, TV3} = {rex*symr1 + rey*symr4 + 
            rez*symr7 + xs1, rex*symr2 + 
            rey*symr5 + rez*symr8 + ys1, 
           rex*symr3 + rey*symr6 + 
            rez*symr9}; {rex, rey, rez} = 
          {TV1, TV2, TV3}; {TV1, TV2, TV3} = {r1*symr1 + r2*symr4 + 
            r3*symr7, r1*symr2 + r2*
             symr5 + r3*symr8, r1*symr3 + 
            r2*symr6 + r3*symr9}; 
         {TV4, TV5, TV6} = {r4*symr1 + r5*symr4 + 
            r6*symr7, r4*symr2 + r5*
             symr5 + r6*symr8, r4*symr3 + 
            r5*symr6 + r6*symr9}; 
         {TV7, TV8, TV9} = {r7*symr1 + r8*symr4 + 
            r9*symr7, r7*symr2 + r8*
             symr5 + r9*symr8, r7*symr3 + 
            r8*symr6 + r9*symr9}; 
         {{r1, r2, r3}, {r4, r5, 
            r6}, {r7, r8, r9}} = 
          {{TV1, TV2, TV3}, {TV4, TV5, TV6}, {TV7, TV8, TV9}}; {rsx, rsy, rsz, 
          rex, rey, rez, r1, r2, 
          r3, r4, r5, r6, r7, r8, 
          r9, wavelength, intensity, intensityscale, 
          ray1, index, opticallength, phaseshift, 
          componentnumber, surfacenumber, surfaceid, ray2, 
          ray3, sid, frontid, sno1, sno2, 
          sno3, sno4, rayslot, generationnumber, 
          intersectionnumber, pol1, pol2, pol3}) & ) /@ 
       {{0., 0., 0., 0., 0., 0., 0.9659258262890683, 0.25881904510252074, 0., -0.25881904510252074, 0.9659258262890683, 0., 0., 0., 1., 0.532, 100., 
         0.09090909090909091, 0., 1.0002694514570802, 0., 0., 1., 0., 0., 0., 0., 16497., 16497., 1., 11., 1., 0., 1., 0., 0., 0., 1., 0.}, 
        {0., 0., 0., 0., 0., 0., 0.9781476007338057, 0.20791169081775934, 0., -0.20791169081775934, 0.9781476007338057, 0., 0., 0., 1., 0.532, 100., 
         0.09090909090909091, 0., 1.0002694514570802, 0., 0., 1., 0., 0., 0., 0., 16497., 16497., 2., 11., 2., 0., 1., 0., 0., 0., 1., 0.}, 
        {0., 0., 0., 0., 0., 0., 0.9876883405951378, 0.15643446504023087, 0., -0.15643446504023087, 0.9876883405951378, 0., 0., 0., 1., 0.532, 100., 
         0.09090909090909091, 0., 1.0002694514570802, 0., 0., 1., 0., 0., 0., 0., 16497., 16497., 3., 11., 3., 0., 1., 0., 0., 0., 1., 0.}, 
        {0., 0., 0., 0., 0., 0., 0.9945218953682733, 0.10452846326765348, 0., -0.10452846326765348, 0.9945218953682733, 0., 0., 0., 1., 0.532, 100., 
         0.09090909090909091, 0., 1.0002694514570802, 0., 0., 1., 0., 0., 0., 0., 16497., 16497., 4., 11., 4., 0., 1., 0., 0., 0., 1., 0.}, 
        {0., 0., 0., 0., 0., 0., 0.9986295347545738, 0.052335956242943835, 0., -0.052335956242943835, 0.9986295347545738, 0., 0., 0., 1., 0.532, 100., 
         0.09090909090909091, 0., 1.0002694514570802, 0., 0., 1., 0., 0., 0., 0., 16497., 16497., 5., 11., 5., 0., 1., 0., 0., 0., 1., 0.}, 
        {0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1., 0.532, 100., 0.09090909090909091, 0., 1.0002694514570802, 0., 0., 1., 0., 0., 0., 0., 16497., 
         16497., 6., 11., 6., 0., 1., 0., 0., 0., 1., 0.}, {0., 0., 0., 0., 0., 0., 0.9986295347545738, -0.052335956242943835, 0., 0.052335956242943835, 
         0.9986295347545738, 0., 0., 0., 1., 0.532, 100., 0.09090909090909091, 0., 1.0002694514570802, 0., 0., 1., 0., 0., 0., 0., 16497., 16497., 7., 11., 7., 0., 
         1., 0., 0., 0., 1., 0.}, {0., 0., 0., 0., 0., 0., 0.9945218953682733, -0.10452846326765348, 0., 0.10452846326765348, 0.9945218953682733, 0., 0., 0., 1., 
         0.532, 100., 0.09090909090909091, 0., 1.0002694514570802, 0., 0., 1., 0., 0., 0., 0., 16497., 16497., 8., 11., 8., 0., 1., 0., 0., 0., 1., 0.}, 
        {0., 0., 0., 0., 0., 0., 0.9876883405951378, -0.15643446504023087, 0., 0.15643446504023087, 0.9876883405951378, 0., 0., 0., 1., 0.532, 100., 
         0.09090909090909091, 0., 1.0002694514570802, 0., 0., 1., 0., 0., 0., 0., 16497., 16497., 9., 11., 9., 0., 1., 0., 0., 0., 1., 0.}, 
        {0., 0., 0., 0., 0., 0., 0.9781476007338057, -0.20791169081775934, 0., 0.20791169081775934, 0.9781476007338057, 0., 0., 0., 1., 0.532, 100., 
         0.09090909090909091, 0., 1.0002694514570802, 0., 0., 1., 0., 0., 0., 0., 16497., 16497., 10., 11., 10., 0., 1., 0., 0., 0., 1., 0.}, 
        {0., 0., 0., 0., 0., 0., 0.9659258262890683, -0.25881904510252074, 0., 0.25881904510252074, 0.9659258262890683, 0., 0., 0., 1., 0.532, 100., 
         0.09090909090909091, 0., 1.0002694514570802, 0., 0., 1., 0., 0., 0., 0., 16497., 16497., 11., 11., 11., 0., 1., 0., 0., 0., 1., 0.}}]; 
   propagatedrays = (({rsx, rsy, rsz, rex, 
         rey, rez, r1, r2, r3, 
         r4, r5, r6, r7, r8, r9, 
         wavelength, intensity, intensityscale, ray1, 
         index, opticallength, phaseshift, componentnumber, 
         surfacenumber, surfaceid, ray2, ray3, sid, 
         frontid, sno1, sno2, sno3, sno4, 
         rayslot, oldgenerationnumber, intersectionnumber, pol1, 
         pol2, pol3} = #1; componentnumber = 1.; frontid = 1.; 
       {rsx, rsy, rsz, rex, rey, 
        rez, r1, r2, r3, r4, r5, 
        r6, r7, r8, r9, wavelength, 
        intensity, intensityscale, ray1, index, 
        opticallength, phaseshift, componentnumber, surfacenumber, 
        surfaceid, ray2, ray3, sid, frontid, 
        sno1, sno2, sno3, sno4, rayslot, 
        generationnumber, intersectionnumber, pol1, pol2, 
        pol3}) & ) /@ propagatedrays; thresholdintensity = 1.; 
   generationlimit = 200.; reportedrays = 
    {{rsx, rsy, rsz, rex, rey, 
      rez, r1, r2, r3, r4, r5, 
      r6, r7, r8, r9, wavelength, 
      intensity, intensityscale, ray1, index, 
      opticallength, phaseshift, componentnumber, surfacenumber, 
      surfaceid, ray2, ray3, sid, frontid, 
      sno1, sno2, sno3, sno4, rayslot, 
      generationnumber, intersectionnumber, pol1, pol2, 
      pol3}}; generationnumber = 0.; (inrays = propagatedrays; 
    outrays = {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}; deflectedrayspresent = True; 
    reportedrayspresent = False; propagate = True; 
    While[deflectedrayspresent && propagate, generationnumber += 1.; 
      (deflecteddummyrayspresent = False; deflectedrayspresent = False; 
       reportedrayspresent = False; reporteddummyrayspresent = True; localdummyoutrayspresent = 
        False; localoutrayspresent = False; {deflectedlocalrays, localoutrays, 
         reportedlocalrays} = (Flatten[#1, 1] & ) /@ 
         Transpose[(({rsx, rsy, rsz, rex, rey, 
               rez, r1, r2, r3, r4, 
               r5, r6, r7, r8, r9, 
               wavelength, intensity, intensityscale, ray1, 
               index, opticallength, phaseshift, componentnumber, 
               surfacenumber, surfaceid, ray2, ray3, 
               sid, frontid, sno1, sno2, sno3, 
               sno4, rayslot, oldgenerationnumber, intersectionnumber, 
               pol1, pol2, pol3} = #1; If[intensity > 
                thresholdintensity && generationnumber < generationlimit, 
              (resonatenumber = 0; resonateraylength = 1.7976931348623157*^308; 
                resonateintersections = {{-1., 0., 0.}}; MapIndexed[If[#1[[1,1]] > -1. && #1[[1,1]] < resonateraylength, 
                   resonateraylength = #1[[1,1]]; resonateintersections = #1; resonatenumber = 
                     First[#2], resonatenumber] & , intersectionslist = {{{-1., 0., 0.}}}; 
                  holeQlist = {0}; ({rsx, rsy, rsz, 
                     rex, rey, rez, r1, r2, 
                     r3, r4, r5, r6, r7, 
                     r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     oldgenerationnumber, intersectionnumber, pol1, pol2, 
                     pol3} = #1; index1 = index; xunit = 
                    r1*Re[Cos[0.017453292519943295*aa]] + r2*Re[Sin[0.017453292519943295*aa]]; 
                   yunit = r2*Re[Cos[0.017453292519943295*aa]] - 1.*r1*
                      Re[Sin[0.017453292519943295*aa]]; zunit = r3*Re[Cos[0.017453292519943295*aa]^2 + 
                       Sin[0.017453292519943295*aa]^2]; xunit = Chop[xunit]; yunit = 
                    Chop[yunit]; zunit = Chop[zunit]; magnitude = 
                    (xunit^2 + yunit^2 + zunit^2)^2^(-1); xunit = 
                    xunit/magnitude; yunit = yunit/
                     magnitude; zunit = zunit/magnitude; 
                   rx = (rex - 1.*Re[xa1])*Re[Cos[0.017453292519943295*aa]] + (rey - 1.*Re[ya1])*
                      Re[Sin[0.017453292519943295*aa]]; ry = (rey - 1.*Re[ya1])*Re[Cos[0.017453292519943295*aa]] - 
                     1.*(rex - 1.*Re[xa1])*Re[Sin[0.017453292519943295*aa]]; rz = rez*
                     Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]; rx = Chop[rx]; 
                   ry = Chop[ry]; rz = Chop[rz]; 
                   terminatedeflection = False; intersections = {{-1., 0., 0.}}; intersections = 
                    Block[{TempVar, TempVar1, TempVar2, TempVar3, RayL, SCx, SCy}, TempVar1 = 1.; TempVar2 = 1.; TempVar3 = 1.; TempVar = {{-1., 0., 0.}}; RayL = 1.; 
                      SCx = 1.; SCy = 1.; TempVar1 = rx^2 + ry^2 + rz^2; TempVar2 = 
                       rx*xunit + ry*yunit + rz*
                         zunit - xunit/ca1; TempVar3 = -TempVar1 + TempVar2^2 + (2*rx)/ca1; 
                      TempVar3 = Re[TempVar3]; If[N[Chop[TempVar3]] >= 0., TempVar3 = Sqrt[TempVar3]; TempVar = {{RayL, SCx, SCy} = {-TempVar2 - TempVar3, 
                             ry - (TempVar2 + TempVar3)*yunit, rz - (TempVar2 + TempVar3)*
                               zunit}; If[N[Chop[RayL]] >= 0. && SCx^2 + SCy^2 <= 0.25*Abs[apa]^2 && Abs[ry + 
                                 RayL*yunit - SCx] + Abs[rz + RayL*zunit - SCy] + 
                               Abs[-ca1^(-1) + rx + RayL*xunit + Sqrt[ca1^(-2) - SCx^2 - SCy^2]/Sign[ca1]] < 1.*^-7, 
                            {RayL, SCx, SCy}, {-1., 0., 0.}], {RayL, SCx, SCy} = {-TempVar2 + TempVar3, ry + (-TempVar2 + TempVar3)*
                               yunit, rz + (-TempVar2 + TempVar3)*zunit}; If[N[Chop[RayL]] >= 0. && 
                             SCx^2 + SCy^2 <= 0.25*Abs[apa]^2 && Abs[ry + RayL*yunit - SCx] + 
                               Abs[rz + RayL*zunit - SCy] + Abs[-ca1^(-1) + rx + 
                                 RayL*xunit + Sqrt[ca1^(-2) - SCx^2 - SCy^2]/Sign[ca1]] < 1.*^-7, {RayL, SCx, SCy}, {-1., 0., 0.}]}; 
                        If[TempVar[[1]] == TempVar[[2]], {TempVar[[1]]}, TempVar], {{-1., 0., 0.}}]]; AppendTo[intersections, {-1., 0., 0.}]; 
                   holeQ = False; If[ !(Chop[surfaceid - 16519.] == 0. || N[Chop[ray1]] < 
                        0.0001) || generationnumber == 1., minresonateraylength = 0., 
                    minresonateraylength = 0.0001]; Block[{RayL, SCx, SCy}, RayL = 1.; SCx = 1.; SCy = 1.; dummymatrix = 
                      Union[Select[intersections, ({RayL, SCx, SCy} = #1; {RayL, SCx, SCy} == {-1., 0., 0.} || N[Chop[RayL]] >= 
                            minresonateraylength) & ]]; intersections = dummymatrix]; 
                   unfilteredintersections = intersections; If[Length[intersections] > 1, 
                    dummymatrix = Take[intersections, {2}]; intersections = 
                      dummymatrix; 1., 1.]; AppendTo[intersectionslist, intersections]; 
                   If[holeQ, AppendTo[holeQlist, 1], AppendTo[holeQlist, 0]]; 
                   {rsx, rsy, rsz, rex, rey, 
                     rez, r1, r2, r3, r4, 
                     r5, r6, r7, r8, r9, 
                     wavelength, intensity, intensityscale, ray1, 
                     index, opticallength, phaseshift, componentnumber, 
                     surfacenumber, surfaceid, ray2, ray3, 
                     sid, frontid, sno1, sno2, sno3, 
                     sno4, rayslot, oldgenerationnumber, intersectionnumber, 
                     pol1, pol2, pol3} = #1; index1 = 
                    index; xunit = r1*Re[Cos[0.017453292519943295*aa]] + r2*
                      Re[Sin[0.017453292519943295*aa]]; yunit = r2*Re[Cos[0.017453292519943295*aa]] - 
                     1.*r1*Re[Sin[0.017453292519943295*aa]]; zunit = r3*
                     Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]; xunit = Chop[xunit]; 
                   yunit = Chop[yunit]; zunit = Chop[zunit]; 
                   magnitude = (xunit^2 + yunit^2 + zunit^2)^2^(-1); 
                   xunit = xunit/magnitude; yunit = 
                    yunit/magnitude; zunit = zunit/
                     magnitude; rx = Re[Cos[0.017453292519943295*aa]]*(rex - 
                       1.*Re[xa1 + ta*Cos[0.017453292519943295*aa]]) + Re[Sin[0.017453292519943295*aa]]*(rey - 
                       1.*Re[ya1 + ta*Sin[0.017453292519943295*aa]]); ry = -1.*(rex - 
                       1.*Re[xa1 + ta*Cos[0.017453292519943295*aa]])*Re[Sin[0.017453292519943295*aa]] + Re[Cos[0.017453292519943295*aa]]*(rey - 
                       1.*Re[ya1 + ta*Sin[0.017453292519943295*aa]]); rz = rez*Re[Cos[0.017453292519943295*aa]^2 + 
                       Sin[0.017453292519943295*aa]^2]; rx = Chop[rx]; ry = 
                    Chop[ry]; rz = Chop[rz]; terminatedeflection = False; 
                   intersections = {{-1., 0., 0.}}; intersections = Block[{TempVar, TempVar1, TempVar2, TempVar3, RayL, 
                      SCx, SCy}, TempVar1 = 1.; TempVar2 = 1.; TempVar3 = 1.; TempVar = {{-1., 0., 0.}}; RayL = 1.; SCx = 1.; SCy = 1.; TempVar1 = 
                       rx^2 + ry^2 + rz^2; TempVar2 = rx*
                         xunit + ry*yunit + rz*zunit - 
                        xunit/ca2; TempVar3 = -TempVar1 + TempVar2^2 + (2*rx)/ca2; TempVar3 = Re[TempVar3]; 
                      If[N[Chop[TempVar3]] >= 0., TempVar3 = Sqrt[TempVar3]; TempVar = {{RayL, SCx, SCy} = {-TempVar2 - TempVar3, ry - 
                              (TempVar2 + TempVar3)*yunit, rz - (TempVar2 + TempVar3)*zunit}; 
                           If[N[Chop[RayL]] >= 0. && SCx^2 + SCy^2 <= 0.25*Abs[apa]^2 && Abs[ry + RayL*yunit - SCx] + 
                               Abs[rz + RayL*zunit - SCy] + Abs[-ca2^(-1) + rx + 
                                 RayL*xunit + Sqrt[ca2^(-2) - SCx^2 - SCy^2]/Sign[ca2]] < 1.*^-7, {RayL, SCx, SCy}, {-1., 0., 0.}], 
                          {RayL, SCx, SCy} = {-TempVar2 + TempVar3, ry + (-TempVar2 + TempVar3)*yunit, 
                             rz + (-TempVar2 + TempVar3)*zunit}; If[N[Chop[RayL]] >= 0. && SCx^2 + SCy^2 <= 
                              0.25*Abs[apa]^2 && Abs[ry + RayL*yunit - SCx] + Abs[rz + 
                                 RayL*zunit - SCy] + Abs[-ca2^(-1) + rx + RayL*xunit + 
                                 Sqrt[ca2^(-2) - SCx^2 - SCy^2]/Sign[ca2]] < 1.*^-7, {RayL, SCx, SCy}, {-1., 0., 0.}]}; If[TempVar[[1]] == TempVar[[2]], 
                         {TempVar[[1]]}, TempVar], {{-1., 0., 0.}}]]; AppendTo[intersections, {-1., 0., 0.}]; holeQ = 
                    False; If[ !(Chop[surfaceid - 16520.] == 0. || N[Chop[ray1]] < 0.0001) || 
                     generationnumber == 1., minresonateraylength = 0., minresonateraylength = 
                     0.0001]; Block[{RayL, SCx, SCy}, RayL = 1.; SCx = 1.; SCy = 1.; dummymatrix = Union[Select[
                        intersections, ({RayL, SCx, SCy} = #1; {RayL, SCx, SCy} == {-1., 0., 0.} || N[Chop[RayL]] >= 
                            minresonateraylength) & ]]; intersections = dummymatrix]; 
                   unfilteredintersections = intersections; If[Length[intersections] > 1, 
                    dummymatrix = Take[intersections, {2}]; intersections = 
                      dummymatrix; 1., 1.]; AppendTo[intersectionslist, intersections]; 
                   If[holeQ, AppendTo[holeQlist, 1], AppendTo[holeQlist, 0]]; 
                   {rsx, rsy, rsz, rex, rey, 
                     rez, r1, r2, r3, r4, 
                     r5, r6, r7, r8, r9, 
                     wavelength, intensity, intensityscale, ray1, 
                     index, opticallength, phaseshift, componentnumber, 
                     surfacenumber, surfaceid, ray2, ray3, 
                     sid, frontid, sno1, sno2, sno3, 
                     sno4, rayslot, oldgenerationnumber, intersectionnumber, 
                     pol1, pol2, pol3} = #1; index1 = 
                    index; xunit = r1*Re[Cos[0.017453292519943295*ab]] + r2*
                      Re[Sin[0.017453292519943295*ab]]; yunit = r2*Re[Cos[0.017453292519943295*ab]] - 
                     1.*r1*Re[Sin[0.017453292519943295*ab]]; zunit = r3*
                     Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]; xunit = Chop[xunit]; 
                   yunit = Chop[yunit]; zunit = Chop[zunit]; 
                   magnitude = (xunit^2 + yunit^2 + zunit^2)^2^(-1); 
                   xunit = xunit/magnitude; yunit = 
                    yunit/magnitude; zunit = zunit/
                     magnitude; rx = (rex - 1.*Re[xb1])*Re[Cos[0.017453292519943295*ab]] + 
                     (rey - 1.*Re[yb1])*Re[Sin[0.017453292519943295*ab]]; ry = (rey - 1.*Re[yb1])*
                      Re[Cos[0.017453292519943295*ab]] - 1.*(rex - 1.*Re[xb1])*Re[Sin[0.017453292519943295*ab]]; rz = 
                    rez*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]; rx = 
                    Chop[rx]; ry = Chop[ry]; rz = 
                    Chop[rz]; terminatedeflection = False; intersections = {{-1., 0., 0.}}; 
                   intersections = Block[{TempVar, TempVar1, TempVar2, TempVar3, RayL, SCx, SCy}, TempVar1 = 1.; TempVar2 = 1.; TempVar3 = 1.; 
                      TempVar = {{-1., 0., 0.}}; RayL = 1.; SCx = 1.; SCy = 1.; TempVar1 = rx^2 + ry^2 + 
                        rz^2; TempVar2 = rx*xunit + ry*
                         yunit + rz*zunit - xunit/cb1; 
                      TempVar3 = -TempVar1 + TempVar2^2 + (2*rx)/cb1; TempVar3 = Re[TempVar3]; If[N[Chop[TempVar3]] >= 0., 
                       TempVar3 = Sqrt[TempVar3]; TempVar = {{RayL, SCx, SCy} = {-TempVar2 - TempVar3, ry - (TempVar2 + TempVar3)*
                               yunit, rz - (TempVar2 + TempVar3)*zunit}; If[N[Chop[RayL]] >= 0. && 
                             SCx^2 + SCy^2 <= 0.25*Abs[apb]^2 && Abs[ry + RayL*yunit - SCx] + 
                               Abs[rz + RayL*zunit - SCy] + Abs[-cb1^(-1) + rx + 
                                 RayL*xunit + Sqrt[cb1^(-2) - SCx^2 - SCy^2]/Sign[cb1]] < 1.*^-7, {RayL, SCx, SCy}, {-1., 0., 0.}], 
                          {RayL, SCx, SCy} = {-TempVar2 + TempVar3, ry + (-TempVar2 + TempVar3)*yunit, 
                             rz + (-TempVar2 + TempVar3)*zunit}; If[N[Chop[RayL]] >= 0. && SCx^2 + SCy^2 <= 
                              0.25*Abs[apb]^2 && Abs[ry + RayL*yunit - SCx] + Abs[rz + 
                                 RayL*zunit - SCy] + Abs[-cb1^(-1) + rx + RayL*xunit + 
                                 Sqrt[cb1^(-2) - SCx^2 - SCy^2]/Sign[cb1]] < 1.*^-7, {RayL, SCx, SCy}, {-1., 0., 0.}]}; If[TempVar[[1]] == TempVar[[2]], 
                         {TempVar[[1]]}, TempVar], {{-1., 0., 0.}}]]; AppendTo[intersections, {-1., 0., 0.}]; holeQ = 
                    False; If[ !(Chop[surfaceid - 16547.] == 0. || N[Chop[ray1]] < 0.0001) || 
                     generationnumber == 1., minresonateraylength = 0., minresonateraylength = 
                     0.0001]; Block[{RayL, SCx, SCy}, RayL = 1.; SCx = 1.; SCy = 1.; dummymatrix = Union[Select[
                        intersections, ({RayL, SCx, SCy} = #1; {RayL, SCx, SCy} == {-1., 0., 0.} || N[Chop[RayL]] >= 
                            minresonateraylength) & ]]; intersections = dummymatrix]; 
                   unfilteredintersections = intersections; If[Length[intersections] > 1, 
                    dummymatrix = Take[intersections, {2}]; intersections = 
                      dummymatrix; 1., 1.]; AppendTo[intersectionslist, intersections]; 
                   If[holeQ, AppendTo[holeQlist, 1], AppendTo[holeQlist, 0]]; 
                   {rsx, rsy, rsz, rex, rey, 
                     rez, r1, r2, r3, r4, 
                     r5, r6, r7, r8, r9, 
                     wavelength, intensity, intensityscale, ray1, 
                     index, opticallength, phaseshift, componentnumber, 
                     surfacenumber, surfaceid, ray2, ray3, 
                     sid, frontid, sno1, sno2, sno3, 
                     sno4, rayslot, oldgenerationnumber, intersectionnumber, 
                     pol1, pol2, pol3} = #1; index1 = 
                    index; xunit = r1*Re[Cos[0.017453292519943295*ab]] + r2*
                      Re[Sin[0.017453292519943295*ab]]; yunit = r2*Re[Cos[0.017453292519943295*ab]] - 
                     1.*r1*Re[Sin[0.017453292519943295*ab]]; zunit = r3*
                     Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]; xunit = Chop[xunit]; 
                   yunit = Chop[yunit]; zunit = Chop[zunit]; 
                   magnitude = (xunit^2 + yunit^2 + zunit^2)^2^(-1); 
                   xunit = xunit/magnitude; yunit = 
                    yunit/magnitude; zunit = zunit/
                     magnitude; rx = Re[Cos[0.017453292519943295*ab]]*(rex - 
                       1.*Re[xb1 + tb*Cos[0.017453292519943295*ab]]) + Re[Sin[0.017453292519943295*ab]]*(rey - 
                       1.*Re[yb1 + tb*Sin[0.017453292519943295*ab]]); ry = -1.*(rex - 
                       1.*Re[xb1 + tb*Cos[0.017453292519943295*ab]])*Re[Sin[0.017453292519943295*ab]] + Re[Cos[0.017453292519943295*ab]]*(rey - 
                       1.*Re[yb1 + tb*Sin[0.017453292519943295*ab]]); rz = rez*Re[Cos[0.017453292519943295*ab]^2 + 
                       Sin[0.017453292519943295*ab]^2]; rx = Chop[rx]; ry = 
                    Chop[ry]; rz = Chop[rz]; terminatedeflection = False; 
                   intersections = {{-1., 0., 0.}}; intersections = Block[{TempVar, TempVar1, TempVar2, TempVar3, RayL, 
                      SCx, SCy}, TempVar1 = 1.; TempVar2 = 1.; TempVar3 = 1.; TempVar = {{-1., 0., 0.}}; RayL = 1.; SCx = 1.; SCy = 1.; TempVar1 = 
                       rx^2 + ry^2 + rz^2; TempVar2 = rx*
                         xunit + ry*yunit + rz*zunit - 
                        xunit/cb2; TempVar3 = -TempVar1 + TempVar2^2 + (2*rx)/cb2; TempVar3 = Re[TempVar3]; 
                      If[N[Chop[TempVar3]] >= 0., TempVar3 = Sqrt[TempVar3]; TempVar = {{RayL, SCx, SCy} = {-TempVar2 - TempVar3, ry - 
                              (TempVar2 + TempVar3)*yunit, rz - (TempVar2 + TempVar3)*zunit}; 
                           If[N[Chop[RayL]] >= 0. && SCx^2 + SCy^2 <= 0.25*Abs[apb]^2 && Abs[ry + RayL*yunit - SCx] + 
                               Abs[rz + RayL*zunit - SCy] + Abs[-cb2^(-1) + rx + 
                                 RayL*xunit + Sqrt[cb2^(-2) - SCx^2 - SCy^2]/Sign[cb2]] < 1.*^-7, {RayL, SCx, SCy}, {-1., 0., 0.}], 
                          {RayL, SCx, SCy} = {-TempVar2 + TempVar3, ry + (-TempVar2 + TempVar3)*yunit, 
                             rz + (-TempVar2 + TempVar3)*zunit}; If[N[Chop[RayL]] >= 0. && SCx^2 + SCy^2 <= 
                              0.25*Abs[apb]^2 && Abs[ry + RayL*yunit - SCx] + Abs[rz + 
                                 RayL*zunit - SCy] + Abs[-cb2^(-1) + rx + RayL*xunit + 
                                 Sqrt[cb2^(-2) - SCx^2 - SCy^2]/Sign[cb2]] < 1.*^-7, {RayL, SCx, SCy}, {-1., 0., 0.}]}; If[TempVar[[1]] == TempVar[[2]], 
                         {TempVar[[1]]}, TempVar], {{-1., 0., 0.}}]]; AppendTo[intersections, {-1., 0., 0.}]; holeQ = 
                    False; If[ !(Chop[surfaceid - 16548.] == 0. || N[Chop[ray1]] < 0.0001) || 
                     generationnumber == 1., minresonateraylength = 0., minresonateraylength = 
                     0.0001]; Block[{RayL, SCx, SCy}, RayL = 1.; SCx = 1.; SCy = 1.; dummymatrix = Union[Select[
                        intersections, ({RayL, SCx, SCy} = #1; {RayL, SCx, SCy} == {-1., 0., 0.} || N[Chop[RayL]] >= 
                            minresonateraylength) & ]]; intersections = dummymatrix]; 
                   unfilteredintersections = intersections; If[Length[intersections] > 1, 
                    dummymatrix = Take[intersections, {2}]; intersections = 
                      dummymatrix; 1., 1.]; AppendTo[intersectionslist, intersections]; 
                   If[holeQ, AppendTo[holeQlist, 1], AppendTo[holeQlist, 0]]); 
                  holeQlist = Drop[holeQlist, 1]; Drop[intersectionslist, 1]]); Which[
                resonatenumber == 1, {rsx, rsy, rsz, 
                   rex, rey, rez, r1, r2, 
                   r3, r4, r5, r6, r7, 
                   r8, r9, wavelength, intensity, 
                   intensityscale, ray1, index, opticallength, 
                   phaseshift, componentnumber, surfacenumber, surfaceid, 
                   ray2, ray3, sid, frontid, sno1, 
                   sno2, sno3, sno4, rayslot, 
                   oldgenerationnumber, intersectionnumber, pol1, pol2, 
                   pol3} = #1; index1 = index; xunit = 
                  r1*Re[Cos[0.017453292519943295*aa]] + r2*Re[Sin[0.017453292519943295*aa]]; 
                 yunit = r2*Re[Cos[0.017453292519943295*aa]] - 1.*r1*
                    Re[Sin[0.017453292519943295*aa]]; zunit = r3*Re[Cos[0.017453292519943295*aa]^2 + 
                     Sin[0.017453292519943295*aa]^2]; xunit = Chop[xunit]; yunit = 
                  Chop[yunit]; zunit = Chop[zunit]; magnitude = 
                  (xunit^2 + yunit^2 + zunit^2)^2^(-1); xunit = 
                  xunit/magnitude; yunit = yunit/
                   magnitude; zunit = zunit/magnitude; 
                 rx = (rex - 1.*Re[xa1])*Re[Cos[0.017453292519943295*aa]] + (rey - 1.*Re[ya1])*
                    Re[Sin[0.017453292519943295*aa]]; ry = (rey - 1.*Re[ya1])*Re[Cos[0.017453292519943295*aa]] - 
                   1.*(rex - 1.*Re[xa1])*Re[Sin[0.017453292519943295*aa]]; rz = rez*
                   Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]; rx = Chop[rx]; 
                 ry = Chop[ry]; rz = Chop[rz]; 
                 terminatedeflection = False; intersections = {{-1., 0., 0.}}; intersections = 
                  resonateintersections; holeQ = holeQlist[[
                    resonatenumber]] == 1; {newcomponentnumber, newsurfacenumber, 
                   newsurfaceid} = {1., 1., 16519.}; {componentnumber, surfacenumber, 
                   surfaceid} = {1., 1., 16519.}; {ray1, ray2, ray3} = 
                  intersections[[1]]; rsx = rex; rsy = 
                  rey; rsz = rez; rxhlp = rx + 
                   ray1*xunit; ryhlp = ry + ray1*
                    yunit; rzhlp = rz + ray1*zunit; 
                 rex = Re[xa1] + rxhlp*Re[Cos[0.017453292519943295*aa]] - 1.*ryhlp*
                    Re[Sin[0.017453292519943295*aa]]; rey = Re[ya1] + ryhlp*Re[Cos[0.017453292519943295*aa]] + 
                   rxhlp*Re[Sin[0.017453292519943295*aa]]; rez = rzhlp*
                   Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]; rex = Chop[rex]; 
                 rey = Chop[rey]; rez = Chop[rez]; 
                 If[rex == 0 && rey == 0 && rez == 0, rex = 
                    rxhlp + Re[xa1]; rey = ryhlp + Re[ya1]; rez = 
                    rzhlp, rex]; intersectionnumber = intersectionnumber + 
                   1.; opticallength = opticallength + ray1*index1; 
                 completedray = {rsx, rsy, rsz, rex, 
                   rey, rez, r1, r2, r3, 
                   r4, r5, r6, r7, r8, 
                   r9, wavelength, intensity, intensityscale, 
                   ray1, index, opticallength, phaseshift, 
                   componentnumber, surfacenumber, surfaceid, ray2, 
                   ray3, sid, frontid, sno1, sno2, 
                   sno3, sno4, rayslot, generationnumber, 
                   intersectionnumber, pol1, pol2, pol3}; 
                 reportedrayspresent = True; intersectionspresent = True; deflectedrayspresent = 
                  True; localdummyoutrayspresent = True; (surfnorm = 
                   1.*Re[{{1., -1.*(ray2/Sign[ca1]/(ca1^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5), 
                       -1.*(ray3/Sign[ca1]/(ca1^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5)}, 
                      {ray2/Sign[ca1]/(ca1^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5, 1., 0.}, 
                      {ray3/Sign[ca1]/(ca1^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5, 0., 1.}}]; 
                  dummymatrix = ({{(Re[Cos[0.017453292519943295*aa]]*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2])/
                         (Re[Cos[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2] + Re[Sin[0.017453292519943295*aa]]^2*
                           Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]), (-1.*Re[Sin[0.017453292519943295*aa]]*
                          Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2])/(Re[Cos[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^
                              2 + Sin[0.017453292519943295*aa]^2] + Re[Sin[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*
                                aa]^2]), 0.}, {(Re[Sin[0.017453292519943295*aa]]*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2])/
                         (Re[Cos[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2] + Re[Sin[0.017453292519943295*aa]]^2*
                           Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]), (Re[Cos[0.017453292519943295*aa]]*Re[Cos[0.017453292519943295*aa]^
                             2 + Sin[0.017453292519943295*aa]^2])/(Re[Cos[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*
                                aa]^2] + Re[Sin[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]), 0.}, 
                       {0., 0., (Re[Cos[0.017453292519943295*aa]]^2 + Re[Sin[0.017453292519943295*aa]]^2)/(Re[Cos[0.017453292519943295*aa]]^2*
                           Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2] + Re[Sin[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^
                              2 + Sin[0.017453292519943295*aa]^2])}} . #1 & ) /@ surfnorm; surfnorm = 
                   (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; surfacerotationmatrix = {{Re[Cos[0.017453292519943295*aa]], 
                     Re[Sin[0.017453292519943295*aa]], 0.}, {-1.*Re[Sin[0.017453292519943295*aa]], Re[Cos[0.017453292519943295*aa]], 0.}, 
                    {0., 0., Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]}}; surfacetranslationvector = 
                   {Re[xa1], Re[ya1], 0.}; activeregion = 1.; completedray = {rsx, 
                    rsy, rsz, rex, rey, rez, 
                    r1, r2, r3, r4, r5, 
                    r6, r7, r8, r9, wavelength, 
                    intensity, intensityscale, ray1, index, 
                    opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; componentnumber = 
                   newcomponentnumber; surfacenumber = newsurfacenumber; 
                  surfaceid = newsurfaceid; originalray = {rsx, 
                    rsy, rsz, rex, rey, rez, 
                    r1, r2, r3, r4, r5, 
                    r6, r7, r8, r9, wavelength, 
                    intensity, intensityscale, ray1, index, 
                    opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; ZAxis = 
                   First[surfnorm]*Sign[First[surfnorm] . {r1, r2, 
                       r3}]; raytiltz = ZAxis . {r1, 
                     r2, r3}; raytiltx = 1 - raytiltz^2; 
                  If[raytiltx < 2.220446049250313*^-16, raytiltx = 0., raytiltx = 
                    Sqrt[raytiltx]]; XAxis = {r1, r2, 
                     r3} - ZAxis*raytiltz; If[Abs[raytiltx] < 1.*^-7, 
                   XAxis = surfnorm[[2]], XAxis = XAxis/
                     raytiltx]; YAxis = {-(ZAxis[[3]]*XAxis[[2]]) + 
                     ZAxis[[2]]*XAxis[[3]], ZAxis[[3]]*XAxis[[1]] - 
                     ZAxis[[1]]*XAxis[[3]], -(ZAxis[[2]]*XAxis[[1]]) + 
                     ZAxis[[1]]*XAxis[[2]]}; YAxis = YAxis/
                    Sqrt[YAxis . YAxis]; PAxis = ZAxis*
                     {r1, r2, r3} . XAxis - XAxis*
                     {r1, r2, r3} . ZAxis; PAxis = 
                   PAxis/Sqrt[PAxis . PAxis]; normalizedpolarization = 
                   {pol1, pol2, pol3}[[{1, 2}]]; normalizedpolarization = 
                   normalizedpolarization/Sqrt[normalizedpolarization . normalizedpolarization]; 
                  allcreatedamplitudes = {{0.*I, 0.*I, 0.*I}}; amplitude = Sqrt[intensity]*
                    (normalizedpolarization[[1]]*{r4, r5, r6} + 
                     normalizedpolarization[[2]]*{r7, r8, r9}*
                      Exp[I*pol3*0.017453292519943295]); perpendicularamplitude = amplitude . 
                    YAxis; parallelamplitude = amplitude . PAxis; 
                  If[activeregion*(surfacerotationmatrix . {r1, r2, 
                        r3})[[1]] > 0, index2 = Sqrt[2.2718929 - 0.010108077*wavelength^2 + 
                       0.010592509/wavelength^2 + 0.00020816965/wavelength^4 - 7.6472538*^-6/
                        wavelength^6 + 4.9240991*^-7/wavelength^8]; index1 = 
                     index; If[Abs[index2 - index1] < 0.01, alternateindex = 
                       Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/wavelength^
                           2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/wavelength^6 + 
                         2.830549017940716*^-12/wavelength^8]; If[Abs[alternateindex - index1] > 
                        Abs[index2 - index1], index1 = index2; 
                        index2 = alternateindex, index1 = alternateindex], 
                     index1 = Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/
                         wavelength^2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/
                         wavelength^6 + 2.830549017940716*^-12/wavelength^8]], index2 = 
                     Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/wavelength^
                         2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/wavelength^6 + 
                       2.830549017940716*^-12/wavelength^8]; index1 = index; 
                    If[Abs[index2 - index1] < 0.01, alternateindex = Sqrt[2.2718929 - 
                         0.010108077*wavelength^2 + 0.010592509/wavelength^2 + 0.00020816965/
                          wavelength^4 - 7.6472538*^-6/wavelength^6 + 4.9240991*^-7/wavelength^8]; 
                      If[Abs[alternateindex - index1] > Abs[index2 - index1], 
                       index1 = index2; index2 = alternateindex, 
                       index1 = alternateindex], index1 = Sqrt[2.2718929 - 0.010108077*
                         wavelength^2 + 0.010592509/wavelength^2 + 0.00020816965/wavelength^4 - 
                        7.6472538*^-6/wavelength^6 + 4.9240991*^-7/wavelength^8]]]; n12 = 
                   index1/index2; index = index2; 
                  If[1 - n12^2*raytiltx^2 >= 0., dummyvector = n12*
                       {r1, r2, r3} + ZAxis*(-(n12*
                          raytiltz) + Sqrt[1 - n12^2*raytiltx^2]); {r1, 
                      r2, r3} = dummyvector; 1., 
                   If[{r1, r2, r3} . surfnorm[[1]] > 0, 
                     phaseshift = phaseshift + 0.5, 1.]; index2 = index1; 
                    index = index1; dummyvector = {r1, 
                       r2, r3} - 2*ZAxis*ZAxis . {r1, 
                         r2, r3}; {r1, r2, r3} = 
                     dummyvector; dummyvector = {pol1, pol2, 
                       pol3} + {0, 0, -180}; {pol1, pol2, pol3} = 
                     dummyvector; 1.]; allcreatedrays = {{rsx, rsy, 
                     rsz, rex, rey, rez, r1, 
                     r2, r3, r4, r5, r6, 
                     r7, r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     generationnumber, intersectionnumber, pol1, pol2, 
                     pol3}}; originalrotation = Partition[Take[originalray, {7, 15}], 3]; 
                  previoustilt = originalrotation[[1]]; transmatrix = 
                   MapIndexed[({rsx, rsy, rsz, rex, 
                        rey, rez, r1, r2, r3, 
                        r4, r5, r6, r7, r8, 
                        r9, wavelength, intensity, intensityscale, 
                        ray1, index, opticallength, phaseshift, 
                        componentnumber, surfacenumber, surfaceid, ray2, 
                        ray3, sid, frontid, sno1, sno2, 
                        sno3, sno4, rayslot, oldgenerationnumber, 
                        intersectionnumber, pol1, pol2, pol3} = #1; 
                      rotationaxis = Cross[{1., 0., 0.}, {r1, r2, r3}]; 
                      axismagnitude = Sqrt[rotationaxis . rotationaxis]; 
                      If[Chop[axismagnitude] != 0., rotationaxis = rotationaxis/
                          axismagnitude; axismatrix = (#1*rotationaxis & ) /@ 
                          rotationaxis; crossmatrix = {{0., -rotationaxis[[3]], 
                           rotationaxis[[2]]}, {rotationaxis[[3]], 0., -rotationaxis[[1]]}, 
                          {-rotationaxis[[2]], rotationaxis[[1]], 0.}}; cosangle = 
                         {1., 0., 0.} . {r1, r2, r3}; If[1. - cosangle < 
                          2.220446049250313*^-16, cosangle = 1.]; sinangle = 1 - cosangle^2; 
                        If[sinangle < 2.220446049250313*^-16, sinangle = 0., sinangle = 
                          Sqrt[sinangle]]; axismatrix + ({{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}} - 
                           axismatrix)*cosangle + crossmatrix*sinangle, 
                       {{r1, r2, r3}, {r4, r5, 
                          r6}, {r7, r8, r9}} = 
                        If[Sign[{1., 0., 0.} . {r1, r2, r3}] === 1, {{1., 0., 0.}, {0., 1., 0.}, 
                          {0., 0., 1.}}, {{-1., 0., 0.}, {0., -1., 0.}, {0., 0., 1.}}]]) & , allcreatedrays]; 
                  allcreatedrays = MapIndexed[({rsx, rsy, rsz, 
                        rex, rey, rez, r1, r2, 
                        r3, r4, r5, r6, r7, 
                        r8, r9, wavelength, intensity, 
                        intensityscale, ray1, index, opticallength, 
                        phaseshift, componentnumber, surfacenumber, 
                        surfaceid, ray2, ray3, sid, 
                        frontid, sno1, sno2, sno3, 
                        sno4, rayslot, oldgenerationnumber, 
                        intersectionnumber, pol1, pol2, pol3} = #1; 
                      {{r1, r2, r3}, {r4, r5, 
                         r6}, {r7, r8, r9}} = 
                       Function[mm, (mm . #1 & ) /@ {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}][
                        transmatrix[[First[#2]]]]; tiltmagnitude = {r1, r2, 
                         r3} . {r1, r2, r3}; 
                      If[Abs[tiltmagnitude - 1.] > 1.*^-6, dummyvector = {r1, 
                           r2, r3}/Sqrt[tiltmagnitude]; {r1, 
                          r2, r3} = dummyvector; dummymatrix = 
                         (#1/Sqrt[#1 . #1] & ) /@ {{r1, r2, r3}, {r4, 
                            r5, r6}, {r7, r8, r9}}; 
                        {{r1, r2, r3}, {r4, r5, 
                           r6}, {r7, r8, r9}} = 
                         dummymatrix; 1., 1.]; If[allcreatedamplitudes != {{0.*I, 0.*I, 0.*I}}, 
                       amplitude = allcreatedamplitudes[[First[#2]]]; hor = 
                         {r4, r5, r6} . amplitude; 
                        ver = {r7, r8, r9} . amplitude; 
                        If[Abs[hor] > 0. || Abs[ver] > 0., {pol1, pol2, 
                            pol3} = {Abs[hor], Abs[ver], 0.}; dummyvector = 
                           {pol1, pol2, pol3}/Sqrt[{pol1, 
                               pol2, pol3} . {pol1, pol2, 
                               pol3}]; {pol1, pol2, pol3} = 
                           dummyvector; dummyvector = ReplacePart[{pol1, 
                             pol2, pol3}, (If[ver == 0., 0., Arg[ver]] - 
                              If[hor == 0., 0., Arg[hor]])/Degree, 3]; {pol1, 
                            pol2, pol3} = dummyvector]; 1., 1.]; {rsx, 
                       rsy, rsz, rex, rey, rez, 
                       r1, r2, r3, r4, r5, 
                       r6, r7, r8, r9, wavelength, 
                       intensity, intensityscale, ray1, index, 
                       opticallength, phaseshift, componentnumber, 
                       surfacenumber, surfaceid, ray2, ray3, 
                       sid, frontid, sno1, sno2, sno3, 
                       sno4, rayslot, generationnumber, intersectionnumber, 
                       pol1, pol2, pol3}) & , allcreatedrays]); 
                 {allcreatedrays, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 
                    0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, {completedray}}, 
                resonatenumber == 2, {rsx, rsy, rsz, 
                   rex, rey, rez, r1, r2, 
                   r3, r4, r5, r6, r7, 
                   r8, r9, wavelength, intensity, 
                   intensityscale, ray1, index, opticallength, 
                   phaseshift, componentnumber, surfacenumber, surfaceid, 
                   ray2, ray3, sid, frontid, sno1, 
                   sno2, sno3, sno4, rayslot, 
                   oldgenerationnumber, intersectionnumber, pol1, pol2, 
                   pol3} = #1; index1 = index; xunit = 
                  r1*Re[Cos[0.017453292519943295*aa]] + r2*Re[Sin[0.017453292519943295*aa]]; 
                 yunit = r2*Re[Cos[0.017453292519943295*aa]] - 1.*r1*
                    Re[Sin[0.017453292519943295*aa]]; zunit = r3*Re[Cos[0.017453292519943295*aa]^2 + 
                     Sin[0.017453292519943295*aa]^2]; xunit = Chop[xunit]; yunit = 
                  Chop[yunit]; zunit = Chop[zunit]; magnitude = 
                  (xunit^2 + yunit^2 + zunit^2)^2^(-1); xunit = 
                  xunit/magnitude; yunit = yunit/
                   magnitude; zunit = zunit/magnitude; 
                 rx = Re[Cos[0.017453292519943295*aa]]*(rex - 1.*Re[xa1 + ta*Cos[0.017453292519943295*aa]]) + 
                   Re[Sin[0.017453292519943295*aa]]*(rey - 1.*Re[ya1 + ta*Sin[0.017453292519943295*aa]]); ry = 
                  -1.*(rex - 1.*Re[xa1 + ta*Cos[0.017453292519943295*aa]])*Re[Sin[0.017453292519943295*aa]] + Re[Cos[0.017453292519943295*aa]]*
                    (rey - 1.*Re[ya1 + ta*Sin[0.017453292519943295*aa]]); rz = rez*
                   Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]; rx = Chop[rx]; 
                 ry = Chop[ry]; rz = Chop[rz]; 
                 terminatedeflection = False; intersections = {{-1., 0., 0.}}; intersections = 
                  resonateintersections; holeQ = holeQlist[[
                    resonatenumber]] == 1; {newcomponentnumber, newsurfacenumber, 
                   newsurfaceid} = {1., 2., 16520.}; {componentnumber, surfacenumber, 
                   surfaceid} = {1., 2., 16520.}; {ray1, ray2, ray3} = 
                  intersections[[1]]; rsx = rex; rsy = 
                  rey; rsz = rez; rxhlp = rx + 
                   ray1*xunit; ryhlp = ry + ray1*
                    yunit; rzhlp = rz + ray1*zunit; 
                 rex = rxhlp*Re[Cos[0.017453292519943295*aa]] + Re[xa1 + ta*Cos[0.017453292519943295*aa]] - 
                   1.*ryhlp*Re[Sin[0.017453292519943295*aa]]; rey = ryhlp*
                    Re[Cos[0.017453292519943295*aa]] + rxhlp*Re[Sin[0.017453292519943295*aa]] + Re[ya1 + ta*Sin[0.017453292519943295*aa]]; 
                 rez = rzhlp*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]; 
                 rex = Chop[rex]; rey = Chop[rey]; 
                 rez = Chop[rez]; If[rex == 0 && rey == 0 && 
                   rez == 0, rex = rxhlp + Re[xa1 + ta*Cos[0.017453292519943295*aa]]; 
                   rey = ryhlp + Re[ya1 + ta*Sin[0.017453292519943295*aa]]; rez = 
                    rzhlp, rex]; intersectionnumber = intersectionnumber + 
                   1.; opticallength = opticallength + ray1*index1; 
                 completedray = {rsx, rsy, rsz, rex, 
                   rey, rez, r1, r2, r3, 
                   r4, r5, r6, r7, r8, 
                   r9, wavelength, intensity, intensityscale, 
                   ray1, index, opticallength, phaseshift, 
                   componentnumber, surfacenumber, surfaceid, ray2, 
                   ray3, sid, frontid, sno1, sno2, 
                   sno3, sno4, rayslot, generationnumber, 
                   intersectionnumber, pol1, pol2, pol3}; 
                 reportedrayspresent = True; intersectionspresent = True; deflectedrayspresent = 
                  True; localdummyoutrayspresent = True; (surfnorm = 
                   1.*Re[{{1., -1.*(ray2/Sign[ca2]/(ca2^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5), 
                       -1.*(ray3/Sign[ca2]/(ca2^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5)}, 
                      {ray2/Sign[ca2]/(ca2^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5, 1., 0.}, 
                      {ray3/Sign[ca2]/(ca2^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5, 0., 1.}}]; 
                  dummymatrix = ({{(Re[Cos[0.017453292519943295*aa]]*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2])/
                         (Re[Cos[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2] + Re[Sin[0.017453292519943295*aa]]^2*
                           Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]), (-1.*Re[Sin[0.017453292519943295*aa]]*
                          Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2])/(Re[Cos[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^
                              2 + Sin[0.017453292519943295*aa]^2] + Re[Sin[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*
                                aa]^2]), 0.}, {(Re[Sin[0.017453292519943295*aa]]*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2])/
                         (Re[Cos[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2] + Re[Sin[0.017453292519943295*aa]]^2*
                           Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]), (Re[Cos[0.017453292519943295*aa]]*Re[Cos[0.017453292519943295*aa]^
                             2 + Sin[0.017453292519943295*aa]^2])/(Re[Cos[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*
                                aa]^2] + Re[Sin[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]), 0.}, 
                       {0., 0., (Re[Cos[0.017453292519943295*aa]]^2 + Re[Sin[0.017453292519943295*aa]]^2)/(Re[Cos[0.017453292519943295*aa]]^2*
                           Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2] + Re[Sin[0.017453292519943295*aa]]^2*Re[Cos[0.017453292519943295*aa]^
                              2 + Sin[0.017453292519943295*aa]^2])}} . #1 & ) /@ surfnorm; surfnorm = 
                   (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; surfacerotationmatrix = {{Re[Cos[0.017453292519943295*aa]], 
                     Re[Sin[0.017453292519943295*aa]], 0.}, {-1.*Re[Sin[0.017453292519943295*aa]], Re[Cos[0.017453292519943295*aa]], 0.}, 
                    {0., 0., Re[Cos[0.017453292519943295*aa]^2 + Sin[0.017453292519943295*aa]^2]}}; surfacetranslationvector = 
                   {Re[xa1 + ta*Cos[0.017453292519943295*aa]], Re[ya1 + ta*Sin[0.017453292519943295*aa]], 0.}; activeregion = -1.; 
                  completedray = {rsx, rsy, rsz, rex, 
                    rey, rez, r1, r2, r3, 
                    r4, r5, r6, r7, r8, 
                    r9, wavelength, intensity, intensityscale, 
                    ray1, index, opticallength, phaseshift, 
                    componentnumber, surfacenumber, surfaceid, ray2, 
                    ray3, sid, frontid, sno1, sno2, 
                    sno3, sno4, rayslot, generationnumber, 
                    intersectionnumber, pol1, pol2, pol3}; 
                  componentnumber = newcomponentnumber; surfacenumber = 
                   newsurfacenumber; surfaceid = newsurfaceid; originalray = 
                   {rsx, rsy, rsz, rex, rey, 
                    rez, r1, r2, r3, r4, 
                    r5, r6, r7, r8, r9, 
                    wavelength, intensity, intensityscale, ray1, 
                    index, opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; ZAxis = 
                   First[surfnorm]*Sign[First[surfnorm] . {r1, r2, 
                       r3}]; raytiltz = ZAxis . {r1, 
                     r2, r3}; raytiltx = 1 - raytiltz^2; 
                  If[raytiltx < 2.220446049250313*^-16, raytiltx = 0., raytiltx = 
                    Sqrt[raytiltx]]; XAxis = {r1, r2, 
                     r3} - ZAxis*raytiltz; If[Abs[raytiltx] < 1.*^-7, 
                   XAxis = surfnorm[[2]], XAxis = XAxis/
                     raytiltx]; YAxis = {-(ZAxis[[3]]*XAxis[[2]]) + 
                     ZAxis[[2]]*XAxis[[3]], ZAxis[[3]]*XAxis[[1]] - 
                     ZAxis[[1]]*XAxis[[3]], -(ZAxis[[2]]*XAxis[[1]]) + 
                     ZAxis[[1]]*XAxis[[2]]}; YAxis = YAxis/
                    Sqrt[YAxis . YAxis]; PAxis = ZAxis*
                     {r1, r2, r3} . XAxis - XAxis*
                     {r1, r2, r3} . ZAxis; PAxis = 
                   PAxis/Sqrt[PAxis . PAxis]; normalizedpolarization = 
                   {pol1, pol2, pol3}[[{1, 2}]]; normalizedpolarization = 
                   normalizedpolarization/Sqrt[normalizedpolarization . normalizedpolarization]; 
                  allcreatedamplitudes = {{0.*I, 0.*I, 0.*I}}; amplitude = Sqrt[intensity]*
                    (normalizedpolarization[[1]]*{r4, r5, r6} + 
                     normalizedpolarization[[2]]*{r7, r8, r9}*
                      Exp[I*pol3*0.017453292519943295]); perpendicularamplitude = amplitude . 
                    YAxis; parallelamplitude = amplitude . PAxis; 
                  If[activeregion*(surfacerotationmatrix . {r1, r2, 
                        r3})[[1]] > 0, index2 = Sqrt[2.2718929 - 0.010108077*wavelength^2 + 
                       0.010592509/wavelength^2 + 0.00020816965/wavelength^4 - 7.6472538*^-6/
                        wavelength^6 + 4.9240991*^-7/wavelength^8]; index1 = 
                     index; If[Abs[index2 - index1] < 0.01, alternateindex = 
                       Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/wavelength^
                           2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/wavelength^6 + 
                         2.830549017940716*^-12/wavelength^8]; If[Abs[alternateindex - index1] > 
                        Abs[index2 - index1], index1 = index2; 
                        index2 = alternateindex, index1 = alternateindex], 
                     index1 = Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/
                         wavelength^2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/
                         wavelength^6 + 2.830549017940716*^-12/wavelength^8]], index2 = 
                     Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/wavelength^
                         2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/wavelength^6 + 
                       2.830549017940716*^-12/wavelength^8]; index1 = index; 
                    If[Abs[index2 - index1] < 0.01, alternateindex = Sqrt[2.2718929 - 
                         0.010108077*wavelength^2 + 0.010592509/wavelength^2 + 0.00020816965/
                          wavelength^4 - 7.6472538*^-6/wavelength^6 + 4.9240991*^-7/wavelength^8]; 
                      If[Abs[alternateindex - index1] > Abs[index2 - index1], 
                       index1 = index2; index2 = alternateindex, 
                       index1 = alternateindex], index1 = Sqrt[2.2718929 - 0.010108077*
                         wavelength^2 + 0.010592509/wavelength^2 + 0.00020816965/wavelength^4 - 
                        7.6472538*^-6/wavelength^6 + 4.9240991*^-7/wavelength^8]]]; n12 = 
                   index1/index2; index = index2; 
                  If[1 - n12^2*raytiltx^2 >= 0., dummyvector = n12*
                       {r1, r2, r3} + ZAxis*(-(n12*
                          raytiltz) + Sqrt[1 - n12^2*raytiltx^2]); {r1, 
                      r2, r3} = dummyvector; 1., 
                   If[{r1, r2, r3} . surfnorm[[1]] > 0, 
                     phaseshift = phaseshift + 0.5, 1.]; index2 = index1; 
                    index = index1; dummyvector = {r1, 
                       r2, r3} - 2*ZAxis*ZAxis . {r1, 
                         r2, r3}; {r1, r2, r3} = 
                     dummyvector; dummyvector = {pol1, pol2, 
                       pol3} + {0, 0, -180}; {pol1, pol2, pol3} = 
                     dummyvector; 1.]; allcreatedrays = {{rsx, rsy, 
                     rsz, rex, rey, rez, r1, 
                     r2, r3, r4, r5, r6, 
                     r7, r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     generationnumber, intersectionnumber, pol1, pol2, 
                     pol3}}; originalrotation = Partition[Take[originalray, {7, 15}], 3]; 
                  previoustilt = originalrotation[[1]]; transmatrix = 
                   MapIndexed[({rsx, rsy, rsz, rex, 
                        rey, rez, r1, r2, r3, 
                        r4, r5, r6, r7, r8, 
                        r9, wavelength, intensity, intensityscale, 
                        ray1, index, opticallength, phaseshift, 
                        componentnumber, surfacenumber, surfaceid, ray2, 
                        ray3, sid, frontid, sno1, sno2, 
                        sno3, sno4, rayslot, oldgenerationnumber, 
                        intersectionnumber, pol1, pol2, pol3} = #1; 
                      rotationaxis = Cross[{1., 0., 0.}, {r1, r2, r3}]; 
                      axismagnitude = Sqrt[rotationaxis . rotationaxis]; 
                      If[Chop[axismagnitude] != 0., rotationaxis = rotationaxis/
                          axismagnitude; axismatrix = (#1*rotationaxis & ) /@ 
                          rotationaxis; crossmatrix = {{0., -rotationaxis[[3]], 
                           rotationaxis[[2]]}, {rotationaxis[[3]], 0., -rotationaxis[[1]]}, 
                          {-rotationaxis[[2]], rotationaxis[[1]], 0.}}; cosangle = 
                         {1., 0., 0.} . {r1, r2, r3}; If[1. - cosangle < 
                          2.220446049250313*^-16, cosangle = 1.]; sinangle = 1 - cosangle^2; 
                        If[sinangle < 2.220446049250313*^-16, sinangle = 0., sinangle = 
                          Sqrt[sinangle]]; axismatrix + ({{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}} - 
                           axismatrix)*cosangle + crossmatrix*sinangle, 
                       {{r1, r2, r3}, {r4, r5, 
                          r6}, {r7, r8, r9}} = 
                        If[Sign[{1., 0., 0.} . {r1, r2, r3}] === 1, {{1., 0., 0.}, {0., 1., 0.}, 
                          {0., 0., 1.}}, {{-1., 0., 0.}, {0., -1., 0.}, {0., 0., 1.}}]]) & , allcreatedrays]; 
                  allcreatedrays = MapIndexed[({rsx, rsy, rsz, 
                        rex, rey, rez, r1, r2, 
                        r3, r4, r5, r6, r7, 
                        r8, r9, wavelength, intensity, 
                        intensityscale, ray1, index, opticallength, 
                        phaseshift, componentnumber, surfacenumber, 
                        surfaceid, ray2, ray3, sid, 
                        frontid, sno1, sno2, sno3, 
                        sno4, rayslot, oldgenerationnumber, 
                        intersectionnumber, pol1, pol2, pol3} = #1; 
                      {{r1, r2, r3}, {r4, r5, 
                         r6}, {r7, r8, r9}} = 
                       Function[mm, (mm . #1 & ) /@ {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}][
                        transmatrix[[First[#2]]]]; tiltmagnitude = {r1, r2, 
                         r3} . {r1, r2, r3}; 
                      If[Abs[tiltmagnitude - 1.] > 1.*^-6, dummyvector = {r1, 
                           r2, r3}/Sqrt[tiltmagnitude]; {r1, 
                          r2, r3} = dummyvector; dummymatrix = 
                         (#1/Sqrt[#1 . #1] & ) /@ {{r1, r2, r3}, {r4, 
                            r5, r6}, {r7, r8, r9}}; 
                        {{r1, r2, r3}, {r4, r5, 
                           r6}, {r7, r8, r9}} = 
                         dummymatrix; 1., 1.]; If[allcreatedamplitudes != {{0.*I, 0.*I, 0.*I}}, 
                       amplitude = allcreatedamplitudes[[First[#2]]]; hor = 
                         {r4, r5, r6} . amplitude; 
                        ver = {r7, r8, r9} . amplitude; 
                        If[Abs[hor] > 0. || Abs[ver] > 0., {pol1, pol2, 
                            pol3} = {Abs[hor], Abs[ver], 0.}; dummyvector = 
                           {pol1, pol2, pol3}/Sqrt[{pol1, 
                               pol2, pol3} . {pol1, pol2, 
                               pol3}]; {pol1, pol2, pol3} = 
                           dummyvector; dummyvector = ReplacePart[{pol1, 
                             pol2, pol3}, (If[ver == 0., 0., Arg[ver]] - 
                              If[hor == 0., 0., Arg[hor]])/Degree, 3]; {pol1, 
                            pol2, pol3} = dummyvector]; 1., 1.]; {rsx, 
                       rsy, rsz, rex, rey, rez, 
                       r1, r2, r3, r4, r5, 
                       r6, r7, r8, r9, wavelength, 
                       intensity, intensityscale, ray1, index, 
                       opticallength, phaseshift, componentnumber, 
                       surfacenumber, surfaceid, ray2, ray3, 
                       sid, frontid, sno1, sno2, sno3, 
                       sno4, rayslot, generationnumber, intersectionnumber, 
                       pol1, pol2, pol3}) & , allcreatedrays]); 
                 {allcreatedrays, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 
                    0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, {completedray}}, 
                resonatenumber == 3, {rsx, rsy, rsz, 
                   rex, rey, rez, r1, r2, 
                   r3, r4, r5, r6, r7, 
                   r8, r9, wavelength, intensity, 
                   intensityscale, ray1, index, opticallength, 
                   phaseshift, componentnumber, surfacenumber, surfaceid, 
                   ray2, ray3, sid, frontid, sno1, 
                   sno2, sno3, sno4, rayslot, 
                   oldgenerationnumber, intersectionnumber, pol1, pol2, 
                   pol3} = #1; index1 = index; xunit = 
                  r1*Re[Cos[0.017453292519943295*ab]] + r2*Re[Sin[0.017453292519943295*ab]]; 
                 yunit = r2*Re[Cos[0.017453292519943295*ab]] - 1.*r1*
                    Re[Sin[0.017453292519943295*ab]]; zunit = r3*Re[Cos[0.017453292519943295*ab]^2 + 
                     Sin[0.017453292519943295*ab]^2]; xunit = Chop[xunit]; yunit = 
                  Chop[yunit]; zunit = Chop[zunit]; magnitude = 
                  (xunit^2 + yunit^2 + zunit^2)^2^(-1); xunit = 
                  xunit/magnitude; yunit = yunit/
                   magnitude; zunit = zunit/magnitude; 
                 rx = (rex - 1.*Re[xb1])*Re[Cos[0.017453292519943295*ab]] + (rey - 1.*Re[yb1])*
                    Re[Sin[0.017453292519943295*ab]]; ry = (rey - 1.*Re[yb1])*Re[Cos[0.017453292519943295*ab]] - 
                   1.*(rex - 1.*Re[xb1])*Re[Sin[0.017453292519943295*ab]]; rz = rez*
                   Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]; rx = Chop[rx]; 
                 ry = Chop[ry]; rz = Chop[rz]; 
                 terminatedeflection = False; intersections = {{-1., 0., 0.}}; intersections = 
                  resonateintersections; holeQ = holeQlist[[
                    resonatenumber]] == 1; {newcomponentnumber, newsurfacenumber, 
                   newsurfaceid} = {1., 3., 16547.}; {componentnumber, surfacenumber, 
                   surfaceid} = {1., 3., 16547.}; {ray1, ray2, ray3} = 
                  intersections[[1]]; rsx = rex; rsy = 
                  rey; rsz = rez; rxhlp = rx + 
                   ray1*xunit; ryhlp = ry + ray1*
                    yunit; rzhlp = rz + ray1*zunit; 
                 rex = Re[xb1] + rxhlp*Re[Cos[0.017453292519943295*ab]] - 1.*ryhlp*
                    Re[Sin[0.017453292519943295*ab]]; rey = Re[yb1] + ryhlp*Re[Cos[0.017453292519943295*ab]] + 
                   rxhlp*Re[Sin[0.017453292519943295*ab]]; rez = rzhlp*
                   Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]; rex = Chop[rex]; 
                 rey = Chop[rey]; rez = Chop[rez]; 
                 If[rex == 0 && rey == 0 && rez == 0, rex = 
                    rxhlp + Re[xb1]; rey = ryhlp + Re[yb1]; rez = 
                    rzhlp, rex]; intersectionnumber = intersectionnumber + 
                   1.; opticallength = opticallength + ray1*index1; 
                 completedray = {rsx, rsy, rsz, rex, 
                   rey, rez, r1, r2, r3, 
                   r4, r5, r6, r7, r8, 
                   r9, wavelength, intensity, intensityscale, 
                   ray1, index, opticallength, phaseshift, 
                   componentnumber, surfacenumber, surfaceid, ray2, 
                   ray3, sid, frontid, sno1, sno2, 
                   sno3, sno4, rayslot, generationnumber, 
                   intersectionnumber, pol1, pol2, pol3}; 
                 reportedrayspresent = True; intersectionspresent = True; deflectedrayspresent = 
                  True; localdummyoutrayspresent = True; (surfnorm = 
                   1.*Re[{{1., -1.*(ray2/Sign[cb1]/(cb1^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5), 
                       -1.*(ray3/Sign[cb1]/(cb1^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5)}, 
                      {ray2/Sign[cb1]/(cb1^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5, 1., 0.}, 
                      {ray3/Sign[cb1]/(cb1^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5, 0., 1.}}]; 
                  dummymatrix = ({{(Re[Cos[0.017453292519943295*ab]]*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2])/
                         (Re[Cos[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2] + Re[Sin[0.017453292519943295*ab]]^2*
                           Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]), (-1.*Re[Sin[0.017453292519943295*ab]]*
                          Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2])/(Re[Cos[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^
                              2 + Sin[0.017453292519943295*ab]^2] + Re[Sin[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*
                                ab]^2]), 0.}, {(Re[Sin[0.017453292519943295*ab]]*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2])/
                         (Re[Cos[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2] + Re[Sin[0.017453292519943295*ab]]^2*
                           Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]), (Re[Cos[0.017453292519943295*ab]]*Re[Cos[0.017453292519943295*ab]^
                             2 + Sin[0.017453292519943295*ab]^2])/(Re[Cos[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*
                                ab]^2] + Re[Sin[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]), 0.}, 
                       {0., 0., (Re[Cos[0.017453292519943295*ab]]^2 + Re[Sin[0.017453292519943295*ab]]^2)/(Re[Cos[0.017453292519943295*ab]]^2*
                           Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2] + Re[Sin[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^
                              2 + Sin[0.017453292519943295*ab]^2])}} . #1 & ) /@ surfnorm; surfnorm = 
                   (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; surfacerotationmatrix = {{Re[Cos[0.017453292519943295*ab]], 
                     Re[Sin[0.017453292519943295*ab]], 0.}, {-1.*Re[Sin[0.017453292519943295*ab]], Re[Cos[0.017453292519943295*ab]], 0.}, 
                    {0., 0., Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]}}; surfacetranslationvector = 
                   {Re[xb1], Re[yb1], 0.}; activeregion = 1.; completedray = {rsx, 
                    rsy, rsz, rex, rey, rez, 
                    r1, r2, r3, r4, r5, 
                    r6, r7, r8, r9, wavelength, 
                    intensity, intensityscale, ray1, index, 
                    opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; componentnumber = 
                   newcomponentnumber; surfacenumber = newsurfacenumber; 
                  surfaceid = newsurfaceid; originalray = {rsx, 
                    rsy, rsz, rex, rey, rez, 
                    r1, r2, r3, r4, r5, 
                    r6, r7, r8, r9, wavelength, 
                    intensity, intensityscale, ray1, index, 
                    opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; ZAxis = 
                   First[surfnorm]*Sign[First[surfnorm] . {r1, r2, 
                       r3}]; raytiltz = ZAxis . {r1, 
                     r2, r3}; raytiltx = 1 - raytiltz^2; 
                  If[raytiltx < 2.220446049250313*^-16, raytiltx = 0., raytiltx = 
                    Sqrt[raytiltx]]; XAxis = {r1, r2, 
                     r3} - ZAxis*raytiltz; If[Abs[raytiltx] < 1.*^-7, 
                   XAxis = surfnorm[[2]], XAxis = XAxis/
                     raytiltx]; YAxis = {-(ZAxis[[3]]*XAxis[[2]]) + 
                     ZAxis[[2]]*XAxis[[3]], ZAxis[[3]]*XAxis[[1]] - 
                     ZAxis[[1]]*XAxis[[3]], -(ZAxis[[2]]*XAxis[[1]]) + 
                     ZAxis[[1]]*XAxis[[2]]}; YAxis = YAxis/
                    Sqrt[YAxis . YAxis]; PAxis = ZAxis*
                     {r1, r2, r3} . XAxis - XAxis*
                     {r1, r2, r3} . ZAxis; PAxis = 
                   PAxis/Sqrt[PAxis . PAxis]; normalizedpolarization = 
                   {pol1, pol2, pol3}[[{1, 2}]]; normalizedpolarization = 
                   normalizedpolarization/Sqrt[normalizedpolarization . normalizedpolarization]; 
                  allcreatedamplitudes = {{0.*I, 0.*I, 0.*I}}; amplitude = Sqrt[intensity]*
                    (normalizedpolarization[[1]]*{r4, r5, r6} + 
                     normalizedpolarization[[2]]*{r7, r8, r9}*
                      Exp[I*pol3*0.017453292519943295]); perpendicularamplitude = amplitude . 
                    YAxis; parallelamplitude = amplitude . PAxis; 
                  If[activeregion*(surfacerotationmatrix . {r1, r2, 
                        r3})[[1]] > 0, index2 = Sqrt[3.0539614 - 0.011580432*wavelength^2 + 
                       0.039199816/wavelength^2 + 0.0029462812/wavelength^4 - 0.00020371019000000004/
                        wavelength^6 + 0.000027633568999999998/wavelength^8]; index1 = 
                     index; If[Abs[index2 - index1] < 0.01, alternateindex = 
                       Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/wavelength^
                           2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/wavelength^6 + 
                         2.830549017940716*^-12/wavelength^8]; If[Abs[alternateindex - index1] > 
                        Abs[index2 - index1], index1 = index2; 
                        index2 = alternateindex, index1 = alternateindex], 
                     index1 = Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/
                         wavelength^2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/
                         wavelength^6 + 2.830549017940716*^-12/wavelength^8]], index2 = 
                     Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/wavelength^
                         2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/wavelength^6 + 
                       2.830549017940716*^-12/wavelength^8]; index1 = index; 
                    If[Abs[index2 - index1] < 0.01, alternateindex = Sqrt[3.0539614 - 
                         0.011580432*wavelength^2 + 0.039199816/wavelength^2 + 0.0029462812/
                          wavelength^4 - 0.00020371019000000004/wavelength^6 + 0.000027633568999999998/
                          wavelength^8]; If[Abs[alternateindex - index1] > 
                        Abs[index2 - index1], index1 = index2; 
                        index2 = alternateindex, index1 = alternateindex], 
                     index1 = Sqrt[3.0539614 - 0.011580432*wavelength^2 + 0.039199816/wavelength^
                          2 + 0.0029462812/wavelength^4 - 0.00020371019000000004/wavelength^6 + 0.000027633568999999998/
                         wavelength^8]]]; n12 = index1/index2; 
                  index = index2; If[1 - n12^2*raytiltx^2 >= 0., 
                   dummyvector = n12*{r1, r2, r3} + 
                      ZAxis*(-(n12*raytiltz) + Sqrt[1 - n12^2*
                           raytiltx^2]); {r1, r2, r3} = 
                     dummyvector; 1., If[{r1, r2, r3} . 
                       surfnorm[[1]] > 0, phaseshift = phaseshift + 0.5, 1.]; 
                    index2 = index1; index = index1; 
                    dummyvector = {r1, r2, r3} - 
                      2*ZAxis*ZAxis . {r1, r2, r3}; 
                    {r1, r2, r3} = dummyvector; 
                    dummyvector = {pol1, pol2, pol3} + {0, 0, -180}; 
                    {pol1, pol2, pol3} = dummyvector; 1.]; 
                  allcreatedrays = {{rsx, rsy, rsz, 
                     rex, rey, rez, r1, r2, 
                     r3, r4, r5, r6, r7, 
                     r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     generationnumber, intersectionnumber, pol1, pol2, 
                     pol3}}; originalrotation = Partition[Take[originalray, {7, 15}], 3]; 
                  previoustilt = originalrotation[[1]]; transmatrix = 
                   MapIndexed[({rsx, rsy, rsz, rex, 
                        rey, rez, r1, r2, r3, 
                        r4, r5, r6, r7, r8, 
                        r9, wavelength, intensity, intensityscale, 
                        ray1, index, opticallength, phaseshift, 
                        componentnumber, surfacenumber, surfaceid, ray2, 
                        ray3, sid, frontid, sno1, sno2, 
                        sno3, sno4, rayslot, oldgenerationnumber, 
                        intersectionnumber, pol1, pol2, pol3} = #1; 
                      rotationaxis = Cross[{1., 0., 0.}, {r1, r2, r3}]; 
                      axismagnitude = Sqrt[rotationaxis . rotationaxis]; 
                      If[Chop[axismagnitude] != 0., rotationaxis = rotationaxis/
                          axismagnitude; axismatrix = (#1*rotationaxis & ) /@ 
                          rotationaxis; crossmatrix = {{0., -rotationaxis[[3]], 
                           rotationaxis[[2]]}, {rotationaxis[[3]], 0., -rotationaxis[[1]]}, 
                          {-rotationaxis[[2]], rotationaxis[[1]], 0.}}; cosangle = 
                         {1., 0., 0.} . {r1, r2, r3}; If[1. - cosangle < 
                          2.220446049250313*^-16, cosangle = 1.]; sinangle = 1 - cosangle^2; 
                        If[sinangle < 2.220446049250313*^-16, sinangle = 0., sinangle = 
                          Sqrt[sinangle]]; axismatrix + ({{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}} - 
                           axismatrix)*cosangle + crossmatrix*sinangle, 
                       {{r1, r2, r3}, {r4, r5, 
                          r6}, {r7, r8, r9}} = 
                        If[Sign[{1., 0., 0.} . {r1, r2, r3}] === 1, {{1., 0., 0.}, {0., 1., 0.}, 
                          {0., 0., 1.}}, {{-1., 0., 0.}, {0., -1., 0.}, {0., 0., 1.}}]]) & , allcreatedrays]; 
                  allcreatedrays = MapIndexed[({rsx, rsy, rsz, 
                        rex, rey, rez, r1, r2, 
                        r3, r4, r5, r6, r7, 
                        r8, r9, wavelength, intensity, 
                        intensityscale, ray1, index, opticallength, 
                        phaseshift, componentnumber, surfacenumber, 
                        surfaceid, ray2, ray3, sid, 
                        frontid, sno1, sno2, sno3, 
                        sno4, rayslot, oldgenerationnumber, 
                        intersectionnumber, pol1, pol2, pol3} = #1; 
                      {{r1, r2, r3}, {r4, r5, 
                         r6}, {r7, r8, r9}} = 
                       Function[mm, (mm . #1 & ) /@ {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}][
                        transmatrix[[First[#2]]]]; tiltmagnitude = {r1, r2, 
                         r3} . {r1, r2, r3}; 
                      If[Abs[tiltmagnitude - 1.] > 1.*^-6, dummyvector = {r1, 
                           r2, r3}/Sqrt[tiltmagnitude]; {r1, 
                          r2, r3} = dummyvector; dummymatrix = 
                         (#1/Sqrt[#1 . #1] & ) /@ {{r1, r2, r3}, {r4, 
                            r5, r6}, {r7, r8, r9}}; 
                        {{r1, r2, r3}, {r4, r5, 
                           r6}, {r7, r8, r9}} = 
                         dummymatrix; 1., 1.]; If[allcreatedamplitudes != {{0.*I, 0.*I, 0.*I}}, 
                       amplitude = allcreatedamplitudes[[First[#2]]]; hor = 
                         {r4, r5, r6} . amplitude; 
                        ver = {r7, r8, r9} . amplitude; 
                        If[Abs[hor] > 0. || Abs[ver] > 0., {pol1, pol2, 
                            pol3} = {Abs[hor], Abs[ver], 0.}; dummyvector = 
                           {pol1, pol2, pol3}/Sqrt[{pol1, 
                               pol2, pol3} . {pol1, pol2, 
                               pol3}]; {pol1, pol2, pol3} = 
                           dummyvector; dummyvector = ReplacePart[{pol1, 
                             pol2, pol3}, (If[ver == 0., 0., Arg[ver]] - 
                              If[hor == 0., 0., Arg[hor]])/Degree, 3]; {pol1, 
                            pol2, pol3} = dummyvector]; 1., 1.]; {rsx, 
                       rsy, rsz, rex, rey, rez, 
                       r1, r2, r3, r4, r5, 
                       r6, r7, r8, r9, wavelength, 
                       intensity, intensityscale, ray1, index, 
                       opticallength, phaseshift, componentnumber, 
                       surfacenumber, surfaceid, ray2, ray3, 
                       sid, frontid, sno1, sno2, sno3, 
                       sno4, rayslot, generationnumber, intersectionnumber, 
                       pol1, pol2, pol3}) & , allcreatedrays]); 
                 {allcreatedrays, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 
                    0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, {completedray}}, 
                resonatenumber == 4, {rsx, rsy, rsz, 
                   rex, rey, rez, r1, r2, 
                   r3, r4, r5, r6, r7, 
                   r8, r9, wavelength, intensity, 
                   intensityscale, ray1, index, opticallength, 
                   phaseshift, componentnumber, surfacenumber, surfaceid, 
                   ray2, ray3, sid, frontid, sno1, 
                   sno2, sno3, sno4, rayslot, 
                   oldgenerationnumber, intersectionnumber, pol1, pol2, 
                   pol3} = #1; index1 = index; xunit = 
                  r1*Re[Cos[0.017453292519943295*ab]] + r2*Re[Sin[0.017453292519943295*ab]]; 
                 yunit = r2*Re[Cos[0.017453292519943295*ab]] - 1.*r1*
                    Re[Sin[0.017453292519943295*ab]]; zunit = r3*Re[Cos[0.017453292519943295*ab]^2 + 
                     Sin[0.017453292519943295*ab]^2]; xunit = Chop[xunit]; yunit = 
                  Chop[yunit]; zunit = Chop[zunit]; magnitude = 
                  (xunit^2 + yunit^2 + zunit^2)^2^(-1); xunit = 
                  xunit/magnitude; yunit = yunit/
                   magnitude; zunit = zunit/magnitude; 
                 rx = Re[Cos[0.017453292519943295*ab]]*(rex - 1.*Re[xb1 + tb*Cos[0.017453292519943295*ab]]) + 
                   Re[Sin[0.017453292519943295*ab]]*(rey - 1.*Re[yb1 + tb*Sin[0.017453292519943295*ab]]); ry = 
                  -1.*(rex - 1.*Re[xb1 + tb*Cos[0.017453292519943295*ab]])*Re[Sin[0.017453292519943295*ab]] + Re[Cos[0.017453292519943295*ab]]*
                    (rey - 1.*Re[yb1 + tb*Sin[0.017453292519943295*ab]]); rz = rez*
                   Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]; rx = Chop[rx]; 
                 ry = Chop[ry]; rz = Chop[rz]; 
                 terminatedeflection = False; intersections = {{-1., 0., 0.}}; intersections = 
                  resonateintersections; holeQ = holeQlist[[
                    resonatenumber]] == 1; {newcomponentnumber, newsurfacenumber, 
                   newsurfaceid} = {1., 4., 16548.}; {componentnumber, surfacenumber, 
                   surfaceid} = {1., 4., 16548.}; {ray1, ray2, ray3} = 
                  intersections[[1]]; rsx = rex; rsy = 
                  rey; rsz = rez; rxhlp = rx + 
                   ray1*xunit; ryhlp = ry + ray1*
                    yunit; rzhlp = rz + ray1*zunit; 
                 rex = rxhlp*Re[Cos[0.017453292519943295*ab]] + Re[xb1 + tb*Cos[0.017453292519943295*ab]] - 
                   1.*ryhlp*Re[Sin[0.017453292519943295*ab]]; rey = ryhlp*
                    Re[Cos[0.017453292519943295*ab]] + rxhlp*Re[Sin[0.017453292519943295*ab]] + Re[yb1 + tb*Sin[0.017453292519943295*ab]]; 
                 rez = rzhlp*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]; 
                 rex = Chop[rex]; rey = Chop[rey]; 
                 rez = Chop[rez]; If[rex == 0 && rey == 0 && 
                   rez == 0, rex = rxhlp + Re[xb1 + tb*Cos[0.017453292519943295*ab]]; 
                   rey = ryhlp + Re[yb1 + tb*Sin[0.017453292519943295*ab]]; rez = 
                    rzhlp, rex]; intersectionnumber = intersectionnumber + 
                   1.; opticallength = opticallength + ray1*index1; 
                 completedray = {rsx, rsy, rsz, rex, 
                   rey, rez, r1, r2, r3, 
                   r4, r5, r6, r7, r8, 
                   r9, wavelength, intensity, intensityscale, 
                   ray1, index, opticallength, phaseshift, 
                   componentnumber, surfacenumber, surfaceid, ray2, 
                   ray3, sid, frontid, sno1, sno2, 
                   sno3, sno4, rayslot, generationnumber, 
                   intersectionnumber, pol1, pol2, pol3}; 
                 reportedrayspresent = True; intersectionspresent = True; deflectedrayspresent = 
                  True; localdummyoutrayspresent = True; (surfnorm = 
                   1.*Re[{{1., -1.*(ray2/Sign[cb2]/(cb2^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5), 
                       -1.*(ray3/Sign[cb2]/(cb2^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5)}, 
                      {ray2/Sign[cb2]/(cb2^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5, 1., 0.}, 
                      {ray3/Sign[cb2]/(cb2^(-2) - 1.*ray2^2 - 1.*ray3^2)^0.5, 0., 1.}}]; 
                  dummymatrix = ({{(Re[Cos[0.017453292519943295*ab]]*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2])/
                         (Re[Cos[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2] + Re[Sin[0.017453292519943295*ab]]^2*
                           Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]), (-1.*Re[Sin[0.017453292519943295*ab]]*
                          Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2])/(Re[Cos[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^
                              2 + Sin[0.017453292519943295*ab]^2] + Re[Sin[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*
                                ab]^2]), 0.}, {(Re[Sin[0.017453292519943295*ab]]*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2])/
                         (Re[Cos[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2] + Re[Sin[0.017453292519943295*ab]]^2*
                           Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]), (Re[Cos[0.017453292519943295*ab]]*Re[Cos[0.017453292519943295*ab]^
                             2 + Sin[0.017453292519943295*ab]^2])/(Re[Cos[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*
                                ab]^2] + Re[Sin[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]), 0.}, 
                       {0., 0., (Re[Cos[0.017453292519943295*ab]]^2 + Re[Sin[0.017453292519943295*ab]]^2)/(Re[Cos[0.017453292519943295*ab]]^2*
                           Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2] + Re[Sin[0.017453292519943295*ab]]^2*Re[Cos[0.017453292519943295*ab]^
                              2 + Sin[0.017453292519943295*ab]^2])}} . #1 & ) /@ surfnorm; surfnorm = 
                   (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; surfacerotationmatrix = {{Re[Cos[0.017453292519943295*ab]], 
                     Re[Sin[0.017453292519943295*ab]], 0.}, {-1.*Re[Sin[0.017453292519943295*ab]], Re[Cos[0.017453292519943295*ab]], 0.}, 
                    {0., 0., Re[Cos[0.017453292519943295*ab]^2 + Sin[0.017453292519943295*ab]^2]}}; surfacetranslationvector = 
                   {Re[xb1 + tb*Cos[0.017453292519943295*ab]], Re[yb1 + tb*Sin[0.017453292519943295*ab]], 0.}; activeregion = -1.; 
                  completedray = {rsx, rsy, rsz, rex, 
                    rey, rez, r1, r2, r3, 
                    r4, r5, r6, r7, r8, 
                    r9, wavelength, intensity, intensityscale, 
                    ray1, index, opticallength, phaseshift, 
                    componentnumber, surfacenumber, surfaceid, ray2, 
                    ray3, sid, frontid, sno1, sno2, 
                    sno3, sno4, rayslot, generationnumber, 
                    intersectionnumber, pol1, pol2, pol3}; 
                  componentnumber = newcomponentnumber; surfacenumber = 
                   newsurfacenumber; surfaceid = newsurfaceid; originalray = 
                   {rsx, rsy, rsz, rex, rey, 
                    rez, r1, r2, r3, r4, 
                    r5, r6, r7, r8, r9, 
                    wavelength, intensity, intensityscale, ray1, 
                    index, opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; ZAxis = 
                   First[surfnorm]*Sign[First[surfnorm] . {r1, r2, 
                       r3}]; raytiltz = ZAxis . {r1, 
                     r2, r3}; raytiltx = 1 - raytiltz^2; 
                  If[raytiltx < 2.220446049250313*^-16, raytiltx = 0., raytiltx = 
                    Sqrt[raytiltx]]; XAxis = {r1, r2, 
                     r3} - ZAxis*raytiltz; If[Abs[raytiltx] < 1.*^-7, 
                   XAxis = surfnorm[[2]], XAxis = XAxis/
                     raytiltx]; YAxis = {-(ZAxis[[3]]*XAxis[[2]]) + 
                     ZAxis[[2]]*XAxis[[3]], ZAxis[[3]]*XAxis[[1]] - 
                     ZAxis[[1]]*XAxis[[3]], -(ZAxis[[2]]*XAxis[[1]]) + 
                     ZAxis[[1]]*XAxis[[2]]}; YAxis = YAxis/
                    Sqrt[YAxis . YAxis]; PAxis = ZAxis*
                     {r1, r2, r3} . XAxis - XAxis*
                     {r1, r2, r3} . ZAxis; PAxis = 
                   PAxis/Sqrt[PAxis . PAxis]; normalizedpolarization = 
                   {pol1, pol2, pol3}[[{1, 2}]]; normalizedpolarization = 
                   normalizedpolarization/Sqrt[normalizedpolarization . normalizedpolarization]; 
                  allcreatedamplitudes = {{0.*I, 0.*I, 0.*I}}; amplitude = Sqrt[intensity]*
                    (normalizedpolarization[[1]]*{r4, r5, r6} + 
                     normalizedpolarization[[2]]*{r7, r8, r9}*
                      Exp[I*pol3*0.017453292519943295]); perpendicularamplitude = amplitude . 
                    YAxis; parallelamplitude = amplitude . PAxis; 
                  If[activeregion*(surfacerotationmatrix . {r1, r2, 
                        r3})[[1]] > 0, index2 = Sqrt[3.0539614 - 0.011580432*wavelength^2 + 
                       0.039199816/wavelength^2 + 0.0029462812/wavelength^4 - 0.00020371019000000004/
                        wavelength^6 + 0.000027633568999999998/wavelength^8]; index1 = 
                     index; If[Abs[index2 - index1] < 0.01, alternateindex = 
                       Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/wavelength^
                           2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/wavelength^6 + 
                         2.830549017940716*^-12/wavelength^8]; If[Abs[alternateindex - index1] > 
                        Abs[index2 - index1], index1 = index2; 
                        index2 = alternateindex, index1 = alternateindex], 
                     index1 = Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/
                         wavelength^2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/
                         wavelength^6 + 2.830549017940716*^-12/wavelength^8]], index2 = 
                     Sqrt[1.0005281111748805 - 2.364805676582992*^-10*wavelength^2 + 2.9781670902982902*^-6/wavelength^
                         2 + 2.6004099469477187*^-8/wavelength^4 + 3.7787839314048715*^-10/wavelength^6 + 
                       2.830549017940716*^-12/wavelength^8]; index1 = index; 
                    If[Abs[index2 - index1] < 0.01, alternateindex = Sqrt[3.0539614 - 
                         0.011580432*wavelength^2 + 0.039199816/wavelength^2 + 0.0029462812/
                          wavelength^4 - 0.00020371019000000004/wavelength^6 + 0.000027633568999999998/
                          wavelength^8]; If[Abs[alternateindex - index1] > 
                        Abs[index2 - index1], index1 = index2; 
                        index2 = alternateindex, index1 = alternateindex], 
                     index1 = Sqrt[3.0539614 - 0.011580432*wavelength^2 + 0.039199816/wavelength^
                          2 + 0.0029462812/wavelength^4 - 0.00020371019000000004/wavelength^6 + 0.000027633568999999998/
                         wavelength^8]]]; n12 = index1/index2; 
                  index = index2; If[1 - n12^2*raytiltx^2 >= 0., 
                   dummyvector = n12*{r1, r2, r3} + 
                      ZAxis*(-(n12*raytiltz) + Sqrt[1 - n12^2*
                           raytiltx^2]); {r1, r2, r3} = 
                     dummyvector; 1., If[{r1, r2, r3} . 
                       surfnorm[[1]] > 0, phaseshift = phaseshift + 0.5, 1.]; 
                    index2 = index1; index = index1; 
                    dummyvector = {r1, r2, r3} - 
                      2*ZAxis*ZAxis . {r1, r2, r3}; 
                    {r1, r2, r3} = dummyvector; 
                    dummyvector = {pol1, pol2, pol3} + {0, 0, -180}; 
                    {pol1, pol2, pol3} = dummyvector; 1.]; 
                  allcreatedrays = {{rsx, rsy, rsz, 
                     rex, rey, rez, r1, r2, 
                     r3, r4, r5, r6, r7, 
                     r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     generationnumber, intersectionnumber, pol1, pol2, 
                     pol3}}; originalrotation = Partition[Take[originalray, {7, 15}], 3]; 
                  previoustilt = originalrotation[[1]]; transmatrix = 
                   MapIndexed[({rsx, rsy, rsz, rex, 
                        rey, rez, r1, r2, r3, 
                        r4, r5, r6, r7, r8, 
                        r9, wavelength, intensity, intensityscale, 
                        ray1, index, opticallength, phaseshift, 
                        componentnumber, surfacenumber, surfaceid, ray2, 
                        ray3, sid, frontid, sno1, sno2, 
                        sno3, sno4, rayslot, oldgenerationnumber, 
                        intersectionnumber, pol1, pol2, pol3} = #1; 
                      rotationaxis = Cross[{1., 0., 0.}, {r1, r2, r3}]; 
                      axismagnitude = Sqrt[rotationaxis . rotationaxis]; 
                      If[Chop[axismagnitude] != 0., rotationaxis = rotationaxis/
                          axismagnitude; axismatrix = (#1*rotationaxis & ) /@ 
                          rotationaxis; crossmatrix = {{0., -rotationaxis[[3]], 
                           rotationaxis[[2]]}, {rotationaxis[[3]], 0., -rotationaxis[[1]]}, 
                          {-rotationaxis[[2]], rotationaxis[[1]], 0.}}; cosangle = 
                         {1., 0., 0.} . {r1, r2, r3}; If[1. - cosangle < 
                          2.220446049250313*^-16, cosangle = 1.]; sinangle = 1 - cosangle^2; 
                        If[sinangle < 2.220446049250313*^-16, sinangle = 0., sinangle = 
                          Sqrt[sinangle]]; axismatrix + ({{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}} - 
                           axismatrix)*cosangle + crossmatrix*sinangle, 
                       {{r1, r2, r3}, {r4, r5, 
                          r6}, {r7, r8, r9}} = 
                        If[Sign[{1., 0., 0.} . {r1, r2, r3}] === 1, {{1., 0., 0.}, {0., 1., 0.}, 
                          {0., 0., 1.}}, {{-1., 0., 0.}, {0., -1., 0.}, {0., 0., 1.}}]]) & , allcreatedrays]; 
                  allcreatedrays = MapIndexed[({rsx, rsy, rsz, 
                        rex, rey, rez, r1, r2, 
                        r3, r4, r5, r6, r7, 
                        r8, r9, wavelength, intensity, 
                        intensityscale, ray1, index, opticallength, 
                        phaseshift, componentnumber, surfacenumber, 
                        surfaceid, ray2, ray3, sid, 
                        frontid, sno1, sno2, sno3, 
                        sno4, rayslot, oldgenerationnumber, 
                        intersectionnumber, pol1, pol2, pol3} = #1; 
                      {{r1, r2, r3}, {r4, r5, 
                         r6}, {r7, r8, r9}} = 
                       Function[mm, (mm . #1 & ) /@ {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}][
                        transmatrix[[First[#2]]]]; tiltmagnitude = {r1, r2, 
                         r3} . {r1, r2, r3}; 
                      If[Abs[tiltmagnitude - 1.] > 1.*^-6, dummyvector = {r1, 
                           r2, r3}/Sqrt[tiltmagnitude]; {r1, 
                          r2, r3} = dummyvector; dummymatrix = 
                         (#1/Sqrt[#1 . #1] & ) /@ {{r1, r2, r3}, {r4, 
                            r5, r6}, {r7, r8, r9}}; 
                        {{r1, r2, r3}, {r4, r5, 
                           r6}, {r7, r8, r9}} = 
                         dummymatrix; 1., 1.]; If[allcreatedamplitudes != {{0.*I, 0.*I, 0.*I}}, 
                       amplitude = allcreatedamplitudes[[First[#2]]]; hor = 
                         {r4, r5, r6} . amplitude; 
                        ver = {r7, r8, r9} . amplitude; 
                        If[Abs[hor] > 0. || Abs[ver] > 0., {pol1, pol2, 
                            pol3} = {Abs[hor], Abs[ver], 0.}; dummyvector = 
                           {pol1, pol2, pol3}/Sqrt[{pol1, 
                               pol2, pol3} . {pol1, pol2, 
                               pol3}]; {pol1, pol2, pol3} = 
                           dummyvector; dummyvector = ReplacePart[{pol1, 
                             pol2, pol3}, (If[ver == 0., 0., Arg[ver]] - 
                              If[hor == 0., 0., Arg[hor]])/Degree, 3]; {pol1, 
                            pol2, pol3} = dummyvector]; 1., 1.]; {rsx, 
                       rsy, rsz, rex, rey, rez, 
                       r1, r2, r3, r4, r5, 
                       r6, r7, r8, r9, wavelength, 
                       intensity, intensityscale, ray1, index, 
                       opticallength, phaseshift, componentnumber, 
                       surfacenumber, surfaceid, ray2, ray3, 
                       sid, frontid, sno1, sno2, sno3, 
                       sno4, rayslot, generationnumber, intersectionnumber, 
                       pol1, pol2, pol3}) & , allcreatedrays]); 
                 {allcreatedrays, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 
                    0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, {completedray}}, True, 
                deflecteddummyrayspresent = True; reporteddummyrayspresent = True; 
                 localoutrayspresent = True; {{{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 
                    1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, {{rsx, rsy, 
                    rsz, rex, rey, rez, r1, 
                    r2, r3, r4, r5, r6, 
                    r7, r8, r9, wavelength, intensity, 
                    intensityscale, ray1, index, opticallength, 
                    phaseshift, componentnumber, surfacenumber, surfaceid, 
                    ray2, ray3, sid, frontid, sno1, 
                    sno2, sno3, sno4, rayslot, 
                    generationnumber, intersectionnumber, pol1, pol2, 
                    pol3}}, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 0., 
                    -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}}], deflecteddummyrayspresent = True; 
               reporteddummyrayspresent = True; localdummyoutrayspresent = True; {{{1., 1., 1., 1., 1., 1., 1., 1., 1., 
                  1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, 
                {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 
                  1., 1., 1., 1., 1., 1., 1., 1., 1.}}, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 
                  0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}}]) & ) /@ inrays]; 
       If[deflectedrayspresent && deflecteddummyrayspresent, deflectedlocalrays = 
         Select[deflectedlocalrays, #1[[23]] > 0. & ]*1., 1.]; If[localoutrayspresent, 
        If[localdummyoutrayspresent, localoutrays = Select[localoutrays, #1[[23]] > 0. & ]; 1., 
          1.]; outrays = Join[outrays, localoutrays]; 1., 1.]; 
       If[reportedrayspresent && reporteddummyrayspresent, 
        reportedlocalrays = Select[reportedlocalrays, #1[[23]] > 0. & ]; 1., 1.]); 
      If[reportedrayspresent, raysreported = True; reportedrays = 
         Join[reportedrays, {{-1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308}}, 
          reportedlocalrays]; 1., 1.]; inrays = deflectedlocalrays; 
      propagate = generationnumber < generationlimit]; 
    propagatedrays = If[Length[outrays] > 1, Drop[outrays, 1], outrays]; 
    If[generationnumber >= generationlimit && deflectedrayspresent, 
     propagatedrays = Join[propagatedrays, deflectedlocalrays]; 1., 1.]); 
   (inrays = propagatedrays; outrays = 
     {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 
       1., 1., 1., 1., 1., 1.}}; deflectedrayspresent = True; reportedrayspresent = False; 
    propagate = True; While[deflectedrayspresent && propagate, 
     generationnumber += 1.; (deflecteddummyrayspresent = False; deflectedrayspresent = False; 
       reportedrayspresent = False; reporteddummyrayspresent = True; localdummyoutrayspresent = 
        False; localoutrayspresent = False; {deflectedlocalrays, localoutrays, 
         reportedlocalrays} = (Flatten[#1, 1] & ) /@ 
         Transpose[(({rsx, rsy, rsz, rex, rey, 
               rez, r1, r2, r3, r4, 
               r5, r6, r7, r8, r9, 
               wavelength, intensity, intensityscale, ray1, 
               index, opticallength, phaseshift, componentnumber, 
               surfacenumber, surfaceid, ray2, ray3, 
               sid, frontid, sno1, sno2, sno3, 
               sno4, rayslot, oldgenerationnumber, intersectionnumber, 
               pol1, pol2, pol3} = #1; If[intensity > 
                thresholdintensity && generationnumber < generationlimit, 
              (resonatenumber = 0; resonateraylength = 1.7976931348623157*^308; 
                resonateintersections = {{-1., 0., 0.}}; MapIndexed[If[#1[[1,1]] > -1. && #1[[1,1]] < resonateraylength, 
                   resonateraylength = #1[[1,1]]; resonateintersections = #1; resonatenumber = 
                     First[#2], resonatenumber] & , intersectionslist = {{{-1., 0., 0.}}}; 
                  holeQlist = {0}; ({rsx, rsy, rsz, 
                     rex, rey, rez, r1, r2, 
                     r3, r4, r5, r6, r7, 
                     r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     oldgenerationnumber, intersectionnumber, pol1, pol2, 
                     pol3} = #1; index1 = index; xunit = 
                    -1.*r1; yunit = -1.*r2; zunit = r3; 
                   xunit = Chop[xunit]; yunit = Chop[yunit]; 
                   zunit = Chop[zunit]; magnitude = (xunit^2 + 
                      yunit^2 + zunit^2)^2^(-1); xunit = xunit/
                     magnitude; yunit = yunit/magnitude; 
                   zunit = zunit/magnitude; rx = 
                    -1.*(-200. + rex); ry = -1.*rey; rz = 
                    rez; rx = Chop[rx]; ry = Chop[ry]; 
                   rz = Chop[rz]; terminatedeflection = False; 
                   intersections = {{-1., 0., 0.}}; intersections = If[Chop[xunit] != 0, 
                     {(If[N[Chop[#1[[1]]]] >= 0. && (#1[[2]] <= 75. && #1[[2]] >= -75. && #1[[3]] >= -75. && #1[[3]] <= 75.), {#1[[1]], #1[[2]], #1[[3]]}, 
                         {-1., 0., 0.}] & )[{-(rx/xunit), ry - (rx*
                           yunit)/xunit, rz - (rx*
                           zunit)/xunit}]}, {{-1., 0., 0.}}]; AppendTo[intersections, 
                    {-1., 0., 0.}]; holeQ = False; If[ !(Chop[surfaceid - 16598.] == 0. || 
                       N[Chop[ray1]] < 0.0001) || generationnumber == 1., minresonateraylength = 
                     0., minresonateraylength = 0.0001]; Block[{RayL, SCx, SCy}, RayL = 1.; SCx = 1.; SCy = 1.; 
                     dummymatrix = Union[Select[intersections, ({RayL, SCx, SCy} = #1; {RayL, SCx, SCy} == 
                            {-1., 0., 0.} || N[Chop[RayL]] >= minresonateraylength) & ]]; intersections = 
                      dummymatrix]; unfilteredintersections = intersections; 
                   If[Length[intersections] > 1, dummymatrix = Take[intersections, {2}]; 
                     intersections = dummymatrix; 1., 1.]; AppendTo[intersectionslist, 
                    intersections]; If[holeQ, AppendTo[holeQlist, 1], 
                    AppendTo[holeQlist, 0]]; {rsx, rsy, rsz, 
                     rex, rey, rez, r1, r2, 
                     r3, r4, r5, r6, r7, 
                     r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     oldgenerationnumber, intersectionnumber, pol1, pol2, 
                     pol3} = #1; index1 = index; xunit = 
                    r2; yunit = -1.*r1; zunit = r3; 
                   xunit = Chop[xunit]; yunit = Chop[yunit]; 
                   zunit = Chop[zunit]; magnitude = (xunit^2 + 
                      yunit^2 + zunit^2)^2^(-1); xunit = xunit/
                     magnitude; yunit = yunit/magnitude; 
                   zunit = zunit/magnitude; rx = 
                    75. + rey; ry = -1.*(-100. + rex); rz = 
                    rez; rx = Chop[rx]; ry = Chop[ry]; 
                   rz = Chop[rz]; terminatedeflection = False; 
                   intersections = {{-1., 0., 0.}}; intersections = If[Chop[xunit] != 0, 
                     {(If[N[Chop[#1[[1]]]] >= 0. && (#1[[2]] <= 100. && #1[[2]] >= -100. && #1[[3]] >= -75. && #1[[3]] <= 75.), {#1[[1]], #1[[2]], #1[[3]]}, 
                         {-1., 0., 0.}] & )[{-(rx/xunit), ry - (rx*
                           yunit)/xunit, rz - (rx*
                           zunit)/xunit}]}, {{-1., 0., 0.}}]; AppendTo[intersections, 
                    {-1., 0., 0.}]; holeQ = False; If[ !(Chop[surfaceid - 16599.] == 0. || 
                       N[Chop[ray1]] < 0.0001) || generationnumber == 1., minresonateraylength = 
                     0., minresonateraylength = 0.0001]; Block[{RayL, SCx, SCy}, RayL = 1.; SCx = 1.; SCy = 1.; 
                     dummymatrix = Union[Select[intersections, ({RayL, SCx, SCy} = #1; {RayL, SCx, SCy} == 
                            {-1., 0., 0.} || N[Chop[RayL]] >= minresonateraylength) & ]]; intersections = 
                      dummymatrix]; unfilteredintersections = intersections; 
                   If[Length[intersections] > 1, dummymatrix = Take[intersections, {2}]; 
                     intersections = dummymatrix; 1., 1.]; AppendTo[intersectionslist, 
                    intersections]; If[holeQ, AppendTo[holeQlist, 1], 
                    AppendTo[holeQlist, 0]]; {rsx, rsy, rsz, 
                     rex, rey, rez, r1, r2, 
                     r3, r4, r5, r6, r7, 
                     r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     oldgenerationnumber, intersectionnumber, pol1, pol2, 
                     pol3} = #1; index1 = index; xunit = 
                    r1; yunit = r2; zunit = r3; 
                   xunit = Chop[xunit]; yunit = Chop[yunit]; 
                   zunit = Chop[zunit]; magnitude = (xunit^2 + 
                      yunit^2 + zunit^2)^2^(-1); xunit = xunit/
                     magnitude; yunit = yunit/magnitude; 
                   zunit = zunit/magnitude; rx = rex; 
                   ry = rey; rz = rez; rx = 
                    Chop[rx]; ry = Chop[ry]; rz = 
                    Chop[rz]; terminatedeflection = False; intersections = {{-1., 0., 0.}}; 
                   intersections = If[Chop[xunit] != 0, {(If[N[Chop[#1[[1]]]] >= 0. && (#1[[2]] <= 75. && 
                           #1[[2]] >= -75. && #1[[3]] >= -75. && #1[[3]] <= 75.), {#1[[1]], #1[[2]], #1[[3]]}, {-1., 0., 0.}] & )[{-(rx/
                          xunit), ry - (rx*yunit)/
                          xunit, rz - (rx*zunit)/
                          xunit}]}, {{-1., 0., 0.}}]; AppendTo[intersections, {-1., 0., 0.}]; 
                   holeQ = False; If[ !(Chop[surfaceid - 16600.] == 0. || N[Chop[ray1]] < 
                        0.0001) || generationnumber == 1., minresonateraylength = 0., 
                    minresonateraylength = 0.0001]; Block[{RayL, SCx, SCy}, RayL = 1.; SCx = 1.; SCy = 1.; dummymatrix = 
                      Union[Select[intersections, ({RayL, SCx, SCy} = #1; {RayL, SCx, SCy} == {-1., 0., 0.} || N[Chop[RayL]] >= 
                            minresonateraylength) & ]]; intersections = dummymatrix]; 
                   unfilteredintersections = intersections; If[Length[intersections] > 1, 
                    dummymatrix = Take[intersections, {2}]; intersections = 
                      dummymatrix; 1., 1.]; AppendTo[intersectionslist, intersections]; 
                   If[holeQ, AppendTo[holeQlist, 1], AppendTo[holeQlist, 0]]; 
                   {rsx, rsy, rsz, rex, rey, 
                     rez, r1, r2, r3, r4, 
                     r5, r6, r7, r8, r9, 
                     wavelength, intensity, intensityscale, ray1, 
                     index, opticallength, phaseshift, componentnumber, 
                     surfacenumber, surfaceid, ray2, ray3, 
                     sid, frontid, sno1, sno2, sno3, 
                     sno4, rayslot, oldgenerationnumber, intersectionnumber, 
                     pol1, pol2, pol3} = #1; index1 = 
                    index; xunit = -1.*r2; yunit = r1; 
                   zunit = r3; xunit = Chop[xunit]; 
                   yunit = Chop[yunit]; zunit = Chop[zunit]; 
                   magnitude = (xunit^2 + yunit^2 + zunit^2)^2^(-1); 
                   xunit = xunit/magnitude; yunit = 
                    yunit/magnitude; zunit = zunit/
                     magnitude; rx = -1.*(-75. + rey); ry = 
                    -100. + rex; rz = rez; rx = 
                    Chop[rx]; ry = Chop[ry]; rz = 
                    Chop[rz]; terminatedeflection = False; intersections = {{-1., 0., 0.}}; 
                   intersections = If[Chop[xunit] != 0, {(If[N[Chop[#1[[1]]]] >= 0. && (#1[[2]] <= 100. && 
                           #1[[2]] >= -100. && #1[[3]] >= -75. && #1[[3]] <= 75.), {#1[[1]], #1[[2]], #1[[3]]}, {-1., 0., 0.}] & )[{-(rx/
                          xunit), ry - (rx*yunit)/
                          xunit, rz - (rx*zunit)/
                          xunit}]}, {{-1., 0., 0.}}]; AppendTo[intersections, {-1., 0., 0.}]; 
                   holeQ = False; If[ !(Chop[surfaceid - 16601.] == 0. || N[Chop[ray1]] < 
                        0.0001) || generationnumber == 1., minresonateraylength = 0., 
                    minresonateraylength = 0.0001]; Block[{RayL, SCx, SCy}, RayL = 1.; SCx = 1.; SCy = 1.; dummymatrix = 
                      Union[Select[intersections, ({RayL, SCx, SCy} = #1; {RayL, SCx, SCy} == {-1., 0., 0.} || N[Chop[RayL]] >= 
                            minresonateraylength) & ]]; intersections = dummymatrix]; 
                   unfilteredintersections = intersections; If[Length[intersections] > 1, 
                    dummymatrix = Take[intersections, {2}]; intersections = 
                      dummymatrix; 1., 1.]; AppendTo[intersectionslist, intersections]; 
                   If[holeQ, AppendTo[holeQlist, 1], AppendTo[holeQlist, 0]]; 
                   {rsx, rsy, rsz, rex, rey, 
                     rez, r1, r2, r3, r4, 
                     r5, r6, r7, r8, r9, 
                     wavelength, intensity, intensityscale, ray1, 
                     index, opticallength, phaseshift, componentnumber, 
                     surfacenumber, surfaceid, ray2, ray3, 
                     sid, frontid, sno1, sno2, sno3, 
                     sno4, rayslot, oldgenerationnumber, intersectionnumber, 
                     pol1, pol2, pol3} = #1; index1 = 
                    index; xunit = r3; yunit = r2; 
                   zunit = -1.*r1; xunit = Chop[xunit]; 
                   yunit = Chop[yunit]; zunit = Chop[zunit]; 
                   magnitude = (xunit^2 + yunit^2 + zunit^2)^2^(-1); 
                   xunit = xunit/magnitude; yunit = 
                    yunit/magnitude; zunit = zunit/
                     magnitude; rx = 75. + rez; ry = 
                    rey; rz = -1.*(-100. + rex); rx = 
                    Chop[rx]; ry = Chop[ry]; rz = 
                    Chop[rz]; terminatedeflection = False; intersections = {{-1., 0., 0.}}; 
                   intersections = If[Chop[xunit] != 0, {(If[N[Chop[#1[[1]]]] >= 0. && (#1[[2]] <= 75. && 
                           #1[[2]] >= -75. && #1[[3]] >= -100. && #1[[3]] <= 100.), {#1[[1]], #1[[2]], #1[[3]]}, {-1., 0., 0.}] & )[{-(rx/
                          xunit), ry - (rx*yunit)/
                          xunit, rz - (rx*zunit)/
                          xunit}]}, {{-1., 0., 0.}}]; AppendTo[intersections, {-1., 0., 0.}]; 
                   holeQ = False; If[ !(Chop[surfaceid - 16602.] == 0. || N[Chop[ray1]] < 
                        0.0001) || generationnumber == 1., minresonateraylength = 0., 
                    minresonateraylength = 0.0001]; Block[{RayL, SCx, SCy}, RayL = 1.; SCx = 1.; SCy = 1.; dummymatrix = 
                      Union[Select[intersections, ({RayL, SCx, SCy} = #1; {RayL, SCx, SCy} == {-1., 0., 0.} || N[Chop[RayL]] >= 
                            minresonateraylength) & ]]; intersections = dummymatrix]; 
                   unfilteredintersections = intersections; If[Length[intersections] > 1, 
                    dummymatrix = Take[intersections, {2}]; intersections = 
                      dummymatrix; 1., 1.]; AppendTo[intersectionslist, intersections]; 
                   If[holeQ, AppendTo[holeQlist, 1], AppendTo[holeQlist, 0]]; 
                   {rsx, rsy, rsz, rex, rey, 
                     rez, r1, r2, r3, r4, 
                     r5, r6, r7, r8, r9, 
                     wavelength, intensity, intensityscale, ray1, 
                     index, opticallength, phaseshift, componentnumber, 
                     surfacenumber, surfaceid, ray2, ray3, 
                     sid, frontid, sno1, sno2, sno3, 
                     sno4, rayslot, oldgenerationnumber, intersectionnumber, 
                     pol1, pol2, pol3} = #1; index1 = 
                    index; xunit = -1.*r3; yunit = 
                    -1.*r2; zunit = -1.*r1; xunit = 
                    Chop[xunit]; yunit = Chop[yunit]; zunit = 
                    Chop[zunit]; magnitude = (xunit^2 + yunit^2 + 
                      zunit^2)^2^(-1); xunit = xunit/magnitude; 
                   yunit = yunit/magnitude; zunit = 
                    zunit/magnitude; rx = -1.*(-75. + rez); 
                   ry = -1.*rey; rz = -1.*(-100. + rex); 
                   rx = Chop[rx]; ry = Chop[ry]; rz = 
                    Chop[rz]; terminatedeflection = False; intersections = {{-1., 0., 0.}}; 
                   intersections = If[Chop[xunit] != 0, {(If[N[Chop[#1[[1]]]] >= 0. && (#1[[2]] <= 75. && 
                           #1[[2]] >= -75. && #1[[3]] >= -100. && #1[[3]] <= 100.), {#1[[1]], #1[[2]], #1[[3]]}, {-1., 0., 0.}] & )[{-(rx/
                          xunit), ry - (rx*yunit)/
                          xunit, rz - (rx*zunit)/
                          xunit}]}, {{-1., 0., 0.}}]; AppendTo[intersections, {-1., 0., 0.}]; 
                   holeQ = False; If[ !(Chop[surfaceid - 16603.] == 0. || N[Chop[ray1]] < 
                        0.0001) || generationnumber == 1., minresonateraylength = 0., 
                    minresonateraylength = 0.0001]; Block[{RayL, SCx, SCy}, RayL = 1.; SCx = 1.; SCy = 1.; dummymatrix = 
                      Union[Select[intersections, ({RayL, SCx, SCy} = #1; {RayL, SCx, SCy} == {-1., 0., 0.} || N[Chop[RayL]] >= 
                            minresonateraylength) & ]]; intersections = dummymatrix]; 
                   unfilteredintersections = intersections; If[Length[intersections] > 1, 
                    dummymatrix = Take[intersections, {2}]; intersections = 
                      dummymatrix; 1., 1.]; AppendTo[intersectionslist, intersections]; 
                   If[holeQ, AppendTo[holeQlist, 1], AppendTo[holeQlist, 0]]); 
                  holeQlist = Drop[holeQlist, 1]; Drop[intersectionslist, 1]]); Which[
                resonatenumber == 1, {rsx, rsy, rsz, 
                   rex, rey, rez, r1, r2, 
                   r3, r4, r5, r6, r7, 
                   r8, r9, wavelength, intensity, 
                   intensityscale, ray1, index, opticallength, 
                   phaseshift, componentnumber, surfacenumber, surfaceid, 
                   ray2, ray3, sid, frontid, sno1, 
                   sno2, sno3, sno4, rayslot, 
                   oldgenerationnumber, intersectionnumber, pol1, pol2, 
                   pol3} = #1; index1 = index; xunit = 
                  -1.*r1; yunit = -1.*r2; zunit = r3; 
                 xunit = Chop[xunit]; yunit = Chop[yunit]; 
                 zunit = Chop[zunit]; magnitude = (xunit^2 + 
                    yunit^2 + zunit^2)^2^(-1); xunit = xunit/
                   magnitude; yunit = yunit/magnitude; 
                 zunit = zunit/magnitude; rx = 
                  -1.*(-200. + rex); ry = -1.*rey; rz = 
                  rez; rx = Chop[rx]; ry = Chop[ry]; 
                 rz = Chop[rz]; terminatedeflection = False; intersections = 
                  {{-1., 0., 0.}}; intersections = resonateintersections; holeQ = 
                  holeQlist[[resonatenumber]] == 1; {newcomponentnumber, 
                   newsurfacenumber, newsurfaceid} = {2., 1., 16598.}; {componentnumber, 
                   surfacenumber, surfaceid} = {2., 1., 16598.}; {ray1, ray2, 
                   ray3} = intersections[[1]]; rsx = rex; 
                 rsy = rey; rsz = rez; rxhlp = 
                  rx + ray1*xunit; ryhlp = ry + 
                   ray1*yunit; rzhlp = rz + ray1*
                    zunit; rex = 200. - 1.*rxhlp; rey = 
                  -1.*ryhlp; rez = rzhlp; rex = 
                  Chop[rex]; rey = Chop[rey]; rez = 
                  Chop[rez]; If[rex == 0 && rey == 0 && rez == 0, 
                  rex = 200. + rxhlp; rey = ryhlp; 
                   rez = rzhlp, rex]; intersectionnumber = 
                  intersectionnumber + 1.; opticallength = opticallength + 
                   ray1*index1; completedray = {rsx, 
                   rsy, rsz, rex, rey, rez, 
                   r1, r2, r3, r4, r5, 
                   r6, r7, r8, r9, wavelength, 
                   intensity, intensityscale, ray1, index, 
                   opticallength, phaseshift, componentnumber, surfacenumber, 
                   surfaceid, ray2, ray3, sid, frontid, 
                   sno1, sno2, sno3, sno4, rayslot, 
                   generationnumber, intersectionnumber, pol1, pol2, 
                   pol3}; reportedrayspresent = True; intersectionspresent = True; 
                 deflectedrayspresent = True; localdummyoutrayspresent = True; (surfnorm = 
                   1.*Re[{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}]; dummymatrix = ({{-1., 0., 0.}, {0., -1., 0.}, {0., 0., 1.}} . #1 & ) /@ 
                    surfnorm; surfnorm = (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; 
                  surfacerotationmatrix = {{-1., 0, 0}, {0, -1., 0}, {0, 0, 1.}}; surfacetranslationvector = 
                   {200., 0, 0}; activeregion = 0.; completedray = {rsx, rsy, 
                    rsz, rex, rey, rez, r1, 
                    r2, r3, r4, r5, r6, 
                    r7, r8, r9, wavelength, intensity, 
                    intensityscale, ray1, index, opticallength, 
                    phaseshift, componentnumber, surfacenumber, surfaceid, 
                    ray2, ray3, sid, frontid, sno1, 
                    sno2, sno3, sno4, rayslot, 
                    generationnumber, intersectionnumber, pol1, pol2, 
                    pol3}; componentnumber = newcomponentnumber; 
                  surfacenumber = newsurfacenumber; surfaceid = newsurfaceid; 
                  originalray = {rsx, rsy, rsz, rex, 
                    rey, rez, r1, r2, r3, 
                    r4, r5, r6, r7, r8, 
                    r9, wavelength, intensity, intensityscale, 
                    ray1, index, opticallength, phaseshift, 
                    componentnumber, surfacenumber, surfaceid, ray2, 
                    ray3, sid, frontid, sno1, sno2, 
                    sno3, sno4, rayslot, generationnumber, 
                    intersectionnumber, pol1, pol2, pol3}; 
                  intensity = 0.; allcreatedrays = {{rsx, rsy, 
                     rsz, rex, rey, rez, r1, 
                     r2, r3, r4, r5, r6, 
                     r7, r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     generationnumber, intersectionnumber, pol1, pol2, 
                     pol3}}); {allcreatedrays, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 
                    -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, 
                  {completedray}}, resonatenumber == 2, {rsx, rsy, 
                   rsz, rex, rey, rez, r1, 
                   r2, r3, r4, r5, r6, 
                   r7, r8, r9, wavelength, intensity, 
                   intensityscale, ray1, index, opticallength, 
                   phaseshift, componentnumber, surfacenumber, surfaceid, 
                   ray2, ray3, sid, frontid, sno1, 
                   sno2, sno3, sno4, rayslot, 
                   oldgenerationnumber, intersectionnumber, pol1, pol2, 
                   pol3} = #1; index1 = index; xunit = 
                  r2; yunit = -1.*r1; zunit = r3; 
                 xunit = Chop[xunit]; yunit = Chop[yunit]; 
                 zunit = Chop[zunit]; magnitude = (xunit^2 + 
                    yunit^2 + zunit^2)^2^(-1); xunit = xunit/
                   magnitude; yunit = yunit/magnitude; 
                 zunit = zunit/magnitude; rx = 
                  75. + rey; ry = -1.*(-100. + rex); rz = 
                  rez; rx = Chop[rx]; ry = Chop[ry]; 
                 rz = Chop[rz]; terminatedeflection = False; intersections = 
                  {{-1., 0., 0.}}; intersections = resonateintersections; holeQ = 
                  holeQlist[[resonatenumber]] == 1; {newcomponentnumber, 
                   newsurfacenumber, newsurfaceid} = {2., 2., 16599.}; {componentnumber, 
                   surfacenumber, surfaceid} = {2., 2., 16599.}; {ray1, ray2, 
                   ray3} = intersections[[1]]; rsx = rex; 
                 rsy = rey; rsz = rez; rxhlp = 
                  rx + ray1*xunit; ryhlp = ry + 
                   ray1*yunit; rzhlp = rz + ray1*
                    zunit; rex = 100. - 1.*ryhlp; rey = 
                  -75. + rxhlp; rez = rzhlp; rex = 
                  Chop[rex]; rey = Chop[rey]; rez = 
                  Chop[rez]; If[rex == 0 && rey == 0 && rez == 0, 
                  rex = 100. + rxhlp; rey = -75. + ryhlp; 
                   rez = rzhlp, rex]; intersectionnumber = 
                  intersectionnumber + 1.; opticallength = opticallength + 
                   ray1*index1; completedray = {rsx, 
                   rsy, rsz, rex, rey, rez, 
                   r1, r2, r3, r4, r5, 
                   r6, r7, r8, r9, wavelength, 
                   intensity, intensityscale, ray1, index, 
                   opticallength, phaseshift, componentnumber, surfacenumber, 
                   surfaceid, ray2, ray3, sid, frontid, 
                   sno1, sno2, sno3, sno4, rayslot, 
                   generationnumber, intersectionnumber, pol1, pol2, 
                   pol3}; reportedrayspresent = True; intersectionspresent = True; 
                 deflectedrayspresent = True; localdummyoutrayspresent = True; (surfnorm = 
                   1.*Re[{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}]; dummymatrix = ({{0., -1., 0.}, {1., 0., 0.}, {0., 0., 1.}} . #1 & ) /@ 
                    surfnorm; surfnorm = (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; 
                  surfacerotationmatrix = {{0, 1., 0}, {-1., 0, 0}, {0, 0, 1.}}; surfacetranslationvector = 
                   {100., -75., 0}; activeregion = 0.; completedray = {rsx, 
                    rsy, rsz, rex, rey, rez, 
                    r1, r2, r3, r4, r5, 
                    r6, r7, r8, r9, wavelength, 
                    intensity, intensityscale, ray1, index, 
                    opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; componentnumber = 
                   newcomponentnumber; surfacenumber = newsurfacenumber; 
                  surfaceid = newsurfaceid; originalray = {rsx, 
                    rsy, rsz, rex, rey, rez, 
                    r1, r2, r3, r4, r5, 
                    r6, r7, r8, r9, wavelength, 
                    intensity, intensityscale, ray1, index, 
                    opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; intensity = 0.; 
                  allcreatedrays = {{rsx, rsy, rsz, 
                     rex, rey, rez, r1, r2, 
                     r3, r4, r5, r6, r7, 
                     r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     generationnumber, intersectionnumber, pol1, pol2, 
                     pol3}}); {allcreatedrays, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 
                    -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, 
                  {completedray}}, resonatenumber == 3, {rsx, rsy, 
                   rsz, rex, rey, rez, r1, 
                   r2, r3, r4, r5, r6, 
                   r7, r8, r9, wavelength, intensity, 
                   intensityscale, ray1, index, opticallength, 
                   phaseshift, componentnumber, surfacenumber, surfaceid, 
                   ray2, ray3, sid, frontid, sno1, 
                   sno2, sno3, sno4, rayslot, 
                   oldgenerationnumber, intersectionnumber, pol1, pol2, 
                   pol3} = #1; index1 = index; xunit = 
                  r1; yunit = r2; zunit = r3; 
                 xunit = Chop[xunit]; yunit = Chop[yunit]; 
                 zunit = Chop[zunit]; magnitude = (xunit^2 + 
                    yunit^2 + zunit^2)^2^(-1); xunit = xunit/
                   magnitude; yunit = yunit/magnitude; 
                 zunit = zunit/magnitude; rx = rex; 
                 ry = rey; rz = rez; rx = 
                  Chop[rx]; ry = Chop[ry]; rz = 
                  Chop[rz]; terminatedeflection = False; intersections = {{-1., 0., 0.}}; 
                 intersections = resonateintersections; holeQ = holeQlist[[
                    resonatenumber]] == 1; {newcomponentnumber, newsurfacenumber, 
                   newsurfaceid} = {2., 3., 16600.}; {componentnumber, surfacenumber, 
                   surfaceid} = {2., 3., 16600.}; {ray1, ray2, ray3} = 
                  intersections[[1]]; rsx = rex; rsy = 
                  rey; rsz = rez; rxhlp = rx + 
                   ray1*xunit; ryhlp = ry + ray1*
                    yunit; rzhlp = rz + ray1*zunit; 
                 rex = rxhlp; rey = ryhlp; rez = 
                  rzhlp; rex = Chop[rex]; rey = 
                  Chop[rey]; rez = Chop[rez]; If[rex == 0 && 
                   rey == 0 && rez == 0, rex = rxhlp; 
                   rey = ryhlp; rez = rzhlp, rex]; 
                 intersectionnumber = intersectionnumber + 1.; opticallength = 
                  opticallength + ray1*index1; completedray = 
                  {rsx, rsy, rsz, rex, rey, 
                   rez, r1, r2, r3, r4, 
                   r5, r6, r7, r8, r9, 
                   wavelength, intensity, intensityscale, ray1, 
                   index, opticallength, phaseshift, componentnumber, 
                   surfacenumber, surfaceid, ray2, ray3, 
                   sid, frontid, sno1, sno2, sno3, 
                   sno4, rayslot, generationnumber, intersectionnumber, 
                   pol1, pol2, pol3}; reportedrayspresent = True; 
                 intersectionspresent = True; deflectedrayspresent = True; 
                 localdummyoutrayspresent = True; (surfnorm = 1.*Re[{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}]; 
                  dummymatrix = (#1 & ) /@ surfnorm; surfnorm = (#1/Sqrt[#1 . #1] & ) /@ 
                    dummymatrix; surfacerotationmatrix = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; 
                  surfacetranslationvector = {0, 0, 0}; activeregion = 0.; completedray = 
                   {rsx, rsy, rsz, rex, rey, 
                    rez, r1, r2, r3, r4, 
                    r5, r6, r7, r8, r9, 
                    wavelength, intensity, intensityscale, ray1, 
                    index, opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; componentnumber = 
                   newcomponentnumber; surfacenumber = newsurfacenumber; 
                  surfaceid = newsurfaceid; originalray = {rsx, 
                    rsy, rsz, rex, rey, rez, 
                    r1, r2, r3, r4, r5, 
                    r6, r7, r8, r9, wavelength, 
                    intensity, intensityscale, ray1, index, 
                    opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; intensity = 0.; 
                  allcreatedrays = {{rsx, rsy, rsz, 
                     rex, rey, rez, r1, r2, 
                     r3, r4, r5, r6, r7, 
                     r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     generationnumber, intersectionnumber, pol1, pol2, 
                     pol3}}); {allcreatedrays, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 
                    -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, 
                  {completedray}}, resonatenumber == 4, {rsx, rsy, 
                   rsz, rex, rey, rez, r1, 
                   r2, r3, r4, r5, r6, 
                   r7, r8, r9, wavelength, intensity, 
                   intensityscale, ray1, index, opticallength, 
                   phaseshift, componentnumber, surfacenumber, surfaceid, 
                   ray2, ray3, sid, frontid, sno1, 
                   sno2, sno3, sno4, rayslot, 
                   oldgenerationnumber, intersectionnumber, pol1, pol2, 
                   pol3} = #1; index1 = index; xunit = 
                  -1.*r2; yunit = r1; zunit = r3; 
                 xunit = Chop[xunit]; yunit = Chop[yunit]; 
                 zunit = Chop[zunit]; magnitude = (xunit^2 + 
                    yunit^2 + zunit^2)^2^(-1); xunit = xunit/
                   magnitude; yunit = yunit/magnitude; 
                 zunit = zunit/magnitude; rx = 
                  -1.*(-75. + rey); ry = -100. + rex; rz = 
                  rez; rx = Chop[rx]; ry = Chop[ry]; 
                 rz = Chop[rz]; terminatedeflection = False; intersections = 
                  {{-1., 0., 0.}}; intersections = resonateintersections; holeQ = 
                  holeQlist[[resonatenumber]] == 1; {newcomponentnumber, 
                   newsurfacenumber, newsurfaceid} = {2., 4., 16601.}; {componentnumber, 
                   surfacenumber, surfaceid} = {2., 4., 16601.}; {ray1, ray2, 
                   ray3} = intersections[[1]]; rsx = rex; 
                 rsy = rey; rsz = rez; rxhlp = 
                  rx + ray1*xunit; ryhlp = ry + 
                   ray1*yunit; rzhlp = rz + ray1*
                    zunit; rex = 100. + ryhlp; rey = 
                  75. - 1.*rxhlp; rez = rzhlp; rex = 
                  Chop[rex]; rey = Chop[rey]; rez = 
                  Chop[rez]; If[rex == 0 && rey == 0 && rez == 0, 
                  rex = 100. + rxhlp; rey = 75. + ryhlp; 
                   rez = rzhlp, rex]; intersectionnumber = 
                  intersectionnumber + 1.; opticallength = opticallength + 
                   ray1*index1; completedray = {rsx, 
                   rsy, rsz, rex, rey, rez, 
                   r1, r2, r3, r4, r5, 
                   r6, r7, r8, r9, wavelength, 
                   intensity, intensityscale, ray1, index, 
                   opticallength, phaseshift, componentnumber, surfacenumber, 
                   surfaceid, ray2, ray3, sid, frontid, 
                   sno1, sno2, sno3, sno4, rayslot, 
                   generationnumber, intersectionnumber, pol1, pol2, 
                   pol3}; reportedrayspresent = True; intersectionspresent = True; 
                 deflectedrayspresent = True; localdummyoutrayspresent = True; (surfnorm = 
                   1.*Re[{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}]; dummymatrix = ({{0., 1., 0.}, {-1., 0., 0.}, {0., 0., 1.}} . #1 & ) /@ 
                    surfnorm; surfnorm = (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; 
                  surfacerotationmatrix = {{0, -1., 0}, {1., 0, 0}, {0, 0, 1.}}; surfacetranslationvector = 
                   {100., 75., 0}; activeregion = 0.; completedray = {rsx, 
                    rsy, rsz, rex, rey, rez, 
                    r1, r2, r3, r4, r5, 
                    r6, r7, r8, r9, wavelength, 
                    intensity, intensityscale, ray1, index, 
                    opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; componentnumber = 
                   newcomponentnumber; surfacenumber = newsurfacenumber; 
                  surfaceid = newsurfaceid; originalray = {rsx, 
                    rsy, rsz, rex, rey, rez, 
                    r1, r2, r3, r4, r5, 
                    r6, r7, r8, r9, wavelength, 
                    intensity, intensityscale, ray1, index, 
                    opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; intensity = 0.; 
                  allcreatedrays = {{rsx, rsy, rsz, 
                     rex, rey, rez, r1, r2, 
                     r3, r4, r5, r6, r7, 
                     r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     generationnumber, intersectionnumber, pol1, pol2, 
                     pol3}}); {allcreatedrays, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 
                    -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, 
                  {completedray}}, resonatenumber == 5, {rsx, rsy, 
                   rsz, rex, rey, rez, r1, 
                   r2, r3, r4, r5, r6, 
                   r7, r8, r9, wavelength, intensity, 
                   intensityscale, ray1, index, opticallength, 
                   phaseshift, componentnumber, surfacenumber, surfaceid, 
                   ray2, ray3, sid, frontid, sno1, 
                   sno2, sno3, sno4, rayslot, 
                   oldgenerationnumber, intersectionnumber, pol1, pol2, 
                   pol3} = #1; index1 = index; xunit = 
                  r3; yunit = r2; zunit = -1.*r1; 
                 xunit = Chop[xunit]; yunit = Chop[yunit]; 
                 zunit = Chop[zunit]; magnitude = (xunit^2 + 
                    yunit^2 + zunit^2)^2^(-1); xunit = xunit/
                   magnitude; yunit = yunit/magnitude; 
                 zunit = zunit/magnitude; rx = 
                  75. + rez; ry = rey; rz = 
                  -1.*(-100. + rex); rx = Chop[rx]; ry = 
                  Chop[ry]; rz = Chop[rz]; terminatedeflection = False; 
                 intersections = {{-1., 0., 0.}}; intersections = resonateintersections; 
                 holeQ = holeQlist[[resonatenumber]] == 1; 
                 {newcomponentnumber, newsurfacenumber, newsurfaceid} = {2., 5., 16602.}; 
                 {componentnumber, surfacenumber, surfaceid} = {2., 5., 16602.}; 
                 {ray1, ray2, ray3} = intersections[[1]]; 
                 rsx = rex; rsy = rey; rsz = 
                  rez; rxhlp = rx + ray1*xunit; 
                 ryhlp = ry + ray1*yunit; rzhlp = 
                  rz + ray1*zunit; rex = 
                  100. - 1.*rzhlp; rey = ryhlp; rez = 
                  -75. + rxhlp; rex = Chop[rex]; rey = 
                  Chop[rey]; rez = Chop[rez]; If[rex == 0 && 
                   rey == 0 && rez == 0, rex = 100. + rxhlp; 
                   rey = ryhlp; rez = -75. + rzhlp, 
                  rex]; intersectionnumber = intersectionnumber + 1.; 
                 opticallength = opticallength + ray1*index1; 
                 completedray = {rsx, rsy, rsz, rex, 
                   rey, rez, r1, r2, r3, 
                   r4, r5, r6, r7, r8, 
                   r9, wavelength, intensity, intensityscale, 
                   ray1, index, opticallength, phaseshift, 
                   componentnumber, surfacenumber, surfaceid, ray2, 
                   ray3, sid, frontid, sno1, sno2, 
                   sno3, sno4, rayslot, generationnumber, 
                   intersectionnumber, pol1, pol2, pol3}; 
                 reportedrayspresent = True; intersectionspresent = True; deflectedrayspresent = 
                  True; localdummyoutrayspresent = True; (surfnorm = 1.*Re[{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}]; 
                  dummymatrix = ({{0., 0., -1.}, {0., 1., 0.}, {1., 0., 0.}} . #1 & ) /@ surfnorm; 
                  surfnorm = (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; surfacerotationmatrix = 
                   {{0, 0, 1.}, {0, 1., 0}, {-1., 0, 0}}; surfacetranslationvector = {100., 0, -75.}; activeregion = 0.; 
                  completedray = {rsx, rsy, rsz, rex, 
                    rey, rez, r1, r2, r3, 
                    r4, r5, r6, r7, r8, 
                    r9, wavelength, intensity, intensityscale, 
                    ray1, index, opticallength, phaseshift, 
                    componentnumber, surfacenumber, surfaceid, ray2, 
                    ray3, sid, frontid, sno1, sno2, 
                    sno3, sno4, rayslot, generationnumber, 
                    intersectionnumber, pol1, pol2, pol3}; 
                  componentnumber = newcomponentnumber; surfacenumber = 
                   newsurfacenumber; surfaceid = newsurfaceid; originalray = 
                   {rsx, rsy, rsz, rex, rey, 
                    rez, r1, r2, r3, r4, 
                    r5, r6, r7, r8, r9, 
                    wavelength, intensity, intensityscale, ray1, 
                    index, opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; intensity = 0.; 
                  allcreatedrays = {{rsx, rsy, rsz, 
                     rex, rey, rez, r1, r2, 
                     r3, r4, r5, r6, r7, 
                     r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     generationnumber, intersectionnumber, pol1, pol2, 
                     pol3}}); {allcreatedrays, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 
                    -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, 
                  {completedray}}, resonatenumber == 6, {rsx, rsy, 
                   rsz, rex, rey, rez, r1, 
                   r2, r3, r4, r5, r6, 
                   r7, r8, r9, wavelength, intensity, 
                   intensityscale, ray1, index, opticallength, 
                   phaseshift, componentnumber, surfacenumber, surfaceid, 
                   ray2, ray3, sid, frontid, sno1, 
                   sno2, sno3, sno4, rayslot, 
                   oldgenerationnumber, intersectionnumber, pol1, pol2, 
                   pol3} = #1; index1 = index; xunit = 
                  -1.*r3; yunit = -1.*r2; zunit = 
                  -1.*r1; xunit = Chop[xunit]; yunit = 
                  Chop[yunit]; zunit = Chop[zunit]; magnitude = 
                  (xunit^2 + yunit^2 + zunit^2)^2^(-1); xunit = 
                  xunit/magnitude; yunit = yunit/
                   magnitude; zunit = zunit/magnitude; 
                 rx = -1.*(-75. + rez); ry = -1.*rey; 
                 rz = -1.*(-100. + rex); rx = Chop[rx]; 
                 ry = Chop[ry]; rz = Chop[rz]; 
                 terminatedeflection = False; intersections = {{-1., 0., 0.}}; intersections = 
                  resonateintersections; holeQ = holeQlist[[
                    resonatenumber]] == 1; {newcomponentnumber, newsurfacenumber, 
                   newsurfaceid} = {2., 6., 16603.}; {componentnumber, surfacenumber, 
                   surfaceid} = {2., 6., 16603.}; {ray1, ray2, ray3} = 
                  intersections[[1]]; rsx = rex; rsy = 
                  rey; rsz = rez; rxhlp = rx + 
                   ray1*xunit; ryhlp = ry + ray1*
                    yunit; rzhlp = rz + ray1*zunit; 
                 rex = 100. - 1.*rzhlp; rey = -1.*ryhlp; 
                 rez = 75. - 1.*rxhlp; rex = Chop[rex]; 
                 rey = Chop[rey]; rez = Chop[rez]; 
                 If[rex == 0 && rey == 0 && rez == 0, rex = 
                    100. + rxhlp; rey = ryhlp; rez = 
                    75. + rzhlp, rex]; intersectionnumber = 
                  intersectionnumber + 1.; opticallength = opticallength + 
                   ray1*index1; completedray = {rsx, 
                   rsy, rsz, rex, rey, rez, 
                   r1, r2, r3, r4, r5, 
                   r6, r7, r8, r9, wavelength, 
                   intensity, intensityscale, ray1, index, 
                   opticallength, phaseshift, componentnumber, surfacenumber, 
                   surfaceid, ray2, ray3, sid, frontid, 
                   sno1, sno2, sno3, sno4, rayslot, 
                   generationnumber, intersectionnumber, pol1, pol2, 
                   pol3}; reportedrayspresent = True; intersectionspresent = True; 
                 deflectedrayspresent = True; localdummyoutrayspresent = True; (surfnorm = 
                   1.*Re[{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}]; dummymatrix = ({{0., 0., -1.}, {0., -1., 0.}, {-1., 0., 0.}} . #1 & ) /@ 
                    surfnorm; surfnorm = (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; 
                  surfacerotationmatrix = {{0, 0, -1.}, {0, -1., 0}, {-1., 0, 0}}; surfacetranslationvector = 
                   {100., 0, 75.}; activeregion = 0.; completedray = {rsx, 
                    rsy, rsz, rex, rey, rez, 
                    r1, r2, r3, r4, r5, 
                    r6, r7, r8, r9, wavelength, 
                    intensity, intensityscale, ray1, index, 
                    opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; componentnumber = 
                   newcomponentnumber; surfacenumber = newsurfacenumber; 
                  surfaceid = newsurfaceid; originalray = {rsx, 
                    rsy, rsz, rex, rey, rez, 
                    r1, r2, r3, r4, r5, 
                    r6, r7, r8, r9, wavelength, 
                    intensity, intensityscale, ray1, index, 
                    opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}; intensity = 0.; 
                  allcreatedrays = {{rsx, rsy, rsz, 
                     rex, rey, rez, r1, r2, 
                     r3, r4, r5, r6, r7, 
                     r8, r9, wavelength, intensity, 
                     intensityscale, ray1, index, opticallength, 
                     phaseshift, componentnumber, surfacenumber, surfaceid, 
                     ray2, ray3, sid, frontid, sno1, 
                     sno2, sno3, sno4, rayslot, 
                     generationnumber, intersectionnumber, pol1, pol2, 
                     pol3}}); {allcreatedrays, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 
                    -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, 
                  {completedray}}, True, deflecteddummyrayspresent = True; 
                 reporteddummyrayspresent = True; localoutrayspresent = True; {{{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 
                    1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, 
                  {{rsx, rsy, rsz, rex, rey, 
                    rez, r1, r2, r3, r4, 
                    r5, r6, r7, r8, r9, 
                    wavelength, intensity, intensityscale, ray1, 
                    index, opticallength, phaseshift, componentnumber, 
                    surfacenumber, surfaceid, ray2, ray3, 
                    sid, frontid, sno1, sno2, sno3, 
                    sno4, rayslot, generationnumber, intersectionnumber, 
                    pol1, pol2, pol3}}, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 
                    1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}}], 
              deflecteddummyrayspresent = True; reporteddummyrayspresent = True; 
               localdummyoutrayspresent = True; {{{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 
                  1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 
                  1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}}, 
                {{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., -1.7976931348623157*^308, 1., 1., 1., 1., 0., -1., -1., -1., 1., 1., 1., 1., 1., 
                  1., 1., 1., 1., 1., 1., 1., 1., 1.}}}]) & ) /@ inrays]; 
       If[deflectedrayspresent && deflecteddummyrayspresent, deflectedlocalrays = 
         Select[deflectedlocalrays, #1[[23]] > 0. & ]*1., 1.]; If[localoutrayspresent, 
        If[localdummyoutrayspresent, localoutrays = Select[localoutrays, #1[[23]] > 0. & ]; 1., 
          1.]; outrays = Join[outrays, localoutrays]; 1., 1.]; 
       If[reportedrayspresent && reporteddummyrayspresent, 
        reportedlocalrays = Select[reportedlocalrays, #1[[23]] > 0. & ]; 1., 1.]); 
      If[reportedrayspresent, raysreported = True; reportedrays = 
         Join[reportedrays, {{-1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, 
            -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308, -1.7976931348623157*^308}}, 
          reportedlocalrays]; 1., 1.]; inrays = deflectedlocalrays; 
      propagate = generationnumber < generationlimit]; 
    propagatedrays = If[Length[outrays] > 1, Drop[outrays, 1], outrays]; 
    If[generationnumber >= generationlimit && deflectedrayspresent, 
     propagatedrays = Join[propagatedrays, deflectedlocalrays]; 1., 1.]); 
   If[Length[reportedrays] > 2, reportedrays = Drop[reportedrays, 2], 
    reportedrays]; ({#1[[1]], #1[[2]], #1[[3]], #1[[4]], #1[[5]], #1[[6]]} & ) /@ reportedrays]];
(*Surfaces*)
rot[pt_,a_,cen_]:=cen+{{Cos[-a],Sin[-a]},{-Sin[-a],Cos[-a]}}.(pt-cen);
surface[r1_,r2_,ap_,x_,t_,y_]:=Join[Table[{r1+x-Sqrt[r1^2-yy^2]*Sign[r1],y+yy} ,{yy,-ap/2,ap/2,ap/20}],Reverse[Table[{r2+x+t-Sqrt[r2^2-yy^2]*Sign[r2],y+yy},{yy,-ap/2,ap/2,ap/20}]]];
system[as1_,as2_,wl1_,wl2_,xs1_,xs2_,ys1_,ys2_,aa_,ab_,apa_,apb_,ca1_,ca2_,cb1_,cb2_,ta_,tb_,xa1_,xb1_,ya1_,yb1_]:={micronsToColor[wl1],
If[(wl1<.380)||(wl1>.780),AbsoluteDashing[{3,3}],Dashing[{}]],Line[Map[{#[[{1,2}]],#[[{4,5}]]}&,rayfunction[as1,wl1,xs1,ys1,aa,ab,apa,apb,ca1,ca2,cb1,cb2,ta,tb,xa1,xb1,ya1,yb1]]],
micronsToColor[wl2],
If[(wl2<.380)||(wl2>.780),AbsoluteDashing[{3,3}],Dashing[{}]],
Line[Map[{#[[{1,2}]],#[[{4,5}]]}&,rayfunction[as2,wl2,xs2,ys2,aa,ab,apa,apb,ca1,ca2,cb1,cb2,ta,tb,xa1,xb1,ya1,yb1]]]};
locatorSize=10;
locatorDot={Lighter[Blue],Disk[{0,0},.5]};
locatorMaster={Lighter[Red],Disk[{0,0},.8]};
locatorLeftRight={Lighter[Blue],AbsoluteThickness[2],Arrowheads[{{.5,1},{-.5,0}}],Arrow[{{-1,0},{1,0}}]};
locatorUpDown={Lighter[Blue],AbsoluteThickness[2],Arrowheads[{{.5,1},{-.5,0}}],Arrow[{{0,-1},{0,1}}]};
locatorAll={Lighter[Blue],AbsoluteThickness[2],Arrowheads[{{.5,1},{-.5,0}}],Arrow[{{-1,0},{1,0}}],Arrow[{{0,-1},{0,1}}]};
locatorLeftRight=locatorUpDown=locatorAll=locatorDot;
lensWithControls[Dynamic[r1_],Dynamic[r2_],Dynamic[a_],Dynamic[ap_],Dynamic[x_],Dynamic[t_],Dynamic[y_]]:={
Dynamic[Rotate[Polygon[surface[1/r1,1/r2,ap,x,t,y]],a,{x,y}]],
Locator[Dynamic[{x,y}],Graphics[locatorMaster,ImageSize->locatorSize]],
Locator[Dynamic[{x,y}+5{Cos[a],Sin[a]},(a=Apply[ArcTan,#-{x,y}])&],Graphics[locatorDot,ImageSize->locatorSize]],

Locator[Dynamic[rot[{x+t,y},a,{x,y}],(t=First[rot[#,-a,{x,y}]]-x)&],Graphics[locatorLeftRight,ImageSize->locatorSize]],
Locator[Dynamic[rot[{((1/r1+x-Sqrt[(1/r1)^2-(ap/2)^2]*Sign[1/r1])+(1/r2+x-Sqrt[(1/r2)^2-(ap/2)^2]*Sign[1/r2]+t))/2,y-ap/2},a,{x,y}],(ap=2(y-Last[rot[#,-a,{x,y}]]))&],Graphics[locatorUpDown,ImageSize->locatorSize]],
Locator[Dynamic[rot[{1/r1+x-Sqrt[(1/r1)^2-(ap/2)^2]*Sign[1/r1],y+ap/2},a,{x,y}],(r1=With[{c=First[rot[#,-a,{x,y}]]},-((8 (x-c))/((ap^2)+4 (x^2)-8 x c+4 (c^2)))])&],Graphics[locatorLeftRight,ImageSize->locatorSize]],
Locator[Dynamic[rot[{1/r2+x-Sqrt[(1/r2)^2-(ap/2)^2]*Sign[1/r2]+t,y+ap/2},a,{x,y}],(With[{c=First[rot[#,-a,{x,y}]]},r2=-((8 (x-(c-t)))/((ap^2)+4 (x^2)-8 x (c-t)+4 (((c-t))^2)))])&],Graphics[locatorLeftRight,ImageSize->locatorSize]]
};
Attributes[LensLayout]={HoldAll};
LensLayout[{wavelength1_,wavelength2_,sourceAngle_,sourceXY1_,sourceXY2_,lensAngle_,thicknessAndAperture1_,thicknessAndAperture2_,curvature1_,curvature2_,a_,b_}]:=Graphics[{
Dynamic[{
system[sourceAngle[[1]]/Degree,sourceAngle[[2]]/Degree,wavelength1,wavelength2,sourceXY1[[1]],sourceXY2[[1]],sourceXY1[[2]],sourceXY2[[2]],lensAngle[[1]]/Degree,lensAngle[[2]]/Degree,thicknessAndAperture1[[2]],thicknessAndAperture2[[2]],curvature1[[1]],curvature1[[2]],curvature2[[1]],curvature2[[2]],thicknessAndAperture1[[1]],thicknessAndAperture2[[1]],a[[1]],b[[1]],a[[2]],b[[2]]],
Arrowheads->Large,Lighter[Blue],
Arrow[{sourceXY1,sourceXY1+10{Cos[sourceAngle[[1]]],Sin[sourceAngle[[1]]]}}],
Arrow[{sourceXY2,sourceXY2+10{Cos[sourceAngle[[2]]],Sin[sourceAngle[[2]]]}}],
Arrowheads->Medium,
Arrow[{a,a+5{Cos[lensAngle[[1]]],Sin[lensAngle[[1]]]}}],
Arrow[{b,b+5{Cos[lensAngle[[2]]],Sin[lensAngle[[2]]]}}]}],
Locator[Dynamic[sourceXY1],Graphics[locatorAll,ImageSize->locatorSize]],
Locator[Dynamic[sourceXY2],Graphics[locatorAll,ImageSize->locatorSize]],
Locator[Dynamic[sourceXY1+10{Cos[sourceAngle[[1]]],Sin[sourceAngle[[1]]]},(sourceAngle[[1]]=Apply[ArcTan,#-sourceXY1])&],Graphics[locatorAll,ImageSize->locatorSize]],
Locator[Dynamic[sourceXY2+10{Cos[sourceAngle[[2]]],Sin[sourceAngle[[2]]]},(sourceAngle[[2]]=Apply[ArcTan,#-sourceXY2])&],Graphics[locatorAll,ImageSize->locatorSize]],
Cyan,Opacity[0.4],
lensWithControls[Dynamic[curvature1[[1]]],Dynamic[curvature1[[2]]],Dynamic[lensAngle[[1]]],Dynamic[thicknessAndAperture1[[2]]],Dynamic[a[[1]]],Dynamic[thicknessAndAperture1[[1]]],Dynamic[a[[2]]]],
RGBColor[.4,.6,.8],
lensWithControls[Dynamic[curvature2[[1]]],Dynamic[curvature2[[2]]],Dynamic[lensAngle[[2]]],Dynamic[thicknessAndAperture2[[2]]],Dynamic[b[[1]]],Dynamic[thicknessAndAperture2[[1]]],Dynamic[b[[2]]]]
},PlotRange->{{-5,205},{-80,80}},ImageSize->500];
)]

CAPTION

This demonstrates an interactive lens system. Move any blue dot to control the location, thickness, and curvature of the lenses, or the starting point or angle of the rays. When you see rays shooting off at odd angles, you are seeing an accurate depiction of what happens when rays either miss the lens or are reflected instead of refracted because of the angle of incidence.

THIS NOTEBOOK IS THE SOURCE CODE FROM

"Designing Optical Systems" from the Wolfram Demonstrations Project
 http://demonstrations.wolfram.com/DesigningOpticalSystems/

Contributed by: Donald Barnhart, developer of Rayica and LensLab

A full-function Wolfram Mathematica system (Version 6 or higher) is required to edit this notebook.
GET WOLFRAM MATHEMATICA »

© Wolfram Demonstrations Project & Contributors   |   Terms of Use   |   Make a new version of this Demonstration »