Graphics Settings

<< Miscellaneous`Units`;
Remove[Gray];
<< Graphics`Master`;
Off[General::spell1];
fs = 7;
fsbig = 9;
width  =  Convert[8.64 Centi Meter,  Point]/Point;
th  =  0.0008;

Photon Trajectory code, needed

Off[General::spell1];
Clear[
        h];


norm[a_]  := (Sqrt[a . a]);

    
myRealDigits[a_Real]  :=  
        Module[{man, exp, ss}, 
{man, exp} = 
            MantissaExponent[Abs[a]]; 
man = 
            ToString[InputForm[man]]; 
ss = 
            StringPosition[man,  "`"]; 
If[
            ss ≠ {},  
            man  =  
              StringDrop[
                man,  {ss[[1, 1]],  
                  StringLength[man]}]]; 
{ToExpression /@ 
              PadRight[Characters[StringDrop[man, 2]], 
                Ceiling[ Precision[a]]],  exp}];

    
mySetPrecision[0 | 0. , _] := 0;
\(\(mySetPrecision[a_, 
          prec_] := \((
rd = 
            myRealDigits[SetPrecision[a, prec + 1]]; 
fd = 
            Take[rd[[1]], 
              prec]; 
If[
            rd[[1, 
                prec + 1]] ≥ 5, 
            fd = ReplacePart[fd, 
                fd[[\(-1\)]] + 1, 
                prec]]; 
ffd = 
            FromDigits[{fd, 
                rd[[2]]}]; 
\(Sign[a]\ N[Numerator[ffd], prec]\)\/N[Denominator[ffd], prec])\);\)

 (*\ 
      not\ necessary\ in\ version\ 5\ *) \)
Unprotect[Round]; 
Round[Interval[{a_,  b_}]]  :=  
        Interval[{Round[a],  Round[b]}];
\(\(Protect[Round];\)\ 
\)
iMin[{Interval[{a_,  b_}],  
            Interval[{c_,  d_}]}]   :=   Interval[{Min[a, c],  
            Min[b, d]}];

    
IntervalToReal[Interval[{a_, b_}]] := ((a + b)/2);
IntervalToReal[a_] := a;
RealToInterval[r_, d_] := 
        Interval[r + {-1, 1} d];

    
diam[Interval[{a_,  b_}]]  :=  b - a;
diam[{i__Interval}]  :=  
        Max[diam  /@  {i}];
h[{a_, b_}] := 
        9 {{(b^2) - (a^2), -2 a b}, {-2 a b, 
              (a^2) - (b^2)}};

    
\(\(\(PhotonTrajectory::usage\  = \ "PhotonTrajectory[p, v, targettime (opts)] gives the set of reflection points and the first and last point. If interval arithmetic is used, the results are all intervals.";\)\(
\)
    \)\)
