/* This applet illustrates the Bose-Einstein statistics of an ideal
   Bose gas in a harmonic trap potential. The free parameters are
   the temperature, the total density of particles in the gas and the
   diameter of the harmonic trap.

   Author: Christopher Moseley */



import java.applet.Applet;
import java.awt.*;

public class BoseStatistics extends Applet {

   double beta;		// the inverse temperature
   double zeta;		// the fugacity
   double omega;	// frequency of the harmonic trap potential
   double ntot;		// total particle density
   int levelNumber;	// number of levels to be shown
   boolean wrongValues;
   BosePanel P;

   public void init() {
	setFont(new Font("Arial(Arabic)",Font.PLAIN,14));
	P = new BosePanel();
	setLayout(new FlowLayout(FlowLayout.RIGHT));
	add(P);
   }

   public void start() {
	double T = (Double.valueOf(P.TempField.getText())).doubleValue();
	ntot = (Double.valueOf(P.DensityField.getText())).doubleValue();
	double D = (Double.valueOf(P.TrapField.getText())).doubleValue();
	if (!(T>=0 & T<=10 & ntot>=0 & ntot<=10 & D>=0 & D<=100)) {
	   wrongValues = (true);
	} else {
	   if (T<0.2) T=0.2;	// if T=0, beta will be infinity
	   // rescaling of the parameters:
	   beta = 1/T*6.5;
	   ntot = ntot*0.2;
	   omega = 1/D*5.0;
	   levelNumber = 1 + (int)(9.5/omega);
	   Fugacity();
	   wrongValues = (false);
	}
	repaint();
   }

   public void paint(Graphics g) {
	g.clearRect(0,0,760,500);
	g.setColor(Color.black);
	g.drawLine(80,450,80,20);	// energy-axis
	g.drawLine(81,450,81,20);
	g.drawLine(80,20,70,30);
	g.drawLine(81,20,91,30);
	g.drawString("Energy",20,50);
	g.drawLine(80,449,740,449);	// density-axis
	g.drawLine(80,450,740,450);
	g.drawLine(740,449,730,439);
	g.drawLine(740,450,730,460);
	g.drawString("Particle density in an energy level",400,470);
	if (wrongValues) {
	   g.setColor(Color.red);
	   g.drawString("VALUES OUT OF RANGE",400,250);
	   return;
	}
	if (beta*beta*beta>1.202/ntot) {	// condition for condensate
	   g.setColor(Color.blue);
	   g.drawString("Condensate exists",550,180);
	} else {
	   g.setColor(Color.red);
	   g.drawString("No condensate",550,180);
	}
	g.setColor(Color.blue);
	// color for ground state
	for (int k=0; k<levelNumber ;k++) {
	   double E = omega*(double)k;
	   double Boltzmann = zeta*Math.exp(-beta*E);	// Boltzmann factor
	   double nk = (k+1)*(k+2)/2*omega*omega*omega*Boltzmann/(1-Boltzmann);
	   /*	=density of particles in the states with energy E.
		The factor (k+1)*(k+2)/2 is the degeneracy. */
	   int xbar = (int)(nk*330+0.5);
	   int ybar = 430 - (int)(omega*40*k);
	   if (xbar!=0) {
		g.drawLine(85,ybar,84+xbar,ybar);
		g.drawLine(85,ybar+1,84+xbar,ybar+1);
	   }
	   g.setColor(Color.red);	// color for exited states
	}
   }

/* The following method determines the fugacity by solving the equation
   ntot = newntot(zeta)
   with given ntot, applying Newton's algorithm. newntot is the total
   particle density as a funcion of zeta. */

   void Fugacity() {
	double newntot;
	double newntotDeriv;
	zeta = ntot/(ntot+omega*omega*omega);	// initial value of zeta
	do {
	   newntot=0;
	   newntotDeriv=0;
	   int k=0;
	   double nk;
	   double nkDeriv;
	   do {
		double E = omega*(double)k;
		double Boltzmann = zeta*Math.exp(-beta*E);
		nk = (k+1)*(k+2)/2*omega*omega*omega*Boltzmann/(1-Boltzmann);
		// same nk as before
		nkDeriv = nk/zeta/(1-Boltzmann);
		newntot = newntot + nk;
		newntotDeriv = newntotDeriv + nkDeriv;
		// the derivative of newntot(zeta) with respect to zeta
		k = k + 1;
	   } while (nk>0.00001); // the infinite sum has to be cut off somewhere
	   zeta = zeta - (newntot-ntot)/newntotDeriv;
	   // an iteration of Newton's algorithm
	} while (newntot-ntot>0.001);	// tolerance for the deviation to ntot
   }

   public boolean action(Event e,Object arg) {
	if (" New ".equals(arg)) {
	   start();
	}
	return true;
   }
}

class BosePanel extends Panel {

   TextField TempField;
   TextField DensityField;
   TextField TrapField;

   BosePanel() {
	Panel P1 = new Panel();	
	P1.setLayout(new GridLayout(4,1));
	P1.add(new Label("temperature (0-10)",Label.LEFT));
	P1.add(new Label("total density (0-10)",Label.LEFT));
	P1.add(new Label("diameter of trap (0-100)",Label.LEFT));
	P1.add(new Label());
	Panel P2 = new Panel();
	P2.setLayout(new GridLayout(4,1));
	P2.add(TempField = new TextField("6"));
	P2.add(DensityField = new TextField("10"));
	P2.add(TrapField = new TextField("2"));
	P2.add(new Button(" New "));
	Panel P3 = new Panel();
	setLayout(new FlowLayout());
	add(P1);
	add(P2);
   }

   public Dimension preferedSize() {
	return new Dimension(150,100);
   }
}