import java.applet.*;
import java.awt.*;
import java.lang.Math;


/**********************************************************************

       Class trojan v1.0

       Written 1997 by W. Schroers / Published under the GPL 1999
       (see the file COPYING)
       Fachbereich Physik
       Bergische Universitaet und Gesamthochschule Wuppertal

       This program simulates the motion of 3 bodies under
       gravitational interaction. For iteration it uses a 5th
       order symplectic integrator which is exactly conserving
       the overall energy if the system.

       The initial configuration has been chosen to reflect the
       situation of the so-called trojan asteroids:
       These are objects in the solar system directly placed at
       the 2 stable Lagrange points of the Sun-Jupiter system.

       This program works with either JDK 1.0 and JDK 1.1

**********************************************************************/


public class trojan extends Applet implements Runnable {

    // Global constants concerning the definition&output of the system
    public static final int dim       =   2,           // we only have 2 dimensions
	                    bodies    =   3;           // and 3 bodies (planet, trojan and sun)
  
    public static final int centerx   =  75,           // center of screen region
	                    centery   =  65;
  
    public static final double cutoff =  0.0001d;      // Short distance cutoff
    public static final double indist = 50.0d;         // Initial sun<->planet separation

    public static final double a[]    = new double[6], // Constants for integrator
	                       b[]    = new double[6];
    
    public static final int SUN       = 2,
	                    TROJAN    = 1,
	                    PLANET    = 0;
                        
    /*
      Physical coordinates, velocities and forces acting on a body
    */
    private double coord[][] = new double[dim][bodies];
    private double velo[][]  = new double[dim][bodies];
    private double force[][] = new double[dim][bodies];

    /*
      Masses of the bodies
    */
    private double masses[]  = new double[bodies];
  
    /*
      Length of timestep
    */
    protected double dt;

    /*
      At which Lagrange point should we start?
    */
    int LagrangeP            = 1;

    /*
      Shall we show the trace of the trojan motion? (yes=1, no=0)
    */
    int    showtrace         = 0;

    /*
      Screen coordinates where the bodies are plotted
    */
    private int plc[][]      = new int[dim][bodies];


