/* program for finding the eigenvalues and eigenfunctions  *
 * of a particle in different potentials,                  *
 * using the Numerov- and Schiess-Algorithms               *
 *                                                         *
 * @author:  Christoph Federrath                           *
 * @version: 17.01.2006                                    */

import java.awt.*;
import java.applet.*;
import java.awt.event.*;

public class Numerov extends Applet implements Runnable, ActionListener, ItemListener
{
    private int W = 500, H = 420, i = -1; //window-size, order i
    private double E = 0.1, dE = 0.1; //eigenvalue E and accuracy
    private static final double min_dE = 1.e-13;
    private static final double h = 0.001; //Numerov stepping value
    private double x, x_m1, x_p1; //independend variable x, x_minus_1, ...
    private double psi, psi_m1, psi_p1; //eigenfunction psi, ...
    private double psi_m1_initial = -h*1.e-20;
	private	double integral = 0., integral_upper_boundary = 10.; //for normalization
	private boolean begin = true;
	private String choice = "Values_1", V = "x^2";
    Thread animator;
    Canvas can, out;
    Button but_next, but_previous;
    Choice Chooser;
    private Graphics g_can, g_out, offGraphics;
	private Image offImage;

    public void init()
    {
		setLayout(null);
		can = new Canvas();
		out = new Canvas();
		but_next = new Button("Next");
		but_previous = new Button("Previous");
		can.setBackground(Color.white);
		can.setBounds(1,1,W,H);
		out.setBackground(new Color(0,200,100));
		out.setBounds(1,H,W,50);

        Panel pan = new Panel();
        pan.setLayout(null);
        pan.setBounds(W+1,1,120,H+49);
        pan.setBackground(new Color(0,200,100));
        pan.setFont(new Font("Verdana",Font.PLAIN,10));
        Choice Chooser = new Choice();
        Chooser.add("Values_1");
        Chooser.add("Values_2");
        Chooser.add("Values_3");
        Chooser.add("Values_4");
        Chooser.add("Values_5");
        Chooser.add("Values_6");
        Chooser.add("Values_7");

        Chooser.setBounds(20,30,80,30);
		but_next.setBounds(20,H/2-35,80,30);
		but_previous.setBounds(20,H/2+5,80,30);

		add(can); add(out); pan.add(Chooser);
		pan.add(but_next); pan.add(but_previous); add(pan);

		Chooser.addItemListener(this);
		but_next.addActionListener(this);
		but_previous.addActionListener(this);

		g_can = can.getGraphics();
		g_out = out.getGraphics();
		g_out.setColor(Color.blue);
		g_out.setFont(new Font("Verdana",Font.BOLD,12));
        offImage = createImage(W,H);
        offGraphics = offImage.getGraphics();
    }//init()

    public void run()
    {
		if(begin)
		{ try{Thread.sleep(100);} catch(InterruptedException exc){};
		  begin = false; }
		
	    drawBackgroundImage(offGraphics);
	    g_can.drawImage(offImage, 0, 0, null); //display Image on screen

		int step = 0; //number of Numerov steps
		integral = 0.; integral_upper_boundary = 10.;


		/* begin of the Numerov-Algorithm */
		x = -10.; psi = 0.; psi_m1 = psi_m1_initial;

		while(x <= 10.) //Numerov-Algorithm, iteration
    	{
    		x_p1 = x + h; x_m1 = x - h;
			psi_p1 = ((2.-5./6.*h*h*k(x))*psi-(1.+h*h/12.*k(x_m1))*psi_m1)
				      /(1.+h*h/12.*k(x_p1));
    		x_m1 = x; x = x_p1;
    		psi_m1 = psi; psi = psi_p1;

			if(step%(int)(1.e-2/h) == 0) { draw(g_can); } //paint not so often

			/* normalization of the wave function psi : integral(|psi|^2*dx)=1 */
            if((x > 0.) && (Math.abs(psi-psi_m1) > 10.*h) && (Math.abs(psi-psi_m1) < 50.*h))
            { integral_upper_boundary = x; }
            if(x < integral_upper_boundary-2.0)
            { integral = integral + h*psi*psi; }
            if((i > -1) && (dE < min_dE) && (integral > 10.))
            { psi_m1_initial = psi_m1_initial*0.1; run(); return; }
			step++;
    	}//while iteration

        if((i > -1) && (dE < min_dE) && (integral > 1.01))
        { psi_m1_initial = psi_m1_initial*0.94; run(); return; }
		if((i > -1) && (dE < min_dE))
		{
        	if(integral < 0.10)
            {psi_m1_initial = psi_m1_initial*100.; run(); return;}
        	if(integral < 0.99)
            {psi_m1_initial = psi_m1_initial*1.14; run(); return;}
		}
		/* normalization end */

		print(); //print values
		if(dE < min_dE)
		{
            System.out.println("order i      = " + i);
            System.out.println("int_up_bound = " + integral_upper_boundary);
            System.out.println("integral     = " + integral);
            System.out.println("psi_m1_init  = " + psi_m1_initial);
            System.out.println();
		}
    }//run()

