(* Testing results of Segev - Heller paper *)
(* Distribution of energies *)
Clear["@"];
(********** INPUT starts **********)
tabdata={ (* x0, y0, omx, omy *)
   {0,0,0.02,0.18}, (* Fig. 1 *)
   {0,0,10.,2.2}, (* Fig. 2 *)
   {0,0,0.45,0.01}, (* Fig. 3 *)
   {0,0,2.,18.}, (* Fig. 4 *)
   {3,0,2.,0.1}, (* Fig. 5 *)
   {3,0,2.,10.}, (* Fig. 6 *)
   {0.4,0.2,0.7,1.5}, (* square pattern *)
   {0.1519,0.2649,0.04689,0.05555}, (* err=38644, the worst *) 
   {0.9740,0.9398,0.5707,0.5647} (* err=30, the best *)
    };
nntab={2,6,12,20,30,40,50}; 
(* Subr. for minimizing W.-function and printing *)
<<D:\sergeev\bgu\math\packages\harm.m;
<<D:\sergeev\purdue\math\packages\printing.m;
(********** INPUT ends **********)
dirtemp="d:\\temp\bilha";
sd=SetDirectory[dirtemp];
If[sd=!=dirtemp,Print["Directory ",dirtemp," does not exist. Exiting..."];
  Abort[]];
psi[n_,x_,om_]=(om/Pi)^(1/4)/(2^(n/2)Sqrt[n!])E^(-1/2 om x^2)*
    HermiteH[n,x Sqrt[om]];
mtabdata=Length[tabdata];
mnn=Length[nntab];
Do[
{x0,y0,omx,omy}=tabdata[[ntabdata]];
Print[StyleForm[If[ntabdata <= 6, StringJoin["Fig. ", ToString[ntabdata], ".   "], ""], 
    "Section"], StyleForm[StringForm[
"\!\(x\_0\) = ``,   \!\(y\_0\) = ``,   \!\(\[Omega]\_x\) = ``,   \!\(\[Omega]\_y\) = ``",
 Sequence @@ (SetPrecision[#1, 3] & ) /@ {x0, y0, omx, omy}], 
    "Section"]]; 
tab={};
Clear[int];
int[n_,om0_,om1_,a_]:=int[n,om0,om1,a]=E^(-a^2/2 om0 om1/(om0+om1))*
   om0^(1/4) om1^(1/4)/(om0+om1)^(n+1/2)*
   Sqrt[2^(n+1) n!]*
      Sum[2^(-2 i)/i! /(n-2 i)! (om1^2-om0^2)^i *
      If[n/2-i==0,1,(a^2 om0^2 om1)^(n/2-i)],{i,0,Floor[n/2]}];
Do[n=nntab[[nn]];
en=n+1;
{px, py, qx, qy, la, w}=harm`harmxy[omx, omy,x0,y0, 1., 1., en];
w=2 w;
{e1,e2}=1/2{px^2+qx^2,py^2+qy^2};
{h1,h2,norm}=Sum[
n1=j;n2=n-j;
(int[n1,omx,1,x0] int[n2,omy,1,y0])^2*
   {n1+1/2,n2+1/2,1}
   ,{j,0,n}];
{en1,en2}={h1,h2}/norm;
{wig,qua}={e1,en1}/en 100;
   tab=Append[tab,
   {n,printing`printF[px,6,2],printing`printF[py,6,2],
   printing`printF[qx,6,2],printing`printF[qy,6,2],
   printing`printF[w,6,2],
   printing`printF[qua,6,2],printing`printF[wig,6,2]}
      ];
,{nn,mnn}];
Print[TableForm[tab, TableHeadings -> 
     {None, {"n", "\!\(\(p\_x\%*\)\)", "\!\(\(p\_y\%*\)\)", "\!\(\(x\^*\)\)", 
       "\!\(\(y\^*\)\)", "-log(\!\(\(\[Rho]\^*\)\))", "\!\(R\_x\)", "\!\(\(R\_x\%*\)\)"}}, 
    TableAlignments -> Center, TableSpacing -> {1, 1}]];
  ,{ntabdata,mtabdata}];