(* This computes the action of proper affine parabolic group on a line and two transversal foliations generated by lines--- Figure3 *)
(* linear part of Phi_t *)
Phi1[t_]:= {{1, t, t^2/2}, {0, 1, t}, {0, 0, 1}}

(* Translational part of Phi_t*)

Phi2[mu1_, t_]:= mu1{t^3/6, t^2/2, t}

(* Inner product for the parabolic coordinates *)

M={{0,0,-1}, {0, 1, 0}, {-1, 0, 0}}

(* See S.3.3.2. The following is denoted by f there *)

y0[r_]:= (3/4)*(r/Sqrt[1-r^2])

(* See S.3.3.2. The following is the line there *)

L[s_, r_]:= { r*s, y0[r],  s* Sqrt[1-r^2]}

(* First foliations are given by leaves given as follows *)

LINE[s_,r_, t_]:= Phi1[t].L[s, r] + Phi2[mu, t]

LINE[s,r, t]

(*This gives a section at x= x0. First we solve for s*)

Solve[r s+((3 r t)/(4 (Sqrt[1-(r^2)])))+(1/2) (Sqrt[1-(r^2)]) s (t^2)+((mu (t^3))/6)==x0, s]

(* Now we plug back s to obtain the actual section at x=x0 *)

Simplify[LINE[((-9 r t-2 mu (Sqrt[1-(r^2)]) (t^3)+12 (Sqrt[1-(r^2)]) x0)/(6 (Sqrt[1-(r^2)]) (2 r+(Sqrt[1-(r^2)]) (t^2)))),r, t]]

(* need condition y0 < mu a/c From the paper*)

(*An animation of the section as x0 changes.*)

Animate[ParametricPlot[{((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{r-> 0.99, mu->1}, {t, -10, 10}], {x0, -10, 10}]

(* Another example of the section with different values *)

Animate[ParametricPlot[{{((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{r-> 0.7, mu->1},{((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{r-> 0.995, mu->1}}, {t, -10, 10}], {x0, -10, 10}]

(* We draw three ruled surfaces using three sections   *)

p1=ParametricPlot3D[{{x0,((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{r-> 0.95, mu->1},{x0,((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{r-> 0.995, mu->1}, {x0,((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{r-> 0.9995, mu->1}},{x0, -50, 50}, {t, -30, 30},  PlotStyle->{Red, Magenta, Pink}, BoxRatios->{1,1,1.3}]

(* Now we look at the second foliation transversal to the first one, now t is fixed and r changes*)

LINE[s, r,t]

Animate[ParametricPlot[{((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{t->1, mu->1}, {r, 0.4, 0.99}], {x0, -10, 10}]

(* The five surface for five different sections *)

p2=ParametricPlot3D[{{x0,((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{t-> 0, mu->1},{x0,((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{t-> 3, mu->1}, {x0,((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{t-> 6, mu->1},{x0,((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{t-> -3, mu->1},{x0,((3 (-3+4 mu) r (Sqrt[1-(r^2)]) (t^2)+2 mu (t^4)+24 t x0-2 (r^2) (-9+mu (t^4)+12 t x0))/(12 (2 r (Sqrt[1-(r^2)])+(t^2)-(r^2) (t^2)))),((3 (-3+4 mu) r t+4 (Sqrt[1-(r^2)]) (mu (t^3)+3 x0))/(6 (2 r+(Sqrt[1-(r^2)]) (t^2))))}/.{t-> -6, mu->1}},{x0, -50, 50}, {r, 0.95, 0.99999}, Lighting->{{"Ambient",White}}, PlotStyle->{Green, LightGreen, Darker[Green], LightBlue, Darker[Blue]}, BoxRatios->{1,1,1.3}]

(* The final figure obtained from p1 and p2. This is Figure 3 in the article *)

Show[p1, p2]