// xydipole.r -- Calculate the x and y components of the dipole
//               matrix element, given a function that computes
//               the eigenvectors.

// This program calculates the double integrals
// int(int(psif*x*psii)dy)dx to obtain the dipole matrix
// elements for a two-dimensional problem.  The wave functions
// need not be normalized because this is taken care of in
// this program.

// Input parameters: nx,ny: number of points to sample in the
//                   x and y directions
//
//                   dx,dy: x and y grid cell sizes
//
//                   qi1,qi2: initial state quantum numbers
//
//                   qf1,qf2: final state quantum numbers
//
//                   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 = xydipole(30,30,0.2,0.2,1,1,1,2,psiho2)
// To extract the results:
//     xres = results.x
//     yres = results.y

xydipole = function(nx,ny,dx,dy,qi1,qi2,qf1,qf2,evcalc)
{
  if (nargs != 9) {
    error("Usage: xydipole(nx,ny,dx,dy,qi1,qi2,qf1,qf2,evcalc)");
  }

// Loop over x-y space, incrementing the integrals
  sumx = 0 + 0i;
  sumy = 0 + 0i;
  normi = 0;
  normf = 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
      psii = evcalc(qi1,qi2,x,y);
      psif = evcalc(qf1,qf2,x,y);

// Increment the sums
      sumx = sumx + conj(psif)*x*psii;
      sumy = sumy + conj(psif)*y*psii;
      normi = normi + real(conj(psii)*psii);
      normf = normf + real(conj(psif)*psif);
    }
  }

// Do the normalizations and return the results in a list
  norm = sqrt(normi*normf);
  xdipole = sumx/norm;
  ydipole = sumy/norm;
  return <<x = xdipole; y = ydipole>>;
}

// 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;
}