ThreeBodyExample.java

package jat.examples.ThreeBodyExample;

import jat.core.cm.TwoBodyAPL;
import jat.core.plot.plot.FrameView;
import jat.core.plot.plot.Plot2DPanel;
import jat.core.plot.plot.PlotPanel;
import jat.core.plot.plot.plots.ScatterPlot;
import jat.coreNOSA.math.MatrixVector.data.VectorN;
import jat.examples.TwoBodyExample.TwoBodyExample;

import java.awt.Color;
import java.text.DecimalFormat;

import javax.swing.JFrame;

import org.apache.commons.lang3.ArrayUtils;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.RealVectorFormat;
import org.apache.commons.math3.ode.FirstOrderIntegrator;
import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;

/*
 * We modified the ThreeBodyProblem by emulating the two body problem, and expanding it to ThreeBody
 * Elements of two body, such as the orbit calculations, were not as heavily focused on
 * We acknowledge that some of our calculations may be a bit inaccurate because some of the math got a bit over our heads
 * 
 * @author Maegan and John
 */
public class ThreeBodyExample {

	public ThreeBodyExample() {
	}

	public static void main(String[] args) {
		// Class that solves a three body problem and plots it
		ThreeBodyExample x = new ThreeBodyExample();

		//initialize variables
		double totalEnergy = 0;
		// set the final time = one orbit period
		double tf = 1.0;

		// create a ThreeBody orbit using three masses and gravitational con
		ThreeBodyAPL sat = new ThreeBodyAPL(9.87, 0.3, 1.0, 2.0);
		//initialize VectorN elements
		double[] y = sat.randv();

		ArrayRealVector v = new ArrayRealVector(y);

		DecimalFormat df2 = new DecimalFormat("#,###,###,##0.00");
		RealVectorFormat format = new RealVectorFormat(df2);
		System.out.println(format.format(v));

		// find out the derivs of the problem
		double[] derivs = sat.derivs(tf, y);

		// set the initial time to zero
		double t0 = 0.0;
		// obtain center of mass
		VectorN cm = sat.center_of_mass(y);

		// obtain current energy
		
		totalEnergy = sat.Energy(y);
		/*
		 * Notes about the calculations:
		 * Dormand Prince is a type of Runge-Kutta diff eq. (another method of checking diff eq)
		 * Errors are calculated looking at 4th and 5th level differential equations
		 * These were most likely chosen by the original creator of the project because they take the error into account fairly accurately
		 * In relation to physics, this calculates the relative orbits for the bodies
		 * The step handler is another type of analyzing the algorithms for each body
		 * These provide fairly accurate results
		 * 
		 * In general, these equations went over our heads, but refer to the NASA site as well as 
		 * look into other resources to get more information on 3-Body calculations
		 * 
		 */
		// propagate the orbit
		FirstOrderIntegrator dp853 = new DormandPrince853Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
		dp853.addStepHandler(sat.stepHandler);
		// double[] y = new double[] { 7000.0, 0, 0, .0, 8, 0 }; // initial
		// state

		dp853.integrate(sat, 0.0, y, 8000, y); // now y contains final state at
		// tf

		Double[] objArray = sat.time.toArray(new Double[sat.time.size()]);
		double[] timeArray = ArrayUtils.toPrimitive(objArray);
		double[] xsolArray = ArrayUtils.toPrimitive(sat.xsol.toArray(new Double[sat.time.size()]));
		double[] ysolArray = ArrayUtils.toPrimitive(sat.ysol.toArray(new Double[sat.time.size()]));

		double[][] XY = new double[timeArray.length][2];

		 int a=0;
		// System.arraycopy(timeArray,0,XY[a],0,timeArray.length);
		// System.arraycopy(ysolArray,0,XY[1],0,ysolArray.length);

		for (int i = 0; i < timeArray.length; i++) {
			XY[i][0] = xsolArray[i];
			XY[i][1] = ysolArray[i];
		}

		Plot2DPanel p = new Plot2DPanel();

		// Plot2DPanel p = new Plot2DPanel(min, max, axesScales, axesLabels);

		ScatterPlot s = new ScatterPlot("orbit", Color.RED, XY);
		// LinePlot l = new LinePlot("sin", Color.RED, XY);
		// l.closed_curve = false;
		// l.draw_dot = true;
		p.addPlot(s);
		p.setLegendOrientation(PlotPanel.SOUTH);
		double plotSize = 10000.;
		double[] min = { -plotSize, -plotSize };
		double[] max = { plotSize, plotSize };

		p.setFixedBounds(min, max);
		
		// test outputs
		// output derivatives
		for (int i = 0; i< derivs.length; i++)
		{
			System.out.println(derivs[i]);
		}

		// output center of mass
		System.out.println(cm);

		// output total energy of the system
		System.out.println(totalEnergy);
		new FrameView(p).setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);

		System.out.println("end");
	}
}