// ho2perturb.r -- Solve the perturbation problem H = (p_x^2 + p_y^2 +
//                 x^2 + y^2)/2 - 0.2*exp(-x^2 - y^2).  (This is
//                 obtained by modifying xydipole.r.)

// This program calculates the double integrals
// int(int(psin*exp(-x^2 -y^2)*psin)dy)dx to obtain the perturbation
// energies for the perturbation potental plus the 2-D harmonic
// oscillator.

// Input parameters: nx,ny: number of points to sample in the
//                   x and y directions
//
//                   dx,dy: x and y grid cell sizes
//
//                   qn1,qn2: quantum numbers of the desired
//                   2-D harmonic oscillator state
//
//                   evcalc: name of function used to calculate
//                   wave functions or eigenvectors --
//                   calling sequence must be
//                   evcalc(q1,q2,x,y) where q1 and q2
//                   are the quantum numbers defining
//                   the state and x and y are the
//                   locations at which the wave function
//                   is to be evaluated

// A sample evcalc function, psiho2, is included at the end of this file.
// This sample calculates the wave function for the 2-D harmonic
// oscillator in Cartesian coordinates, given the x and y oscillator
// quantum numbers, taken to be qx = 1,2,3... and qy = 1,2,3....

// Sample calling sequence:
//     results = ho2perturb(30,30,0.2,0.2,1,2,psiho2)

ho2perturb = function(nx,ny,dx,dy,qn1,qn2,evcalc)
{
  if (nargs != 7) {
    error("Usage: ho2perturb(nx,ny,dx,dy,qn1,qn2,evcalc)");
  }

// Loop over x-y space, incrementing the integrals
  energyp = 0 + 0i;
  norm = 0;
  for (ix in 0:nx) {
    for (iy in 0:ny) {

// Compute the sampling locations for the integrals
// We assume a symmetic integration domain about the origin
      x = (ix - nx/2)*dx;
      y = (iy - ny/2)*dy;

// Compute the initial and final wave functions at each sampled point
      psin = evcalc(qn1,qn2,x,y);

// Increment the sums
      energyp = energyp + conj(psin)*(-0.2*exp(-x^2 - y^2))*psin;
      norm = norm + real(conj(psin)*psin);
    }
  }

// Do the normalization and return the results
  energyp = energyp/norm;
  return energyp;
}

// psiho2.r -- Compute the wave functions for the 2-D harmonic oscillator
psiho2 = function(qx,qy,x,y)
{

// Compute polynomial part of x wave function
  if (qx == 0) {
    xpoly = 1;
  else if (qx == 1) {
    xpoly = x;
  else if (qx == 2) {
    xpoly = 2*x^2 - 1;
  else if (qx == 3) {
    xpoly = 2*x^3 - 3*x;
  else
    error("psiho2: out of range on qx");
  }}}}

// Compute polynomial part of y wave function
  if (qy == 0) {
    ypoly = 1;
  else if (qy == 1) {
    ypoly = y;
  else if (qy == 2) {
    ypoly = 2*y^2 - 1;
  else if (qy == 3) {
    ypoly = 2*y^3 - 3*y;
  else
    error("psiho2: out of range on qy");
  }}}}

// Compute the full wave function and return the value
  wavefn = xpoly*ypoly*exp(-0.5*(x^2 + y^2));
  return wavefn;
}