(* Plotting potentials *)

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

(* Should be defined:
ifdisplay - whether to display graphs
ifoverwrite - whether to overwrite existing files
*)

potplots := Module[
   {file1,file2,file3,iw1,iw2,iw3,rplot,xmin,xmax,ymin,ymax,
   plt,prange,wx,wy,theta,mtraj,plt1,zmin,zmax,dz,tmax,dy,plt2,
   yn,sol,t1,xt,yt,endpotplots},
file1 = ToFileName[subdirhtml, "pot.gif"];
file2 = ToFileName[subdirhtml, "poten.gif"];
file3 = ToFileName[subdirhtml, "potent.gif"];
{iw1,iw2,iw3} = {(file=file1;ifwrite),(file=file2;ifwrite),(file=file3;ifwrite)};
If[!iw1 && !iw2 && !iw3 && !ifdisplay, Goto[endpotplots]];
(* Finding r plot-max *)
rplot = upper[Function[r, w[r]], 10.^-2];
{xmin, xmax} = {ymin, ymax} = rplot {-1, 1};
plt = Plot[{w[Abs[x]], en}, {x, xmin, xmax},
      PlotPoints -> 99,
      AspectRatio -> 0.4,
      DisplayFunction -> Identity];
prange = PlotRange /. AbsoluteOptions[plt];
(*If[NumberQ[rpmax], prange[[1]] = rpmax{-1, 1}];*)
If[NumberQ[zpmin], prange[[2, 1]] = zpmin];
If[NumberQ[zpmax], prange[[2, 2]] = zpmax];
(* Small plot *)
If[iw1 || ifdisplay,
plt = Plot[{w[Abs[x]], en}, {x, xmin, xmax}, PlotPoints -> 99,
      PlotStyle -> {{RGBColor[0, 0, 1],Thickness[0.02]}, {RGBColor[1, 0, 0], 
            Dashing[{0.03, 0.05}],Thickness[0.02]}},
      PlotRange -> prange,
      Frame -> True,
      FrameStyle -> GrayLevel[0.8],
      Axes -> True,
      Ticks -> None,
      AxesOrigin -> {0, 0},
      FrameTicks -> None,
      AspectRatio -> 0.4,
      ImageSize -> 150,
      DisplayFunction -> display];
];
If[iw1, Export[file1, plt]];
(* Larger plot *)
If[iw2 || ifdisplay,
plt = Plot[{w[Abs[x]], en}, {x, xmin, xmax}, PlotPoints -> 99,
      PlotStyle -> {{RGBColor[0, 0, 1],Thickness[0.01]}, {RGBColor[1, 0, 0], 
            Dashing[{0.01, 0.02}],Thickness[0.01]}},
      TextStyle -> {FontFamily -> "Times", FontSize -> 14},
      PlotRange -> prange,
      Frame -> True,
      FrameStyle -> GrayLevel[0.8],
      Axes -> True,
      FrameLabel -> {"x", "w(x,0)"},
      AxesOrigin -> {0, 0},
      FrameTicks -> None,
      AspectRatio -> 0.6,
      ImageSize -> 400,
      DisplayFunction -> display];
];
If[iw2, Export[file2, plt]];
(* 3D - plot *)
If[iw3 || ifdisplay,
wx[x_, y_] = D[w[x, y], x];
wy[x_, y_] = D[w[x, y], y];
theta[x_] := (1 + Sign[x])/2;
mtraj = 50;
plt1 = Plot3D[w[x, y], {x, xmin, xmax}, {y, ymin, ymax},
      TextStyle -> {FontFamily -> "Times", FontSize -> 14},
      AxesLabel -> {"x", "y", "w"},
      PlotPoints -> 3 mtraj,
      Mesh -> False,
      PlotRange -> All,
      DisplayFunction -> Identity];
{zmin, zmax} = Last[PlotRange /. AbsoluteOptions[plt1]];
dz = 0.01 (zmax - zmin);
tmax = 999;
dy = (ymax-ymin)/mtraj;
plt2 = Table[
      yn = ymin + (ntraj - 1/2)dy;
      sol = NDSolve[{
              qx[0] == xmin, px[0] == p[xmin, 0],
              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, tmax},
            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]];
      ParametricPlot3D[
        {(xt = qx[t] /. sol), (yt = qy[t] /. sol), w[xt, yt] + dz,
          {RGBColor[0, 0, 0.7], Thickness[0.002]}}, {t, 0, t1},
      TextStyle -> {FontFamily -> "Times", FontSize -> 14},
      AxesLabel -> {"x", "y", "w"},
        PlotRange -> All,
        DisplayFunction -> Identity], {ntraj, mtraj}];
plt = Show[plt1, plt2,
      DisplayFunction -> display,
      ImageSize -> 450];
];
If[iw3, Export[file3, plt]];
Label[endpotplots];
];