/**
 * Schroedinger equation for the hydrogen atom:
 * bisection algorithm for finding eigenvalues
 * Feynman algorithm for the radial wavefunction
 * optional radial probability density plot.
 * The distance is given in units of a_0, the Bohr radius. Energy is
 * given in units of E_0 = 13,6 eV. Writing the radial wavefunction
 * as R(r) = r^l u(r), we get the equation:
 * d^2u/dr^2 = -Eu - 2*u/r - 2*(l+1)/r*du/dr
 * This equation is solved using Feynman method (Euler method +
 * du/dr calculated at midpoints in the mesh).
 * 
 * The HTML file should be:
 * <applet code= Hatom.class width= 561 height= 281></applet>
 * 
 * C.A. Bertulani, 09/01/2000
 **/

import java.applet.Applet;
import java.awt.*;
import java.awt.event.ActionListener;
import java.awt.event.ActionEvent;
public class Hatom extends Applet
                   implements ActionListener {
  double ymax, zmax;  // 
  int n, l;  // atomic quantum numbers
  int kk=0;   // k=1  ---> action   
  // define buttons
  String b1s="1s"; Button b1=new Button(b1s);
  String b2s="2s"; Button b2=new Button(b2s);
  String b3s="2p"; Button b3=new Button(b3s);
  String b4s="3s"; Button b4=new Button(b4s);
  String b5s="3p"; Button b5=new Button(b5s);
  String b6s="3d"; Button b6=new Button(b6s);
  String b7s="4s"; Button b7=new Button(b7s);
  String b8s="4p"; Button b8=new Button(b8s);
  String b9s="4d"; Button b9=new Button(b9s);
  String b10s="4f"; Button b10=new Button(b10s);
  // check boxes
  TextField tx1=new TextField(20);
  CheckboxGroup chg=new CheckboxGroup();
  Checkbox[] ch=new Checkbox[2];
  // text
  String Caption1, Caption2, Caption3;

  public void init() {
	// ad buttons  
    add(tx1);
    add(b1); add(b2); add(b3); add(b4); add(b5);
    add(b6); add(b7); add(b8); add(b9); add(b10);
	// add check boxes
    add(ch[0]=new Checkbox("Eigen", chg, true));
    add(ch[1]=new Checkbox("Pdens", chg, false));
	// define captions
	Caption1="Numerical solution of the hydrogen atom"; 
    Caption2="Energy, eigenfunction and probability density."
		+ "    Abcissa = distance/Bohr radius";
	// keep aware of actions
    b1.addActionListener(this);
    b2.addActionListener(this);
    b3.addActionListener(this);
    b4.addActionListener(this);
    b5.addActionListener(this);
    b6.addActionListener(this);
    b7.addActionListener(this);
    b8.addActionListener(this);
    b9.addActionListener(this);
    b10.addActionListener(this);
	// background
    setBackground(Color.white);
  }

  public void pin(Graphics g) {
    if (kk==0) {  // Screen layout
	  // Set the size and location relative to the parent's upper left corner, 
	  // in pixels. The four int arguments specify x, y, width, and height, in 
	  // that order. Use these methods to position and size a component when you 
	  // aren't using layout manager
      tx1.setBounds(425, 2, 129, 20);
      ch[0].setBounds(425, 25, 65, 20);
      ch[1].setBounds(490, 25, 65, 20);
      b1.setBounds(425, 117, 30, 20);
      b2.setBounds(425, 94, 30, 20);
      b3.setBounds(458, 94, 30, 20);
      b4.setBounds(425, 71, 30, 20);
      b5.setBounds(458, 71, 30, 20);
      b6.setBounds(491, 71, 30, 20);
      b7.setBounds(425, 48, 30, 20);
      b8.setBounds(458, 48, 30, 20);
      b9.setBounds(491, 48, 30, 20);
      b10.setBounds(524, 48, 30, 20);
      g.setFont(new Font("TimesRoman", Font.PLAIN, 14));
    }
    g.setColor(Color.black);
	// draw grid lines
    for (int i=0; i<=4; i=i+1) {
      g.drawLine(140*i, 0, 140*i, 280);
    }
    for (int i=0; i<=2;i=i+1) {
      g.drawLine(0, 140*i, 560, 140*i);
    }
  }

  public void paint(Graphics g) {    // Main program
    double a, b, e, de;
    int sa, sb, se;
    pin(g);
    g.setColor(Color.red);
	// draw caption
	g.drawString(Caption1,5,255);
	g.drawString(Caption2,5,270);
    if (kk==1) {                                       
      g.setColor(Color.blue);
      e=-1./n/n; zmax=1.05*zmax; // e= analytical value of energy
      if (ch[0].getState()==true) {
        de=.00001; a=e*1.1; b=e/1.1;     // Bisection Algorithm
        sa=fey(a, g); sb=fey(b, g);   // check nodes
        if (sa!=sb) {
          while ((b-a)>de) { 
            e=(a+b)/2; se=fey(e, g);  // bisect and recalculate 
            if (sa==se) {
              a=e;
            }
            else {
              b=e;
            }
          }
          g.setColor(Color.red);  // plot best wavefunction
          e=(a+b)/2; se=fey(e, g);
        }
      }
      else {
        se=fey(e, g);
      }
	  // Success. Got numerical energy.
      g.setColor(Color.black); tx1.setText("E = "+(float)e);
    }
    kk=0;
  }
  int fey(double e, Graphics g) {
    double r, u, v, dv, dr, z;
    int test=0;
    int x1, y1, x2, y2;
    r=0; u=1; v=-u/(l+1); dr=n*n*.002;   // Feynman Algorithm
    dv=-e*u*dr; v=v+dv/2;
    if (ch[0].getState()==true) {
      if (l==0) {
        x1=0; y1=(int)(140*(1-1/ymax));
      }
      else {
        x1=0; y1=0;
      }
      do {        // Plot function R
        r=r+dr; u=u+v*dr;
        z=u*Math.pow(r, l);
        x2=(int)(140*r/n/n); y2=(int)(140*(1-z/ymax));      
        g.drawLine(x1, y1, x2, y2); x1=x2; y1=y2;
        dv=-(e*u+2/r*(u+(l+1)*(v+dv/2)))*dr;  // Schroedinger step
        v=v+dv;
      } while ((z>-ymax)&&(z<ymax));
    }
    else {
      x1=0; y1=0;
      do {         // Plot density (rR)^2
        r=r+dr; u=u+v*dr;
        z=u*Math.pow(r, l+1);
        x2=(int)(140*r/n/n); y2=(int)(140*(1-Math.pow(z/zmax, 2))); 
        g.drawLine(x1, y1, x2, y2); x1=x2; y1=y2;
        dv=-(e*u+2/r*(u+(l+1)*(v+dv/2)))*dr;   // Schroedinger step
        v=v+dv;
      } while (r<(4*n*n));
    }
    if (u>0) {
      test=1;
    }
    else {
      test=0;
    }
    return test;
  }

  public void actionPerformed(ActionEvent e) {         // Buttons
    String tst;
    tst=e.getActionCommand();
	// get actions and initialize values
    if (b1s.equals(tst)) { n=1; l=0; kk=1; ymax=1.1; zmax=.368; }
    if (b2s.equals(tst)) { n=2; l=0; kk=1; ymax=1.1; zmax=.618; }
    if (b3s.equals(tst)) { n=2; l=1; kk=1; ymax=1; zmax=2.165; }
    if (b4s.equals(tst)) { n=3; l=0; kk=1; ymax=1.1; zmax=.828; }
    if (b5s.equals(tst)) { n=3; l=1; kk=1; ymax=1; zmax=2.64; }
    if (b6s.equals(tst)) { n=3; l=2; kk=1; ymax=5; zmax=36.3; }
    if (b7s.equals(tst)) { n=4; l=0; kk=1; ymax=1.1; zmax=1.1; }
    if (b8s.equals(tst)) { n=4; l=1; kk=1; ymax=1; zmax=3.15; }
    if (b9s.equals(tst)) { n=4; l=2; kk=1; ymax=5; zmax=36.5; }
    if (b10s.equals(tst)) { n=4; l=3; kk=1; ymax=100; zmax=1200; }
    repaint();
  }
}