/*
 * Trojan.java by Richard J. Davies
 * from `Introductory Java for Scientists and Engineers'
 * chapter: `Physical Modelling'
 * section: `Example - Planetary Motion'
 *
 * This program simulates a three body system of the Sun, Jupiter and a
 * small asteroid. We integrate this system's motion under Newton's laws,
 * and transform our results into rotation coordinates in which the Sun
 * and Jupiter are stationary. This allows us to see that certain positions
 * of the asteroid relative to the other bodies (the `Trojan points' of
 * Jupiter) are stable.
 * 
 * We use a single 12-component vector to store the state of the system.
 * Taken in pairs, its contents are: the position of the sun, the velocity of
 * the sun, the position of Jupiter, the velocity of Jupiter, the position
 * of the asteroid, the velocity of the asteroid.
 */
import uk.co.jscieng.*;
import java.awt.*;

public class Trojan
{
  // Statistics concerning the Sun's orbit.
  public static final double SUNMASS = 1.99e30;
  
  // Statistics concerning Jupiter's orbit.
  public static final double JUPMASS = 1.90e27;
  public static final double JUPDIST = 7.78e11;
  public static final double JUPPERIOD = 3.74e8;
  public static final double JUPSPEED =
      2 * Math.PI * JUPDIST / JUPPERIOD;
  
  // The initial data for the asteroid.
  public static final double TRJMASS = 6.4e10;
  public static final double TRJANGLE = -1.0;
  public static final double TRJDIST = 8.2e11;
  public static final double TRJPERIOD = 4.2e8;
  public static final double TRJSPEED =
      2 * Math.PI * TRJDIST / TRJPERIOD;
  
  // The constant for gravitational attraction
  public static final double G = 6.67e-11;

  // The size of the steps in our integration
  public static final double STEP = 1e6;

  // The length of the trails that we leave
  // to indicate the asteroid's motion.
  public static final int HISTORY = 1000;
  
  
  // Utility method to add two vectors. Do not
  // use JNL, we cannot assume this knowledge.
  public static double[] add(double[] a, double[] b)
  {
    double[] c = new double[a.length];
    
    for (int i=0; i<a.length; i++)
      c[i] = a[i] + b[i];
    
    return c;
  }
  
  
  // Utility method to multiply a vector
  // by a scalar. Do not use JNL, we cannot assume
  // this knowledge.
  public static double[] mult(double[] a, double b)
  {
    double[] c = new double[a.length];
    
    for (int i=0; i<a.length; i++)
      c[i] = a[i] * b;
    
    return c;
  }


  // Computation of gravitational force on body a
  // from body b. Body a is at (ax,ay) and has mass
  // am. Body b is at (bx,by) and has mass bm.
  public static double[] force(double ax, double ay,
                               double am,
                               double bx, double by,
                               double bm)
  {
    // Calculate the displacement, its magnitude
    // squared and magnitude.
    double dispx = bx - ax;
    double dispy = by - ay;
    double dsqr = dispx*dispx + dispy*dispy;
    double d = Math.sqrt(dsqr);
    
    // Calculate the scalar size of the force.
    double fmag = G * am * bm / dsqr;

    // Create a vector of this size pointing
    // from body a to b.
    double[] f = { dispx * fmag / d,
                   dispy * fmag / d };
    
    return f;
  }

  
  // Compute the derivative of the current state of
  // the system.
  public static double[] deriv(double[] state)
  {
    // The force on Jupiter from the Sun.
    double[] fjs = force(state[4], state[5], JUPMASS,
                         state[0], state[1], SUNMASS);
    
    // The force on the asteroid from the Sun.
    double[] fts = force(state[8], state[9], TRJMASS,
                         state[0], state[1], SUNMASS);
    
    // The force on the asteroid from Jupiter.
    double[] ftj = force(state[8], state[9], TRJMASS,
                         state[4], state[5], JUPMASS);
    
    // The derivative of the state. Positions of objects
    // have that object's velocity as their derivative,
    // Velocities have acceleration (calculated from force)
    // as their derivative.
    double[] ans = { // Velocity of the sun.
                     state[2], state[3],
                     // Acceleration of the sun.
                     (-fjs[0] - fts[0])/SUNMASS,
                     (-fjs[1] - fts[1])/SUNMASS,
                     // Velocity of Jupiter.
                     state[6], state[7],
                     // Acceleration of Jupiter.
                     (fjs[0] - ftj[0])/JUPMASS,
                     (fjs[1] - ftj[1])/JUPMASS,
                     // Velocity of the asteroid.
                     state[10], state[11],
                     // Acceleration of the asteroid.
                     (fts[0] + ftj[0])/TRJMASS,
                     (fts[1] + ftj[1])/TRJMASS};
    return ans;
  }
  

