Needs["Graphics`ImplicitPlot`"];

exact := Module[
{(*filesave,lfile1,lfile2,lfile3,lfile3a,file1,file2,file3,file3a,iw1,iw2,iw3,iw3a,
 pnts,mqmax,mser,eps,epsdelta,epspsi,preclost,rmax1,mplots,imgsz,
 plotdelta,plt1,plt2,plt,psipart,epsm,delta,psixy,theta,res,rtab,
 mr,wr,wrser,psi,zero,zser,asol,psi1,psi2,dif1,dif2,dpsi,ddif1,ddif2,
 psilist,xtab,mx,psir, difr1,difr2, dpsir, ddifr1, ddifr2, psilistr,
 relerr,drelerr,err,clist,clost,r0,psir0,dpsir0,mc1,mcm,mser1,r1,dsol,
 dpsi1,psij,psin,psijr,dpsijr,psinr,dpsinr,csol,cjr,cnr,cr,cx,cy,
 mq0,r0f,crf,deltaf,psi0f,psi1f,psi2f,rminf,rmin,delta0,delta1,mqplot,
 deltadata,xyplot,rplot,dxy,time,rstep,steps,datarad,dt,mqr,
 data1,data,zmax,cline,pictfile,c,tiffiles*)},
filesave = ToFileName[subdirtemp, "psi.mx"];
lfile1 = "delta.gif";
lfile2 = "exactabs.gif";
lfile3 = "exactre.gif";
lfile3a = "exactrea.gif";
lfile4 = "exctabs.gif";
lfile5 = "exctre.gif";
lfile5a = "exctrea.gif";
file1 = ToFileName[subdirhtml, lfile1];
file2 = ToFileName[subdirhtml, lfile2];
file3 = ToFileName[subdirhtml, lfile3];
file3a = ToFileName[subdirhtml, lfile3a];
file4 = ToFileName[subdirhtml, lfile4];
file5 = ToFileName[subdirhtml, lfile5];
file5a = ToFileName[subdirhtml, lfile5a];
{iw1,iw2,iw3,iw3a,iw4,iw5,iw5a} = {(file=file1;ifwrite),(file=file2;ifwrite),(file=file3;ifwrite),(file=file3a;ifwrite),
(file=file4;ifwrite),(file=file5;ifwrite),(file=file5a;ifwrite)};
If[!iw1 && !iw2 && !iw3 && !iw3a && !iw4 && !iw5 && !iw5a && !ifdisplay, Goto[endexact]];

If[iftest === True,
   pnts = 100;
   mqmax = 222,
   pnts = 1000;
   mqmax = 300];

mser = 111;
eps = 10.^-32;
epsdelta = -1;
epspsi = 10.^-8;
preclost = 1000;
rmax1 = If[rpmax>rmax, rpmax, rmax];
mplots = 12;
imgsz = "500x500";
If[OddQ[pnts], Print["Error: pnts should be even! Exiting..."]; Exit[]];

plotdelta := If[iw1 || ifdisplay,
plt1 = ListPlot[deltadata,
      PlotJoined -> True,
      PlotStyle -> {RGBColor[0.8, 0.8, 1]},
      DisplayFunction -> Identity];
plt2 = ListPlot[deltadata,
      PlotJoined -> False,
      PlotStyle -> {RGBColor[0, 0, 1]},
      DisplayFunction -> Identity];
plt = Show[plt1, plt2,
      PlotRange -> All,
      TextStyle -> {FontFamily -> "Times", FontSize -> 14},
      Frame -> True,
      FrameStyle -> GrayLevel[0.8],
      Axes -> True,
      FrameLabel -> {"m", "\[Delta]"},
      AxesOrigin -> {0, 0},
      AspectRatio -> 0.6,
      ImageSize -> 400,
      DisplayFunction -> display];
If[iw1, Export[file1, plt]];
];

If[FileType[filesave] =!= File, (* start calculations *)
(* See Adhikari, Am.J.Phys.54 (4), 362 (1986) *)
psipart := 
  If[r < rminf[mq], 0,
      epsm = If[mq == 0, 1, 2];
      delta = If[mq > mq0, 0, deltaf[mq]];
      epsm*I^mq *If[mq > mq0, BesselJ[mq, p0 r],
          If[r < r0f[mq], psi0f[mq], If[r < rmax1, psi1f[mq], psi2f[mq]]]/crf[mq]] Exp[-I delta]
    ];
psixy[x_Real, y_Real] := (
      r = Sqrt[x^2 + y^2];
      theta = If[r == 0, 0, ArcCos[x/r]];
      If[y < 0, theta = -theta];
      res = Sum[psipart* Cos[mq theta], {mq, 0, mq1}];
      r =.; theta =.;
      res);
rtab = samples[0, 1000, 0.01, 1.01];
mr = Length[rtab];
Do[
    ClearAll[a, r];
    wr = 2 m(en - w[r]) - mq^2/r^2;
    wrser = Series[wr, {r, 0, mser}];
    psi = r^mq(1 + Sum[a[i]r^i, {i, mser}] + O[r]^(mser + 1));
    zero = Chop[1/r D[r D[psi, r], r] + wrser psi]/r^(mq - 1);
    zser = CoefficientList[zero, r];
    asol = Solve[zser == Table[0, {mser}]][[1]];
    psi = Normal[psi /. asol];
    psi1 = Normal[psi + O[r]^(mq + mser - 1)];
    psi2 = Normal[psi + O[r]^(mq + mser - 4)];
    dif1 = psi1 - psi;
    dif2 = psi2 - psi;
    dpsi = D[psi, r];
    ddif1 = D[dif1, r];
    ddif2 = D[dif2, r];
    psilist = List @@ psi;
    (* Finding r0 *)
    xtab = samples[0, 1000, 0.01, 1.01];
    mx = Length[xtab];
    Do[r = xtab[[nx]];
      {psir, difr1, 
      difr2, dpsir, ddifr1, ddifr2, psilistr} = {psi, dif1, dif2, dpsi, 
    ddif1, ddif2, psilist};
      relerr = Abs[{difr1, difr2}/psir];
      drelerr = Abs[{ddifr1, ddifr2}/dpsir];
      err = Max[relerr, drelerr];
      clist = Abs[psilistr/psir];
      clost = Max[clist];
      If[err > eps || clost > preclost, Break[]];
      r0 = xtab[[nx]], {nx, mx}];
    If[r0 > rmax1,
      print[0, "Decreasing r0 to 0.99 rmax1", {mq, r0, rmax1}];
      r0 = 0.99 rmax1;
      ];
    r = r0;
    psir0 = psi;
    dpsir0 = dpsi;
    r =.;
    mcl = clist // Length;
    mcm = mc1;
    Do[mcm = n;
      If[clist[[n]] > eps, Break[]];
      , {n, mcl, 1, -1}];
    mser1 = psilist[[mcm]] /. c_ r^n_ -> n;
    print[6, "mq = ", mq, ",   r0 = ", r0, ",   psi(r0) = ", psir0, ",   eff. 
    order = ", mser1 - mq, "   {err,clost} = ", {err, clost}];
    psi0 = psi;
    (*Continuing integration from r0 to r1 = rmax1*)
    r1 = rmax1;
ClearAll[psi, r];
dsol = NDSolve[{1/r D[r D[psi[r], r], r] + wr psi[r] == 0,
            psi[r0] == psir0, psi'[r0] == dpsir0}, psi, {r, r0, r1}, 
          AccuracyGoal -> 12, PrecisionGoal -> 12,
          MaxStepFraction -> 1/1000,
          (*StartingStepSize -> 0.00000001,*)
          MaxSteps -> 999999][[1]];
psi1 = psi[r] /. dsol;
dpsi1 = D[psi1, r];
    r = rmax1;
{psir, dpsir} = {psi1, dpsi1};
    r =.;
psij = BesselJ[mq, p0 r];
psin = BesselY[mq, p0 r];
dpsij = D[psij, r];
dpsin = D[psin, r];
r = rmax1;
{psijr, dpsijr} = {psij, dpsij};
{psinr, dpsinr} = {psin, dpsin};
r =.;
    Check[
csol = Solve[{cj psijr + cn psinr == 
        psir, cj dpsijr + cn dpsinr == dpsir}, {cj, cn}][[1]];
{cjr, cnr} = {cj, cn} /. csol,
      {cjr, cnr} = {psir/psijr, 0}];
cr = Sqrt[cjr^2 + cnr^2];
{cx, cy} = {cjr, cnr}/cr;
delta = ArcCos[cx];
If[cy < 0, delta = -delta];
psi2 = cjr psij + cnr psin;
    {mq0, r0f[mq], crf[mq], deltaf[mq], psi0f[
    mq], psi1f[mq], psi2f[mq], rminf[mq]} = {mq, r0, cr, delta, psi0, psi1, psi2, 0};
    If[mq > 10 && Abs[delta] < epsdelta && Abs[deltaf[mq - 1]] < epsdelta &&
       Abs[deltaf[mq - 2]] < epsdelta, Break[]];
(* Finding Rmin *)
rmin = 0;
Do[r = rtab[[nr]];
   If[Abs[psipart] > epspsi, Break[]];
   rmin = rtab[[nr]], {nr, mr}];
r=.;    
rminf[mq] = rmin;      
print[1, "mq = ",mq,"  ->  delta = ", delta, "   {rmin, r0} = ", {rmin, r0}];
    , {mq, 0, mqmax}];
print[0, "mq-max(delta=0) = ", mq0];
delta0 = delta1 = 0;
mqplot =.;
deltadata = Table[
        deltaest = 2 delta0 - delta1;
        delta = deltaf[mq];
        If[! NumberQ[mqplot] && Abs[delta] > 10.^-2, mqplot = mq];
        n = Round[(delta - deltaest)/2/Pi];
        delta1 = delta0;
        delta0 = delta - 2 Pi n;
        {mq, delta0}
        , {mq, mq0, 0, -1}] // Reverse;
deltadata = Take[deltadata, mqplot];
plotdelta;
ClerAll[psiex, psixy, rr];
xyplot = If[NumberQ[rpmax], rpmax, rmax]/Sqrt[2]; (* may be set larger, 
    but computations will be much slower! *)
rplot = xyplot Sqrt[2];
(* Finding m - max *)
mq1 = mqmax;
Do[
    If[rminf[mq] < 1.1 rplot, Break[]];
    mq1 = mq;
    , {mq, mqmax, 0, -1}];
If[mq1 == mqmax, Print["Error: number of waves is too small to give 
    necessary accuracy!"]];
print[0, "Number of partial waves: ", mq1];
dxy = 2 xyplot/pnts;
time = Timing[
        data = Table[Table[
                xn = -xyplot + (nx - 1/2)dxy;
                yn = -xyplot + (ny - 1/2)dxy;
                psixy[xn, yn], {nx, pnts}], {ny, pnts/2, pnts}];
        ][[1]];
print[0, "Exact wave function calc. time: ", time];
(* Wavefunction in polar coordinates *)
rstep = rplot/300;
steps = Table[10^n{1, 2, 5}, {n, -9, 9}] // Flatten;
rstep = Sort[steps, Abs[#1 - rstep] < Abs[#2 - rstep] &][[1]];
datarad = Table[
      r = N[r1];
      dt = Table[Chop[psipart], {mq, 0, mq0}];
      mqr = 0;
      Do[
        If[dt[[mq + 1]] =!= 0,
          mqr = mq;
          Break[]], {mq, mq0, 0, -1}];
      Take[dt, mqr + 1], {r1, 0, rplot, rstep}];
r =.;
DumpSave[filesave, {mqmax, mq0, mq1, mqplot, deltadata, mser, eps, 
                epsdelta, epspsi, rmax, rmax1, rplot, xyplot, pnts, data, rstep, datarad}],
(* if filesave exists *)
Get[filesave];
plotdelta;
];             

data1 = Reverse[data];
data = Join[data1, data];
zmax = Max[Abs[data]];
data = data/zmax;
If[iw2 || ifdisplay,
plt = ListDensityPlot[Abs[data],
      Mesh -> False,
      ColorFunction -> (GrayLevel[Sqrt[#]] &),
      ColorFunctionScaling -> False,
      MeshRange -> xyplot{{-1, 1}, {-1, 1}},
      AspectRatio -> 1,
      ImageSize -> 433,
      DisplayFunction -> display];
If[iw2,
   SetDirectory[ToFileName[{dirhtml, handle}]];
   Export["tmp.tif", plt, ImageResolution -> 183];
   Run["convert -scale " <> imgsz <> " tmp.tif " <> lfile2];
   DeleteFile["tmp.tif"];
   UnsetDirectory[]; 
 ];
];
If[iw3 || ifdisplay,
plt = ListDensityPlot[Re[data],
      Mesh -> False,
      ColorFunction -> (Hue[If[# > 0, 0, 0.65], 1/4, Sqrt[Abs[#]]] &),
      ColorFunctionScaling -> False,
      MeshRange -> xyplot{{-1, 1}, {-1, 1}},
      AspectRatio -> 1,
      ImageSize -> 433,
      DisplayFunction -> display];
If[iw3,
   SetDirectory[ToFileName[{dirhtml, handle}]];
   Export["tmp.tif", plt, ImageResolution -> 183];
   Run["convert -scale " <> imgsz <> " tmp.tif " <> lfile3];
   DeleteFile["tmp.tif"];
   UnsetDirectory[]; 
 ];
];
If[iw3a,
dt = Round[100/mplots];
cline = "convert -scale " <> imgsz <> " -delay " <> ToString[dt];
SetDirectory[ToFileName[{dirhtml, handle}]];
Do[
    print[2, "Plotting in progress: ", np, "/", mplots];
    cline = cline <> " " <> ToString[np] <> ".tif";
    pictfile = ToString[np] <> ".tif";
    c = Exp[-(np - 1)/mplots 2 Pi I];
    plt = ListDensityPlot[Re[c data],
        Mesh -> False,
        ColorFunction -> (Hue[If[# > 0, 0, 0.65], 1/4, Sqrt[Abs[#]]] &),
        ColorFunctionScaling -> False,
        MeshRange -> xyplot{{-1, 1}, {-1, 1}},
        AspectRatio -> 1,
        ImageSize -> 433,
        DisplayFunction -> Identity];
    Export[pictfile, plt, ImageResolution -> 183];
    , {np, mplots}];
cline = cline <> " -loop 9999 " <> lfile3a;
Print["Running:   ", cline];
Run[cline];
tiffiles = Table[ToString[np] <> ".tif", {np, mplots}];
DeleteFile /@ tiffiles;
ResetDirectory[];
];
(* Simplified plots *)
If[iw4 || iw5 || iw5a,
   dxy = 2 xyplot/pnts;
   dataxy = Table[
      ny1 = Abs[ny - pnts/2] + 1;
      Table[
        xn = -xyplot + (nx - 1/2)dxy;
        yn = -xyplot + (ny - 1/2)dxy;
        {xn, yn, data[[ny, nx]]}, {nx, pnts}], {ny, pnts}];
   dataxy = Flatten[dataxy, 1];
   psixyi = Interpolation[dataxy];
   ClearAll[psixy, data, dataxy,datarad];
  ];
npts = pnts;
Off[InterpolatingFunction::dmval];
If[iw4,
plt = ContourPlot[
      Abs[zmax psixyi[x, y]],
      {x, -xyplot, xyplot}, {y, -xyplot, xyplot},
      PlotRange -> {0, 3},
      Contours -> 2,
      ColorFunctionScaling -> False,
      ColorFunction -> (If[# < 1, GrayLevel[0.9], 
              If[# < 2, GrayLevel[1], RGBColor[1, 0, 0]]] &),
      PlotPoints -> npts,
      ImageSize -> 433,
      Frame -> True,
      Axes -> False,
      DisplayFunction -> display];
SetDirectory[ToFileName[{dirhtml, handle}]];
Export["tmp.tif", plt, ImageResolution -> 256];
Run["convert -scale 700x700 tmp.tif " <> lfile4];
DeleteFile["tmp.tif"];
UnsetDirectory[];
];
(*
If[iw5,
Off[InterpolatingFunction::"dmval"];
pltre = ImplicitPlot[
      Re[psixyi[x, y]] == 0,
      {x, -xyplot, xyplot}, {y, -xyplot, xyplot},
      PlotStyle -> RGBColor[1, 0, 0],
      PlotPoints -> npts,
      Frame -> True,
      DisplayFunction -> Identity];
pltim = ImplicitPlot[
      Im[psixyi[x, y]] == 0,
      {x, -xyplot, xyplot}, {y, -xyplot, xyplot},
      PlotStyle -> RGBColor[0, 0, 1],
      PlotPoints -> npts,
      Frame -> True,
      DisplayFunction -> Identity];
plt = Show[pltre, pltim,
      ImageSize -> 433,
      Frame -> True,
      Axes -> False,
      DisplayFunction -> display];
SetDirectory[ToFileName[{dirhtml, handle}]];
Export["tmp.tif", plt, ImageResolution -> 256];
Run["convert -scale 700x700 tmp.tif " <> lfile5];
DeleteFile["tmp.tif"];
UnsetDirectory[];
];
If[iw5a,
dt = Round[100/mplots];
cline = "convert -scale 700x700 -delay " <> ToString[dt];
SetDirectory[ToFileName[{dirhtml, handle}]];
Do[
    print[2, "Plotting in progress: ", np, "/", mplots];
    cline = cline <> " " <> ToString[np] <> ".tif";
    pictfile = ToString[np] <> ".tif";
    c = Exp[-(np - 1)/mplots Pi I];
pltre = ImplicitPlot[
      Re[c psixyi[x, y]] == 0,
      {x, -xyplot, xyplot}, {y, -xyplot, xyplot},
      PlotStyle -> RGBColor[1, 0, 0],
      PlotPoints -> npts,
      Frame -> True,
      PlotRange -> {{-xyplot, xyplot}, {-xyplot, xyplot}, Automatic},
      DisplayFunction -> Identity];
pltim = ImplicitPlot[
      Im[c psixyi[x, y]] == 0,
      {x, -xyplot, xyplot}, {y, -xyplot, xyplot},
      PlotStyle -> RGBColor[0, 0, 1],
      PlotPoints -> npts,
      Frame -> True,
      PlotRange -> {{-xyplot, xyplot}, {-xyplot, xyplot}, Automatic},
      DisplayFunction -> Identity];
plt = Show[pltre, pltim,
      ImageSize -> 433,
      Frame -> True,
      Axes -> False,
      DisplayFunction -> Identity];
    Export[pictfile, plt, ImageResolution -> 256];
    , {np, mplots}];
cline = cline <> " -loop 9999 " <> lfile5a;
Print["Running:   ", cline];
Run[cline];
tiffiles = Table[ToString[np] <> ".tif", {np, mplots}];
DeleteFile /@ tiffiles;
ResetDirectory[];
];
*)
If[iw5,
Off[InterpolatingFunction::"dmval"];
plt = ContourPlot[
      Re[zmax psixyi[x, y]],
      {x, -xyplot, xyplot}, {y, -xyplot, xyplot},
      PlotRange -> {-4, 4},
      Contours -> 3,
      ColorFunctionScaling -> False,
      ColorFunction -> (Hue[If[# > 0, 0, 0.65], If[Abs[#] < 2, 1/4, 1], 1] &),
      PlotPoints -> npts,
      ImageSize -> 433,
      Frame -> True,
      Axes -> False,
      DisplayFunction -> display];
SetDirectory[ToFileName[{dirhtml, handle}]];
Export["tmp.tif", plt, ImageResolution -> 233];
Run["convert -scale 700x700 tmp.tif " <> lfile5];
DeleteFile["tmp.tif"];
UnsetDirectory[];
];
If[iw5a,
dt = Round[100/mplots];
cline = "convert -scale 700x700 -delay " <> ToString[dt];
SetDirectory[ToFileName[{dirhtml, handle}]];
Do[
    print[2, "Plotting in progress: ", np, "/", mplots];
    cline = cline <> " " <> ToString[np] <> ".tif";
    pictfile = ToString[np] <> ".tif";
    c = Exp[-(np - 1)/mplots 2 Pi I];
plt = ContourPlot[
      Re[zmax c psixyi[x, y]],
      {x, -xyplot, xyplot}, {y, -xyplot, xyplot},
      PlotRange -> {-4, 4},
      Contours -> 3,
      ColorFunctionScaling -> False,
      ColorFunction -> (Hue[If[# > 0, 0, 0.65], If[Abs[#] < 2, 1/4, 1], 1] &),
      PlotPoints -> npts,
      ImageSize -> 433,
      Frame -> True,
      Axes -> False,
      DisplayFunction -> Identity];
    Export[pictfile, plt, ImageResolution -> 233];
    , {np, mplots}];
cline = cline <> " -loop 9999 " <> lfile5a;
Print["Running:   ", cline];
Run[cline];
tiffiles = Table[ToString[np] <> ".tif", {np, mplots}];
DeleteFile /@ tiffiles;
ResetDirectory[];
];
On[InterpolatingFunction::dmval];
];