// solid.r -- Periodic potential to emulate a solid material consisting
//            of a number of atoms in a row

// input arguments -- n -- number of grid points (200 is a good
//                         starting point)
//                    natoms -- number of atoms (try 2-10)
//                    strength -- strength of atomic potential (try
//                    2-6)

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

solid = function(n,natoms,strength) {

// argument check
  if (nargs != 3) {
    return "Usage: solid(n,natoms,strength)";
  }

// 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);
  wavelength = n*dx/natoms;
  for (i in 1:n) {
    u[i] = -strength*sin(3.14159*(x[i] - x[1])/wavelength)^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);
  }

// make for periodic solutions
  psq[1;n] = -1/(dx*dx);
  psq[n;1] = -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",-8,1,10,"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);
	pgpause();
	pgclose();
      }
  }

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