(* Testing results of Segev - Heller paper *)
(* Phase space results (constrained minimium of the Wigner function) *)

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

(* Load a package to find constraint minimum for harmonic oscillators *)
<<harm.txt;

engap=nmax; (* Note that the energy of quantum zero vibrations (1/2+1/2) is omitted *) 
{px, py, qx, qy, la, w} = harmxy[omx, omy, x0, y0, 1., 1., engap];
If[engap!=0,
   {e1, e2} = 1/2{px^2 + qx^2, py^2 + qy^2};
   r1cl=(e1+1/2)/(nmax+1) 100;
   lnicl=-2 w - 2 Log[Pi];
   Print["<PRE> In phase space: "<>printF[r1cl,6,2]<>
      "% of energy goes to x-mode.  ln ro(q*,p*) ="<>printF[lnicl,8,2]<>"."];
  ];
{f11, f22, f33, f44} = {omx - la, omy - la, 1/omx - la, 1/omy - la};
gradh = Sqrt[qx^2 + qy^2 + px^2 + py^2];
a = (qx^2 f22 f33 f44 + qy^2 f11 f33 f44 + px^2 f11 f22 f44 + 
          py^2 f11 f22 f33)/gradh^2;
If[a>0,
   lnidelt = -2 w - 1/2 Log[Pi a] - Log[gradh];
   Print["                                              "<>
      "ln I (Delta-appr.) ="<>printF[lnidelt,8,2]<>"."],
   Print["                                              "<>
      "ln I (Delta-appr.) = <FONT COLOR=\"#00FF00\">Undefined</FONT> (alpha=0)."]
   ];
Print["   Minimum of W: (qx, px) = ("<>
   printF[qx,5,2]<>", "<>printF[px,5,2]<>")"];
Print["                 (qy, py) = ("<>
   printF[qy,5,2]<>", "<>printF[py,5,2]<>")"];
Print["</PRE>"];
ax = Sqrt[qx^2 + px^2];
phix = If[ax==0,0,ArcCos[qx/ax]];
If[px < 0, phix = -phix];
ay = Sqrt[qy^2 + py^2];
phiy = If[ay==0,0,ArcCos[qy/ay]];
If[py < 0, phiy = -phiy];