(* Testing results of Segev - Heller paper *)
(* Graphs of potentials in three dimensions *)

(* 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;
*)

xym = If[nmax<=20,7.,Sqrt[2.5 engap]];
{xmin, xmax} = {-xym, xym};
{ymin, ymax} = {-xym, xym};
nholes = 12;
vi = omx^2 (x - x0)^2/2 + omy^2(y - y0)^2/2 + engap;
vf = x^2/2 + y^2/2;
npts = 3 nholes + 2;
hole[nx_, ny_] := Mod[nx, 3] == 0 && Mod[ny, 3] == 0;
Off[Plot3D::"plnc", Plot3D::"gval", Graphics3D::"nlist3"];
plt1 = Plot3D[{vf, RGBColor[0.8, 0.4, 0.4]}, {x, xmin, xmax}, {y, ymin, ymax},
   MeshStyle -> {RGBColor[0.3, 0, 0]},
   DisplayFunction -> Identity, 
   PlotPoints -> nholes];
plt2 = Plot3D[(nx = Floor[(x - xmin)/(xmax - xmin)(npts - 1) + 0.5];
      ny = Floor[(y - ymin)/(ymax - ymin)(npts - 1) + 0.5];
      If[hole[nx, ny], {Null, RGBColor[0, 0, 0]},
         {vi, RGBColor[0.7, 1, 1]}]),
   {x, xmin, xmax}, {y, ymin, ymax}, 
   MeshStyle -> {RGBColor[0, 0, 1]},
   DisplayFunction -> Identity, 
   PlotPoints -> npts];
enmaxplt = 1.2 xym^2;
plt3 = ParametricPlot3D[{(xt = ax Cos[t + phix]), (yt = ay Cos[t + phiy]),
     (vf/.{x -> xt, y -> yt}) + 0.5,
     {RGBColor[0, 1, 0], Thickness[0.01]}},
   {t, 0, 2 Pi}, 
   DisplayFunction -> Identity];
ei = engap + 1;
ef = (vf /. {x -> qx, y -> qy}) + 1;
pi = {x0, y0, ei};
pf = {qx, qy, ef};
pnti = {PointSize[0.03], RGBColor[1, 1, 0], Point[pi]};
pntf = {PointSize[0.03], RGBColor[1, 1, 0], Point[pf]};
lineif = {Thickness[0.02], RGBColor[1, 1, 0], Line[{pi, pf}]};
plt4 = Graphics3D[{pnti, pntf, lineif}];
plt = Show[plt1, plt2, plt3, plt4,
   AxesLabel -> {"x", "y", "V"},
   Background -> RGBColor[1, 1, 0.8],
   ImageSize -> 72 5,
   AspectRatio -> 292/360,
   DisplayFunction -> Identity,
   PlotRange -> {0, enmaxplt}];
If[Head[plt] === Graphics3D,
(* Saving picture *)
If[$MachineType === "PC",
  pictfile = "D:\\temp\\harmpots.gif";
  Display[pictfile,plt,"GIF",
     ImageSize -> 220 {1,rat},
     ImageResolution -> 150],
   pictfile = $HomeDirectory<>"/.web/system/temp/harmpots.gif";
   Display[pictfile,plt,"GIF"];
   Run["chmod 755 "<>pictfile]
  ];
(* 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=\"#FFFFCD\"><IMG SRC=\"../system/temp/harmpots.gif?"<>rnd<>"\" WIDTH=360 HEIGHT=292 ALT=\"Donor and acceptor potentials\"></TD>"];
Print["<TD ALIGN=\"LEFT\" BGCOLOR=\"#FFFFFF\"><B>Plots of the <FONT COLOR=\"#0000FF\">donor potential</FONT> (blue porous surface) and the <FONT COLOR=\"#AA0000\">acceptor potential</FONT> (brown solid surface) as functions of two coordinates, <I>x</I> and <I>y</I>."];
Print["Yellow line (sometimes short or hidden) marks the <FONT COLOR=\"#DDDD00\">jumping path</FONT>."];
Print["Green line on the acceptor surface is the <FONT COLOR=\"#00DD00\">classical trajectory</FONT> starting from the leakage point.</B></TD></TR></TABLE></TD></TR></TABLE></P>"];
Print["<P>"];
];