(* Testing results of Segev - Heller paper *)
(* Quantum results *)

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

 (* Eigenfunction of harmonic oscillator *)
psi[n_, x_, om_] := (om/Pi)^(1/4)/(2^(n/2)Sqrt[n!])E^(-1/2 om x^2)*
   HermiteH[n, x Sqrt[om]];
 (* Supporting functions *)
cint[om0_, om1_, a_] := cint[om0, om1, a] =
   E^(-a^2/2 om0 om1/(om0 + om1))*
   om0^(1/4) om1^(1/4);
cnint[n_, om_] := cnint[n, om] =
   1/om^(n + 1/2)*Sqrt[2^(n + 1) n!];
(* Overlap integral between wavefunctions psi[0, x-a, om0] and psi[n, x, om1] *)
int[n_, om0_, om1_, a_] := int[n, om0, om1, a] =
(
   s=If[om0==om1,If[a==0,If[n==0,1,0],(a om0)^n om1^(n/2)/n!],
      nh=Floor[n/2];
      If[a==0,If[2*nh==n,2^(-n)/nh! (om1^2 - om0^2)^nh,0],
         aom=a om0 Sqrt[om1];
         sm=si=aom^n/n!;
         c=(om1^2 - om0^2)/4/aom^2;
         sm+Sum[si=si*c(n-2 i+1)(n-2 i+2)/i,{i, 1, Floor[n/2]}]
        ]
       ];
   cint[om0,om1,a] cnint[n,om0+om1]*s
);
(* Two-dimensional overlap integrals between
 wavefunctions psi[0, x-x0, omx] psi[0, x-y0, omy]
 and psi[nx, x, om=1] psi[nmax-nx, y, om=1], nx = 0, 1, ..., nmax *)
intxy = Table[int[nx, omx, 1, x0] int[nmax - nx, omy, 1, y0], {nx, 0, nmax}];
{h1, h2, norm} = Sum[intxy[[j + 1]]^2 {j + 1/2, nmax - j + 1/2, 1}, {j, 0, nmax}];
If[norm==0,
   Print["<PRE><FONT COLOR=\"#FF0000\">The final function is identically zero</FONT>"];
   Print["Nothing was computed.\n</PRE>"];
   Exit[] ];
r1=h1/(h1+h2) 100;
lni=Log[norm];
Print["<PRE>Quantum results: "<>printF[r1,6,2]<>
   "% of energy goes to x-mode.          ln I ="<>printF[lni,8,2]<>".</PRE>"];