This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Numerics Volume, Chapter 1

Consider the following 2D scattering system: A light ray is (potentially repeatedly) reflected in the interior region formed by three circles centered at the vertices of a regular triangle [✶], [✶], [✶]. The ray can escape through the spaces between the circles. For (at least) 200 rays starting at equidistant points along the line connecting the (centers of) two circles, in (at least) 200 discrete directions, visualize the dependence of the three possible exit points on the starting points and the starting angles. Use circles of 9/10 and 99/100 of their maximal size.

Off[General::spell]; Off[General::spell1];

threeCircleConfiguration[ρ_] := 
With[{𝓈 = N[Sqrt[3], Precision[ρ]]},
{Circle[{𝓈, -1}/2, ρ], Circle[{-𝓈, -1}/2, ρ], Circle[{0, 1}, ρ], 
 Line[{{𝓈 - ρ, -1 + 𝓈 ρ}, {ρ, 2 - 𝓈 ρ}}/2], 
 Line[{{-𝓈 + ρ, -1 + 𝓈 ρ}, {-ρ, 2 - 𝓈 ρ}}/2], 
 Line[{{𝓈 - 2ρ, -1}, {-𝓈 + 2ρ, -1}}/2]}]

Show[Graphics[threeCircleConfiguration[0.5]],
     PlotRange -> All, AspectRatio -> Automatic, Frame -> True];

Given a ray Ray[p, d] starting at the point p and moves in (normalized) direction d, the function nearestIntersection calculates the nearest intersection of this ray with a circle, or with a line segment joining the circles.