    /*
      Method to initialize our applet, draw the background and the Sun
    */
    public void init() {

	double mass_sum;                  // Sum of body masses
	double velnorm;                   // Length of trojan velocity vector
	double vtx = 0.0d, 
	       vty = 0.0d;                // Initial velocities of trojan

	int i;

	// We have a black sky, of course ;-)
	this.setBackground(Color.black);

	// Get the startup parameters from
	// the HTML file
	try {
	    // Iteration step length (default: 0.1d)
	    String timestep = this.getParameter("dt");
	    try {
		dt = Double.valueOf(timestep).doubleValue();
	    } catch (NumberFormatException e) { dt = 0.1d; }
	    // Initial Lagrange point (4-5)
	    // NOTE: velocity is adjusted automatically to reflect correct orbit
	    String lagrpoint = this.getParameter("lp");
	    try {
		LagrangeP = Integer.valueOf(lagrpoint).intValue();
	    } catch (NumberFormatException e) { LagrangeP = 4; }
	    if ((LagrangeP < 4) || (LagrangeP > 5)) LagrangeP = 4;
	    // Initial velocity of trojan in x direction (default: 0.0d)
	    // (Deviations from the velocity of the Lagrange point)
	    String vTrojan_x = this.getParameter("vtx");
	    try {
		vtx = Double.valueOf(vTrojan_x).doubleValue();
	    } catch (NumberFormatException e) { vtx = 0.0d; }
	    // Initial velocity of trojan in y direction (default: 0.0d)
	    String vTrojan_y = this.getParameter("vty");
	    try {
		vty = Double.valueOf(vTrojan_y).doubleValue();
	    } catch (NumberFormatException e) { vty = 0.0d; }
	    // Show the trace of the trojan motion (default: no=0)
	    String traceDisplay = this.getParameter("showtrace");
	    try {
		showtrace = Integer.valueOf(traceDisplay).intValue();
	    } catch (NumberFormatException e) { showtrace = 0; }
	} catch (NullPointerException e0) { 
	    dt  = 0.1d;
	    vtx = 0.0d;
	    vty = 0.0d; 
	}

	// Set the masses
	masses[PLANET]   =   30.0d; // the planet is comparatively light
	masses[TROJAN]   =   0.01d; // the trojan does not pull the others
	masses[SUN]      = 1000.0d; // the sun is very heavy
	mass_sum         =    0.0d;
	for (i=0;i<bodies;i++) mass_sum += masses[i];

	// Initial positions of planet and the sun
	coord[0][PLANET] =    indist;	// The planet
	coord[1][PLANET] =    0.0d;
	coord[0][SUN]    =    0.0d;	// The sun sits in the middle
	coord[1][SUN]    =    0.0d;
                            
	// Initial velocities of the bodies
	velo[0][PLANET]  =    0.0d; // The planet
	velo[1][PLANET]  = masses[SUN]/Math.sqrt(mass_sum*indist);
	velo[0][SUN]     =    0.0d; // The sun
	velo[1][SUN]     = -masses[PLANET]/Math.sqrt(mass_sum*indist);
      
	// Treat the trojan seperately based on the choice of the Lagrange point
	// We have to be very careful in adjusting initial position&velocities!
	switch (LagrangeP) {
	case 4: 
	    // The lower stable point
	    coord[0][TROJAN] = indist/2.0d;
	    coord[1][TROJAN] = Math.sqrt(3.0d)/2.0d*indist;
	    velo[0][TROJAN]  = -Math.sqrt(3.0d)/2.0d*indist;
	    velo[1][TROJAN]  = indist*(masses[SUN]/mass_sum-0.5d);
	    velnorm = Math.sqrt((masses[SUN]*masses[SUN]+masses[SUN]*masses[PLANET]+
				 masses[PLANET]*masses[PLANET])/(mass_sum*indist)/
				(velo[0][TROJAN]*velo[0][TROJAN]+
				 velo[1][TROJAN]*velo[1][TROJAN]));
	    break;
	default: 
	case 5: 
	    // The upper stable point
	    coord[0][TROJAN] = indist/2.0d;
	    coord[1][TROJAN] = -Math.sqrt(3.0d)/2.0d*indist;
	    velo[0][TROJAN]  = Math.sqrt(3.0d)/2.0d*indist;
	    velo[1][TROJAN]  = indist*(masses[SUN]/mass_sum-0.5d);
	    velnorm = Math.sqrt((masses[SUN]*masses[SUN]+masses[SUN]*masses[PLANET]+
				 masses[PLANET]*masses[PLANET])/(mass_sum*indist)/
				(velo[0][TROJAN]*velo[0][TROJAN]+
				 velo[1][TROJAN]*velo[1][TROJAN]));
	    break;
	}
	velo[0][TROJAN] *= velnorm;
	velo[1][TROJAN] *= velnorm;

	// Before the first redraw, positions must be known
	getScreencoords();

	// Set the constants for the symplectic integrator
	a[0] = 7.0d/24.0d;
	a[1] = 3.0d/4.0d;
	a[2] = -1.0d/24.0d;
	b[0] = 2.0d/3.0d;
	b[1] = -2.0d/3.0d;
	b[2] = 1.0d;    

	// Add the initial velocity to the velocity of the Lagrange point
	velo[0][TROJAN] += vtx;
	velo[1][TROJAN] += vty;
    }


    /*
      Convert the physical coordinates (the double array) to screen
      coordinates (the integer array)
    */
    public void getScreencoords() {
	plc[0][0] = (int)coord[0][0];
	plc[1][0] = (int)coord[1][0];
	plc[0][1] = (int)coord[0][1];
	plc[1][1] = (int)coord[1][1];
	plc[0][2] = (int)coord[0][2];
	plc[1][2] = (int)coord[1][2];
    }

	
    /*
      Compute the forces acting on the bodies
    */
    public void getForces() {
	double dist,dx,dy;
	int    i,j;
    
	for (i=0;i<bodies;i++) {
	    force[0][i] = 0.0d;
	    force[1][i] = 0.0d;
	    for (j=0;j<bodies;j++) if (i!=j) {
		dx = coord[0][i] - coord[0][j];
		dy = coord[1][i] - coord[1][j];
		dist = Math.sqrt(cutoff+dx*dx+dy*dy);
		dist = masses[j]/(dist*dist*dist);
		force[0][i] -= dx*dist;
		force[1][i] -= dy*dist;
	    }
	}
    }
  

