(* Input:  nqx,nqy - Oscillator quantum numbers
           omx,omy - Oscillator frequencies
           morder - Order of perturbation theory
           ndigits - Working precision
Example of input line: nqx=9;nqy=1;omx=1;omy=11/10;morder=12;ndigits=64;
*)
(* Testing that nqx, nqy lie in a legitimate range *)
If[IntegerQ[nqx]&&nqx>=0&&IntegerQ[nqy]&&nqy>=0,
Print["Oscillator quantum numbers: nx = ",nqx,",  ny = ",nqy],
Print["Oscillator quantum numbers must be non-negative integers!"];
Exit[]];
(* Testing that omx, omy lie in a legitimate range *)
If[Positive[omx]&&Positive[omy],
Print["Oscillator frequences: wx = ",omx//InputForm,",  wy = ",omy//InputForm],
Print["Oscillator frequences must be positive real numbers!"];
Exit[]];
(* Testing that morder lies in a legitimate range *)
If[IntegerQ[morder]&&morder>0,Print["Order of perturbation theory: N = ",morder],
Print["Order of perturbation theory must be a positive integer!"];
Exit[]];
(* Testing that ndigits lies in a legitimate range *)
If[IntegerQ[ndigits]&&morder>=0,Print["Working precision: ",ndigits,"\n"],
Print["Working precision must be a non-negative integer!"];
Exit[]];
If[ndigits==0,ndigits=Infinity];
omx=SetPrecision[omx//Rationalize,ndigits];
omy=SetPrecision[omy//Rationalize,ndigits];
ncoef=morder+1;
Off[Power::infy];
Off[Infinity::indet];
Print["Printing PT coefficients for the Barbanis potential wx^2 x^2/2 + wy^2 y^2/2 + g^(1/2) x y^2 ...\n"];
energy[0]=(nqx+1/2)*omx+(nqy+1/2)*omy;
Print["E[0] = ",energy[0]//InputForm," (Unperturbed harmonic-oscillator energy)"];
   "(* PT coefficients for Barbanis potential *)";
mdouble=morder*2;
mbasisx=nqx+1*morder;
mbasisy=nqy+2*morder;
kax=nqx;
kay=nqy;
kayh=Floor[kay/2];
kymin=0;
If[kayh*2!=kay,kymin=1];
enxy=Table[0,{m,1,morder}];
xn=Table[Null,{kx,0,mbasisx},{k,-1,1,2}];
yn=Table[Null,{ky,kymin,mbasisy,2},{k,-2,2,2}];
kxm=0;
If[Floor[mbasisx/2]*2!=mbasisx,kxm=1];
psi=Table[Null,{n,0,mdouble},{kx,kxm,mbasisx,2},{ky,kymin,mbasisy,2}];
(* Calculating xn[i,(-1,1)]=<i|x|(i-1,i+1> and yn[i,(-2,0,2)]=<i|y^2|(i-2,i,i+2)> *)
Do[i1=i+1;
   xn[[i1,1]]=1;xn[[i1,2]]=(i+1)/(2*omx),{i,0,mbasisx}];
Do[ih=Floor[i/2]+1;
   yn[[ih,1]]=1;yn[[ih,2]]=(2*i+1)/(2*omy);
   yn[[ih,3]]=(i+1)*(i+2)/(4*omy^2),{i,kymin,mbasisy,2}];
(*     Harmonic oscillator wave function (zero order) *)
      psi[[1,Floor[nqx/2]+1,Floor[nqy/2]+1]]=1;
(*     Successive calculation of the expansion coefficients *)
kbxmin0=kax;
kbxmax0=kax;
kbymin0=kay;
kbymax0=kay;
Do[n9=n-1;
   nh=Floor[n/2];
   kanx=kax+n;
   kanxh=Floor[kanx/2];
   kxmin=0;
   If[kanxh*2!=kanx,kxmin=1];
(*     Calculating of wavefunctions at n-th step *)
   If[nh*2==n,enxy[[nh]]=0];
   mdx=Min[n,mdouble-n];
   kbxmin=Max[kxmin,kax-mdx];
   kbxmax=kax+mdx;
   mdy=2*mdx;
   kbymin=Max[kymin,kay-mdy];
   kbymax=kay+mdy;
   Do[kbxh=Floor[kbx/2]+1;Do[a=0;kbyh=Floor[kby/2]+1;
ncmin=Max[Abs[kbx-kax],Ceiling[Abs[kby-kay]/2]];
nbcmax=n-ncmin;
nbcmaxh=Floor[nbcmax/2];
If[nbcmaxh*2!=nbcmax,ncmin=ncmin+1];
If[ncmin<=n9,
Do[nbc=n-nc;
   nbch=nbc/2;
   a=a+psi[[nc+1,kbxh,kbyh]]*enxy[[nbch]],{nc,ncmin,n9,2}]];
kcxmin=Max[kbxmin0,kbx-1];
kcxmax=Min[kbxmax0,kbx+1];
kcymin=Max[kbymin0,kby-2];
kcymax=Min[kbymax0,kby+2];
Do[kcxh=Floor[kcx/2]+1;Do[kcyh=Floor[kcy/2]+1;
         a=a-psi[[n,kcxh,kcyh]]*xn[[kcx+1,(kbx-kcx+1)/2+1]]*yn[[kcyh,(kby-kcy)/2+2]],
         {kcy,kcymin,kcymax,2}],
      {kcx,kcxmin,kcxmax,2}];
If[kbx!=kax||kby!=kay,
        a=a/((kbx-kax)*omx+(kby-kay)*omy)];
psi[[n+1,kbxh,kbyh]]=a,{kby,kbymin,kbymax,2}],
      {kbx,kbxmin,kbxmax,2}];
   If[nh*2==n,
   enxy[[nh]]=-psi[[n+1,Floor[kax/2]+1,Floor[kay/2]+1]];
   energy[nh]=enxy[[nh]];
   Print["E[",nh,"] = ",energy[nh]//InputForm];
   If[enxy[[nh]]==Indeterminate,
      Print[];Print["Fermi resonance was encountered or insufficient precision. Higher order coefficients are indeterminate!"];Exit[]];
   psi[[n+1,Floor[kax/2]+1,Floor[kay/2]+1]]=0];
kbxmin0=kbxmin;
kbxmax0=kbxmax;
kbymin0=kbymin;
kbymax0=kbymax,{n,1,mdouble}];