\(\(zx = {#\_\(\([[\)\(3\)\(]]\)\), #\_\(\([[\)\(1\)\(]]\)\)} &;\)\)
\(\(yx = {#\_\(\([[\)\(2\)\(]]\)\), #\_\(\([[\)\(1\)\(]]\)\)} &;\)\)

<< Graphics`Shapes`

\(pp3DtoLine[g_] := 
      Line[Append[
          Map[#\_\(\([[\)\(1, 1\)\(]]\)\) &, 
            g\_\(\([[\)\(1\)\(]]\)\)], 
          g\_\(\([[\)\(1, \(-1\), 1, 2\)\(]]\)\)]]\)

Spherical

H[r_, k_, 
        b_] := (1/2) 
          k . k - (1/2) ((2 Norm[r] - 
                1)/Norm)[r]^4 (((r . k))^2) - (1/2) (((((b (Norm[r] - 1))/(Norm[r] (b - 1)))))^2)
dHdk[r_, k_, b_] := 
      k - ((2 Norm[r] - 1)/Norm)[r]^4 (r . 
              k) r 
dHdr[r_, k_, 
        b_] := (1/Norm)[r]^4 (-(2 Norm[r] - 1) (r . k) 
              k + ((3 Norm[r] - 2)/Norm)[r]^2 (((r . k))^2) 
              r - ((((b/(b - 1))))^2) (Norm[r] - 1) 
              r)

ref[ r_, k0_, b_] := 
      Block[{k, kx, ky, kz, n, 
          kL}, 
k = {kx, ky, kz}; 
n = 
          (r/Norm)[r]; 
kL = 
          Map[Chop[k /. #, (10^(-9))] &, 
            Solve[{ H[r, k, b] == 0 , 
                Drop[(k - (k . n) n) - (k0 - (k0 . 
                                n) n), -1] == {0, 0}}, 
              k]]; 
kL  =  
          Select[kL, (Im[#] == {0, 0, 
                    0}) &]; 
kL  =  
          Select[kL, (dHdk[r, #, b] . n < 
                  0) &]; 
(Subscript[kL, [[1]]])
]

eqs[r0_, k0_, 
        b_] := {
D[r[τ], τ] == 
          dHdk[r[τ], k[τ], b], 
D[
            k[τ], τ] == -dHdr[r[τ], k[τ], 
              b], 
r[0] == r0, 
        k[0] == k0
}

fromSpherical[{r_,  θ_,  ϕ_} ] :=  
      r {Sin[θ] Cos[ϕ], Sin[θ] Sin[ϕ], 
          Cos[θ]}

b = 2;
r0L = 
        Table[fromSpherical[{b,  π + ArcSin[x],  
              0} ], {x, - .975,  .975,  .05}];
Timing[rfL = 
          Map[ (r[τ] /. (Subscript[(NDSolve[eqs[#, ref[#, {0, 0, 1}, b], b], {r[τ], k[τ]}, {τ, 0, 100}, StoppingTest -> (Norm[r[τ]] > 1.01 b)]), [[1]]])) &, r0L];]
\(\(tL = 
        Map[τ /. 
              FindRoot[
                Norm[#] == 
                  b, {τ,  .9\ #\_\(\([[\)\(0, 1, 1, 2\)\(]]\)\), #\_\(\([[\)\(0, 1, 1, 2\)\(]]\)\)}] &, rfL];\)\)
gRays = 
        MapThread[
          ParametricPlot[ zx[#1], {τ, 0, #2}, 
              DisplayFunction -> Identity] &, {rfL, 
            tL}];
gCircles = 
        Graphics[ {Circle[{0, 0}, b], 
            Circle[{0, 0}, 1]}];
Show[gRays, gCircles, DisplayFunction -> $DisplayFunction, 
      AspectRatio -> Automatic, 
      PlotRange -> (b +  .01) {{-1, 1}, {-1, 
              1}}]

b = 2;
r0L = 
        Table[fromSpherical[{b,  π + ArcSin[x],  
              0 } ], {x, - .975,  .975,  .1}];
Timing[rfL = 
          Map[ (r[τ] /. (Subscript[(NDSolve[eqs[#, ref[#, fromSpherical[{1,  - .08,  π/2} ], b], b], {r[τ], k[τ]}, {τ, 0, 100}, StoppingTest -> (Norm[r[τ]] > 1.01 b)]), [[1]]])) &, 
            r0L];]
\(\(tL = 
        Map[τ /. 
              FindRoot[
                Norm[#] == 
                  b, {τ,  .9\ #\_\(\([[\)\(0, 1, 1, 2\)\(]]\)\), #\_\(\([[\)\(0, 1, 1, 2\)\(]]\)\)}] &, rfL];\)\)
gRays = 
        MapThread[
          ParametricPlot3D[ #1, {τ, 0, #2}, 
              DisplayFunction -> Identity] &, {rfL, 
            tL}];
Show[gRays, Graphics3D[Prepend[Sphere[1, 50, 50], EdgeForm[]]], 
      DisplayFunction -> $DisplayFunction, 
      PlotRange -> 
        2 {{-1, 1}, {-1, 1}, {-1, 1}}]

3D Sphere Figure

n = 20;
m = 80;
θ = 3. ;
zmax = 2.5;
b = 2;
r0L = 
        Table[fromSpherical[{b,  θ,  ϕ } ], {ϕ, 0, 
            2 π (1 - 1/n), 
            2 π/n}];
Timing[rfL = 
          Map[ (r[τ] /. (Subscript[(NDSolve[eqs[#, ref[#, {0, 0, 1}, b], b], {r[τ], k[τ]}, {τ, 0, 100}, StoppingTest -> (Norm[r[τ]] > 1.01 b)]), [[1]]])) &, r0L];]
\(\(tL = 
        Map[τ /. 
              FindRoot[
                Norm[#] == 
                  b, {τ,  .9\ #\_\(\([[\)\(0, 1, 1, 2\)\(]]\)\), #\_\(\([[\)\(0, 1, 1, 2\)\(]]\)\)}] &, rfL];\)\)
gRays = 
        MapThread[
          ParametricPlot3D[ #1, {τ, 0, #2}, 
              DisplayFunction -> Identity] &, {rfL, 
            tL}];
hemi = (Subscript[(ParametricPlot3D[b {Sin[θ] Cos[ϕ], Sin[θ] Sin[ϕ], Cos[θ]}, {θ, 0, π}, {ϕ, -π/2, π/2}, PlotPoints -> m, DisplayFunction -> Identity]), [[1]]]);
lines  =  
        Join[
Map[
            Line[{ReplacePart[#, -zmax, 3], #}] &, 
            r0L], 
Map[
            Line[{ReplacePart[#, b Cos[π - θ], 3], 
                  ReplacePart[#, zmax, 3]}] &, 
            r0L]
];
ψ = 4.1;
g = 
      Show[gRays, 
        Graphics3D[Join[{EdgeForm[]}, hemi, lines, Sphere[1, m, m]]], 
        DisplayFunction -> $DisplayFunction, 
        PlotRange -> {2 {-1, 1}, b {-1, 1}, zmax {-1, 1}}, 
        ViewVertical -> {0, 1, 0}, 
        ViewPoint -> {2 Cos[ψ], 1, 2 Sin[ψ]}, 
        Boxed -> False, Axes -> False]

anim = 
        Table[
gPhotons = 
            Graphics3D[
              Join[{RGBColor[0, 1,  .5], AbsoluteThickness[5]}, 
                Map[pp3DtoLine[
                      ParametricPlot3D[ #, {τ, τ1, τ1 +  .3}, 
                        DisplayFunction -> Identity]] &, 
                  rfL]]]; 
Show[gRays, gPhotons, 
            Graphics3D[Join[{EdgeForm[]}, hemi, lines, Sphere[1, m, m]]], 
            DisplayFunction -> $DisplayFunction, 
            PlotRange -> {2 {-1, 1}, b {-1, 1}, 
                zmax {-1, 1}}, ViewVertical -> {0, 1, 0}, 
            ViewPoint -> {2 Cos[ψ], 1, 2 Sin[ψ]}, 
            Boxed -> False, 
            Axes -> False]
, {τ1, 0, 
            (Subscript[tL, [[1]]]) -  .3, (((Subscript[tL, [[1]]]) -  .3)/40)}];

Export["photon.gif", anim, "GIF", ImageSize -> 800, 
        ConversionOptions -> {"AnimationDisplayTime" -> 
              0.05, "Loop" -> True}];

Cylindrical

T=(({{1, 0, 0}, {0, 1, 0}, {0, 0, 0}}));
Z=(({{0, 0, 0}, {0, 0, 0}, {0, 0, 1}}));

H[ρ_, k_, 
        b_] := (1/2) 
          k . T . k + (1/2) (((((b (Norm[ρ] - 1))/(Norm[ρ] (b - 1)))))^2) (k . Z . k - 
              1) - (1/2) ((2 Norm[ρ] - 
                1)/Norm)[ρ]^4 (((ρ . k))^2)
dHdk[ρ_, k_, b_] :=  
      T . k + (((((b (Norm[ρ] - 1))/(Norm[ρ] (b - 1)))))^2) 
          Z . k - ((2 Norm[ρ] - 1)/Norm)[ρ]^4 (ρ . 
              k) ρ
dHdr[ρ_, k_, 
        b_] := (1/Norm)[ρ]^4 (((((b/(b - 1))))^2) (Norm[ρ] - 1) (k . Z . k - 
                  1) ρ + ((3  Norm[ρ] - 
                    2)/Norm)[ρ]^2 (((ρ . 
                      k))^2) ρ - (2 Norm[ρ] - 
                  1) (ρ . k) T . k)

eqs[r0_, k0_, 
        b_] := {
D[r[τ], τ] == 
          dHdk[T . r[τ], k[τ], b], 
D[
            k[τ], τ] == -dHdr[T . r[τ], k[τ], 
              b], 
r[0] == r0, 
        k[0] == k0
}

ref[ r_, k0_, b_] := 
      Block[{k, kx, ky, kz, n, 
          kL}, 
k = {kx, ky, kz}; 
n = 
          T . (r/Norm)[T . r]; 
kL = 
          Map[Chop[k /. #, (10^(-9))] &, 
            Solve[{ H[T . r, k, b] == 0 , 
                Drop[(k - (k . n) n) - (k0 - (k0 . 
                                n) n), {2}] == {0, 0}}, 
              k]]; 
kL  =  
          Select[kL, (Im[#] == {0, 0, 
                    0}) &]; 
kL  =  
          Select[kL, (dHdk[T . r, #, b] . n < 
                  0) &]; 
(Subscript[kL, [[1]]])
]

refOut[ r_, k0_] := 
      Block[{k, kx, ky, kz, n, 
          kL}, 
k = {kx, ky, kz}; 
n = 
          T . (r/Norm)[T . r]; 
kL = 
          Map[Chop[k /. #, (10^(-6))] &, 
            Solve[{ k . k == 1 , 
                Drop[(k - (k . n) n) - (k0 - (k0 . 
                                n) n), {2}] == {0, 0}}, 
              k]]; 
kL  =  
          Select[kL, (Im[#] == {0, 0, 
                    0}) &]; 
kL  =  
          Select[kL, (# . n > 
                  0) &]; 
(Subscript[kL, [[1]]])
]

fromCylindrical[{ρ_, ϕ_,  
          z_} ] :=  {ρ Cos[ϕ], ρ Sin[ϕ], 
        z}

b = 2;
r0L = 
        Table[fromCylindrical[{b,  -ArcCos[x],  
              0} ], {x, - .95,  .95,  .1}];
k0L  =  Map[ {0, 1, 0} &, r0L];
{rfL, kfL} = 
        Transpose[
          MapThread[ ({r[τ], 
                    k[τ]} /. (Subscript[(NDSolve[eqs[#1, ref[#1, #2, b], b], {r[τ], k[τ]}, {τ, 0, 100}, StoppingTest -> (Norm[T . r[τ]] > 1.01 b)]), [[1]]])) &, {r0L, k0L}]];
\(\(tL = 
        Map[τ /. 
              FindRoot[
                Norm[T . #] == 
                  b, {τ,  .9\ #\_\(\([[\)\(0, 1, 1, 2\)\(]]\)\), #\_\(\([[\)\(0, 1, 1, 2\)\(]]\)\)}] &, rfL];\)\)
r1L = 
        MapThread[(#1 /. τ -> #2) &, {rfL, 
            tL}];
k1L = 
        MapThread[
          refOut[#1 /. τ -> #3, #2 /. τ -> #3] &, {rfL, 
            kfL, tL}];
gLines = 
        Graphics[
          Join[ MapThread[Line[{yx[#1], yx[#1 - 3  #2]}] &, {r0L,  k0L}], 
            MapThread[
              Line[{yx[#1], yx[#1 + 3  #2]}] &, {r1L,  
                k1L}]]];
gCircles = 
        Graphics[ {Circle[{0, 0}, b], 
            Circle[{0, 0}, 1]}];
gRays = 
        MapThread[
          ParametricPlot[ yx[#1], {τ, 0, #2}, 
              DisplayFunction -> Identity] &, {rfL, 
            tL}];
Show[gLines, gRays, gCircles, DisplayFunction -> $DisplayFunction, 
      AspectRatio -> Automatic, 
      PlotRange -> (b +  .5) {{-1, 1}, {-1, 1}}]

b = 2; ψ =  .7;
r0L = 
        Table[b {x, -(Sqrt[1 - (x^2)]), -Tan[ψ] (Sqrt[1 - (x^2)])}, {x, - .975,  .975,  .1}];
k0L  =  
        Map[ {0, Cos[ψ], Sin[ψ]} &, 
          r0L];
{rfL, kfL} = 
        Transpose[
          MapThread[ ({r[τ], 
                    k[τ]} /. (Subscript[(NDSolve[eqs[#1, ref[#1, #2, b], b], {r[τ], k[τ]}, {τ, 0, 100}, StoppingTest -> (Norm[T . r[τ]] > 1.01 b)]), [[1]]])) &, {r0L, k0L}]];
\(tL = 
      Map[τ /. 
            FindRoot[
              Norm[T . #] == 
                b, {τ,  .9\ #\_\(\([[\)\(0, 1, 1, 2\)\(]]\)\), #\_\(\([[\)\(0, 1, 1, 2\)\(]]\)\)}] &, rfL]; Max[tL]\)
r1L = 
        MapThread[(#1 /. τ -> #2) &, {rfL, 
            tL}];
k1L = 
        MapThread[
          refOut[#1 /. τ -> #3, #2 /. τ -> #3] &, {rfL, 
            kfL, tL}];
n = 50; l = 5;
hemicylinder = 
        Table[Polygon[{{b Cos[θ], 
                b Sin[θ], -l/2}, {b Cos[θ], 
                b Sin[θ], l/2}, {b Cos[θ + π/n], 
                b Sin[θ + π/n], 
                l/2}, {b Cos[θ + π/n], 
                b Sin[θ + π/n], -l/2}, {b Cos[θ], 
                b Sin[θ], -l/2}}], {θ, π/2, 
            3 π/2 - π/n, π/n}];
gLines = 
      Graphics3D[
        Join[ MapThread[Line[{#1, #1 - 3  #2}] &, {r0L,  k0L}], 
          MapThread[
            Line[{#1, #1 + 3  #2}] &, {r1L,  
              k1L}]]]
gRays = 
        MapThread[
          ParametricPlot3D[ #1, {τ, 0, #2}, 
              DisplayFunction -> Identity] &, {rfL, 
            tL}];
g1 = 
      Show[gLines, gRays, 
        Graphics3D[Join[{EdgeForm[]}, Cylinder[1, 2.5, 100], hemicylinder]], 
        DisplayFunction -> $DisplayFunction, 
        PlotRange -> 2.5 {{-1, 1}, {-1, 1}, {-1, 1}}, 
        ViewPoint -> {9,  0,  3}]
g2 = 
      Show[gLines, gRays, 
        Graphics3D[Join[{EdgeForm[]}, Cylinder[1, 2.5, 100], hemicylinder]], 
        DisplayFunction -> $DisplayFunction, 
        PlotRange -> 2.5 {{-1, 1}, {-1, 1}, {-1, 1}}, 
        ViewPoint -> {1,  0,  9}]