    /*
      Compute the total energy of the configuration
    */
    public double energy() {
	double epot, ekin, dx, dy;
	int    i, j;
    
	epot = 0.0d; ekin = 0.0d;
	for (i=0;i<bodies;i++) {
	    ekin += 0.5d*masses[i]*(velo[0][i]*velo[0][i]+velo[1][i]*velo[1][i]);
	    for (j=0;j<bodies;j++) if (i!=j) {
		dx = coord[0][i] - coord[0][j];
		dy = coord[1][i] - coord[1][j];
		epot -= masses[j]*masses[i]/Math.sqrt(dx*dx+dy*dy);
	    }
	}
	return ekin+epot;
    }
  
      
    /*
      Iterate one time step (symplectic 3rd order integrator)
    */
    public void iterate() {
	int k,i,d;
    
	for (k=0;k<3;k++) {
	    // First do the position update
	    for (i=0;i<bodies;i++)
		for (d=0;d<dim;d++)
		    coord[d][i] += a[k]*dt*velo[d][i];
	    // Next compute the forces from the current positions
	    getForces();
	    // Finally update the velocities
	    for (i=0;i<bodies;i++)
		for (d=0;d<dim;d++)
		    velo[d][i]  += b[k]*dt*force[d][i];
	}
    }
	
	
    /*
      Draw the bodies
      by first removing the old ones, getting the new screen coordinates
      and then plotting them on the new positions
    */	
    public void paint(Graphics g) {
	// First erase the bodies in the old positions
	g.setColor(Color.black);
	g.fillOval(centerx+plc[0][PLANET],centery+plc[1][0],5,5);
	if (showtrace == 0)
	    g.fillOval(centerx+plc[0][TROJAN],centery+plc[1][1],3,3);
	g.fillOval(centerx+plc[0][SUN]-9,centery+plc[1][2]-9,18,18);
	// Then compute the new positions
	iterate();
	getScreencoords();
	// And finally display the bodies at their new locations
	g.setColor(Color.red);
	g.fillOval(centerx+plc[0][PLANET],centery+plc[1][PLANET],5,5);
	g.setColor(Color.white);
	g.fillOval(centerx+plc[0][TROJAN],centery+plc[1][TROJAN],3,3);
	// The sun get's some really nice graphics (shouldn't slow things
	// down tooo much ;^)
	g.setColor(Color.red);
	g.fillOval(centerx-4+plc[0][2],centery-4+plc[1][2],8,8);
	g.setColor(Color.yellow);
	g.fillOval(centerx-3+plc[0][2],centery-3+plc[1][2],6,6);
	g.setColor(Color.white);
	g.fillOval(centerx-1+plc[0][2],centery+plc[1][2]-1,2,2);
	g.drawLine(centerx-8+plc[0][2],centery+plc[1][2],centerx+8+plc[0][2],centery+plc[1][2]);
	g.drawLine(centerx+plc[0][2],centery+plc[1][2]-8,centerx+plc[0][2],centery+plc[1][2]+8);
	g.drawLine(centerx-6+plc[0][2],centery+plc[1][2]-6,centerx+6+plc[0][2],centery+plc[1][2]+6);
	g.drawLine(centerx+6+plc[0][2],centery+plc[1][2]-6,centerx-6+plc[0][2],centery+plc[1][2]+6);
    }


    /*
      Thread handling
    */
    private Thread gravthread = null;


    /*
      Launch the primary animation thread
    */
    public void start() {
	if (gravthread == null) {
	    gravthread = new Thread(this);
	    gravthread.start();
	}
    }


    /*
      Dispose the animation thread
    */
    public void stop() {
	if ((gravthread != null) && (gravthread.isAlive()))
	    gravthread.stop();
	gravthread = null;
    }


    /*
      The thread code
    */
    public void run() {
	while (gravthread != null) {
	    try {
		Thread.sleep(50);
	    }
	    catch(InterruptedException e){} repaint();
	}
	gravthread = null;
    }
  
    public void update(Graphics g) {
	paint(g);
    }

}