    public double V(double x) //potential V(x) can be changed by the user
	{
        if(choice == "Values_1") //harmonic potential
        { V = "x^2"; return x*x; }
        if(choice == "Values_2") //anharmonic potential
        { V = "0.2*x^4+x^2"; return x*x + 0.2*x*x*x*x; }
        if(choice == "Values_3") //Doppelmuldenpotential
        { V = "0.02*x^4-x^2+13"; return - x*x + 0.02*x*x*x*x + 13.; }
        if(choice == "Values_4") //Morse potential 1
        { V = "10*(1-exp[-0.3*(x+5)])^2";
          return 10.*Math.pow((1.-Math.exp(-0.3*(x+5.))),2.); }
        if(choice == "Values_5") //Morse potential 2
        { V = "10*(1-exp[-0.6*(x+2)])^2";
          return 10.*Math.pow((1.-Math.exp(-0.6*(x+2.))),2.) + 1.; }
        if(choice == "Values_6") //test potential
        { V = "(x-1)^2+1"; return (x-1.)*(x-1.) + 1.; }
        if(choice == "Values_7") //Kasten potential
        { V = "Kasten";
          if(x < -3 || x > 3) return 1.e2;
          else return 0.0;
        }
		else return x*x;
	}

	public double k(double x) //function k(x) in the diff.-equ. y''(x)+k(x)*y(x) = 0
	{ return (2*E - V(x)); }//k(x)

	public void drawBackgroundImage(Graphics gr)
	{
        gr.setColor(Color.white);
		gr.fillRect(0,0,W,H);
        gr.setColor(Color.black);
		gr.drawRect(0,0,W-1,H-1);
		gr.drawLine(W/2,20,W/2,H);
		for(int z = 0; z <= 10; z++) //draw ticks
		{gr.drawLine((int)(z*W/10.),H, (int)(z*W/10.),H-5);}
		gr.setFont(new Font("Verdana",Font.BOLD,8));
		gr.drawString("-8",(int)(W/10.-4.),H-6);
		gr.drawString("+8",(int)(9.*W/10.-6.),H-6);

		gr.setFont(new Font("Verdana",Font.BOLD,10));
		gr.setColor(Color.red);
		gr.drawLine(6,15,24,15);
		gr.drawString("psi",30,20);
		gr.setColor(new Color(0,150,0));
		gr.drawLine(6,35,24,35);
		gr.drawString("|psi|^2",30,40);
		gr.setColor(Color.blue);
		gr.drawLine(6,55,24,55);
		gr.drawString("potential = " + V,30,60);

        /* draw potential */
		int px, pyV;
        for(int z = 0; z <= 3000; z++)
        {
        	px = (int)(z/3000.*W);
            pyV = H - (int)(V(2.*z/300.-10.)*H/20.);// 2....-10....../20.........
            gr.setColor(Color.blue);
            gr.drawLine(px,pyV,px,pyV);
        }
        if(choice == "Values_7")
        { gr.drawLine((int)(7./20.*W),H,(int)(7./20.*W),0);
		  gr.drawLine((int)(13./20.*W),H,(int)(13./20.*W),0);
        }

		gr.setColor(Color.black);
		gr.setFont(new Font("Verdana",Font.BOLD,12));
		if(dE > min_dE)
		{ gr.drawString("Schiess-Algorithm",W/2-58,14); }
		if((dE < min_dE) && ((integral < 0.99) || (integral > 1.01)))
		{ gr.drawString("Normalization",W/2-45,14); }
	}

