(* Testing results of Segev - Heller paper *)
(* Density diagram of the wave function of the final state *)
(* Optimized for speed *)

(* Input:  x0, y0 - displacements
           omx, omy - frequencies
           nmax - quantum number of the acceptor state
Example of input line: {x0, y0, omx, omy} = {1.1519, 1.2649, 0.4689, 0.5555};nmax=10;
*)

If[!(nmax<=20),
   Print["<P><FONT COLOR=\"#0000EE\">Density diagram of the wave function in the final state is not drawn because the quantum number is too large.</FONT></P>"],
xym = 7.;
{xmin, xmax} = {-xym, xym};
{ymin, ymax} = {-xym, xym};
chopdata=10.^(-5);
(* Rough estimation *)
npoints=10;
dx = (xmax - xmin)/npoints;
dy = (ymax - ymin)/npoints;
psix = psiy = 
   Table[psi[n, xmin + nx dx, 1.], {n, 0, nmax}, {nx, 0, npoints}];
azmax = 0;
psidata = Table[
   xn = xmin + nx dx;
   yn = ymin + ny dy;
   zn = Sum[psix[[j + 1, nx + 1]] psiy[[nmax - j + 1, ny + 1]] intxy[[j + 1]]
      , {j, 0, nmax}];
   azn=Abs[zn];
   If[azn > azmax, azmax = azn ];
   {xn, yn, zn}, {nx, 0, npoints}, {ny, 0, npoints}];
psidata = Flatten[psidata, 1];
Psi0 = Interpolation[psidata];
abschop = chopdata azmax;
(* Fine estimation *)
npoints=40;
dx = (xmax - xmin)/npoints;
dy = (ymax - ymin)/npoints;
psiy = psix = 
   Table[psi[n, xmin + nx dx, 1.], {n, 0, nmax}, {nx, 0, npoints}];
zmin00 = 10.^99;
zmax00 = -zmin00;
psidata = Table[
   xn = xmin + nx dx;
   yn = ymin + ny dy;
   zn0=Psi0[xn,yn];
   zn=0;
   If[Abs[zn0]>abschop,
     zn = Sum[psix[[j + 1, nx + 1]] psiy[[nmax - j + 1, ny + 1]] intxy[[j + 1]]
      , {j, 0, nmax}] ];
   If[zn < zmin00, {xmin00, ymin00, zmin00} = {xn, yn, zn} ];
   If[zn > zmax00, {xmax00, ymax00, zmax00} = {xn, yn, zn} ];
   {xn, yn, zn}, {nx, 0, npoints}, {ny, 0, npoints}];
psidata = Flatten[psidata, 1];
Psi1 = Interpolation[psidata];
If[zmin00>=0,zmin0=0,
{zmin0, fmin} =
   FindMinimum[Psi1[x, y], {x, xmin00 - dx/100, 
      xmin00 + dx/100}, {y, ymin00 - dy/100, ymin00 + dy/100}];
{xmin0, ymin0} = {x, y} /. fmin; ];
If[zmax00<=0,zmax0=0,
{zmax0, fmax} = FindMinimum[-Psi1[x, y], {x, xmax00 - dx/100, 
   xmax00 + dx/100}, {y, ymax00 - dy/100, ymax00 + dy/100}];
{xmax0, ymax0} = {x, y} /. fmax;
zmax0 = -zmax0; ];
(*If[zmin0>=zmax0,Print["<FONT COLOR=\"#FF0000\">The final function seems to be zero</FONT>"];Exit[]];*)
zscale = Max[Abs[zmin0], Abs[zmax0]];
(* Drawing figures *)
npoints = 80/4;
zmmm = Max[Abs[zmin0], Abs[zmax0]];
cf = (z = zmin0 + #(zmax0 - zmin0);
   Hue[If[z > 0, 0, 0.65], 1/4, z1 = Abs[z/zscale]*1.1; 
   If[z1 <= 1, z1, 1]]) &;
t2 = Timing[
plt1 =ContourPlot[Psi1[x, y]/zscale, {x, xmin, xmax}, {y, ymin, ymax}, 
   PlotPoints -> npoints + 1,
   Contours -> 60,
   PlotRange -> {zmin0, zmax0}/zscale,
   AspectRatio -> 1/1,
   FrameLabel -> {x, y},
   ContourLines -> False,
   Background -> RGBColor[1, 1, 0.9],
   DisplayFunction -> Identity,
   ImageSize -> 72 {5, 5},
   ColorFunction -> cf];
   ] // First;
plt2 = ParametricPlot[{ax Cos[t + phix], ay Cos[t + phiy]}, {t, 0, 2 Pi},
   PlotRange -> {{xmin, xmax}, {ymin, ymax}},
   AspectRatio -> 1/1,
   FrameLabel -> {x, y},
   DisplayFunction -> Identity,
   ImageSize -> 72 {5, 5},
   PlotStyle -> {Thickness[0.006], RGBColor[0, 1, 0]}];
(*
plt3 = ListPlot[{{qx,qy}},
   SymbolStyle -> {RGBColor[1, 1, 0],Thickness[0.002]},
   SymbolShape -> MakeSymbol[RegularPolygon[22, 6, {0, 0}, 0, 7]],
   DisplayFunction -> Identity];
*)
plt3 = ListPlot[{{qx,qy}},
   PlotMarkers -> {"\!\(\*\nStyleBox[\"\[CircleTimes]\",\nFontColor->RGBColor[1, 1, 0]]\)", 20}];
plt=Show[{plt1, plt2, plt3},
  PlotRange -> 1.04{{xmin, xmax}, {ymin, ymax}},
  AspectRatio -> 1/1,
  FrameLabel -> {x, y},
  DisplayFunction -> Identity,
  TextStyle -> {12},
  ImageSize -> {400, 400}
  ];
(* Print[t2]; *)
If[Head[plt] === Graphics,
(* Saving picture *)
pictfile = "harmwfun.jpg";
SetDirectory[ToFileName[{"..","..","..","system","temp"}]];
Export[pictfile,plt,"JPEG"];
ResetDirectory[];
(* Reference to the picture *)
rnd=ToString[Floor[10^8 Random[]]];
Print["<P><TABLE BORDER=1 CELLSPACING=0 CELLPADDING=0><TR><TD><TABLE BORDER=\"0\" CELLPADDING=\"8\" CELLSPACING=\"0\"><TR ALIGN=\"CENTER\">"];
Print["<TD BGCOLOR=\"#FFFFEE\"><IMG SRC=\"http://" <>
   httpHost <> "/system/temp/"<>pictfile<>"?"<>rnd<>"\" WIDTH=400 HEIGHT=400 ALT=\"Density plot of the final wave function\"></TD>"];
Print["<TD ALIGN=\"LEFT\" BGCOLOR=\"#FFFFFF\"><B>Density plot of the final wave function. Yellow dot marks the <FONT COLOR=\"#DDDD00\">leakage point</FONT> (<I>x</I><SUB>*</SUB>, <I>y</I><SUB>*</SUB>)."];
Print["Green line is the <FONT COLOR=\"#00DD00\">classical trajectory</FONT> starting from the leakage point.</B>"];
If[!(nmax<=15),Print["<P><FONT SIZE=\"-2\" COLOR=\"#DD0000\">Warning. Resolution of plotting may be insufficient for accurate display of rapid oscillations of the wavefunction with so large quantum number.</FONT></P>"];];
Print["</TD></TR></TABLE></TD></TR></TABLE></P>"];
Print["<P>"];
]; (* if there is graphics *)
]; (* if nmax<=20 *)