(* nearest intersection (if any) of a ray with a circle *)
nearestIntersection[Ray[p_, d_], Circle[q_, ρ_]] := 
Module[{eqs = (p - q + t d).(p - q + t d) - ρ^2, sol},
  (* next intersection if "forward" *)
  sol = Select[t /. Solve[eqs == 0, t], (Im[#] == 0 && # > 0)&];
  If[sol === {}, {}, {t, p + t d} /. t -> Min[sol]]]

(* nearest intersection (if any) of a ray with a line segment *)
nearestIntersection[Ray[p_, d_], Line[{p1_, p2_}]] := 
Module[{eqs = Thread[p1 + s (p2 - p1) == p + t d], solPre, sol},
  (* the ray is potentially parallel to the line segment *)
  solPre = Solve[eqs, {s, t}]; 
  If[solPre === {}, {}, (* crossing happens inside segment? *)
     sol = Select[{s, t} /. solPre, (0 <= #[[1]] <= 1 && 0 < #[[2]])&];
     If[sol === {}, {}, {t, p + t d} /. t -> Min[Last /@ sol]]]]

To carry out the reflection at a point on the circle, we use the function reflect.

(* reflection of a ray at the point s at a circle *)
reflect[Ray[_, d_], s_, Circle[q_, _]] := 
Module[{n = #/Sqrt[#.#]&[s - q]}, Ray[s, d - 2d.n n]]

To quickly decide if the intersection of a ray with a circle happens with the "inside" part of the circle, we define the angle ranges for the three circles.

Evaluate[angleRange /@ {1, 2, 3}] = 
         Interval /@ (Pi/3 {{2, 3}, {0, 1}, {-2, -1}});

Now, we have all functions defined that are needed to carry out one reflection (or the final escape step) for a ray. The function step does this. In general, we will have to use a few hundred digits to carry out many reflections and obtain a correct result. But to find the next individual scatterer, we can use much lower precision to speed up the calculations.

step[ρ_, {rayHP:Ray[p_, d_], lastCircle___}] :=
Module[{tCC, ray, possibleCircles, iCs, iCs1, α, result, 
        nextEvent, possibleLines, isL, isL1, tLMins, posL},
 tCC = threeCircleConfiguration[ρ];
 (* use reduced precision to find next event *)
 ray = N[rayHP, 25];
 (* next reflection is not on last circle *)
 possibleCircles = Delete[Transpose[{{1, 2, 3}, Take[tCC, 3]}],
                          {lastCircle}];
 (* potential intersection with remaining circles *)
 iCs = {#1, nearestIntersection[ray, #2]}& @@@ possibleCircles;    
 (* reflection happens on the inside circle segment *)
 iCs1 = Select[iCs, (#[[2]] =!= {} && 
            IntervalMemberQ[angleRange[#[[1]]],
                     ArcTan @@ (#[[2, 2]] - tCC[[#[[1]]]][[1]])])&];
 (* is next event a reflection? *)
 If[iCs1 =!= {}, (* next event is a reflection *)
    nextEvent = "Reflection";
    α = iCs1[[Position[#, Min[#]]&[#[[2, 1]]& /@ iCs1][[1, 1]], 1]];   
    {circleNumber, circle} = {α, tCC[[α]]}];
(* escape *)
If[nextEvent =!= "Reflection", nextEvent = "Escape";  
   (* happens only once, so test all three cases *)
   possibleLines = Transpose[{{1, 2, 3}, Take[tCC, -3]}];
   isL = {#1, nearestIntersection[ray, #2]}& @@@ possibleLines;
   isL1 = DeleteCases[isL, {_, {}}];
   tLMins = #[[2, 1]]& /@ isL1;                                 
   posL = Position[tLMins, Min[tLMins]][[1, 1]];
   result = {isL1[[posL, 1]], isL1[[posL, 2, 2]]}];
(* return result *)
If[nextEvent === "Reflection",
   (* high-precision reflection *)
    𝓈 = nearestIntersection[rayHP, circle][[2]];
   {reflect[rayHP, 𝓈, tCC[[circleNumber]]], circleNumber}, 
   result]]

The full path of a ray until it escapes (or until the ray data do not have enough precision anymore to follow the ray) can be calculated using the function rayPath.

rayPath[ρ_, ray:Ray[p_, d_]] :=
NestWhileList[step[ρ, #]&, {ray}, 
              (* stop at exit or insufficient precision *)
              ((Head[First[#]] === Ray) && (Precision[#] > 20))&]

To visualize the path of a ray, we implement the function rayPathGraphic. It colors the path from red to blue and draws a gray dot at the escape point. To concentrate on the ray path, we show only the relevant parts of the circles.

rayPathGraphic[ρ_, rayPath_, opts___] :=
With[{(* reflection points *)
      points = If[Head[#[[1]]] === Ray, #[[1, 1]], #[[2]]]& /@ rayPath,
      (* exit point *)
      exitPoint = If[Head[rayPath[[-1, 1]]] =!= Ray, 
                     {rayPath[[-1, 2]]}, {}],
       xMin = -Sqrt[3]/2 + ρ/2, xMax = Sqrt[3]/2 - ρ/2, 
       yMin = -1/2, yMax = 1 - Sqrt[3]/2 ρ, δ = (Sqrt[3] - ρ)/20}, 
Show[Graphics[
{(* the three scattering circles *)
 {Thickness[0.006], GrayLevel[0], 
  MapIndexed[Function[{x, p}, Insert[x, angleRange[p[[1]]][[1]], 3]],
           Take[threeCircleConfiguration[ρ], 3]],
  Take[threeCircleConfiguration[ρ], -3]},
 (* ray path *)
 {Thickness[0.0001], 
  MapIndexed[{Hue[0.78 (#2[[1]] - 1)/(Length[points] - 1)], Line[#1]}&, 
             Partition[points, 2, 1]]},
 (* exit point *)
 {PointSize[0.025], GrayLevel[1/2], Point /@ exitPoint}}],
     PlotRange -> {{xMin - δ, xMax + δ}, {yMin - δ, yMax + δ}}, 
     opts, AspectRatio -> Automatic, Frame -> True, FrameTicks -> False]]

Here are four example paths. The parameter for the rightmost path are chosen to give 10 internal reflections.

Show[GraphicsArray[
rayPathGraphic[1/2, rayPath[1/2, Ray[{0, -1/2}, #] // N[#, 100]&],
               DisplayFunction -> Identity]& /@
        {{0, 1}, {11, 1}, {1, 1}, (* long ray path *)
         {65690846499036600543, 33168393638474293644}}]]

Here is an animation showing the dependence of the paths from the starting angle. All paths start at the center between two circles.

With[{ε = 10^-6}, 
 Do[rayPathGraphic[4/5, rayPath[4/5, Ray[{0, -1/2}, 
                                ExpandAll[{Cos[φ], Sin[φ]}]] // N[#, 250]&]],
    {φ, ε, Pi - ε, (Pi - 2ε)/90}]]

For an easier specification of the starting points of the rays for arbitrary circle radii, we define a scaled radius (Subscript[ρ, s])=2 (3^(-1/2))ρ. So, (Subscript[ρ, s]) varies between 0 and 1. The initial position is measured with respect to the length of the line segment and ranges over the interval (-1, 1).

scaledρ[ρRel_] := Sqrt[3]/2 ρRel
      
scaledRayData[ρRel_, {ξ0_, φ0_}] :=
Module[{ρ = scaledρ[ρRel]},  
      {ρ, Ray[{ξ0 (Sqrt[3]/2 - ρ), -1/2}, {Cos[φ0], Sin[φ0]}]}]

Here are two examples of paths for different circle radii.

Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
Function[{ρRel, ξ0, φ0, prec}, rayPathGraphic[Sqrt[3]/2 ρRel,
         rayPath @@ N[scaledRayData[ρRel, {ξ0, φ0}], prec]]] @@@
         {{9/10, -1/2, Pi/6, 300}, {99/100, 0, Pi/6, 300}}]]];

Now, we analyze the path length on a grid in the ξ0,φ0-plane. The function exitResult returns a two element list of the form {pathLength,exitNumber}. We use circles of relative size 9/10.

exitResult[ρRel_, {ξ0_, φ0_}, prec_] :=
Module[{rP}, 
  rP = rayPath @@ N[scaledRayData[ρRel, {ξ0, φ0}], prec]; 
  {Length[rP], (* return 0 if no exit was reached *)
   If[MemberQ[Last[rP], _Ray, Infinity], 0, rP[[-1, 1]]]}]

To generate the exit number data for a dense grid in the ξ0,φ0-plane, we make use of the mirror symmetry of the circle configuration along the y-axis. The function makeξφExitData calculates the exit numbers for angles 0<φ≤π/2 explicitly and calculates the exit numbers for angles π/2≤φ<π/2 using symmetry.

makeξφExitData[er_, ρRel_, pp_, prec_] :=
Module[{φM = Pi If[EvenQ[pp], 1/2 , 1/2 - (1 - 2/pp)/pp/2], t},
 (* join data and double mirrored data *)
 Join[#, Map[Reverse, Map[{{#[[1, 1]], Pi - #[[1, 2]]}, 
             {#[[2, 1]], #[[2, 2]] /. {1 -> 2, 2 -> 1}}}&, 
            If[EvenQ[pp], Drop[#, -1], #], {2}], {0, 1}]]&[
 (* half of the data *)
 Table[{{ξ0, φ0}, er[ρRel, {ξ0, φ0}, prec]},
        {φ0, Pi/pp, φM, (φM - Pi/pp)/If[EvenQ[pp], pp/2, (pp - 1)/2]},
        {ξ0, -1 + 1/pp, 1 - 1/pp, (2 - 2/pp)/pp}]]]

Here are the path length and exit numbers for the circles of 90% their possible maximum size.

data[9/10] = makeξφExitData[exitResult, 9/10, 200, 250];

We display the path length as a density plot and the exit number we indicate by red (exit 1), green (exit 2), and blue (exit 3) squares. In case too much precision was lost in the calculation and the exit number could not be determined, we draw a white square.

exitResultGraphics[exitData_] := 
Show[GraphicsArray[
Block[{$DisplayFunction = Identity,
       opts = Sequence[Mesh -> False, PlotRange -> All,
       MeshRange -> {{-1 + 1/#, 1 - 1/#}&[Length[exitData[[1]]] - 1],
                     Pi {1/#, 1 - 1/#}&[Length[exitData] - 1]}]},
      {(* path length graphic *)
       ListDensityPlot[Log @ Map[#[[2, 1]]&, exitData, {2}], opts, 
                      ColorFunction -> GrayLevel],
       (* exit number graphic *)
       ListDensityPlot[Map[#[[2, 2]]&, exitData, {2}], opts, 
                      ColorFunction -> 
                      (If[# == 0, (* no exit found *) RGBColor[1, 1, 1],
                          MapAt[1&, RGBColor[0, 0, 0], #]]&),
                      ColorFunctionScaling -> False]}]]]

For the 9/10 circle, we see bands of equal exit numbers.

exitResultGraphics[data[9/10]]

The next two graphics show one of the longest ray paths and the cumulative distribution of all path lengths.

(* extract longest path parameters *)
extractLongestPath[data_] := Expand[Part[data, ##]& @@@ 
        Position[#, Max[#]]&[Map[#[[2, 1]]&, data, {2}]]]
        
maxData = extractLongestPath[data[9/10]]

(* function to visualize longest path and 
   cumulative path length distribution *)
showMaxDataResults[ρRel_, data_, prec_] := 
Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
{(* the path with the most reflections *)
 Module[{maxData = extractLongestPath[data], ξ0, φ0}, 
         {ξ0, φ0} = maxData[[1, 1]];
     rayPathGraphic[Sqrt[3]/2 ρRel,
                    rayPath @@ N[scaledRayData[ρRel, {ξ0, φ0}], prec]]],
 (* cumulative distribution of path lengths *)
 ListPlot[MapIndexed[{#2[[1]], Log[#1]}&, 
            Sort[Flatten[Map[#[[2, 1]]&, data, {2}]]]],
         PlotRange -> All]}]]];

showMaxDataResults[9/10, data[9/10], 300]

For circles of relative size 99/100, we have potentially very long paths and so need many digits to calculate the exit number. To not slow down the calculation of the exit number of shorter paths by using excessively many digits, we define a function adapativeExitResult that starts with precision prec. In case no exit can be found, adapativeExitResult iteratively doubles the precision until an exit is found.

adapativeExitResult[ρRel_, {ξ0_, φ0_}, precInit_] := 
Module[{er, prec = precInit},
       While[er = exitResult[ρRel, {ξ0, φ0}, prec]; 
             Last[er] === 0, prec = 2 prec]; er]

The next input shows an example where exitResult fails, but adapativeExitResult succeeds.

{exitResult[##], adapativeExitResult[##]}& @@
           {99/100, {121/147, 59/441 Pi}, 200}

Using a 2 GHz computer, we can now in about two hours calculate the exits on a 200×200 grid. (Because all rays are independent, this problem is, of course, ideal for parallelization.)

(data[99/100] = makeξφExitData[adapativeExitResult, 
                               99/100, 200, 250]); // Timing

Now, the bands from above have nearly disappeared and the points of a given color have nearby points of the other two colors [✶].

exitResultGraphics[data[99/100]]

showMaxDataResults[99/100, data[99/100], 800]