Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Needs["Graphics`ImplicitPlot`"];
If[DownValues[Transform2DPlot`Transform2DGraphicsPlot] === {}, 
   Get["C:\\sergeev\\tulane\\math\\packages\\transform2dplot.m"]];
(* Calc. of of mapping (y,t) -> (qx,qy) *)
mapping := Module[
   {file1,file2,file3,iw1,iw2,iw3,rplot,
xmin,xmax,ymin,ymax,theta,mtraj,tm,dy,tmax,yn,y0,ndig,sol,t1,xx,yy,
dt,plt,plt1,plt2,ddq,ddp,pltq,pltp,pltq1,pltp1,gr,endmapping},
file1 = ToFileName[subdirhtml, "map.gif"];
file2 = ToFileName[subdirhtml, "caustics.gif"];
file3 = ToFileName[subdirhtml, "caustyt.gif"];
{iw1,iw2,iw3} = {(file=file1;ifwrite),(file=file2;ifwrite),(file=file3;ifwrite)};
If[!iw1 && !iw2 && !iw3 && !ifdisplay, Goto[endmapping]];
(* Finding r plot-max *)
rplot = If[NumberQ[rpmax], rpmax/Sqrt[2], upper[Function[r, w[r]], 10.^-4]];
{xmin, xmax} = {ymin, ymax} = rplot {-1, 1};
theta[x_] = (1 + Sign[x])/2;
mtraj = 200;
tm = 9999;
dy = ymax/mtraj;
tmax = 0;
Do[yn = (ntraj - 1/2)dy;
    sol = 
      NDSolve[{qx[0] == xmin, px[0] == p0, qy[0] == yn, py[0] == 0, 
            qx'[t] == px[t]/m, qy'[t] == py[t]/m, px'[t] == -wx[qx[t], qy[t]],
             py'[t] == -wy[qx[t], qy[t]]}, {qx, qy, px, py}, {t, 0, tm},
          MaxStepSize -> rplot/100,
          Method -> {EventLocator, 
              "Event" -> 
                theta[xmax - qx[t]]theta[qx[t] - (xmin - 0.001)]theta[
                    ymax - qy[t]]theta[qy[t] - ymin]}][[1]];
    t1 = InterpolatingFunctionDomain[qx /. sol][[1, -1]];
    If[t1 > tmax, tmax = t1];
    , {ntraj, mtraj}];
print[9, "tmax = ", tmax];
tmax = 1.1 tmax;
If[handle === "modgauss", tmax = 180.];(* special abnormal case *)
y0 = -ymax/10000;
ndig = 12;
If[handle === "yukawa", y0 = 0.01; ndig = 8];(* to avoid singularity at y=0 for Yukawa potential *)
sol = NDSolve[{
          qx[y, 0] == xmin, px[y, 0] == p0, qy[y, 0] == y, py[y, 0] == 0, 
          s[y, 0] == 0,
          Derivative[0, 1][qx][y, t] == px[y, t]/m,
          Derivative[0, 1][qy][y, t] == py[y, t]/m,
          Derivative[0, 1][px][y, t] == -wx[qx[y, t], qy[y, t]],
          Derivative[0, 1][py][y, t] == -wy[qx[y, t], qy[y, t]],
          Derivative[0, 1][s][y, t] == (px[y, t]^2 + py[y, t]^2)/m
          }, {qx, qy, px, py, s}, {y, y0, ymax}, {t, 0, tmax},
        Sequence @@ If[iftest === True, {MaxStepFraction -> 1/500},
            {AccuracyGoal -> ndig,
              PrecisionGoal -> ndig,
              MaxStepFraction -> 1/1000(* 2000 - maximum *)}]
        ][[1]];
xx[y_, t_] = qx[Abs[y], t] /. sol;
yy[y_, t_] = Sign[y] qy[Abs[y], t] /. sol;
print[1, "Field of trajectories was calculated"];
If[iw1 || ifdisplay,
mtraj = 50;
dy = (ymax - ymin)/mtraj;
dt = dy/(p0/m);
plt1 = Table[
      yn = ymin + (ntraj - 1/2)dy;
      ParametricPlot[{xx[yn,t], yy[yn,t]}, {t, 0, tmax},
        PlotStyle -> {{RGBColor[0, 0, 1], Thickness[0.002]}},
        TextStyle -> {FontFamily -> "Times", FontSize -> 14},
        FrameLabel -> {"x","y"},
        Frame -> True,
        Release[Sequence @@ If[iftest === True, {PlotPoints -> 99},
              {PlotPoints -> 999,
                PlotDivision -> 99,
                MaxBend -> 1}]],
        DisplayFunction -> Identity], {ntraj, mtraj}];
plt2 = Table[
      ParametricPlot[{xx[y,tn], yy[y,tn]}, {y, ymin, ymax},
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[0.002]}},
        TextStyle -> {FontFamily -> "Times", FontSize -> 14},
        FrameLabel -> {"x","y"},
        Frame -> True,
        Release[Sequence @@ If[iftest === True, {PlotPoints -> 99},
              {PlotPoints -> 999,
                PlotDivision -> 99,
                MaxBend -> 1}]],
        DisplayFunction -> Identity], {tn, 0, tmax, dt}];
plt = Show[plt2, plt1, DisplayFunction -> display,
      ImageSize -> 500,
      Axes -> None,
      PlotRange -> {{xmin, xmax}, {ymin, ymax}},
      Background -> RGBColor[1, 1, 1],
      AspectRatio -> 1];
If[iw1,
    Export[file1, plt, ImageResolution -> 8 72];
    Run["mogrify -scale 500x500 " <> file1];
   ];
];
If[!iw2 && !iw3 && !ifdisplay, Goto[endmapping]];
(*
    Caustics 
*)
ddq[y_, t_] = Det[D[{qx[y, t], qy[y, t]}, {{y, t}}]] /. sol;
ddp[y_, t_] = Det[D[{px[y, t], py[y, t]}, {{y, t}}]] /. sol;
If[iw2 || ifdisplay,
pltq = ImplicitPlot[ddq[y // Abs, t] == 0, {y, ymin, ymax}, {t, 0, tmax},
        PlotPoints -> If[iftest === True, 111, 999], PlotRange -> All,
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[0.005]}},
        AspectRatio -> 1, DisplayFunction -> Identity] // Graphics;
pltp = ImplicitPlot[ddp[y // Abs, t] == 0, {y, ymin, ymax}, {t, 0, tmax},
        PlotPoints -> If[iftest === True, 111, 999], PlotRange -> All,
        PlotStyle -> {{RGBColor[0, 0, 1], Thickness[0.005]}},
        AspectRatio -> 1, DisplayFunction -> Identity] // Graphics;
pltq1 = Transform2DGraphicsPlot[pltq // Graphics, Function[{y, t}, {xx[y, 
          t], yy[y, t]}], ResolutionLength -> 0.1, AspectRatio -> 
    Automatic, Axes -> False,
      TextStyle -> {FontFamily -> "Times", FontSize -> 14},
      FrameLabel -> {"x","y"},
      PlotRange -> {{xmin, xmax}, {ymin, ymax}},
      Frame -> True, DisplayFunction -> Identity];
pltp1 = Transform2DGraphicsPlot[pltp // Graphics, 
          Function[{y, t}, {xx[y, t], yy[y, t]}], 
      ResolutionLength -> 0.1, AspectRatio -> Automatic, Axes -> False,
      TextStyle -> {FontFamily -> "Times", FontSize -> 14},
      FrameLabel -> {"x","y"},
      PlotRange -> {{xmin, xmax}, {ymin, ymax}},
      Frame -> True, DisplayFunction -> Identity];
gr = GraphicsArray[{pltq1, pltp1}];
plt = Show[gr, ImageSize -> 800, DisplayFunction -> display];
If[iw2, Export[file2, plt]];
];
(* in (y,t) coordinates *)
If[iw3 || ifdisplay,
pltq = ContourPlot[ddq[y // Abs, t], {t, 0, tmax}, {y, ymin, ymax},
      PlotPoints -> If[iftest === True, 111, 999], ColorFunction -> (If[# < 0, Hue[0, 0, 1], Hue[0.65, 
      0.05, 0.9]] &),
      ColorFunctionScaling -> False,
      Contours -> 1,
      PlotRange -> {-1, 1},
      FrameLabel -> {"t", "y"},
      TextStyle -> {FontFamily -> "Times", FontSize -> 14},
      ContourStyle -> {RGBColor[1, 0, 0], Thickness[0.005]},
      DisplayFunction -> Identity,
      AspectRatio -> 1, ImageSize -> 400];
pltp = ContourPlot[-ddp[y // Abs, t], {t, 0, tmax}, {y, ymin, ymax},
      PlotPoints -> If[iftest === True, 111, 999], ColorFunction -> (If[# < 0, Hue[0, 0, 1], 
            Hue[0.65, 0.05, 0.9]] &),
      ColorFunctionScaling -> False,
      Contours -> 1,
      PlotRange -> {-1, 1},
      FrameLabel -> {"t", "y"},
      TextStyle -> {FontFamily -> "Times", FontSize -> 14},
      ContourStyle -> {RGBColor[0, 0, 1], Thickness[0.005]},
      DisplayFunction -> Identity,
      AspectRatio -> 1, ImageSize -> 400];
gr = GraphicsArray[{pltq, pltp}];
plt = Show[gr, ImageSize -> 800, DisplayFunction -> display];
If[iw3, Export[file3, plt]];
];
Label[endmapping];
];