StartIntervalPrecision::usage = "StartIntervalPrecision is an \
option to PhotonTrajectory that sets the starting size of the intervals. A \
setting of automatic means that the accuracy goal is used for the initial \
size.";
IntervalPrecisionStep::usage = "IntervalPrecisionStep is an option \
to PhotonTrajectory that sets the number of powers of 10 by which the \
interval size is decreased after a failure.";
IntervalArithmetic::usage = "IntervalArithmetic is a possible \
choice to the Method option of PhotonTrajectory.";
FixedPrecision::usage = "FixedPrecision is a possible choice to the \
Method option of PhotonTrajectory.";
SignificanceArithmetic::usage  =  "SignificanceArithmetic is a \
possible choice to the Method option of PhotonTrajectory. When using this \
method it is essential that the initial conditions be approximate real \
numbers.";
FixedWorkingPrecision::usage  =  "FixedWorkingPrecision is an \
option to PhotonTrajectory that sets the fixed precision to be used when the \
FixedPrecision method us chosen.";

    
\(\(Options[
          PhotonTrajectory]\  := \ {StartIntervalPrecision\  -> 
            Automatic, AccuracyGoal\  -> \ 12, \ 
          IntervalPrecisionStep\  -> \ 1, \ 
          FixedWorkingPrecision -> 16, \ Trace\  -> False, \ 
          Method\  -> \ 
            SignificanceArithmetic};\)
\ \)
\(\(PhotonTrajectory[p_, v_, targettime_, \ opts___Rule] := 
        Module[{ans, mat, pint, t, a, 
            b}, 
{ag, stint, \ ips, met, \ fp, \ 
              trQ}\  = \({AccuracyGoal, \ StartIntervalPrecision, \ 
                  IntervalPrecisionStep, Method, \ \ FixedWorkingPrecision, \ 
                  Trace}\  /. \ {opts}\)\  /. \ 
              Options[PhotonTrajectory]; 

If[stint === Automatic, \ 
            stint\  = \ ag]; 

Which[
            met\  === \ 
              FixedPrecision, \ 

sp[
                a_]\  := \ a /. \ 
                z_Real\  :> \ mySetPrecision[z, fp]; 
p1 = 
              sp[N[p, fp]]; v1 = sp[N[v, fp]]; pdata = {sp[p1]}; 

            resttime = targettime; 

            While[resttime > 0, 
              m = Round[p1 + \(2\ v1\)\/3]; 
s = 
                t /. Solve[
                    norm[p1 + t\ v1 - m]^2 == 1/9, 
                    t]; 
t0 = 
                Min[Cases[s, _?Positive]]; 
If[
                t0 < resttime, 
p1 = sp[p1 + t0\ v1]; 
                v1 = sp[h[p1 - m] . v1], \ \ 
\ \ t0 = 
                  Min[resttime, 2\/3]; 
                p1 = sp[p1 + t0\ v1]]; 
resttime = 
                sp[resttime - t0]; 
\ \ \ AppendTo[pdata, 
                p1]]; 
\ 
            pdata, 


met\  === \ IntervalArithmetic, 
lastpoint = 
              Interval[{\(-∞\), ∞}]; intsize = stint; 
            wp = Max[17, 2 + intsize]; 
While[
              diam[lastpoint] > 10\^\(-ag\), 
pint = 
                N[\((RealToInterval[#, 10\^\(-intsize\)] &)\) /@ p, 
                  wp]; 
vint = 
                N[\((RealToInterval[#, 10\^\(-intsize\)] &)\) /@ v, 
                  wp]; 
pathdata = {N[pint, 
                    wp]}; 
resttime = 
                Interval[{targettime, targettime}]; 
While[
                resttime > 0, 
                m = Round[pint + \(2\ vint\)\/3]; 
S = 
                  t /. \((Solve[\((pint + t\ vint - 
                                m)\) . \((pint + t\ vint - m)\) == 1/9, 
                        t])\); 
If[
                  FreeQ[S, 
                    Power[Interval[{_?Negative, \ _?
                            Positive}], _]], 
\ \ \ \ \ \ \ \ \ \ \ T = Min[Cases[S, _?Positive]], 
                  Break[]]; 
Which[
T ≤ resttime, 
\ \ pint += T\ vint; 
                  vint = h[pint - m] . vint; \ 
                  resttime -= 
                    T, 
T > 
                      resttime\  && \ \ resttime ≥ 2/3, 
                  pint += 2/3\ vint; \ 
                  resttime -= 2/3, 
T > resttime\  && \ 
                    resttime < 2/3, pint += resttime\ vint; 
                  resttime = 0, 
True, \ 
                  Break[]]; 
AppendTo[pathdata, 
                  pint]; 
If[
                  Precision\ [\ {resttime, \ pint, \ vint, \ T}] < \ ag, \ 
                  Break[]]]; 
intsize += ips; \ 
              wp\  = Max[17, intsize + 2]; 
lastpoint = 
                pint + Table[
                    RealToInterval[\(-Max[Abs[resttime]]\), 
                      Max[Abs[resttime]]], {2}]]; 
If[
              trQ, \ Print[
                StringForm[\*""Initial condition interval radius is (10^(-``))."", 
                  intsize]]]; 
pathdata, 

met\  === \ 
              SignificanceArithmetic, 
p1 = p; 
            v1 = v; 
\ \ \ \ pathdata = {p1}; 
            resttime = targettime; 
\ \ \ \ \ While[
              resttime > 0, c = Round[p1 + \(2\/3\) v1]; 
              time = Min[
                  Cases[t /. Solve[norm[p1 + t\ v1 - c]^2 == 1/9, t], _?
                      Positive]]; 
If[
                time < resttime, 
\ \ \ p1 += time\ v1; 
                v1 = h[p1 - c] . v1, 
time = 
                  Min[resttime, 2/3]; 
                p1 += time\ v1]; 
resttime -= 
                time; 
\ \ AppendTo[pathdata, 
                p1]]; 
pathdata]];\)\)

Figure 2.1

variation[
        d_] := 
(centerdata  =  
          norm[PhotonTrajectory[N[{(1/2), (1/10)}], {1, 0}, 10, 
                Method -> 
                  SignificanceArithmetic][[-1]]]; 
        data = Table[
            norm[PhotonTrajectory[
                  N[{(1/2), 
                        (1/10)} + {Cos[θ], 
                          Sin[θ]} (10^(-d))], {1, 0}, 10, 
                  Method -> 
                    SignificanceArithmetic][[-1]]], {θ, 0. , 
              2 π - ((2 π)/10), ((2 π)/10)}]; 
Max[
          Table[Abs[
              centerdata - 
                data[[i]]], {i, 
              Length[data]}]])

vardata = Table[{ii, Log[10, variation[ii]]},  {ii, 2, 17}]

vardata = {{2, 0.38881987698466264`}, {3, 
            0.35387969990650703`}, {4, 0.25879882807752136`}, {5, 
            0.6527475238845092`}, {6, 0.006718562813029293`}, {7, 
            0.22322151216116007`}, {8, 
            0.32439927840600546`}, {9, -0.03531558296993449`}, {10, -0.1428239035814522`}, {11, -0.1974518247016603`}, {12, -1.0993143839546142`}, {13, -2.1837196898179645`}, {14, -3.2222717027838192`}, {15, -4.2183344038535076`}, {16, -5.0048852098643`}, {17, -5.825222291795251`}};

    
ListPlot[vardata, Frame -> ! True, Axes -> Automatic,  
        AxesOrigin -> {0, -6}, 
        Ticks -> {Join[

({#, 
                      10^StringForm["-``", #]} &) /@ 
                Range[2, 17, 3],  ({#, ""} &) /@ 
                Complement[Range[2, 17],  Range[2, 17, 3]]], 
            Prepend[({#1, (10^StringForm)["-``", -#]} &) /@ 
                Range[-1, -6, -1], {0, 1}]},  
        PlotStyle -> PointSize[ .02],  
        PlotRange -> {{0, 17.5}, {-6, 1}},  
        TextStyle -> FontSize -> fs,  
        AxesLabel ->  {DisplayForm[
              StyleBox["  Initial\n\
  variation", FontSize -> fsbig,  
                FontFamily  ->  "Times"]], 
            DisplayForm[
              StyleBox["Final\nvariation",  FontSize -> fsbig,  
                FontFamily  ->  "Times"]]},  
        ImageSize -> width];

Figure 2.2

h[{a_, b_}] := 
        9 {{(b^2) - (a^2), -2 a b}, {-2 a b, 
              (a^2) - (b^2)}};
p  =  {0.35, -0.3};
gra = 
        GrayLevel[ .2];
v =  {Cos[1.2],  
        Sin[1.2]}; 
    top = p + 
        3.3  v;
gr  =  { 
          Table[Line[ {{-1/2, i},  {5/2,  i}}],  {i,  -1/2,  
              5/2}], 
Table[
            Line[{{i,  -1/2},  {i,  5/2}}],  {i,  -1/2,  
              5/2}], 
 
          Table[ Circle[{i, j},  1/3],  {i,  0,  2},  {j,  0,  
              2}],  
{ 
            Table[Point[{i, j}],  {i,  0,  2},  {j,  0,  
                2}]}, 
ds =  .01; 
          ds1 =  .02; 
Line[{p,  p + 3.2 v}]; 
          Thickness[ .008], 
{Black,  
            Line[{p,  
                p  +  (t /. Solve[(p  +  
                                    t v )[[1]] ==  .5, 
                            t][[1]]) 
                    v}], 
gra, {gra, Dashing[{ds, ds1}],  
              Line[{p  +  (t /. Solve[(p  +  
                                      t v )[[1]] ==  .5, 
                              t][[1]]) v,  
                  p  +  (t /. Solve[(p  +  
                                      t v )[[2]] ==  .5, 
                              t][[1]]) 
                      v}]}, 
Black,  
            Line[{p  +  (t /. Solve[(p  +  
                                    t v )[[2]] ==  .5, 
                            t][[1]]) v,  
                p  +  (t /. Solve[(p  +  
                                    t v )[[2]] == 1.5, 
                            t][[1]]) v}], 
{gra, 
              Dashing[{ds, ds1}],  
              Line[{p  +  (t /. Solve[(p  +  
                                      t v )[[2]] == 1.5, 
                              t][[1]]) v,  
                  p  +  (t /. Solve[(p  +  
                                      t v )[[2]] == 2.5, 
                              t][[1]]) 
                      v}]}, 

Black,  
            Line[{top,  
                p  +  (t /. Solve[(p  +  
                                    t v )[[2]] == 2.5, 
                            t][[1]]) v}]}, 
{Arrow[
              p + 3.1 v,  top,  
              HeadWidth ->   .8]}, 
{PointSize[
              0.04], {GrayLevel[ .9], 
              Disk[p,   .04]}, 
Thickness[ .00001],  
            Circle[p,   .04]}, 
{Hue[0]; 
            PointSize[ .01],  
            Point[qqq = 
                p  +  (t /. Solve[(Sqrt[(p  +  t v - {1, 1} ) . (p  +  t v - {1, 1})]) == 1/3, t][[1]]) 
                    v], 
{Thickness[ .0004], 
              Line[{qqq, 
                  qqq + 1.35 
                       h[qqq - {1, 
                              1}] . (v)}]}}};
im1=Show[Graphics[{gr, Text["P", p - { .11,  .095},TextStyle -> FontSize -> 
                            fsbig - 1]}], AspectRatio -> Automatic, Axes -> True, AxesOrigin -> {-1/2,  -1/2},TextStyle -> {FontFamily  ->  "Times", 
                  FontSize -> fs},PlotRange ->  {{- .6, 2.6}, {- .6, 2.8}},
Ticks  ->  {Range[0,  
                    2.5,   .5],  Range[0,  2.5,   .5]}, DisplayFunction -> Identity];

pp = {(-1/3 (Sqrt[2])/2), (1/
                3 (Sqrt[2])/2)};
\(im2 = 
          Show[Graphics[{Circle[{0, 0}, \ 1/3], PointSize[ .03], \ 
                Point[{0, 0}], 
Point[\@2/2 {\(-1\), 1}/3], 
                Line[{{0, 0}, \ 
                    pp, \ {\(- .5\), \  .5}}], 
Text["1/3", \ {\ \(- .048\), \  .1}, \ {0, \ \(-1\)}], 
Text[\*""((Sqrt[2])/2)-(1/3)"", \ {\(- .28\); \(- .26\), \  .42}, \ {0, \ 
                    0}], 
Text[\*""(Sqrt[2])-(2/3)"", \ { .29, \ \(- .31\)}, \ {0, \ 
                    0}], 
Line[{{\(- .5\), \(- .5\)}, \ { .5, \(- .5\)}, \ { .5,  .5}, \ {\(- .5\),  .5}, {\(- .5\), \(- .5\)}}], 
Line[{{\(- .5\), \ 
                      1\/6\ \((\(-3\) + 2\ \@2)\)}, \ {1\/6\ \((3 - 
                            2\ \@2)\), \  .5}}], 
Line[{{ .5, \ 1\/6\ \((3 - 2\ \@2)\)}, \ {1\/6\ \((\(-3\) + 
                            2\ \@2)\), \(-\  .5\)}}], 
PointSize[ .02]}], \ AspectRatio -> Automatic, \ 
            DisplayFunction -> Identity, 
            TextStyle -> {FontFamily\  -> \ "Times", 
                FontSize -> fsbig - 2}];\)

Show[GraphicsArray[{im1, im2}], 
        ImageSize -> 1.4 width];

Figure 2.3

path = PhotonTrajectory[N[{(1/2), (1/10)}], {1, 0}, 20];
path1 = 
        PhotonTrajectory[N[{(1/2), (1/10)}, 40], {1, 0}, 20];
point = 
        path1[[-1]];

xmax = 1.9; xmin  =  -1. ; ymax = 1.8;
Show[
        Graphics[{PointSize[ .02], 
{Thickness[ .0015],  
              Dashing[{ .01,  .01}], Point  /@  (pts = path1); 
              Line[pts], 
Point[
                pts[[-1]]]}, 
 
            Point  /@  (pts =  path); 
            Line[pts], 
Point[
              point], 
Point[pts[[1]]], 
            Point[pts[[-1]]], 
Table[{{GrayLevel[
                    0.85], Disk[{i, j}, 1/3]}, 
                Circle[{i, j},  1/3]},  {i, -1, 2},  {j, -2, 
                2}], 
Line[{{xmin,  -2},  {xmax,  -2},  {xmax,  ymax},  {xmin,  ymax},  {xmin,  -2}}]}],  
        AspectRatio -> Automatic,  
        Frame -> ! True,  
PlotRange -> {{xmin, 
              xmax +  .01}, {-2.01, ymax}}, FrameTicks -> None,  
        ImageSize -> 1.05 width];

Figure 2.4

answerInterval[d_,  opts___]  :=  
        norm[PhotonTrajectory[{(1/2), (1/10)}, {1, 0}, 10, 
              AccuracyGoal -> d, StartIntervalPrecision -> d + 28, 
              opts,  Method -> 
                IntervalArithmetic][[-1]]];

    
answerFixedPrecision[d_,  opts___]  :=  
        norm[PhotonTrajectory[{(1/2), (1/10)}, {1, 0}, 10, 
              Method -> FixedPrecision,  
              FixedWorkingPrecision -> 
                d][[-1]]];

    
answerSigArith[d_]  :=  
        norm[PhotonTrajectory[N[{(1/2), (1/10)}, d + 40], {1, 0}, 
              10][[-1]]];

timedata = 
      Table[{n, 
Developer`ClearCache[]; 
 Timing[{n, answerInterval[n]}][[1, 
            1]]}, {n, 500, 30000, 500}]

Cached times

timedata = {{500, 0.10000000000000142`}, {1000, 
            0.13000000000000256`}, {1500, 0.18999999999999773`}, {2000, 
            0.240000000000002`}, {2500, 0.3200000000000003`}, {3000, 
            0.3999999999999986`}, {3500, 0.46999999999999886`}, {4000, 
            0.5500000000000043`}, {4500, 0.6699999999999946`}, {5000, 
            0.7600000000000051`}, {5500, 0.8499999999999943`}, {6000, 
            0.980000000000004`}, {6500, 1.0599999999999952`}, {7000, 
            1.1900000000000048`}, {7500, 1.3299999999999983`}, {8000, 
            1.4500000000000028`}, {8500, 1.5799999999999983`}, {9000, 
            1.6999999999999957`}, {9500, 1.8400000000000034`}, {10000, 
            1.9699999999999989`}, {10500, 2.1099999999999994`}, {11000, 
            2.25`}, {11500, 2.4199999999999946`}, {12000, 
            2.5800000000000054`}, {12500, 2.759999999999991`}, {13000, 
            2.8900000000000006`}, {13500, 3.1400000000000006`}, {14000, 
            3.3400000000000034`}, {14500, 3.510000000000005`}, {15000, 
            3.759999999999991`}, {15500, 3.940000000000012`}, {16000, 
            4.099999999999994`}, {16500, 4.269999999999996`}, {17000, 
            4.469999999999999`}, {17500, 4.640000000000001`}, {18000, 
            4.810000000000002`}, {18500, 5.079999999999998`}, {19000, 
            5.240000000000009`}, {19500, 5.279999999999987`}, {20000, 
            5.59999999999998`}, {20500, 5.77000000000001`}, {21000, 
            6.050000000000011`}, {21500, 6.25`}, {22000, 
            6.430000000000007`}, {22500, 6.6299999999999955`}, {23000, 
            6.810000000000002`}, {23500, 7.090000000000003`}, {24000, 
            7.259999999999991`}, {24500, 7.490000000000009`}, {25000, 
            7.699999999999989`}, {25500, 7.97999999999999`}, {26000, 
            8.160000000000025`}, {26500, 8.349999999999994`}, {27000, 
            8.48999999999998`}, {27500, 8.840000000000003`}, {28000, 
            9.120000000000005`}, {28500, 9.329999999999984`}, {29000, 
            9.45999999999998`}, {29500, 10.25`}, {30000, 
            10.550000000000011`}};

Plot

pl = 
        ListPlot[
          timedata,   
Ticks  ->  {Range[10000,  
                30000,  10000],  
              Range[2, 14, 2]},   
AxesLabel -> {None,  
              DisplayForm[
                StyleBox["CPU time [sec]", 
                  FontFamily  ->  "Times",  
                  FontSize -> fsbig]]}, PlotJoined -> True,  
          PlotStyle -> Thickness[th], 
          Epilog  -> {PointSize[ .01], Point  /@  timedata,  
              Text["Digits",  {30000,  .2},  {1, -1},  
                TextStyle -> {FontFamily  ->  "Times", 
                    FontSize -> fsbig}]},  
          ImageSize ->  
            width,   TextStyle -> {FontFamily  ->  "Times", 
              FontSize ->  fs}];

Figure 2.5 (required for Fig 2.6 too)

iMin[{Interval[{a_,  b_}],  
            Interval[{c_,  d_}]}]   :=   Interval[{Min[a, c],  
            Min[b, d]}];
iMin[{}]   :=  ∞;

    
RealToInterval[r_, d_] := 
        Interval[r + {-1, 1} d];

    
diam[Interval[{a_,  b_}]]  :=  b - a;
diam[{i__Interval}]  :=  
        Max[diam  /@  {i}];
h[{a_, b_}] := 
        9 {{(b^2) - (a^2), -2 a b}, {-2 a b, 
              (a^2) - (b^2)}};

    
ReliableTrajectory::usage  =  "ReliableTrajectory[p, v, \
targettime, (opts)] gives a set of {time, position} intervals that give the \
reflection points and the times they are there. The last position is an \
interval of diameter less than 10^-ag, where ag is the AccuracyGoal \
setting.";

    
StartIntervalPrecision::usage = "StartIntervalPrecision is an \
option to ReliableTrajectory that sets the starting size of the intervals. A \
setting of automatic means that the accuracy goal is used for the initial \
size.";

    
IntervalPrecisionStep::usage = "IntervalPrecisionStep is an \
option to ReliableTrajectory that sets the number of powers of 10 by which \
the interval size is decreased after a failure.";

    
ShowSize::usage  =  "ShowSize is an option to \
ReliableTrajectory that asks that the initial condition interval size that \
works be printed.";

    
Options[
          ReliableTrajectory]  :=  {StartIntervalPrecision  -> 
            Automatic, AccuracyGoal  ->  12,  
          IntervalPrecisionStep  ->  1, 
          ShowSize -> False};
 
\(\(ReliableTrajectory[p_, v_, targettime_, \ opts___Rule] := 
        Module[{lastpoint = Interval[{\(-∞\), ∞}], 
            stint, \ ips, \ trQ, pint, intsize, t, resttime, pathdata, S, 
            T}, 
{ag, stint, \ ips, \ 
              trQ}\  = \({AccuracyGoal, \ StartIntervalPrecision, \ 
                  IntervalPrecisionStep, \ ShowSize} /. {opts}\) /. 
              Options[ReliableTrajectory]; 
If[
            stint === Automatic, \ 
            stint\  = \ ag]; 
intsize = stint; 
          wp = Max[17, 2 + intsize]; 
While[
            diam[lastpoint] > 10\^\(-ag\)\ , 
pint = 
              N[\((RealToInterval[#, 10\^\(-intsize\)] &)\) /@ p, 
                wp]; 
vint = 
              N[\((RealToInterval[#, 10\^\(-intsize\)] &)\) /@ v, 
                wp]; 
pathdata = {N[pint, 
                  wp]}; 
resttime = 
              Interval[{targettime, 
                  targettime}]; 

While[
              resttime > 0, 
              m = Round[pint + 2/3\ vint]; 
S = 
                t /. \((Solve[\((pint + t\ vint - 
                              m)\) . \((pint + t\ vint - m)\) == 1/9, 
                      t])\); 
If[
                FreeQ[S, 
                  Power[Interval[{_?Negative, \ _?
                          Positive}], _]], 
\ \ \ \ \ \ \ \ \ \ \ T = iMin[Cases[S, _?Positive]], 
                Break[]]; 
Which[
T ≤ resttime, 
\ \ pint += T\ vint; 
                vint = h[pint - m] . vint; \ 
                resttime -= 
                  T, 
T > 
                    resttime\  && \ \ resttime ≥ 2/3, 
                pint += 2/3\ vint; \ 
                resttime -= 2/3, 
T > resttime\  && \ 
                  resttime < 2/3, pint += resttime\ vint; 
                resttime = 0, 
True, \ 
                Break[]]; 
AppendTo[pathdata, 
                pint]; 
If[
                Precision\ [\ {resttime, \ pint, \ vint, \ T}] < \ ag, \ 
                Break[]]]; 

intsize += 
              ips; \ wp\  = 
              Max[17, intsize + 2]; 
lastpoint = 
              pint + Table[
                  RealToInterval[\(-Max[Abs[resttime]]\), 
                    Max[Abs[resttime]]], {2}]]; 
If[trQ, \ 
            Print[StringForm[\*""Initial condition interval radius is (10^(-``))."", 
                intsize]]]; 
pathdata];\)\)

tt = 2000;
(rr = 
          ReliableTrajectory[{1/2, 1/10}, {1, 0}, tt,  
            AccuracyGoal -> 12, StartIntervalPrecision -> 5459, 
            IntervalPrecisionStep -> 1, 
            ShowSize -> True]);

path  = 
        N[ rr  /.  x_Interval :>  IntervalToReal[x], 
          16];

path // Length

Show[
        Graphics[{Line[path], 
PointSize[ .015],  
            Point[path[[-1]]]}],  Frame -> True,  
        FrameTicks -> {Range[0, 60, 20], Range[0, 40, 20],  None,  
            None}, AspectRatio -> Automatic,  
        TextStyle -> {FontFamily  ->  "Times", 
            FontSize -> fs}];

Figure 2.6

good[p_]  := 
        Abs[ Norm[p  -  Round[p]] - 1/3]  <  .0001;
pathfixed = Select[path, good];
dists = 
        Table[Norm[ (pathfixed[[i + 1]] - 
                pathfixed[[i]])],  {i, 
            Length[pathfixed] - 1}];
sorted  =  Sort[dists];

tr = 
        Transpose[{sorted, 
            Range[Length[dists]]/
              Length[dists]}];
\(\(ListPlot[tr, 
        PlotStyle -> 
          PointSize[ .006], 
Ticks -> {{1/3, 
              3, {\@2 - 2\/3, \ \@2 - "2"\/3}, {\@5 - 
                  2/3, \ \@5 - "2"\/3}, {\@10 - 
                  2/3, \ \@10 - "2"\/3}, 4\ }, \ {0,  .25, \  .5,  .75, 
              1}}, 
PlotRange -> {{ .25, \ 
              4.0001}, \ {\(- .02\), 
              1.1}}, 
\ \ TextStyle -> {FontFamily\  -> \ "Times", FontSize -> fs}, 
        AxesLabel -> {None, 
            DisplayForm[
              StyleBox[\ "Probability", 
                FontFamily\  -> \ "Times", \ 
                FontSize -> 
                  fsbig]]}, 
AxesOrigin -> { .25, \(- .02\)}, 
Prolog -> {{GrayLevel[0.3], \ 
              Thickness[th/2], 
              Dashing[{ .02, \  .02}], 
\ \(Line[{{#, \(- .02\)}, {#,  .92}}] &\)\  /@ \ {\@2 - 2/3, \@5 - 2/3, \@10 - 2/3, \@14/
                    3, \@41/3}}}, 
        Epilog\  -> {Line\ [{{ .25, \ 
                  1.03}, { .25, \ \(- .03\)}, \ {4.1, \(- .03\)}}]; \ 
            Text[\@14/3, {\@14/3, 1}], Text[\@41/3, {\@41/3, 1}], 
            Text["Length", {4. , \ \(- .0\)}, \ {1, \(-1\)}, \ 
              TextStyle -> {FontFamily\  -> \ "Times", 
                  FontSize -> fsbig}]}];\)\)