  // A method to change coordinate systems. Given
  // the current state of the system, it returns a
  // translated and rotated state in which the Sun
  // has been moved to the origin and Jupiter
  // rotated on to the x-axis.
  public static double[] chngcoord(double[] state)
  {
    int i;
    
    // Create storage to store translated state.
    double[] tmp = new double[state.length];

    // Copy all elements across unchanged
    // This is only correct for velocities.
    for (i=0; i<state.length; i++)
      tmp[i] = state[i];
    
    // Then go back and copy positions across again
    // translating them as we do so.
    for (i=0; i<state.length; i+=4)
      tmp[i] = state[i] - state[0];
    
    for (i=1; i<state.length; i+=4)
      tmp[i] = state[i] - state[1];
    
    // Create storage for translated and rotated
    // state.
    double[] ans = new double[state.length];
    
    // Take dot products with position of Jupiter to
    // get new x-coordinates.
    for (i=0; i<tmp.length; i+=2)
      ans[i] = (tmp[i]*tmp[4] + tmp[i+1]*tmp[5])/JUPDIST;
    
    // Take dot products with an orthogonal vector to
    // get new y-coordinates.
    for (i=1; i<state.length; i+=2)
      ans[i] = (-tmp[i-1]*tmp[5] + tmp[i]*tmp[4])/JUPDIST;
    
    return ans;
  }
  
  
  // The main method actually performs the integration.
  public static void main(String[] argv)
  {
    // Set up the initial state of the system.
    double[] state = { // Position of the Sun.
                       0, 0,
                       // Velocity of the Sun.
                       0, 0,
                       // Position of Jupiter.
                       JUPDIST, 0,
                       // Velocity of Jupiter.
                       0, JUPSPEED,
                       // Position of the asteroid.
                       TRJDIST * Math.cos(TRJANGLE),
                       TRJDIST * Math.sin(TRJANGLE),
                       // Velocity of the asteroid.
                       -TRJSPEED * Math.sin(TRJANGLE),
                       TRJSPEED * Math.cos(TRJANGLE) };

    // Rotate this to nice coordinates.
    double[] statechg = chngcoord(state);
    
    // Set up the array for storing trails showing
    // motion of asteroid.
    double[][] pastpos = new double[HISTORY][4];
    int i = 0;
    
    // Show a graphics window of the right size.
    SciGraph.showGraph(-3*JUPDIST, 3*JUPDIST,
                       -3*JUPDIST, 3*JUPDIST);
    
    // Keep looping until the user stops us.
    while (true)
    {
      // Perform Runge-Kutta integration.
      double[] d = deriv(state);
      
      double[] k1 = mult(deriv(state), STEP);
      double[] k2 = mult(deriv(add(state, mult(k1, 0.5))),
                         STEP);
      double[] k3 = mult(deriv(add(state, mult(k2, 0.5))),
                         STEP);
      double[] k4 = mult(deriv(add(state, k3)), STEP);
      
      double[] newstate = add(state,
                  add(mult(k1, 1.0/6),
                      add(mult(k2, 1.0/3),
                          add(mult(k3, 1.0/3), mult(k4, 1.0/6))
                         )
                     )
                 );
      
      // Rotate the result into nice coordinates.
      double[] newstatechg = chngcoord(newstate);
      
      // Draw the position of the Sun in yellow.
      SciGraph.line(statechg[0], statechg[1],
                    newstatechg[0], newstatechg[1],
                    Color.yellow);
      
      // Draw the position of Jupiter in blue.
      SciGraph.line(statechg[4], statechg[5],
                    newstatechg[4], newstatechg[5],
                    Color.blue);
      
      // Cover the part of the asteroid's trail that
      // has expired in white.
      SciGraph.line(pastpos[i][0], pastpos[i][1],
                    pastpos[i][2], pastpos[i][3],
                    Color.white);
 
      // Record the new part of the asteroid's trail
      // in its place.
      pastpos[i][0] = statechg[8];
      pastpos[i][1] = statechg[9];
      pastpos[i][2] = newstatechg[8];
      pastpos[i][3] = newstatechg[9];

      // Draw this new part of the trail in red.
      SciGraph.line(pastpos[i][0], pastpos[i][1],
                    pastpos[i][2], pastpos[i][3],
                    Color.red);
      
      // Update the state of the system.
      state = newstate;
      statechg = newstatechg;
      i = (i+1) % HISTORY;
    }
  }
}