Lens Doublets

Manipulate[
 Column[
  {Show[{Graphics[{{Opacity[0.8], GeometricTransformation[
         {AbsoluteThickness[1], AbsolutePointSize[1], EdgeForm[], 
          RGBColor[0.4, 0.6, 0.8], Polygon[Dynamic[
            (Function[{element}, Take[#1, element]] /@ {{1, 27}, 
                {28, 54}} & )[polygons11188572008312134018969656[r11, 
              r12, r13]]]]}, Dynamic[BesselJ[1., 3.]; 
           {{{Cos[Degree*θ13], -Sin[Degree*θ13]}, 
             {Sin[Degree*θ13], Cos[Degree*θ13]}}, {x13, 0}}]], 
        GeometricTransformation[{AbsoluteThickness[1], 
          AbsolutePointSize[1], EdgeForm[], RGBColor[0.4, 0.6, 0.8], 
          Polygon[Dynamic[(Function[{element}, Take[#1, element]] /@ {
                {1, 27}, {28, 54}} & )[
             polygons12188572008312134018987715[r21, r22, r23]]]]}, 
         Dynamic[BesselJ[1., 3.]; {{{Cos[Degree*θ23], 
              -Sin[Degree*θ23]}, {Sin[Degree*θ23], Cos[Degree*
                θ23]}}, {x23, 0}}]], 
        staticpolygons13188572008312134018989499, Opacity[1], 
        RGBColor[0.038747429704798886, 0.9161617552424038, 
         0.03984721226918142], AbsoluteThickness[Medium], 
        Line[Dynamic[({#1[[{1, 2}]], #1[[{4, 5}]]} & ) /@ 
           Select[raytracefunction1188572008312134018995817[r11, r12, 
             r13, r21, r22, r23, x13, x23, θ13, θ23], 
            #1[[1]] > -9.87654321*^8 & ], SynchronousUpdating -> 
           False]], Arrowheads[Large], Thickness[Large], 
        Dashing[{Small, Small}], RGBColor[1, 1, 0], 
        Arrow[Dynamic[If[laststartposθ13 != {x13, 0} || 
             lastθ13 =!= θ13, laststartposθ13 = {x13, 0}; 
             lastθ13 = θ13; {xθ13, yθ13} = laststartposθ13 + 
               radiusθ13*{Cos[θ13*Degree], Sin[θ13*Degree]}]; 
           {laststartposθ13, {xθ13, yθ13}}, 
          ({laststartposθ13, {xθ13, yθ13}} = #1) & ]], 
        Tooltip[Locator[Dynamic[BesselJ[1., 3.]; 
            (If[{xθ13, yθ13} != lastlocposθ13 || lastθ13 != 
                θ13, lastlocposθ13 = {xθ13, yθ13}; lastθ13 = 
                θ13; {xθ13, yθ13} = laststartposθ13 + 
                 radiusθ13*{Cos[θ13*Degree], Sin[θ13*Degree]}]; 
             {xθ13, yθ13}), ({xθ13, yθ13} = #1; 
             xoffθ13 = xθ13 - x13; yoffθ13 = yθ13 + 0; 
             radiusθ13 = Sqrt[xoffθ13^2 + yoffθ13^2]; 
             θ13 = Chop[ArcTan[xoffθ13, yoffθ13]/Degree]; 
             {xθ13, yθ13}) & ], Graphics[{RGBColor[1/3, 1/3, 1], 
            Disk[{0, 0}, 0.5]}, ImageSize -> 10]], 
         Dynamic[StringJoin["θ13 = ", ToString[θ13]]]], 
        Arrowheads[Large], Thickness[Large], Dashing[{Small, Small}], 
        RGBColor[1, 1, 0], Arrow[Dynamic[
          If[laststartposθ23 != {x23, 0} || lastθ23 =!= θ23, 
            laststartposθ23 = {x23, 0}; lastθ23 = θ23; 
             {xθ23, yθ23} = laststartposθ23 + radiusθ23*
                {Cos[θ23*Degree], Sin[θ23*Degree]}]; 
           {laststartposθ23, {xθ23, yθ23}}, 
          ({laststartposθ23, {xθ23, yθ23}} = #1) & ]], 
        Tooltip[Locator[Dynamic[BesselJ[1., 3.]; 
            (If[{xθ23, yθ23} != lastlocposθ23 || lastθ23 != 
                θ23, lastlocposθ23 = {xθ23, yθ23}; lastθ23 = 
                θ23; {xθ23, yθ23} = laststartposθ23 + 
                 radiusθ23*{Cos[θ23*Degree], Sin[θ23*Degree]}]; 
             {xθ23, yθ23}), ({xθ23, yθ23} = #1; 
             xoffθ23 = xθ23 - x23; yoffθ23 = yθ23 + 0; 
             radiusθ23 = Sqrt[xoffθ23^2 + yoffθ23^2]; 
             θ23 = Chop[ArcTan[xoffθ23, yoffθ23]/Degree]; 
             {xθ23, yθ23}) & ], Graphics[{RGBColor[1/3, 1/3, 1], 
            Disk[{0, 0}, 0.5]}, ImageSize -> 10]], 
         Dynamic[StringJoin["θ23 = ", ToString[θ23]]]], 
        Tooltip[Locator[Dynamic[BesselJ[1., 3.]; {x13, 0}, 
           (x13 = Chop[#1[[1]]]) & ], Graphics[{RGBColor[1, 1/3, 1/3], 
            Disk[{0, 0}, 0.8]}, ImageSize -> 10]], 
         Dynamic[StringJoin["x13 = ", ToString[x13]]]], 
        Tooltip[Locator[Dynamic[BesselJ[1., 3.]; {x23, 0}, 
           (x23 = Chop[#1[[1]]]) & ], Graphics[{RGBColor[1, 1/3, 1/3], 
            Disk[{0, 0}, 0.8]}, ImageSize -> 10]], 
         Dynamic[StringJoin["x23 = ", ToString[x23]]]]}}]}, 
    {PlotLabel -> "lens doublets", PlotRange -> {{0, 200}, {-15, 15}}, 
     ImageSize -> {500, 200}, Axes -> True, ImageSize -> 800, 
     AspectRatio -> Automatic}]}, 
  Center],
 {{r11, 71.89}, 50, 80, ControlType -> None}, 
 {{r12, 26.81}, 20, 50, ControlType -> None}, 
 {{r13, -38.35}, -50, -30, ControlType -> None}, 
 {{r21, 38.35}, 30, 50, ControlType -> None}, 
 {{r22, -26.81}, -50, -20, ControlType -> None}, 
 {{r23, -71.89}, -80, -50, ControlType -> None}, 
 {{x13, 50}, 20, 150, ControlType -> None}, 
 {{x23, 100}, 20, 150, ControlType -> None}, 
 {{θ13, 0}, 0, 360, ControlType -> None}, 
 {{θ23, 0}, 0, 360, ControlType -> None}, 
 Style[Grid[{{"r11", Manipulator[Dynamic[r11], {50, 80}], 
     Dynamic[r11]}}]], 
 Style[Grid[{{"r12", Manipulator[Dynamic[r12], {20, 50}], 
     Dynamic[r12]}}]], 
 Style[Grid[{{"r13", Manipulator[Dynamic[r13], {-50, -30}], 
     Dynamic[r13]}}]], 
 Style[Grid[{{"r21", Manipulator[Dynamic[r21], {30, 50}], 
     Dynamic[r21]}}]], 
 Style[Grid[{{"r22", Manipulator[Dynamic[r22], {-50, -20}], 
     Dynamic[r22]}}]], 
 Style[Grid[{{"r23", Manipulator[Dynamic[r23], {-80, -50}], 
     Dynamic[r23]}}]], {{xθ13, 55.}, ControlType -> None}, 
 {{yθ13, 0}, ControlType -> None}, {{radiusθ13, 5.}, 
  ControlType -> None}, {{xoffθ13, 5.}, ControlType -> None}, 
 {{yoffθ13, 0}, ControlType -> None}, {{laststartposθ13, {50, 0}}, 
  ControlType -> None}, {{lastlocposθ13, {55., 0}}, 
  ControlType -> None}, {{lastθ13, 0}, ControlType -> None}, 
 {{xθ23, 105.}, ControlType -> None}, 
 {{yθ23, 0}, ControlType -> None}, {{radiusθ23, 5.}, 
  ControlType -> None}, {{xoffθ23, 5.}, ControlType -> None}, 
 {{yoffθ23, 0}, ControlType -> None}, {{laststartposθ23, {100, 0}}, 
  ControlType -> None}, {{lastlocposθ23, {105., 0}}, 
  ControlType -> None}, {{lastθ23, 0}, ControlType -> None}, 
 {{staticpolygons13188572008312134018989499, 
   {GrayLevel[0], AbsoluteThickness[0.5], 
    Line[{{{200., 100.}, {200., -100.}, {200., -100.}, {200., 100.}, 
       {200., 100.}}, {{200., -100.}, {0., -100.}, {0., -100.}, 
       {200., -100.}, {200., -100.}}, {{0., -100.}, {0., 100.}, 
       {0., 100.}, {0., -100.}, {0., -100.}}, {{0., 100.}, {200., 100.}, 
       {200., 100.}, {0., 100.}, {0., 100.}}, {{200., -100.}, 
       {200., 100.}, {0., 100.}, {0., -100.}, {200., -100.}}, 
      {{200., 100.}, {200., -100.}, {0., -100.}, {0., 100.}, 
       {200., 100.}}}]}}, ControlType -> None}, 
 Initialization :> (
   raytracefunction1188572008312134018995817 = 
    Compile[{{r11, _Real}, {r12, _Real}, {r13, _Real}, {r21, _Real}, 
      {r22, _Real}, {r23, _Real}, {x13, _Real}, {x23, _Real}, 
      {θ13, _Real}, {θ23, _Real}}, 
     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, LocRayL, rayslot, raysreported, 
       raytiltx, raytiltz, reporteddummyrayspresent, reportedlocalrays, 
       reportedrays, reportedrayspresent, resonateintersections, 
       resonatenumber, resonateraylength, rex, rey, rez, rotationaxis, 
       rsx, LocRSx, rsy, LocRSy, rsz, LocRSz, RTx, LocRTy, LocRTz, rx, 
       rxhlp, ry, ryhlp, rz, rzhlp, LocSCx, LocSCy, 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, LocTV1, LocTV10, LocTV2, LocTV3, 
       LocTV4, LocTV5, LocTV6, LocTV7, LocTV8, LocTV9, 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.; LocRayL = 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.; 
       LocRSx = 1.; rsy = 1.; LocRSy = 1.; rsz = 1.; LocRSz = 1.; 
       RTx = 1.; LocRTy = 1.; LocRTz = 1.; rx = 1.; rxhlp = 1.; ry = 1.; 
       ryhlp = 1.; rz = 1.; rzhlp = 1.; LocSCx = 1.; LocSCy = 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.}}}; 
       LocTV1 = 1.; LocTV10 = 1.; LocTV2 = 1.; LocTV3 = 1.; LocTV4 = 1.; 
       LocTV5 = 1.; LocTV6 = 1.; LocTV7 = 1.; LocTV8 = 1.; LocTV9 = 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 = {{0., 0., 0., 0., 0., 0., 
          0.9961946980917455, 0.08715574274765817, 0., 
          -0.08715574274765817, 0.9961946980917455, 0., 0., 0., 1., 
          0.532, 100., 0.16666666666666669, 0., 1.00026945145708, 0., 
          0., 1., 0., 0., 0., 0., 92010., 92010., 1., 6., 1., 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.16666666666666669, 0., 1.00026945145708, 0., 
          0., 1., 0., 0., 0., 0., 92010., 92010., 2., 6., 2., 0., 1., 
          0., 0., 0., 1., 0.}, {0., 0., 0., 0., 0., 0., 
          0.9998476951563913, 0.01745240643728351, 0., 
          -0.01745240643728351, 0.9998476951563913, 0., 0., 0., 1., 
          0.532, 100., 0.16666666666666669, 0., 1.00026945145708, 0., 
          0., 1., 0., 0., 0., 0., 92010., 92010., 3., 6., 3., 0., 1., 
          0., 0., 0., 1., 0.}, {0., 0., 0., 0., 0., 0., 
          0.9998476951563913, -0.01745240643728351, 0., 
          0.01745240643728351, 0.9998476951563913, 0., 0., 0., 1., 
          0.532, 100., 0.16666666666666669, 0., 1.00026945145708, 0., 
          0., 1., 0., 0., 0., 0., 92010., 92010., 4., 6., 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.16666666666666669, 0., 1.00026945145708, 0., 
          0., 1., 0., 0., 0., 0., 92010., 92010., 5., 6., 5., 0., 1., 
          0., 0., 0., 1., 0.}, {0., 0., 0., 0., 0., 0., 
          0.9961946980917455, -0.08715574274765817, 0., 
          0.08715574274765817, 0.9961946980917455, 0., 0., 0., 1., 
          0.532, 100., 0.16666666666666669, 0., 1.00026945145708, 0., 
          0., 1., 0., 0., 0., 0., 92010., 92010., 6., 6., 6., 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., -9.87654321*^8, 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 = 9.87654321*^8; 
                    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*θ13]] + r2*Re[Sin[
                        0.017453292519943295*θ13]]; yunit = 
                        r2*Re[Cos[0.017453292519943295*θ13]] - 
                        1.*r1*Re[Sin[0.017453292519943295*θ13]]; 
                       zunit = r3*Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^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[x13])*Re[Cos[
                        0.017453292519943295*θ13]] + rey*Re[Sin[
                        0.017453292519943295*θ13]]; ry = 
                        rey*Re[Cos[0.017453292519943295*θ13]] - 
                        1.*(rex - 1.*Re[x13])*Re[Sin[
                        0.017453292519943295*θ13]]; rz = rez*
                        Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^2]; rx = 
                        Chop[rx]; ry = Chop[ry]; rz = Chop[rz]; 
                       terminatedeflection = False; intersections = 
                        {{-1., 0., 0.}}; intersections = Block[{TempVar, 
                        TempVar1, TempVar2, TempVar3, LocRayL, LocSCx, 
                        LocSCy}, TempVar1 = 1.; TempVar2 = 1.; 
                        TempVar3 = 1.; TempVar = {{-1., 0., 0.}}; 
                        LocRayL = 1.; LocSCx = 1.; LocSCy = 1.; 
                        TempVar1 = rx^2 + ry^2 + rz^2; TempVar2 = 
                        rx*xunit + ry*yunit + rz*zunit - r11*xunit; 
                        TempVar3 = -TempVar1 + TempVar2^2 + 2*rx*r11; 
                        TempVar3 = Re[TempVar3]; If[N[Chop[TempVar3]] >= 
                        0., TempVar3 = Sqrt[TempVar3]; TempVar = 
                        {{LocRayL, LocSCx, LocSCy} = {-TempVar2 - 
                        TempVar3, ry - (TempVar2 + TempVar3)*yunit, 
                        rz - (TempVar2 + TempVar3)*zunit}; If[
                        N[Chop[LocRayL]] >= 0. && LocSCx^2 + LocSCy^2 <= 
                        161.29 && Abs[ry + LocRayL*yunit - LocSCx] + 
                        Abs[rz + LocRayL*zunit - LocSCy] + Abs[-r11 + 
                        rx + LocRayL*xunit + Sqrt[r11^2 - LocSCx^2 - 
                        LocSCy^2]*Sign[r11]] < 1.*^-7, {LocRayL, LocSCx, 
                        LocSCy}, {-1., 0., 0.}], {LocRayL, LocSCx, 
                        LocSCy} = {-TempVar2 + TempVar3, ry + 
                        (-TempVar2 + TempVar3)*yunit, rz + (-TempVar2 + 
                        TempVar3)*zunit}; If[N[Chop[LocRayL]] >= 0. && 
                        LocSCx^2 + LocSCy^2 <= 161.29 && Abs[ry + 
                        LocRayL*yunit - LocSCx] + Abs[rz + LocRayL*
                        zunit - LocSCy] + Abs[-r11 + rx + LocRayL*
                        xunit + Sqrt[r11^2 - LocSCx^2 - LocSCy^2]*Sign[
                        r11]] < 1.*^-7, {LocRayL, LocSCx, LocSCy}, 
                        {-1., 0., 0.}]}; If[TempVar[[1]] == TempVar[[
                        2]], {TempVar[[1]]}, TempVar], {{-1., 0., 
                        0.}}]]; AppendTo[intersections, {-1., 0., 0.}]; 
                       holeQ = False; If[ !(Chop[surfaceid - 92155.] == 
                        0. || N[Chop[ray1]] < 0.0001) || 
                        generationnumber == 1., minresonateraylength = 
                        0., minresonateraylength = 0.0001]; Block[
                        {LocRayL, LocSCx, LocSCy}, LocRayL = 1.; 
                        LocSCx = 1.; LocSCy = 1.; dummymatrix = Union[
                        Select[intersections, ({LocRayL, LocSCx, 
                        LocSCy} = #1; {LocRayL, LocSCx, LocSCy} == 
                        {-1., 0., 0.} || N[Chop[LocRayL]] >= 
                        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*θ13]] + r2*Re[Sin[
                        0.017453292519943295*θ13]]; yunit = 
                        r2*Re[Cos[0.017453292519943295*θ13]] - 
                        1.*r1*Re[Sin[0.017453292519943295*θ13]]; 
                       zunit = r3*Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^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*θ13]]*(rex - 
                        1.*Re[x13 + 2.1*Cos[0.017453292519943295*
                        θ13]]) + (rey - 2.1*Re[Sin[
                        0.017453292519943295*θ13]])*Re[Sin[
                        0.017453292519943295*θ13]]; ry = 
                        Re[Cos[0.017453292519943295*θ13]]*(rey - 
                        2.1*Re[Sin[0.017453292519943295*θ13]]) - 
                        1.*(rex - 1.*Re[x13 + 2.1*Cos[
                        0.017453292519943295*θ13]])*Re[Sin[
                        0.017453292519943295*θ13]]; rz = rez*
                        Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^2]; rx = 
                        Chop[rx]; ry = Chop[ry]; rz = Chop[rz]; 
                       terminatedeflection = False; intersections = 
                        {{-1., 0., 0.}}; intersections = Block[{TempVar, 
                        TempVar1, TempVar2, TempVar3, LocRayL, LocSCx, 
                        LocSCy}, TempVar1 = 1.; TempVar2 = 1.; 
                        TempVar3 = 1.; TempVar = {{-1., 0., 0.}}; 
                        LocRayL = 1.; LocSCx = 1.; LocSCy = 1.; 
                        TempVar1 = rx^2 + ry^2 + rz^2; TempVar2 = 
                        rx*xunit + ry*yunit + rz*zunit - r12*xunit; 
                        TempVar3 = -TempVar1 + TempVar2^2 + 2*rx*r12; 
                        TempVar3 = Re[TempVar3]; If[N[Chop[TempVar3]] >= 
                        0., TempVar3 = Sqrt[TempVar3]; TempVar = 
                        {{LocRayL, LocSCx, LocSCy} = {-TempVar2 - 
                        TempVar3, ry - (TempVar2 + TempVar3)*yunit, 
                        rz - (TempVar2 + TempVar3)*zunit}; If[
                        N[Chop[LocRayL]] >= 0. && LocSCx^2 + LocSCy^2 <= 
                        161.29 && Abs[ry + LocRayL*yunit - LocSCx] + 
                        Abs[rz + LocRayL*zunit - LocSCy] + Abs[-r12 + 
                        rx + LocRayL*xunit + Sqrt[r12^2 - LocSCx^2 - 
                        LocSCy^2]*Sign[r12]] < 1.*^-7, {LocRayL, LocSCx, 
                        LocSCy}, {-1., 0., 0.}], {LocRayL, LocSCx, 
                        LocSCy} = {-TempVar2 + TempVar3, ry + 
                        (-TempVar2 + TempVar3)*yunit, rz + (-TempVar2 + 
                        TempVar3)*zunit}; If[N[Chop[LocRayL]] >= 0. && 
                        LocSCx^2 + LocSCy^2 <= 161.29 && Abs[ry + 
                        LocRayL*yunit - LocSCx] + Abs[rz + LocRayL*
                        zunit - LocSCy] + Abs[-r12 + rx + LocRayL*
                        xunit + Sqrt[r12^2 - LocSCx^2 - LocSCy^2]*Sign[
                        r12]] < 1.*^-7, {LocRayL, LocSCx, LocSCy}, 
                        {-1., 0., 0.}]}; If[TempVar[[1]] == TempVar[[
                        2]], {TempVar[[1]]}, TempVar], {{-1., 0., 
                        0.}}]]; AppendTo[intersections, {-1., 0., 0.}]; 
                       holeQ = False; If[ !(Chop[surfaceid - 92156.] == 
                        0. || N[Chop[ray1]] < 0.0001) || 
                        generationnumber == 1., minresonateraylength = 
                        0., minresonateraylength = 0.0001]; Block[
                        {LocRayL, LocSCx, LocSCy}, LocRayL = 1.; 
                        LocSCx = 1.; LocSCy = 1.; dummymatrix = Union[
                        Select[intersections, ({LocRayL, LocSCx, 
                        LocSCy} = #1; {LocRayL, LocSCx, LocSCy} == 
                        {-1., 0., 0.} || N[Chop[LocRayL]] >= 
                        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*θ13]] + r2*Re[Sin[
                        0.017453292519943295*θ13]]; yunit = 
                        r2*Re[Cos[0.017453292519943295*θ13]] - 
                        1.*r1*Re[Sin[0.017453292519943295*θ13]]; 
                       zunit = r3*Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^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*θ13]]*(rex - 
                        1.*Re[x13 + 9.6*Cos[0.017453292519943295*
                        θ13]]) + (rey - 9.6*Re[Sin[
                        0.017453292519943295*θ13]])*Re[Sin[
                        0.017453292519943295*θ13]]; ry = 
                        Re[Cos[0.017453292519943295*θ13]]*(rey - 
                        9.6*Re[Sin[0.017453292519943295*θ13]]) - 
                        1.*(rex - 1.*Re[x13 + 9.6*Cos[
                        0.017453292519943295*θ13]])*Re[Sin[
                        0.017453292519943295*θ13]]; rz = rez*
                        Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^2]; rx = 
                        Chop[rx]; ry = Chop[ry]; rz = Chop[rz]; 
                       terminatedeflection = False; intersections = 
                        {{-1., 0., 0.}}; intersections = Block[{TempVar, 
                        TempVar1, TempVar2, TempVar3, LocRayL, LocSCx, 
                        LocSCy}, TempVar1 = 1.; TempVar2 = 1.; 
                        TempVar3 = 1.; TempVar = {{-1., 0., 0.}}; 
                        LocRayL = 1.; LocSCx = 1.; LocSCy = 1.; 
                        TempVar1 = rx^2 + ry^2 + rz^2; TempVar2 = 
                        rx*xunit + ry*yunit + rz*zunit - r13*xunit; 
                        TempVar3 = -TempVar1 + TempVar2^2 + 2*rx*r13; 
                        TempVar3 = Re[TempVar3]; If[N[Chop[TempVar3]] >= 
                        0., TempVar3 = Sqrt[TempVar3]; TempVar = 
                        {{LocRayL, LocSCx, LocSCy} = {-TempVar2 - 
                        TempVar3, ry - (TempVar2 + TempVar3)*yunit, 
                        rz - (TempVar2 + TempVar3)*zunit}; If[
                        N[Chop[LocRayL]] >= 0. && LocSCx^2 + LocSCy^2 <= 
                        161.29 && Abs[ry + LocRayL*yunit - LocSCx] + 
                        Abs[rz + LocRayL*zunit - LocSCy] + Abs[-r13 + 
                        rx + LocRayL*xunit + Sqrt[r13^2 - LocSCx^2 - 
                        LocSCy^2]*Sign[r13]] < 1.*^-7, {LocRayL, LocSCx, 
                        LocSCy}, {-1., 0., 0.}], {LocRayL, LocSCx, 
                        LocSCy} = {-TempVar2 + TempVar3, ry + 
                        (-TempVar2 + TempVar3)*yunit, rz + (-TempVar2 + 
                        TempVar3)*zunit}; If[N[Chop[LocRayL]] >= 0. && 
                        LocSCx^2 + LocSCy^2 <= 161.29 && Abs[ry + 
                        LocRayL*yunit - LocSCx] + Abs[rz + LocRayL*
                        zunit - LocSCy] + Abs[-r13 + rx + LocRayL*
                        xunit + Sqrt[r13^2 - LocSCx^2 - LocSCy^2]*Sign[
                        r13]] < 1.*^-7, {LocRayL, LocSCx, LocSCy}, 
                        {-1., 0., 0.}]}; If[TempVar[[1]] == TempVar[[
                        2]], {TempVar[[1]]}, TempVar], {{-1., 0., 
                        0.}}]]; AppendTo[intersections, {-1., 0., 0.}]; 
                       holeQ = False; If[ !(Chop[surfaceid - 92157.] == 
                        0. || N[Chop[ray1]] < 0.0001) || 
                        generationnumber == 1., minresonateraylength = 
                        0., minresonateraylength = 0.0001]; Block[
                        {LocRayL, LocSCx, LocSCy}, LocRayL = 1.; 
                        LocSCx = 1.; LocSCy = 1.; dummymatrix = Union[
                        Select[intersections, ({LocRayL, LocSCx, 
                        LocSCy} = #1; {LocRayL, LocSCx, LocSCy} == 
                        {-1., 0., 0.} || N[Chop[LocRayL]] >= 
                        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*θ23]] + r2*Re[Sin[
                        0.017453292519943295*θ23]]; yunit = 
                        r2*Re[Cos[0.017453292519943295*θ23]] - 
                        1.*r1*Re[Sin[0.017453292519943295*θ23]]; 
                       zunit = r3*Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^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[x23])*Re[Cos[
                        0.017453292519943295*θ23]] + rey*Re[Sin[
                        0.017453292519943295*θ23]]; ry = 
                        rey*Re[Cos[0.017453292519943295*θ23]] - 
                        1.*(rex - 1.*Re[x23])*Re[Sin[
                        0.017453292519943295*θ23]]; rz = rez*
                        Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^2]; rx = 
                        Chop[rx]; ry = Chop[ry]; rz = Chop[rz]; 
                       terminatedeflection = False; intersections = 
                        {{-1., 0., 0.}}; intersections = Block[{TempVar, 
                        TempVar1, TempVar2, TempVar3, LocRayL, LocSCx, 
                        LocSCy}, TempVar1 = 1.; TempVar2 = 1.; 
                        TempVar3 = 1.; TempVar = {{-1., 0., 0.}}; 
                        LocRayL = 1.; LocSCx = 1.; LocSCy = 1.; 
                        TempVar1 = rx^2 + ry^2 + rz^2; TempVar2 = 
                        rx*xunit + ry*yunit + rz*zunit - r21*xunit; 
                        TempVar3 = -TempVar1 + TempVar2^2 + 2*rx*r21; 
                        TempVar3 = Re[TempVar3]; If[N[Chop[TempVar3]] >= 
                        0., TempVar3 = Sqrt[TempVar3]; TempVar = 
                        {{LocRayL, LocSCx, LocSCy} = {-TempVar2 - 
                        TempVar3, ry - (TempVar2 + TempVar3)*yunit, 
                        rz - (TempVar2 + TempVar3)*zunit}; If[
                        N[Chop[LocRayL]] >= 0. && LocSCx^2 + LocSCy^2 <= 
                        161.29 && Abs[ry + LocRayL*yunit - LocSCx] + 
                        Abs[rz + LocRayL*zunit - LocSCy] + Abs[-r21 + 
                        rx + LocRayL*xunit + Sqrt[r21^2 - LocSCx^2 - 
                        LocSCy^2]*Sign[r21]] < 1.*^-7, {LocRayL, LocSCx, 
                        LocSCy}, {-1., 0., 0.}], {LocRayL, LocSCx, 
                        LocSCy} = {-TempVar2 + TempVar3, ry + 
                        (-TempVar2 + TempVar3)*yunit, rz + (-TempVar2 + 
                        TempVar3)*zunit}; If[N[Chop[LocRayL]] >= 0. && 
                        LocSCx^2 + LocSCy^2 <= 161.29 && Abs[ry + 
                        LocRayL*yunit - LocSCx] + Abs[rz + LocRayL*
                        zunit - LocSCy] + Abs[-r21 + rx + LocRayL*
                        xunit + Sqrt[r21^2 - LocSCx^2 - LocSCy^2]*Sign[
                        r21]] < 1.*^-7, {LocRayL, LocSCx, LocSCy}, 
                        {-1., 0., 0.}]}; If[TempVar[[1]] == TempVar[[
                        2]], {TempVar[[1]]}, TempVar], {{-1., 0., 
                        0.}}]]; AppendTo[intersections, {-1., 0., 0.}]; 
                       holeQ = False; If[ !(Chop[surfaceid - 92310.] == 
                        0. || N[Chop[ray1]] < 0.0001) || 
                        generationnumber == 1., minresonateraylength = 
                        0., minresonateraylength = 0.0001]; Block[
                        {LocRayL, LocSCx, LocSCy}, LocRayL = 1.; 
                        LocSCx = 1.; LocSCy = 1.; dummymatrix = Union[
                        Select[intersections, ({LocRayL, LocSCx, 
                        LocSCy} = #1; {LocRayL, LocSCx, LocSCy} == 
                        {-1., 0., 0.} || N[Chop[LocRayL]] >= 
                        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*θ23]] + r2*Re[Sin[
                        0.017453292519943295*θ23]]; yunit = 
                        r2*Re[Cos[0.017453292519943295*θ23]] - 
                        1.*r1*Re[Sin[0.017453292519943295*θ23]]; 
                       zunit = r3*Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^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*θ23]]*(rex - 
                        1.*Re[x23 + 7.5*Cos[0.017453292519943295*
                        θ23]]) + (rey - 7.5*Re[Sin[
                        0.017453292519943295*θ23]])*Re[Sin[
                        0.017453292519943295*θ23]]; ry = 
                        Re[Cos[0.017453292519943295*θ23]]*(rey - 
                        7.5*Re[Sin[0.017453292519943295*θ23]]) - 
                        1.*(rex - 1.*Re[x23 + 7.5*Cos[
                        0.017453292519943295*θ23]])*Re[Sin[
                        0.017453292519943295*θ23]]; rz = rez*
                        Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^2]; rx = 
                        Chop[rx]; ry = Chop[ry]; rz = Chop[rz]; 
                       terminatedeflection = False; intersections = 
                        {{-1., 0., 0.}}; intersections = Block[{TempVar, 
                        TempVar1, TempVar2, TempVar3, LocRayL, LocSCx, 
                        LocSCy}, TempVar1 = 1.; TempVar2 = 1.; 
                        TempVar3 = 1.; TempVar = {{-1., 0., 0.}}; 
                        LocRayL = 1.; LocSCx = 1.; LocSCy = 1.; 
                        TempVar1 = rx^2 + ry^2 + rz^2; TempVar2 = 
                        rx*xunit + ry*yunit + rz*zunit - r22*xunit; 
                        TempVar3 = -TempVar1 + TempVar2^2 + 2*rx*r22; 
                        TempVar3 = Re[TempVar3]; If[N[Chop[TempVar3]] >= 
                        0., TempVar3 = Sqrt[TempVar3]; TempVar = 
                        {{LocRayL, LocSCx, LocSCy} = {-TempVar2 - 
                        TempVar3, ry - (TempVar2 + TempVar3)*yunit, 
                        rz - (TempVar2 + TempVar3)*zunit}; If[
                        N[Chop[LocRayL]] >= 0. && LocSCx^2 + LocSCy^2 <= 
                        161.29 && Abs[ry + LocRayL*yunit - LocSCx] + 
                        Abs[rz + LocRayL*zunit - LocSCy] + Abs[-r22 + 
                        rx + LocRayL*xunit + Sqrt[r22^2 - LocSCx^2 - 
                        LocSCy^2]*Sign[r22]] < 1.*^-7, {LocRayL, LocSCx, 
                        LocSCy}, {-1., 0., 0.}], {LocRayL, LocSCx, 
                        LocSCy} = {-TempVar2 + TempVar3, ry + 
                        (-TempVar2 + TempVar3)*yunit, rz + (-TempVar2 + 
                        TempVar3)*zunit}; If[N[Chop[LocRayL]] >= 0. && 
                        LocSCx^2 + LocSCy^2 <= 161.29 && Abs[ry + 
                        LocRayL*yunit - LocSCx] + Abs[rz + LocRayL*
                        zunit - LocSCy] + Abs[-r22 + rx + LocRayL*
                        xunit + Sqrt[r22^2 - LocSCx^2 - LocSCy^2]*Sign[
                        r22]] < 1.*^-7, {LocRayL, LocSCx, LocSCy}, 
                        {-1., 0., 0.}]}; If[TempVar[[1]] == TempVar[[
                        2]], {TempVar[[1]]}, TempVar], {{-1., 0., 
                        0.}}]]; AppendTo[intersections, {-1., 0., 0.}]; 
                       holeQ = False; If[ !(Chop[surfaceid - 92311.] == 
                        0. || N[Chop[ray1]] < 0.0001) || 
                        generationnumber == 1., minresonateraylength = 
                        0., minresonateraylength = 0.0001]; Block[
                        {LocRayL, LocSCx, LocSCy}, LocRayL = 1.; 
                        LocSCx = 1.; LocSCy = 1.; dummymatrix = Union[
                        Select[intersections, ({LocRayL, LocSCx, 
                        LocSCy} = #1; {LocRayL, LocSCx, LocSCy} == 
                        {-1., 0., 0.} || N[Chop[LocRayL]] >= 
                        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*θ23]] + r2*Re[Sin[
                        0.017453292519943295*θ23]]; yunit = 
                        r2*Re[Cos[0.017453292519943295*θ23]] - 
                        1.*r1*Re[Sin[0.017453292519943295*θ23]]; 
                       zunit = r3*Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^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*θ23]]*(rex - 
                        1.*Re[x23 + 9.6*Cos[0.017453292519943295*
                        θ23]]) + (rey - 9.6*Re[Sin[
                        0.017453292519943295*θ23]])*Re[Sin[
                        0.017453292519943295*θ23]]; ry = 
                        Re[Cos[0.017453292519943295*θ23]]*(rey - 
                        9.6*Re[Sin[0.017453292519943295*θ23]]) - 
                        1.*(rex - 1.*Re[x23 + 9.6*Cos[
                        0.017453292519943295*θ23]])*Re[Sin[
                        0.017453292519943295*θ23]]; rz = rez*
                        Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^2]; rx = 
                        Chop[rx]; ry = Chop[ry]; rz = Chop[rz]; 
                       terminatedeflection = False; intersections = 
                        {{-1., 0., 0.}}; intersections = Block[{TempVar, 
                        TempVar1, TempVar2, TempVar3, LocRayL, LocSCx, 
                        LocSCy}, TempVar1 = 1.; TempVar2 = 1.; 
                        TempVar3 = 1.; TempVar = {{-1., 0., 0.}}; 
                        LocRayL = 1.; LocSCx = 1.; LocSCy = 1.; 
                        TempVar1 = rx^2 + ry^2 + rz^2; TempVar2 = 
                        rx*xunit + ry*yunit + rz*zunit - r23*xunit; 
                        TempVar3 = -TempVar1 + TempVar2^2 + 2*rx*r23; 
                        TempVar3 = Re[TempVar3]; If[N[Chop[TempVar3]] >= 
                        0., TempVar3 = Sqrt[TempVar3]; TempVar = 
                        {{LocRayL, LocSCx, LocSCy} = {-TempVar2 - 
                        TempVar3, ry - (TempVar2 + TempVar3)*yunit, 
                        rz - (TempVar2 + TempVar3)*zunit}; If[
                        N[Chop[LocRayL]] >= 0. && LocSCx^2 + LocSCy^2 <= 
                        161.29 && Abs[ry + LocRayL*yunit - LocSCx] + 
                        Abs[rz + LocRayL*zunit - LocSCy] + Abs[-r23 + 
                        rx + LocRayL*xunit + Sqrt[r23^2 - LocSCx^2 - 
                        LocSCy^2]*Sign[r23]] < 1.*^-7, {LocRayL, LocSCx, 
                        LocSCy}, {-1., 0., 0.}], {LocRayL, LocSCx, 
                        LocSCy} = {-TempVar2 + TempVar3, ry + 
                        (-TempVar2 + TempVar3)*yunit, rz + (-TempVar2 + 
                        TempVar3)*zunit}; If[N[Chop[LocRayL]] >= 0. && 
                        LocSCx^2 + LocSCy^2 <= 161.29 && Abs[ry + 
                        LocRayL*yunit - LocSCx] + Abs[rz + LocRayL*
                        zunit - LocSCy] + Abs[-r23 + rx + LocRayL*
                        xunit + Sqrt[r23^2 - LocSCx^2 - LocSCy^2]*Sign[
                        r23]] < 1.*^-7, {LocRayL, LocSCx, LocSCy}, 
                        {-1., 0., 0.}]}; If[TempVar[[1]] == TempVar[[
                        2]], {TempVar[[1]]}, TempVar], {{-1., 0., 
                        0.}}]]; AppendTo[intersections, {-1., 0., 0.}]; 
                       holeQ = False; If[ !(Chop[surfaceid - 92312.] == 
                        0. || N[Chop[ray1]] < 0.0001) || 
                        generationnumber == 1., minresonateraylength = 
                        0., minresonateraylength = 0.0001]; Block[
                        {LocRayL, LocSCx, LocSCy}, LocRayL = 1.; 
                        LocSCx = 1.; LocSCy = 1.; dummymatrix = Union[
                        Select[intersections, ({LocRayL, LocSCx, 
                        LocSCy} = #1; {LocRayL, LocSCx, LocSCy} == 
                        {-1., 0., 0.} || N[Chop[LocRayL]] >= 
                        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*θ13]] + r2*Re[Sin[
                        0.017453292519943295*θ13]]; yunit = 
                      r2*Re[Cos[0.017453292519943295*θ13]] - 1.*r1*
                        Re[Sin[0.017453292519943295*θ13]]; zunit = 
                      r3*Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^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[x13])*Re[Cos[0.017453292519943295*
                        θ13]] + rey*Re[Sin[0.017453292519943295*
                        θ13]]; ry = rey*Re[Cos[0.017453292519943295*
                        θ13]] - 1.*(rex - 1.*Re[x13])*Re[Sin[
                        0.017453292519943295*θ13]]; rz = rez*
                       Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^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., 
                       92155.}; {componentnumber, surfacenumber, 
                       surfaceid} = {1., 1., 92155.}; {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[x13] + rxhlp*Re[Cos[0.017453292519943295*
                        θ13]] - 1.*ryhlp*Re[Sin[0.017453292519943295*
                        θ13]]; rey = ryhlp*Re[Cos[0.017453292519943295*
                        θ13]] + rxhlp*Re[Sin[0.017453292519943295*
                        θ13]]; rez = rzhlp*Re[Cos[0.017453292519943295*
                        θ13]^2 + Sin[0.017453292519943295*θ13]^2]; 
                     rex = Chop[rex]; rey = Chop[rey]; rez = Chop[rez]; 
                     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.*((Sign[r11]*ray2)/
                        (r11^2 - 1.*ray2^2 - 1.*ray3^2)^0.5), -1.*
                        ((Sign[r11]*ray3)/(r11^2 - 1.*ray2^2 - 1.*ray3^
                        2)^0.5)}, {(Sign[r11]*ray2)/(r11^2 - 1.*ray2^2 - 
                        1.*ray3^2)^0.5, 1., 0.}, {(Sign[r11]*ray3)/
                        (r11^2 - 1.*ray2^2 - 1.*ray3^2)^0.5, 0., 1.}}]; 
                      dummymatrix = ({{Re[Cos[0.017453292519943295*
                        θ13]], -1.*Re[Sin[0.017453292519943295*
                        θ13]], 0.}, {Re[Sin[0.017453292519943295*
                        θ13]], Re[Cos[0.017453292519943295*θ13]], 
                        0.}, {0., 0., 1.}} . #1 & ) /@ surfnorm; 
                      surfnorm = (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; 
                      surfacerotationmatrix = {{Re[Cos[
                        0.017453292519943295*θ13]], Re[Sin[
                        0.017453292519943295*θ13]], 0.}, 
                        {-1.*Re[Sin[0.017453292519943295*θ13]], 
                        Re[Cos[0.017453292519943295*θ13]], 0.}, 
                        {0., 0., 1.}}; surfacetranslationvector = 
                       {Re[x13], 0., 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 = 
                       {{I*0., I*0., I*0.}}; 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.7105646 - 
                        0.009121199600000001*wavelength^2 + 0.027760538/
                        wavelength^2 + 0.0012739656/wavelength^4 - 
                        0.000047889342/wavelength^6 + 8.0562028*^-6/
                        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.7105646 - 
                        0.009121199600000001*wavelength^2 + 0.027760538/
                        wavelength^2 + 0.0012739656/wavelength^4 - 
                        0.000047889342/wavelength^6 + 8.0562028*^-6/
                        wavelength^8]; If[Abs[alternateindex - index1] > 
                        Abs[index2 - index1], index1 = index2; index2 = 
                        alternateindex, index1 = alternateindex], 
                        index1 = Sqrt[2.7105646 - 0.009121199600000001*
                        wavelength^2 + 0.027760538/wavelength^2 + 
                        0.0012739656/wavelength^4 - 0.000047889342/
                        wavelength^6 + 8.0562028*^-6/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[previoustilt, {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 = 
                        previoustilt . {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[previoustilt . {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 & ) /@ 
                        originalrotation][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 != {{I*0., I*0., I*0.}}, 
                        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., -9.87654321*^8, 
                        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*θ13]] + r2*Re[Sin[
                        0.017453292519943295*θ13]]; yunit = 
                      r2*Re[Cos[0.017453292519943295*θ13]] - 1.*r1*
                        Re[Sin[0.017453292519943295*θ13]]; zunit = 
                      r3*Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^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*θ13]]*(rex - 
                        1.*Re[x13 + 2.1*Cos[0.017453292519943295*
                        θ13]]) + (rey - 2.1*Re[Sin[
                        0.017453292519943295*θ13]])*Re[Sin[
                        0.017453292519943295*θ13]]; ry = 
                      Re[Cos[0.017453292519943295*θ13]]*(rey - 
                        2.1*Re[Sin[0.017453292519943295*θ13]]) - 
                       1.*(rex - 1.*Re[x13 + 2.1*Cos[
                        0.017453292519943295*θ13]])*Re[Sin[
                        0.017453292519943295*θ13]]; rz = rez*
                       Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^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., 
                       92156.}; {componentnumber, surfacenumber, 
                       surfaceid} = {1., 2., 92156.}; {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*θ13]] + 
                       Re[x13 + 2.1*Cos[0.017453292519943295*θ13]] - 
                       1.*ryhlp*Re[Sin[0.017453292519943295*θ13]]; 
                     rey = ryhlp*Re[Cos[0.017453292519943295*θ13]] + 
                       2.1*Re[Sin[0.017453292519943295*θ13]] + 
                       rxhlp*Re[Sin[0.017453292519943295*θ13]]; 
                     rez = rzhlp*Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^2]; rex = 
                      Chop[rex]; rey = Chop[rey]; rez = Chop[rez]; 
                     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.*((Sign[r12]*ray2)/
                        (r12^2 - 1.*ray2^2 - 1.*ray3^2)^0.5), -1.*
                        ((Sign[r12]*ray3)/(r12^2 - 1.*ray2^2 - 1.*ray3^
                        2)^0.5)}, {(Sign[r12]*ray2)/(r12^2 - 1.*ray2^2 - 
                        1.*ray3^2)^0.5, 1., 0.}, {(Sign[r12]*ray3)/
                        (r12^2 - 1.*ray2^2 - 1.*ray3^2)^0.5, 0., 1.}}]; 
                      dummymatrix = ({{Re[Cos[0.017453292519943295*
                        θ13]], -1.*Re[Sin[0.017453292519943295*
                        θ13]], 0.}, {Re[Sin[0.017453292519943295*
                        θ13]], Re[Cos[0.017453292519943295*θ13]], 
                        0.}, {0., 0., 1.}} . #1 & ) /@ surfnorm; 
                      surfnorm = (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; 
                      surfacerotationmatrix = {{Re[Cos[
                        0.017453292519943295*θ13]], Re[Sin[
                        0.017453292519943295*θ13]], 0.}, 
                        {-1.*Re[Sin[0.017453292519943295*θ13]], 
                        Re[Cos[0.017453292519943295*θ13]], 0.}, 
                        {0., 0., 1.}}; surfacetranslationvector = 
                       {Re[x13 + 2.1*Cos[0.017453292519943295*θ13]], 
                        2.1*Re[Sin[0.017453292519943295*θ13]], 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 = 
                       {{I*0., I*0., I*0.}}; 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[
                        2.7105646 - 0.009121199600000001*wavelength^2 + 
                        0.027760538/wavelength^2 + 0.0012739656/
                        wavelength^4 - 0.000047889342/wavelength^6 + 
                        8.0562028*^-6/wavelength^8]; If[Abs[
                        alternateindex - index1] > Abs[index2 - index1], 
                        index1 = index2; index2 = alternateindex, 
                        index1 = alternateindex], index1 = Sqrt[
                        2.7105646 - 0.009121199600000001*wavelength^2 + 
                        0.027760538/wavelength^2 + 0.0012739656/
                        wavelength^4 - 0.000047889342/wavelength^6 + 
                        8.0562028*^-6/wavelength^8]], index2 = Sqrt[
                        2.7105646 - 0.009121199600000001*wavelength^2 + 
                        0.027760538/wavelength^2 + 0.0012739656/
                        wavelength^4 - 0.000047889342/wavelength^6 + 
                        8.0562028*^-6/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[previoustilt, {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 = 
                        previoustilt . {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[previoustilt . {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 & ) /@ 
                        originalrotation][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 != {{I*0., I*0., I*0.}}, 
                        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., -9.87654321*^8, 
                        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*θ13]] + r2*Re[Sin[
                        0.017453292519943295*θ13]]; yunit = 
                      r2*Re[Cos[0.017453292519943295*θ13]] - 1.*r1*
                        Re[Sin[0.017453292519943295*θ13]]; zunit = 
                      r3*Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^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*θ13]]*(rex - 
                        1.*Re[x13 + 9.6*Cos[0.017453292519943295*
                        θ13]]) + (rey - 9.6*Re[Sin[
                        0.017453292519943295*θ13]])*Re[Sin[
                        0.017453292519943295*θ13]]; ry = 
                      Re[Cos[0.017453292519943295*θ13]]*(rey - 
                        9.6*Re[Sin[0.017453292519943295*θ13]]) - 
                       1.*(rex - 1.*Re[x13 + 9.6*Cos[
                        0.017453292519943295*θ13]])*Re[Sin[
                        0.017453292519943295*θ13]]; rz = rez*
                       Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^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., 
                       92157.}; {componentnumber, surfacenumber, 
                       surfaceid} = {1., 3., 92157.}; {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*θ13]] + 
                       Re[x13 + 9.6*Cos[0.017453292519943295*θ13]] - 
                       1.*ryhlp*Re[Sin[0.017453292519943295*θ13]]; 
                     rey = ryhlp*Re[Cos[0.017453292519943295*θ13]] + 
                       9.6*Re[Sin[0.017453292519943295*θ13]] + 
                       rxhlp*Re[Sin[0.017453292519943295*θ13]]; 
                     rez = rzhlp*Re[Cos[0.017453292519943295*θ13]^2 + 
                        Sin[0.017453292519943295*θ13]^2]; rex = 
                      Chop[rex]; rey = Chop[rey]; rez = Chop[rez]; 
                     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.*((Sign[r13]*ray2)/
                        (r13^2 - 1.*ray2^2 - 1.*ray3^2)^0.5), -1.*
                        ((Sign[r13]*ray3)/(r13^2 - 1.*ray2^2 - 1.*ray3^
                        2)^0.5)}, {(Sign[r13]*ray2)/(r13^2 - 1.*ray2^2 - 
                        1.*ray3^2)^0.5, 1., 0.}, {(Sign[r13]*ray3)/
                        (r13^2 - 1.*ray2^2 - 1.*ray3^2)^0.5, 0., 1.}}]; 
                      dummymatrix = ({{Re[Cos[0.017453292519943295*
                        θ13]], -1.*Re[Sin[0.017453292519943295*
                        θ13]], 0.}, {Re[Sin[0.017453292519943295*
                        θ13]], Re[Cos[0.017453292519943295*θ13]], 
                        0.}, {0., 0., 1.}} . #1 & ) /@ surfnorm; 
                      surfnorm = (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; 
                      surfacerotationmatrix = {{Re[Cos[
                        0.017453292519943295*θ13]], Re[Sin[
                        0.017453292519943295*θ13]], 0.}, 
                        {-1.*Re[Sin[0.017453292519943295*θ13]], 
                        Re[Cos[0.017453292519943295*θ13]], 0.}, 
                        {0., 0., 1.}}; surfacetranslationvector = 
                       {Re[x13 + 9.6*Cos[0.017453292519943295*θ13]], 
                        9.6*Re[Sin[0.017453292519943295*θ13]], 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 = 
                       {{I*0., I*0., I*0.}}; 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[previoustilt, {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 = 
                        previoustilt . {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[previoustilt . {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 & ) /@ 
                        originalrotation][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 != {{I*0., I*0., I*0.}}, 
                        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., -9.87654321*^8, 
                        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*θ23]] + r2*Re[Sin[
                        0.017453292519943295*θ23]]; yunit = 
                      r2*Re[Cos[0.017453292519943295*θ23]] - 1.*r1*
                        Re[Sin[0.017453292519943295*θ23]]; zunit = 
                      r3*Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^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[x23])*Re[Cos[0.017453292519943295*
                        θ23]] + rey*Re[Sin[0.017453292519943295*
                        θ23]]; ry = rey*Re[Cos[0.017453292519943295*
                        θ23]] - 1.*(rex - 1.*Re[x23])*Re[Sin[
                        0.017453292519943295*θ23]]; rz = rez*
                       Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^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., 
                       92310.}; {componentnumber, surfacenumber, 
                       surfaceid} = {1., 4., 92310.}; {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[x23] + rxhlp*Re[Cos[0.017453292519943295*
                        θ23]] - 1.*ryhlp*Re[Sin[0.017453292519943295*
                        θ23]]; rey = ryhlp*Re[Cos[0.017453292519943295*
                        θ23]] + rxhlp*Re[Sin[0.017453292519943295*
                        θ23]]; rez = rzhlp*Re[Cos[0.017453292519943295*
                        θ23]^2 + Sin[0.017453292519943295*θ23]^2]; 
                     rex = Chop[rex]; rey = Chop[rey]; rez = Chop[rez]; 
                     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.*((Sign[r21]*ray2)/
                        (r21^2 - 1.*ray2^2 - 1.*ray3^2)^0.5), -1.*
                        ((Sign[r21]*ray3)/(r21^2 - 1.*ray2^2 - 1.*ray3^
                        2)^0.5)}, {(Sign[r21]*ray2)/(r21^2 - 1.*ray2^2 - 
                        1.*ray3^2)^0.5, 1., 0.}, {(Sign[r21]*ray3)/
                        (r21^2 - 1.*ray2^2 - 1.*ray3^2)^0.5, 0., 1.}}]; 
                      dummymatrix = ({{Re[Cos[0.017453292519943295*
                        θ23]], -1.*Re[Sin[0.017453292519943295*
                        θ23]], 0.}, {Re[Sin[0.017453292519943295*
                        θ23]], Re[Cos[0.017453292519943295*θ23]], 
                        0.}, {0., 0., 1.}} . #1 & ) /@ surfnorm; 
                      surfnorm = (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; 
                      surfacerotationmatrix = {{Re[Cos[
                        0.017453292519943295*θ23]], Re[Sin[
                        0.017453292519943295*θ23]], 0.}, 
                        {-1.*Re[Sin[0.017453292519943295*θ23]], 
                        Re[Cos[0.017453292519943295*θ23]], 0.}, 
                        {0., 0., 1.}}; surfacetranslationvector = 
                       {Re[x23], 0., 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 = 
                       {{I*0., I*0., I*0.}}; 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[previoustilt, {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 = 
                        previoustilt . {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[previoustilt . {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 & ) /@ 
                        originalrotation][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 != {{I*0., I*0., I*0.}}, 
                        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., -9.87654321*^8, 
                        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 = r1*Re[Cos[
                        0.017453292519943295*θ23]] + r2*Re[Sin[
                        0.017453292519943295*θ23]]; yunit = 
                      r2*Re[Cos[0.017453292519943295*θ23]] - 1.*r1*
                        Re[Sin[0.017453292519943295*θ23]]; zunit = 
                      r3*Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^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*θ23]]*(rex - 
                        1.*Re[x23 + 7.5*Cos[0.017453292519943295*
                        θ23]]) + (rey - 7.5*Re[Sin[
                        0.017453292519943295*θ23]])*Re[Sin[
                        0.017453292519943295*θ23]]; ry = 
                      Re[Cos[0.017453292519943295*θ23]]*(rey - 
                        7.5*Re[Sin[0.017453292519943295*θ23]]) - 
                       1.*(rex - 1.*Re[x23 + 7.5*Cos[
                        0.017453292519943295*θ23]])*Re[Sin[
                        0.017453292519943295*θ23]]; rz = rez*
                       Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^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., 5., 
                       92311.}; {componentnumber, surfacenumber, 
                       surfaceid} = {1., 5., 92311.}; {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*θ23]] + 
                       Re[x23 + 7.5*Cos[0.017453292519943295*θ23]] - 
                       1.*ryhlp*Re[Sin[0.017453292519943295*θ23]]; 
                     rey = ryhlp*Re[Cos[0.017453292519943295*θ23]] + 
                       7.5*Re[Sin[0.017453292519943295*θ23]] + 
                       rxhlp*Re[Sin[0.017453292519943295*θ23]]; 
                     rez = rzhlp*Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^2]; rex = 
                      Chop[rex]; rey = Chop[rey]; rez = Chop[rez]; 
                     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.*((Sign[r22]*ray2)/
                        (r22^2 - 1.*ray2^2 - 1.*ray3^2)^0.5), -1.*
                        ((Sign[r22]*ray3)/(r22^2 - 1.*ray2^2 - 1.*ray3^
                        2)^0.5)}, {(Sign[r22]*ray2)/(r22^2 - 1.*ray2^2 - 
                        1.*ray3^2)^0.5, 1., 0.}, {(Sign[r22]*ray3)/
                        (r22^2 - 1.*ray2^2 - 1.*ray3^2)^0.5, 0., 1.}}]; 
                      dummymatrix = ({{Re[Cos[0.017453292519943295*
                        θ23]], -1.*Re[Sin[0.017453292519943295*
                        θ23]], 0.}, {Re[Sin[0.017453292519943295*
                        θ23]], Re[Cos[0.017453292519943295*θ23]], 
                        0.}, {0., 0., 1.}} . #1 & ) /@ surfnorm; 
                      surfnorm = (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; 
                      surfacerotationmatrix = {{Re[Cos[
                        0.017453292519943295*θ23]], Re[Sin[
                        0.017453292519943295*θ23]], 0.}, 
                        {-1.*Re[Sin[0.017453292519943295*θ23]], 
                        Re[Cos[0.017453292519943295*θ23]], 0.}, 
                        {0., 0., 1.}}; surfacetranslationvector = 
                       {Re[x23 + 7.5*Cos[0.017453292519943295*θ23]], 
                        7.5*Re[Sin[0.017453292519943295*θ23]], 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 = 
                       {{I*0., I*0., I*0.}}; 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.7105646 - 
                        0.009121199600000001*wavelength^2 + 0.027760538/
                        wavelength^2 + 0.0012739656/wavelength^4 - 
                        0.000047889342/wavelength^6 + 8.0562028*^-6/
                        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]], 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[2.7105646 - 0.009121199600000001*wavelength^
                        2 + 0.027760538/wavelength^2 + 0.0012739656/
                        wavelength^4 - 0.000047889342/wavelength^6 + 
                        8.0562028*^-6/wavelength^8]; If[Abs[
                        alternateindex - index1] > Abs[index2 - index1], 
                        index1 = index2; index2 = alternateindex, 
                        index1 = alternateindex], index1 = Sqrt[
                        2.7105646 - 0.009121199600000001*wavelength^2 + 
                        0.027760538/wavelength^2 + 0.0012739656/
                        wavelength^4 - 0.000047889342/wavelength^6 + 
                        8.0562028*^-6/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[previoustilt, {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 = 
                        previoustilt . {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[previoustilt . {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 & ) /@ 
                        originalrotation][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 != {{I*0., I*0., I*0.}}, 
                        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., -9.87654321*^8, 
                        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 = r1*Re[Cos[
                        0.017453292519943295*θ23]] + r2*Re[Sin[
                        0.017453292519943295*θ23]]; yunit = 
                      r2*Re[Cos[0.017453292519943295*θ23]] - 1.*r1*
                        Re[Sin[0.017453292519943295*θ23]]; zunit = 
                      r3*Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^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*θ23]]*(rex - 
                        1.*Re[x23 + 9.6*Cos[0.017453292519943295*
                        θ23]]) + (rey - 9.6*Re[Sin[
                        0.017453292519943295*θ23]])*Re[Sin[
                        0.017453292519943295*θ23]]; ry = 
                      Re[Cos[0.017453292519943295*θ23]]*(rey - 
                        9.6*Re[Sin[0.017453292519943295*θ23]]) - 
                       1.*(rex - 1.*Re[x23 + 9.6*Cos[
                        0.017453292519943295*θ23]])*Re[Sin[
                        0.017453292519943295*θ23]]; rz = rez*
                       Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^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., 6., 
                       92312.}; {componentnumber, surfacenumber, 
                       surfaceid} = {1., 6., 92312.}; {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*θ23]] + 
                       Re[x23 + 9.6*Cos[0.017453292519943295*θ23]] - 
                       1.*ryhlp*Re[Sin[0.017453292519943295*θ23]]; 
                     rey = ryhlp*Re[Cos[0.017453292519943295*θ23]] + 
                       9.6*Re[Sin[0.017453292519943295*θ23]] + 
                       rxhlp*Re[Sin[0.017453292519943295*θ23]]; 
                     rez = rzhlp*Re[Cos[0.017453292519943295*θ23]^2 + 
                        Sin[0.017453292519943295*θ23]^2]; rex = 
                      Chop[rex]; rey = Chop[rey]; rez = Chop[rez]; 
                     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.*((Sign[r23]*ray2)/
                        (r23^2 - 1.*ray2^2 - 1.*ray3^2)^0.5), -1.*
                        ((Sign[r23]*ray3)/(r23^2 - 1.*ray2^2 - 1.*ray3^
                        2)^0.5)}, {(Sign[r23]*ray2)/(r23^2 - 1.*ray2^2 - 
                        1.*ray3^2)^0.5, 1., 0.}, {(Sign[r23]*ray3)/
                        (r23^2 - 1.*ray2^2 - 1.*ray3^2)^0.5, 0., 1.}}]; 
                      dummymatrix = ({{Re[Cos[0.017453292519943295*
                        θ23]], -1.*Re[Sin[0.017453292519943295*
                        θ23]], 0.}, {Re[Sin[0.017453292519943295*
                        θ23]], Re[Cos[0.017453292519943295*θ23]], 
                        0.}, {0., 0., 1.}} . #1 & ) /@ surfnorm; 
                      surfnorm = (#1/Sqrt[#1 . #1] & ) /@ dummymatrix; 
                      surfacerotationmatrix = {{Re[Cos[
                        0.017453292519943295*θ23]], Re[Sin[
                        0.017453292519943295*θ23]], 0.}, 
                        {-1.*Re[Sin[0.017453292519943295*θ23]], 
                        Re[Cos[0.017453292519943295*θ23]], 0.}, 
                        {0., 0., 1.}}; surfacetranslationvector = 
                       {Re[x23 + 9.6*Cos[0.017453292519943295*θ23]], 
                        9.6*Re[Sin[0.017453292519943295*θ23]], 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 = 
                       {{I*0., I*0., I*0.}}; 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.7105646 - 
                        0.009121199600000001*wavelength^2 + 0.027760538/
                        wavelength^2 + 0.0012739656/wavelength^4 - 
                        0.000047889342/wavelength^6 + 8.0562028*^-6/
                        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.7105646 - 
                        0.009121199600000001*wavelength^2 + 0.027760538/
                        wavelength^2 + 0.0012739656/wavelength^4 - 
                        0.000047889342/wavelength^6 + 8.0562028*^-6/
                        wavelength^8]; If[Abs[alternateindex - index1] > 
                        Abs[index2 - index1], index1 = index2; index2 = 
                        alternateindex, index1 = alternateindex], 
                        index1 = Sqrt[2.7105646 - 0.009121199600000001*
                        wavelength^2 + 0.027760538/wavelength^2 + 
                        0.0012739656/wavelength^4 - 0.000047889342/
                        wavelength^6 + 8.0562028*^-6/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[previoustilt, {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 = 
                        previoustilt . {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[previoustilt . {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 & ) /@ 
                        originalrotation][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 != {{I*0., I*0., I*0.}}, 
                        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., -9.87654321*^8, 
                        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., 
                        -9.87654321*^8, 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., -9.87654321*^8, 
                        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., 
                      -9.87654321*^8, 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., -9.87654321*^8, 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., -9.87654321*^8, 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, 
              {{-9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8}}, 
              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., -9.87654321*^8, 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 = 9.87654321*^8; 
                    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]] <= 
                        100. && #1[[2]] >= -100. && #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 - 92440.] == 
                        0. || N[Chop[ray1]] < 0.0001) || 
                        generationnumber == 1., minresonateraylength = 
                        0., minresonateraylength = 0.0001]; Block[
                        {LocRayL, LocSCx, LocSCy}, LocRayL = 1.; 
                        LocSCx = 1.; LocSCy = 1.; dummymatrix = Union[
                        Select[intersections, ({LocRayL, LocSCx, 
                        LocSCy} = #1; {LocRayL, LocSCx, LocSCy} == 
                        {-1., 0., 0.} || N[Chop[LocRayL]] >= 
                        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 = 100. + 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]] >= -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 - 92441.] == 0. || 
                        N[Chop[ray1]] < 0.0001) || generationnumber == 
                        1., minresonateraylength = 0., 
                        minresonateraylength = 0.0001]; Block[{LocRayL, 
                        LocSCx, LocSCy}, LocRayL = 1.; LocSCx = 1.; 
                        LocSCy = 1.; dummymatrix = Union[Select[
                        intersections, ({LocRayL, LocSCx, LocSCy} = #1; 
                        {LocRayL, LocSCx, LocSCy} == {-1., 0., 0.} || 
                        N[Chop[LocRayL]] >= 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]] <= 100. && #1[[2]] >= -100. && 
                        #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 - 92442.] == 0. || 
                        N[Chop[ray1]] < 0.0001) || generationnumber == 
                        1., minresonateraylength = 0., 
                        minresonateraylength = 0.0001]; Block[{LocRayL, 
                        LocSCx, LocSCy}, LocRayL = 1.; LocSCx = 1.; 
                        LocSCy = 1.; dummymatrix = Union[Select[
                        intersections, ({LocRayL, LocSCx, LocSCy} = #1; 
                        {LocRayL, LocSCx, LocSCy} == {-1., 0., 0.} || 
                        N[Chop[LocRayL]] >= 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.*(-100. + 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]] >= -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 - 92443.] == 0. || 
                        N[Chop[ray1]] < 0.0001) || generationnumber == 
                        1., minresonateraylength = 0., 
                        minresonateraylength = 0.0001]; Block[{LocRayL, 
                        LocSCx, LocSCy}, LocRayL = 1.; LocSCx = 1.; 
                        LocSCy = 1.; dummymatrix = Union[Select[
                        intersections, ({LocRayL, LocSCx, LocSCy} = #1; 
                        {LocRayL, LocSCx, LocSCy} == {-1., 0., 0.} || 
                        N[Chop[LocRayL]] >= 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 = 100. + 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]] <= 100. && #1[[2]] >= -100. && 
                        #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 - 92444.] == 0. || 
                        N[Chop[ray1]] < 0.0001) || generationnumber == 
                        1., minresonateraylength = 0., 
                        minresonateraylength = 0.0001]; Block[{LocRayL, 
                        LocSCx, LocSCy}, LocRayL = 1.; LocSCx = 1.; 
                        LocSCy = 1.; dummymatrix = Union[Select[
                        intersections, ({LocRayL, LocSCx, LocSCy} = #1; 
                        {LocRayL, LocSCx, LocSCy} == {-1., 0., 0.} || 
                        N[Chop[LocRayL]] >= 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.*(-100. + 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]] <= 100. && #1[[2]] >= -100. && 
                        #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 - 92445.] == 0. || 
                        N[Chop[ray1]] < 0.0001) || generationnumber == 
                        1., minresonateraylength = 0., 
                        minresonateraylength = 0.0001]; Block[{LocRayL, 
                        LocSCx, LocSCy}, LocRayL = 1.; LocSCx = 1.; 
                        LocSCy = 1.; dummymatrix = Union[Select[
                        intersections, ({LocRayL, LocSCx, LocSCy} = #1; 
                        {LocRayL, LocSCx, LocSCy} == {-1., 0., 0.} || 
                        N[Chop[LocRayL]] >= 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., 92440.}; 
                     {componentnumber, surfacenumber, surfaceid} = 
                      {2., 1., 92440.}; {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]; 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., -9.87654321*^8, 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 = 100. + 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., 92441.}; 
                     {componentnumber, surfacenumber, surfaceid} = 
                      {2., 2., 92441.}; {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 = -100. + rxhlp; 
                     rez = rzhlp; rex = Chop[rex]; rey = Chop[rey]; 
                     rez = Chop[rez]; 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., 
                        -100., 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., -9.87654321*^8, 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., 
                       92442.}; {componentnumber, surfacenumber, 
                       surfaceid} = {2., 3., 92442.}; {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]; 
                     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., -9.87654321*^8, 
                        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.*(-100. + 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., 92443.}; 
                     {componentnumber, surfacenumber, surfaceid} = 
                      {2., 4., 92443.}; {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 = 100. - 1.*rxhlp; 
                     rez = rzhlp; rex = Chop[rex]; rey = Chop[rey]; 
                     rez = Chop[rez]; 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., 
                        100., 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., -9.87654321*^8, 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 = 100. + 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., 92444.}; 
                     {componentnumber, surfacenumber, surfaceid} = 
                      {2., 5., 92444.}; {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 = -100. + 
                       rxhlp; rex = Chop[rex]; rey = Chop[rey]; 
                     rez = Chop[rez]; 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, -100.}; 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., -9.87654321*^8, 
                        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.*(-100. + 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., 
                       92445.}; {componentnumber, surfacenumber, 
                       surfaceid} = {2., 6., 92445.}; {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 = 
                      100. - 1.*rxhlp; rex = Chop[rex]; rey = Chop[rey]; 
                     rez = Chop[rez]; 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, 100.}; 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., -9.87654321*^8, 
                        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., 
                        -9.87654321*^8, 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., -9.87654321*^8, 
                        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., 
                      -9.87654321*^8, 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., -9.87654321*^8, 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., -9.87654321*^8, 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, 
              {{-9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8, 
                -9.87654321*^8, -9.87654321*^8, -9.87654321*^8}}, 
              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]]; polygons11188572008312134018969656 = 
    Compile[{r11, r12, r13}, {{r11 - 1.*Sqrt[-161.29 + r11^2]*Sign[r11], 
       -12.7}, {r11 - 1.*Sqrt[-112.00694444444441 + r11^2]*Sign[r11], 
       -10.583333333333332}, {r11 - 1.*Sqrt[-71.68444444444445 + r11^2]*
         Sign[r11], -8.466666666666667}, 
      {r11 - 1.*Sqrt[-40.322500000000005 + r11^2]*Sign[r11], 
       -6.3500000000000005}, {r11 - 1.*Sqrt[-17.92111111111112 + r11^2]*
         Sign[r11], -4.233333333333334}, 
      {r11 - 1.*Sqrt[-4.48027777777778 + r11^2]*Sign[r11], 
       -2.116666666666667}, 
      {r11 - 1.*Sqrt[-3.1554436208840472*^-30 + r11^2]*Sign[r11], 
       -1.7763568394002505*^-15}, 
      {r11 - 1.*Sqrt[-4.480277777777765 + r11^2]*Sign[r11], 
       2.1166666666666636}, {r11 - 1.*Sqrt[-17.921111111111088 + r11^2]*
         Sign[r11], 4.233333333333331}, 
      {r11 - 1.*Sqrt[-40.32249999999997 + r11^2]*Sign[r11], 
       6.349999999999998}, {r11 - 1.*Sqrt[-71.68444444444441 + r11^2]*
         Sign[r11], 8.466666666666665}, 
      {r11 - 1.*Sqrt[-112.00694444444434 + r11^2]*Sign[r11], 
       10.583333333333329}, {r11 - 1.*Sqrt[-161.28999999999988 + r11^2]*
         Sign[r11], 12.699999999999996}, 
      {2.1 + r12 - 1.*Sqrt[-161.29 + r12^2]*Sign[r12], 12.7}, 
      {2.1 + r12 - 1.*Sqrt[-112.00694444444441 + r12^2]*Sign[r12], 
       10.583333333333332}, 
      {2.1 + r12 - 1.*Sqrt[-71.68444444444445 + r12^2]*Sign[r12], 
       8.466666666666667}, 
      {2.1 + r12 - 1.*Sqrt[-40.322500000000005 + r12^2]*Sign[r12], 
       6.3500000000000005}, 
      {2.1 + r12 - 1.*Sqrt[-17.92111111111112 + r12^2]*Sign[r12], 
       4.233333333333334}, 
      {2.1 + r12 - 1.*Sqrt[-4.48027777777778 + r12^2]*Sign[r12], 
       2.116666666666667}, 
      {2.1 + r12 - 1.*Sqrt[-3.1554436208840472*^-30 + r12^2]*Sign[r12], 
       1.7763568394002505*^-15}, 
      {2.1 + r12 - 1.*Sqrt[-4.480277777777765 + r12^2]*Sign[r12], 
       -2.1166666666666636}, 
      {2.1 + r12 - 1.*Sqrt[-17.921111111111088 + r12^2]*Sign[r12], 
       -4.233333333333331}, 
      {2.1 + r12 - 1.*Sqrt[-40.32249999999997 + r12^2]*Sign[r12], 
       -6.349999999999998}, 
      {2.1 + r12 - 1.*Sqrt[-71.68444444444441 + r12^2]*Sign[r12], 
       -8.466666666666665}, 
      {2.1 + r12 - 1.*Sqrt[-112.00694444444434 + r12^2]*Sign[r12], 
       -10.583333333333329}, 
      {2.1 + r12 - 1.*Sqrt[-161.28999999999988 + r12^2]*Sign[r12], 
       -12.699999999999996}, {r11 - 1.*Sqrt[-161.29 + r11^2]*Sign[r11], 
       -12.7}, {2.1 + r12 - 1.*Sqrt[-161.29 + r12^2]*Sign[r12], -12.7}, 
      {2.1 + r12 - 1.*Sqrt[-112.00694444444441 + r12^2]*Sign[r12], 
       -10.583333333333332}, 
      {2.1 + r12 - 1.*Sqrt[-71.68444444444445 + r12^2]*Sign[r12], 
       -8.466666666666667}, 
      {2.1 + r12 - 1.*Sqrt[-40.322500000000005 + r12^2]*Sign[r12], 
       -6.3500000000000005}, 
      {2.1 + r12 - 1.*Sqrt[-17.92111111111112 + r12^2]*Sign[r12], 
       -4.233333333333334}, 
      {2.1 + r12 - 1.*Sqrt[-4.48027777777778 + r12^2]*Sign[r12], 
       -2.116666666666667}, 
      {2.1 + r12 - 1.*Sqrt[-3.1554436208840472*^-30 + r12^2]*Sign[r12], 
       -1.7763568394002505*^-15}, 
      {2.1 + r12 - 1.*Sqrt[-4.480277777777765 + r12^2]*Sign[r12], 
       2.1166666666666636}, 
      {2.1 + r12 - 1.*Sqrt[-17.921111111111088 + r12^2]*Sign[r12], 
       4.233333333333331}, 
      {2.1 + r12 - 1.*Sqrt[-40.32249999999997 + r12^2]*Sign[r12], 
       6.349999999999998}, 
      {2.1 + r12 - 1.*Sqrt[-71.68444444444441 + r12^2]*Sign[r12], 
       8.466666666666665}, 
      {2.1 + r12 - 1.*Sqrt[-112.00694444444434 + r12^2]*Sign[r12], 
       10.583333333333329}, 
      {2.1 + r12 - 1.*Sqrt[-161.28999999999988 + r12^2]*Sign[r12], 
       12.699999999999996}, {9.6 + r13 - 1.*Sqrt[-161.29 + r13^2]*
         Sign[r13], 12.7}, 
      {9.6 + r13 - 1.*Sqrt[-112.00694444444441 + r13^2]*Sign[r13], 
       10.583333333333332}, 
      {9.6 + r13 - 1.*Sqrt[-71.68444444444445 + r13^2]*Sign[r13], 
       8.466666666666667}, 
      {9.6 + r13 - 1.*Sqrt[-40.322500000000005 + r13^2]*Sign[r13], 
       6.3500000000000005}, 
      {9.6 + r13 - 1.*Sqrt[-17.92111111111112 + r13^2]*Sign[r13], 
       4.233333333333334}, 
      {9.6 + r13 - 1.*Sqrt[-4.48027777777778 + r13^2]*Sign[r13], 
       2.116666666666667}, 
      {9.6 + r13 - 1.*Sqrt[-3.1554436208840472*^-30 + r13^2]*Sign[r13], 
       1.7763568394002505*^-15}, 
      {9.6 + r13 - 1.*Sqrt[-4.480277777777765 + r13^2]*Sign[r13], 
       -2.1166666666666636}, 
      {9.6 + r13 - 1.*Sqrt[-17.921111111111088 + r13^2]*Sign[r13], 
       -4.233333333333331}, 
      {9.6 + r13 - 1.*Sqrt[-40.32249999999997 + r13^2]*Sign[r13], 
       -6.349999999999998}, 
      {9.6 + r13 - 1.*Sqrt[-71.68444444444441 + r13^2]*Sign[r13], 
       -8.466666666666665}, 
      {9.6 + r13 - 1.*Sqrt[-112.00694444444434 + r13^2]*Sign[r13], 
       -10.583333333333329}, 
      {9.6 + r13 - 1.*Sqrt[-161.28999999999988 + r13^2]*Sign[r13], 
       -12.699999999999996}, {2.1 + r12 - 1.*Sqrt[-161.29 + r12^2]*
         Sign[r12], -12.7}}]; polygons12188572008312134018987715 = 
    Compile[{r21, r22, r23}, {{r21 - 1.*Sqrt[-161.29 + r21^2]*Sign[r21], 
       -12.7}, {r21 - 1.*Sqrt[-112.00694444444441 + r21^2]*Sign[r21], 
       -10.583333333333332}, {r21 - 1.*Sqrt[-71.68444444444445 + r21^2]*
         Sign[r21], -8.466666666666667}, 
      {r21 - 1.*Sqrt[-40.322500000000005 + r21^2]*Sign[r21], 
       -6.3500000000000005}, {r21 - 1.*Sqrt[-17.92111111111112 + r21^2]*
         Sign[r21], -4.233333333333334}, 
      {r21 - 1.*Sqrt[-4.48027777777778 + r21^2]*Sign[r21], 
       -2.116666666666667}, 
      {r21 - 1.*Sqrt[-3.1554436208840472*^-30 + r21^2]*Sign[r21], 
       -1.7763568394002505*^-15}, 
      {r21 - 1.*Sqrt[-4.480277777777765 + r21^2]*Sign[r21], 
       2.1166666666666636}, {r21 - 1.*Sqrt[-17.921111111111088 + r21^2]*
         Sign[r21], 4.233333333333331}, 
      {r21 - 1.*Sqrt[-40.32249999999997 + r21^2]*Sign[r21], 
       6.349999999999998}, {r21 - 1.*Sqrt[-71.68444444444441 + r21^2]*
         Sign[r21], 8.466666666666665}, 
      {r21 - 1.*Sqrt[-112.00694444444434 + r21^2]*Sign[r21], 
       10.583333333333329}, {r21 - 1.*Sqrt[-161.28999999999988 + r21^2]*
         Sign[r21], 12.699999999999996}, 
      {7.5 + r22 - 1.*Sqrt[-161.29 + r22^2]*Sign[r22], 12.7}, 
      {7.5 + r22 - 1.*Sqrt[-112.00694444444441 + r22^2]*Sign[r22], 
       10.583333333333332}, 
      {7.5 + r22 - 1.*Sqrt[-71.68444444444445 + r22^2]*Sign[r22], 
       8.466666666666667}, 
      {7.5 + r22 - 1.*Sqrt[-40.322500000000005 + r22^2]*Sign[r22], 
       6.3500000000000005}, 
      {7.5 + r22 - 1.*Sqrt[-17.92111111111112 + r22^2]*Sign[r22], 
       4.233333333333334}, 
      {7.5 + r22 - 1.*Sqrt[-4.48027777777778 + r22^2]*Sign[r22], 
       2.116666666666667}, 
      {7.5 + r22 - 1.*Sqrt[-3.1554436208840472*^-30 + r22^2]*Sign[r22], 
       1.7763568394002505*^-15}, 
      {7.5 + r22 - 1.*Sqrt[-4.480277777777765 + r22^2]*Sign[r22], 
       -2.1166666666666636}, 
      {7.5 + r22 - 1.*Sqrt[-17.921111111111088 + r22^2]*Sign[r22], 
       -4.233333333333331}, 
      {7.5 + r22 - 1.*Sqrt[-40.32249999999997 + r22^2]*Sign[r22], 
       -6.349999999999998}, 
      {7.5 + r22 - 1.*Sqrt[-71.68444444444441 + r22^2]*Sign[r22], 
       -8.466666666666665}, 
      {7.5 + r22 - 1.*Sqrt[-112.00694444444434 + r22^2]*Sign[r22], 
       -10.583333333333329}, 
      {7.5 + r22 - 1.*Sqrt[-161.28999999999988 + r22^2]*Sign[r22], 
       -12.699999999999996}, {r21 - 1.*Sqrt[-161.29 + r21^2]*Sign[r21], 
       -12.7}, {7.5 + r22 - 1.*Sqrt[-161.29 + r22^2]*Sign[r22], -12.7}, 
      {7.5 + r22 - 1.*Sqrt[-112.00694444444441 + r22^2]*Sign[r22], 
       -10.583333333333332}, 
      {7.5 + r22 - 1.*Sqrt[-71.68444444444445 + r22^2]*Sign[r22], 
       -8.466666666666667}, 
      {7.5 + r22 - 1.*Sqrt[-40.322500000000005 + r22^2]*Sign[r22], 
       -6.3500000000000005}, 
      {7.5 + r22 - 1.*Sqrt[-17.92111111111112 + r22^2]*Sign[r22], 
       -4.233333333333334}, 
      {7.5 + r22 - 1.*Sqrt[-4.48027777777778 + r22^2]*Sign[r22], 
       -2.116666666666667}, 
      {7.5 + r22 - 1.*Sqrt[-3.1554436208840472*^-30 + r22^2]*Sign[r22], 
       -1.7763568394002505*^-15}, 
      {7.5 + r22 - 1.*Sqrt[-4.480277777777765 + r22^2]*Sign[r22], 
       2.1166666666666636}, 
      {7.5 + r22 - 1.*Sqrt[-17.921111111111088 + r22^2]*Sign[r22], 
       4.233333333333331}, 
      {7.5 + r22 - 1.*Sqrt[-40.32249999999997 + r22^2]*Sign[r22], 
       6.349999999999998}, 
      {7.5 + r22 - 1.*Sqrt[-71.68444444444441 + r22^2]*Sign[r22], 
       8.466666666666665}, 
      {7.5 + r22 - 1.*Sqrt[-112.00694444444434 + r22^2]*Sign[r22], 
       10.583333333333329}, 
      {7.5 + r22 - 1.*Sqrt[-161.28999999999988 + r22^2]*Sign[r22], 
       12.699999999999996}, {9.6 + r23 - 1.*Sqrt[-161.29 + r23^2]*
         Sign[r23], 12.7}, 
      {9.6 + r23 - 1.*Sqrt[-112.00694444444441 + r23^2]*Sign[r23], 
       10.583333333333332}, 
      {9.6 + r23 - 1.*Sqrt[-71.68444444444445 + r23^2]*Sign[r23], 
       8.466666666666667}, 
      {9.6 + r23 - 1.*Sqrt[-40.322500000000005 + r23^2]*Sign[r23], 
       6.3500000000000005}, 
      {9.6 + r23 - 1.*Sqrt[-17.92111111111112 + r23^2]*Sign[r23], 
       4.233333333333334}, 
      {9.6 + r23 - 1.*Sqrt[-4.48027777777778 + r23^2]*Sign[r23], 
       2.116666666666667}, 
      {9.6 + r23 - 1.*Sqrt[-3.1554436208840472*^-30 + r23^2]*Sign[r23], 
       1.7763568394002505*^-15}, 
      {9.6 + r23 - 1.*Sqrt[-4.480277777777765 + r23^2]*Sign[r23], 
       -2.1166666666666636}, 
      {9.6 + r23 - 1.*Sqrt[-17.921111111111088 + r23^2]*Sign[r23], 
       -4.233333333333331}, 
      {9.6 + r23 - 1.*Sqrt[-40.32249999999997 + r23^2]*Sign[r23], 
       -6.349999999999998}, 
      {9.6 + r23 - 1.*Sqrt[-71.68444444444441 + r23^2]*Sign[r23], 
       -8.466666666666665}, 
      {9.6 + r23 - 1.*Sqrt[-112.00694444444434 + r23^2]*Sign[r23], 
       -10.583333333333329}, 
      {9.6 + r23 - 1.*Sqrt[-161.28999999999988 + r23^2]*Sign[r23], 
       -12.699999999999996}, {7.5 + r22 - 1.*Sqrt[-161.29 + r22^2]*
         Sign[r22], -12.7}}]), Deployed -> False, 
AutorunSequencing -> {1, 2, 3, 7, 8, 9, 10},
 SynchronousInitialization -> False, SynchronousUpdating -> True]

CAPTION

A lens doublet designates a lens having three spherical surfaces. The radii of curvatures are specified by r11, r12, and r13 for the first lens, and r21, r22, and r23 for the second lens. In addition, a lens doublet is specified by two thicknesses: t1, t2 set to 2.1 and 7.5 mm. The circular aperture is set to 25.5 mm and the material is composed of SF5 and BK7 glass.

To change the  lens positions along the line, drag the red locators. The blue locators rotate the objects. Use the sliders to change the radii of curvatures of the lenses.

DETAILS

This sample optics Demonstration was created using Optica 3, an application package available for Mathematica. It was created by selecting component elements and setting symbolic variables from the GUI popup menus. The variables that users select as symbolic become controllable with sliders and locators.
See http://www.opticasoftware.com for details about additional ray-tracing products.

THIS NOTEBOOK IS THE SOURCE CODE FROM

"Lens Doublets" from the Wolfram Demonstrations Project
 http://demonstrations.wolfram.com/LensDoublets/

Contributed by: Ann Williamson and Donald Barnhart

(Optica Software, a division of Barnhart Optical Research)

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 »