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]];\)\)

Digits Agree Count

\(\(\(
\)\(\(DigitsAgreeCount[a_, 
          b_] := \((prec = 
            Ceiling@Min[Precision /@ {a, b}]; {{ad, ae}, {bd, 
                be}} = \(RealDigits[#, 10, prec] &\) /@ {a, b}; 
          If[ae ≠ be ∨ a\ b ≤ 0, Return[0]]; 
          If[ad == bd, Return@Length[ad]]; {{com}} = 
            Position[MapThread[Equal, {ad, bd}], False, 1, 1] - 1; 
          com)\);\)


    \(DigitsAgreeCount[Interval[{a_, b_}]] := 
        DigitsAgreeCount[a, b];\)\)\)\)

Table 2.1 (uses PhotonTrajectory code)

true = 
        Norm[PhotonTrajectory[{(1/2), (1/10)}, {1, 0}, 
              10, 
Method -> FixedPrecision, 
              FixedWorkingPrecision -> 
                200][[-1]]];

\(DisplayForm[
      StyleBox[GridBox[
          Prepend[Table[
{prec, 
aa = 
                  Norm[\(PhotonTrajectory[{1/2, 1/10}, \ {1, 0}, \ 10, 
                        Method -> FixedPrecision, \ 
                        FixedWorkingPrecision -> 
                          prec]\)[\([\(-1\)]\)]]; \ 
NumberForm[aa, prec], 
DigitsAgreeCount[
                  aa, true]}, \ {prec, \ 5, \ 
                30}], {"Precision", \ "Computed Distance", \ "Correct
Digits"}], \ GridFrame -> 1, \ 
          ColumnAlignments -> {Right, \ Left, Right}, \ 
          RowLines -> {1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 
              1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1}], \ 
        Background -> GrayLevel[0.9], \ FontFamily -> Times, \ 
        FontSize -> 9]]\)