// ident.r -- Plot the probability distribution for joint square well
// eigenfunctions with two particles.  |psi(x1,x2)|^2 is plotted versus
// x1 and x2, the positions of the first and second particles.
// The first plot shows the distribution for two non-identical particles.
// The second plot shows the distribution for two identical Fermions.
// The third plot shows the distribution for two identical Bosons.

// Input parameters are two integers representing the energy levels
// occupied by the two particles.  1 = ground state, 2 = first excited
// state, etc.

ident = function(n1,n2) {

// check to see if we have the correct number of arguments
  if (nargs != 2) {
    error("Usage: ident(n1,n2)");
  }

// compute the x positions and define pi
  x = [0:30]/30;
  pi = 3.14159;

// create the infinite square well eigenfunctions
  psi1 = sin(n1*pi*x);
  psi2 = sin(n2*pi*x);

// create the joint eigenfunctions
  psi = psi1'*psi2/1.414;
  psia = (psi1'*psi2 - psi2'*psi1)/2;
  psis = (psi1'*psi2 + psi2'*psi1)/2;

// compute the probability distributions
  prob = psi.^2;
  proba = psia.^2;
  probs = psis.^2;

// set up the plots
  pgopen();
  pgclear();
  pgwindow(1,0,0,5,5);
  pgwindow(2,5,0,10,5);
  pgwindow(3,10,0,15,5);
  pgaxes(1,0,1,3,"x1",0,1,3,"x2");
  pgaxes(2,0,1,3,"x1",0,1,3,"x2");
  pgaxes(3,0,1,3,"x1",0,1,3,"x2");

// two non-identical particles
  pgconto(1,0.1,5,prob,x,x');
  pgfill(1,-.3,.3,9,prob,x,x');

// two identical Fermions
  pgconto(2,0.1,5,proba,x,x');
  pgfill(2,-.3,.3,9,proba,x,x');

// two identical Bosons
  pgconto(3,0.1,5,probs,x,x');
  pgfill(3,-.3,.3,9,probs,x,x');

// finish the plot
  pgpause();
  pgclose();

// return the probability distributions and x as a list
  return <<x = x; prob = prob; proba = proba; probs = probs>>;
}