// adwell_demo.r -- This rfile demonstrates the solution to an asymmetric
//                  double well potential.  The parameter s allows the
//                  distance between the two wells to be adjusted from
//                  exact overlap to complete separation.

// input arguments -- n -- number of grid points (200 is a good
//                         starting point)
//                    s -- separation of well centers (try values
//                         from 0 to 6)
//                    a -- degree of asymmetry of wells -- 1 is
//                         symmetric

// graphics -- the program interactively plots requested eigenvectors
//             as a function of x before returning -- pgplot is assumed

// output -- a list containing the elements --
//           x -- a column vector containing the x values
//           eigvals -- a row vector containing the eigenvalues
//           eigvecs -- a matrix containing columns with the
//                      corresponding eigenvectors
//           psq -- the square of the momentum operator

adwell_demo = function(n,s,a) {

// argument check
  if (nargs != 3) {
    return "Usage: adwell_demo(n,s,a)";
  }

// define x and the potential energy operator uop
  dx = 1/sqrt(n);
  x = [1 - n/2: n/2]'*dx;
  u = x;
  elevel = x;
  uop = zeros(n,n);
  for (i in 1:n) {
    u[i] = -6.*(exp(-(x[i] - s/2)^2)/a + a*exp(-(x[i] + s/2)^2));
    uop[i;i] = u[i];
  }

// define the momentum squared operator
  psq = zeros(n,n);

// do the diagonal terms
  for (i in 1: n) {
    psq[i;i] = 2/(dx*dx);
  }

// do the off-diagonal terms
  for (i in 1: n - 1) {
    psq[i;i + 1] = -1/(dx*dx);
    psq[i + 1;i] = -1/(dx*dx);
  }

// define the hamiltonian
  ham = 0.5*psq + uop;

// find eigenvalues and eigenvectors using rlab's symmetric solver
  results = eigs(ham);

// plot eigenvectors and potential energy function
  m = 1;
  xmax = n*dx/2;
  while (m > 0) {
      printf("Type desired eigenstate (1 - n; 0 terminates): ");
      inlist = getline("stdin");
      m = inlist.[1];
      if (m > 0 && m <= n) {
        for (i in 1:n) {
          elevel[i] = results.val[m];
        }
        pgopen();
	pgclear();
	pgwindow(1,0,0,7.5,10);
	pgwindow(2,7.5,0,15,10);
	pgaxes(1,-xmax,xmax,6,"x",-.3,.3,7,"psi");
	pgaxes(2,-xmax,xmax,6,"x",-10,2,7,"potential");
	pgline(1,5,x,real(results.vec[;m]));
	pgline(2,9,x,u);
	pgline(2,13,x,elevel);
        sprintf(title,"Eigenstate %d, energy = %f",m,results.val[m]);
	pglabel(1,-xmax,.25,title);
	sprintf(title,"s = %.1f; a = %.8f",s,a);
	pglabel(2,-xmax,1,title);
	pgpause();
	pgclose();
      }
  }

// return the results */
  return <<x = x; eigvals = results.val; eigvecs = results.vec; psq = psq>>;
}