(* Testing results of Segev - Heller paper *)
(* Distribution of the energy over two modes *)

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

Module[{om, x0, en, a, b, w},
    wharm[om_, x0_, en_] := (
        (* p = 0 *)
        a = Sqrt[x0^2/2/ en];
        w = om en (1 - a)^2;
        b = 1 - 1/om^2;
        If[a < b, Min[w, en/om(1 - a^2/b)], w]
        )
    ];
Clear[nx];
Off[NIntegrate::"ncvb"];
normcl = NIntegrate[
   E^(2(w - wharm[omx, x0, nx] - wharm[omy, y0, nmax - nx])),
   {nx, 0.001, nmax - 0.001}];
nxcl = (nmax+1) r1cl/100-1/2;
pr1cl = E^(2(w - wharm[omx, x0, nxcl] - wharm[omy, y0, nmax-nxcl]))/normcl;
pltdata = Table[{(nx + 1/2)/(nmax + 1) 100,
   intxy[[nx + 1]]^2/norm}, {nx, 0, nmax}];
ymax = 1.1 Max[Transpose[pltdata][[2]]];
ycut = ymax/300;
{xmin, xmax} = {0, 100};
mplt = Length[pltdata];
Do[{xn, yn} = pltdata[[n]]; If[yn > ycut, Break[]]; xmin = xn, {n, mplt}];
Do[{xn, yn} = pltdata[[n]]; If[yn > ycut, Break[]]; 
    xmax = xn, {n, mplt, 1, -1}];
xmean = (xmin + xmax)/2;
{xmin, xmax} = {xmin, xmax} - xmean;
{xmin, xmax} = 1.1 {xmin, xmax};
{xmin, xmax} = {xmin, xmax} + xmean;
{xmin0,xmax0} = 100{(-0.01+1/2)/(nmax+1),(nmax+1/2+0.01)/(nmax+1)};
xmin = Max[xmin,xmin0];
xmax = Min[xmax,xmax0];
  (* Plotting *)
p0 = ListPlot[pltdata, PlotRange -> {{xmin, xmax}, {0, ymax}},
   PlotStyle -> {RGBColor[0, 0, 0.3], AbsolutePointSize[5.]},
   DisplayFunction -> Identity];
pltdata1 = 
   Join[{{First[pltdata][[1]], 0}}, pltdata, {{Last[pltdata][[1]], 0}}];
p1 = Graphics[{RGBColor[0.6, 0.6, 0.9], Polygon[pltdata1]},
   PlotRange -> {{xmin, xmax}, {0, ymax}}];
p2 = Graphics[{RGBColor[0, 0, 0.2], Line[pltdata]},
   PlotRange -> {{xmin, xmax}, {0, ymax}}];
p3 = Graphics[{{Thickness[.02], RGBColor[0.3, 0, 0.6], 
      Line[{{r1, -ymax/50}, {r1, ymax}}]},
   {AbsolutePointSize[8.], RGBColor[0, 1, 0], 
      Point[{r1cl, pr1cl}]}},
   PlotRange -> {{xmin, xmax}, {0, ymax}}];
Clear[nx];
p4 = Plot[(nx = (nmax+1) r/100-1/2;
   nx = Max[nx, 0.001];
   nx = Min[nx, nmax - 0.001];
   E^(2(w - wharm[omx, x0, nx] - wharm[omy, y0, nmax-nx]))/normcl),
   {r, xmin, xmax},
   PlotRange -> {0, ymax},
   PlotStyle -> {RGBColor[0, 1, 0]},
   DisplayFunction -> Identity
];
rat=1/GoldenRatio//N;
plt = Show[p1, p0, p2, p4, p3,
    Axes -> True,
    PlotRange -> {{xmin, xmax}, {0, ymax}},
	AxesOrigin -> {xmin, -ymax/50},
   Background -> RGBColor[1, 1, 0.8],
   ImageSize -> 72 5,
   AspectRatio -> rat,
    DisplayFunction -> Identity];
If[Head[plt] === Graphics,
(* Saving picture *)
If[$MachineType === "PC",
  pictfile = "D:\\temp\\harmener.gif";
  Display[pictfile,plt,"GIF",
     ImageSize -> 220 {1,rat},
     ImageResolution -> 150],
   pictfile = $HomeDirectory<>"/.web/system/temp/harmener.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/harmener.gif?"<>rnd<>"\" WIDTH=360 HEIGHT=222 ALT=\"Probability that definite percentage of the total energy goes to x-mode\">
<BR><FONT SIZE=\"-1\">Percentage of the total energy shared by <I>x</I>-mode</FONT></TD>"];
Print["<TD ALIGN=\"LEFT\" BGCOLOR=\"#FFFFFF\"><B>Probability that definite percentage of the total energy goes to <I>x</I>-mode."];
Print["Blue diagram shows the <FONT COLOR=\"#0000FF\">quantum results</FONT>."];
Print["Vertical bar marks the everage quantum share."];
Print["Green line shows results of maximization of the <FONT COLOR=\"#00DD00\">Wigner function</FONT> in phase space."];
Print["Green dot marks absolute maximum of the Wigner function for the most favorable <I>x</I>-share of the energy.</B></TD></TR></TABLE></TD></TR></TABLE></P>"];
Print["<P>"];
];