	public void draw(Graphics gr)
	{
		int px, pyWF, pyWFsquare; //pixelcoords
        px = W/2 + (int)(x/10.*W/2.); /// /10. ............

        /* draw wave-function psi */
        pyWF = H - (int)((psi*5. + E)*H/20.);///.20........
        gr.setColor(Color.red);
        gr.drawLine(px,pyWF,px,pyWF);

        /* draw |wave-function|^2 */
        pyWFsquare = H - (int)((psi*psi*5. + E)*H/20.);//..20.........
        gr.setColor(new Color(0,150,0));
        gr.drawLine(px,pyWFsquare,px,pyWFsquare);
	}

	public void print()
    {
		g_out.clearRect(0,0,W,50);
		g_out.setColor(Color.black);
		g_out.drawString("order:",60,20);
		g_out.drawString("i = "+String.valueOf(i),60,40);
		g_out.drawString("energy:",W/2-60,20);
		g_out.drawString("E = "+String.valueOf(Math.round(1.e9*E)/1.e9),W/2-60,40);
		g_out.drawString("accuracy:",W/2+110,20);
		g_out.drawString("dE = "+String.valueOf((float)dE),W/2+110,40);
    }//print()

    public void actionPerformed(ActionEvent evt)
    {
		if(evt.getActionCommand()==but_next.getActionCommand())
		{
                dE = 0.1;
                E = E+dE;
                run(); //sets a new value for psi
                i++;

                if(i%2 == 0) //i == odd
                {
                    while(dE > min_dE) //Schiess-Algorithm
                    {
                        if(psi > 0)
                        {
                            E = E+dE;
                            run();
                        }
                        else
                        {
                            E = E-dE;
                            dE = dE/10;
                            E = E+dE;
                            run();
                        }
                    }//while
                }//if i == odd

                if(i%2 == 1) //i == even
                {
                    while(dE > min_dE) //Schiess-Algorithm
                    {
                        if(psi < 0)
                        {
                            E = E+dE;
                            run();
                        }
                        else
                        {
                            E = E-dE;
                            dE = dE/10;
                            E = E+dE;
                            run();
                        }
                    }//while
                }//if i == even
		}//but_next

    	if(evt.getActionCommand()==but_previous.getActionCommand())
		{
			if(i < 1) {} //only positive values for i are allowed
			else
			{
                dE = 0.1;
                E = E-dE;
                run(); //sets a new value for psi
                i--;

                if(i%2 == 0) //i == odd
                {
                    while(dE > min_dE) //Schiess-Algorithm
                    {
                        if(psi < 0)
                        {
                            E = E-dE;
                            run();
                        }
                        else
                        {
                            E = E+dE;
                            dE = dE/10;
                            E = E-dE;
                            run();
                        }
                    }//while
                }//if i == odd

                if(i%2 == 1) //i == even
                {
                    while(dE > min_dE) //Schiess-Algorithm
                    {
                        if(psi > 0)
                        {
                            E = E-dE;
                            run();
                        }
                        else
                        {
                            E = E+dE;
                            dE = dE/10;
                            E = E-dE;
                            run();
                        }
                    }//while
                }//if i == even
			}//else i > 1
    	}//but_previous
	}//action

	public void itemStateChanged(ItemEvent evt) //Choice object
	{
        if(evt.getItem() == "Values_1")
        {
	        stop();
	        choice = "Values_1";
			i = -1; E = 0.1; dE = 0.1;
			V(0); start();
        }
        if(evt.getItem() == "Values_2")
        {
	        stop();
	        choice = "Values_2";
			i = -1; E = 0.1; dE = 0.1;
			V(0); start();
        }
        if(evt.getItem() == "Values_3")
        {
	        stop();
	        choice = "Values_3";
			i = -1; E = 0.1; dE = 0.1;
			V(0); start();
        }
        if(evt.getItem() == "Values_4")
        {
	        stop();
	        choice = "Values_4";
			i = -1; E = 0.1; dE = 0.1;
			V(0); start();
        }
        if(evt.getItem() == "Values_5")
        {
	        stop();
	        choice = "Values_5";
			i = -1; E = 0.1; dE = 0.1;
			V(0); start();
        }
        if(evt.getItem() == "Values_6")
        {
	        stop();
	        choice = "Values_6";
			i = -1; E = 0.1; dE = 0.1;
			V(0); start();
        }
        if(evt.getItem() == "Values_7")
        {
	        stop();
	        choice = "Values_7";
			i = -1; E = 0.1; dE = 0.1;
			V(0); start();
        }
	} //itemStateChanged(evt)

    public void start()
    {
		animator = new Thread(this);
    	animator.start();
	}//start()

    public void stop()
    {
		animator = null;
    }//stop()
}