// dipole_demo.r -- calculate the amplitudes of the dipole moment for
//                  for a set of bound states -- the dipole moment is
//                  defined as X_if = psi_f'*xop*psi_i, where X_if
//                  is the dipole moment associated with transitions
//                  between states psi_i and psi_f and xop is the
//                  position (x) operator

// input parameters -- istate -- number of initial state
//                     maxstates -- maximum number of states to be
//                                  considered
//                     x -- column vector containing x values
//                     vecs -- matrix containing eigenvectors in
//                             columns

// obtain x and vecs from demos that compute eigenvalues and
// eigenvectors of bound states, such as ho_demo.r or anho_demo.r

// graphics -- plots the dipole moment against the final state number

// output -- returns a list containing the initial state number,
//           a column vector containing the final state numbers,
//           and a column vector containing the computed dipole
//           moments

dipole_demo = function(istate,maxstates,x,vecs) {
  if (nargs != 4) {
    return "Usage: dipole_demo(istate,maxstates,x,vecs)"
  }

// compute the x operator
  n = length(x);
  xop = zeros(n,n);
  for (i in 1:n) {
    xop[i;i] = x[i];
  }

// compute the matrix elements -- they are real for bound states
  index = [1:maxstates]';
  matrixels = index;
  psi_i = vecs[;istate];
  maxval = 0;
  for (i in 1:maxstates) {
    psi_f = vecs[;i];
    matrixels[i] = real(psi_i'*xop*psi_f);
    if (abs(matrixels[i]) > maxval) {
      maxval = abs(matrixels[i]);
    }
  }

// do plotting
  pgplot([index,matrixels],1,maxstates,maxstates,"state number",-2,2,5,"dipole moment");

// return stuff
  return <<istate = istate; index = index; matrixels = matrixels>>;
}