From 38236f37a560e993ec4eb88c510d2a92fe1d793f Mon Sep 17 00:00:00 2001 From: mxd12001 Date: Fri, 4 Dec 2015 18:40:20 -0500 Subject: [PATCH 1/7] This will emulate the jat.examplesTwoBodyExample. We are now in the process of changing our derivatives. We will be looking into the DormandPrince853Integrator, FirstOrderIntegrator, and FirstOrderDifferentialEquations --- .../ThreeBodyExample/ThreeBodyAPL.java | 12 ++++++++-- .../ThreeBodyExample/ThreeBodyExample.java | 22 +++++++++---------- 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java index 28104e0..f3265db 100644 --- a/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java +++ b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java @@ -1,14 +1,22 @@ package jat.examples.ThreeBodyExample; +import java.util.ArrayList; + import org.apache.commons.math3.ode.FirstOrderDifferentialEquations; import jat.coreNOSA.cm.ThreeBody; import jat.coreNOSA.math.MatrixVector.data.VectorN; public class ThreeBodyAPL extends ThreeBody implements FirstOrderDifferentialEquations{ + double initial_ta; + public ArrayList time = new ArrayList(); + public ArrayList xsol = new ArrayList(); + public ArrayList ysol = new ArrayList(); + public ArrayList zsol = new ArrayList(); - public ThreeBodyAPL(double G, double m1, double m2, double m3) { - super(G, m1, m2, m3); + public ThreeBodyAPL(double mu, VectorN r, VectorN v) { + super(mu, r, v); + initial_ta = ta; } @Override diff --git a/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java b/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java index 485a326..339acb3 100644 --- a/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java +++ b/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java @@ -5,8 +5,6 @@ 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; @@ -20,29 +18,30 @@ import org.apache.commons.math3.ode.FirstOrderIntegrator; import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator; public class ThreeBodyExample { + public ThreeBodyExample() { } - + public static void main(String[] args) { - // TODO Auto-generated method stub + ThreeBodyExample x = new ThreeBodyExample(); - // create a ThreeBody orbit using three masses and gravitational con - ThreeBodyAPL sat = new ThreeBodyAPL(10, 0.3, 1.0, 2.0); - //initialize VectorN elements + // create a TwoBody orbit using orbit elements + ThreeBodyAPL sat = new ThreeBodyAPL(7000.0, 0.3, 0.0, 0.0, 0.0, 0.0); + 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(2.0, y); + // find out the period of the orbit + double period = sat.period(); // set the final time = one orbit period - double tf = 2.0; + double tf = period; // set the initial time to zero double t0 = 0.0; @@ -92,5 +91,4 @@ public class ThreeBodyExample { System.out.println("end"); } - } From 9428b969723ad2acde3f53f9d058d5b2d47ac0da Mon Sep 17 00:00:00 2001 From: mxd12001 Date: Fri, 4 Dec 2015 19:12:58 -0500 Subject: [PATCH 2/7] Edits to ThreeBodyAPL. Still needs to be updated with better calculations for the relative gravity and motion --- src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java index f3265db..28104e0 100644 --- a/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java +++ b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java @@ -1,22 +1,14 @@ package jat.examples.ThreeBodyExample; -import java.util.ArrayList; - import org.apache.commons.math3.ode.FirstOrderDifferentialEquations; import jat.coreNOSA.cm.ThreeBody; import jat.coreNOSA.math.MatrixVector.data.VectorN; public class ThreeBodyAPL extends ThreeBody implements FirstOrderDifferentialEquations{ - double initial_ta; - public ArrayList time = new ArrayList(); - public ArrayList xsol = new ArrayList(); - public ArrayList ysol = new ArrayList(); - public ArrayList zsol = new ArrayList(); - public ThreeBodyAPL(double mu, VectorN r, VectorN v) { - super(mu, r, v); - initial_ta = ta; + public ThreeBodyAPL(double G, double m1, double m2, double m3) { + super(G, m1, m2, m3); } @Override From 03fe0c5a395d5879bd16776fc4fe38734c8e6611 Mon Sep 17 00:00:00 2001 From: mxd12001 Date: Fri, 4 Dec 2015 19:38:52 -0500 Subject: [PATCH 3/7] Our merge failed. I have no clue what I am commiting at this point but I need to to try pulling from John's brach and rewriting what I did horribly wrong --- src/jat/core/cm/ThreeBodyAPL.java | 274 ++++++++++++++++++ .../ThreeBodyExample/ThreeBodyAPL.java | 47 +++ .../ThreeBodyExample/ThreeBodyExample.java | 94 ++++++ 3 files changed, 415 insertions(+) create mode 100644 src/jat/core/cm/ThreeBodyAPL.java create mode 100644 src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java create mode 100644 src/jat/examples/ThreeBodyExample/ThreeBodyExample.java diff --git a/src/jat/core/cm/ThreeBodyAPL.java b/src/jat/core/cm/ThreeBodyAPL.java new file mode 100644 index 0000000..18f186f --- /dev/null +++ b/src/jat/core/cm/ThreeBodyAPL.java @@ -0,0 +1,274 @@ +/* JAT: Java Astrodynamics Toolkit + * + Copyright 2012 Tobias Berthold + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + */ + +package jat.core.cm; + +import jat.coreNOSA.algorithm.integrators.Printable; +import jat.coreNOSA.cm.Constants; +import jat.coreNOSA.cm.TwoBody; +import jat.coreNOSA.math.MatrixVector.data.Matrix; +import jat.coreNOSA.math.MatrixVector.data.VectorN; + +import java.util.ArrayList; + +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.commons.math3.ode.FirstOrderDifferentialEquations; +import org.apache.commons.math3.ode.sampling.StepHandler; +import org.apache.commons.math3.ode.sampling.StepInterpolator; + +public class ThreeBodyAPL extends ThreeBody implements FirstOrderDifferentialEquations { + double initial_ta; + public ArrayList time = new ArrayList(); + public ArrayList xsol = new ArrayList(); + public ArrayList ysol = new ArrayList(); + public ArrayList zsol = new ArrayList(); + + public ThreeBodyAPL(double mu, VectorN r, VectorN v) { + super(mu, r, v); + initial_ta = ta; + } + + public ThreeBodyAPL(double a, double e, double i, double raan, double w, double ta) { + super(a, e, i, raan, w, ta); + } + + @Override + public void computeDerivatives(double t, double[] y, double[] yDot) { + + Vector3D r = new Vector3D(y[0], y[1], y[2]); + double rnorm = r.getNorm(); + double r3 = rnorm * rnorm * rnorm; + double k = -1. * this.mu / r3; + yDot[0] = y[3]; + yDot[1] = y[4]; + yDot[2] = y[5]; + yDot[3] = k * y[0]; + yDot[4] = k * y[1]; + yDot[5] = k * y[2]; + } + + public StepHandler stepHandler = new StepHandler() { + public void init(double t0, double[] y0, double t) { + } + + public void handleStep(StepInterpolator interpolator, boolean isLast) { + double t = interpolator.getCurrentTime(); + double[] y = interpolator.getInterpolatedState(); + // System.out.println(t + " " + y[0] + " " + y[1]+ " " + y[2]); + time.add(t); + xsol.add(y[0]); + ysol.add(y[1]); + zsol.add(y[2]); + } + }; + + @Override + public int getDimension() { + return 6; + } + + public VectorN position(double t) { + + double[] temp = new double[6]; + + // Determine step size + double n = this.meanMotion(); + + // determine initial E and M + double sqrome2 = Math.sqrt(1.0 - this.e * this.e); + double cta = Math.cos(this.ta); + double sta = Math.sin(this.ta); + double sine0 = (sqrome2 * sta) / (1.0 + this.e * cta); + double cose0 = (this.e + cta) / (1.0 + this.e * cta); + double e0 = Math.atan2(sine0, cose0); + + double ma = e0 - this.e * Math.sin(e0); + + ma = ma + n * t; + double ea = solveKepler(ma, this.e); + + double sinE = Math.sin(ea); + double cosE = Math.cos(ea); + double den = 1.0 - this.e * cosE; + + double sinv = (sqrome2 * sinE) / den; + double cosv = (cosE - this.e) / den; + + this.ta = Math.atan2(sinv, cosv); + if (this.ta < 0.0) { + this.ta = this.ta + 2.0 * Constants.pi; + } + + temp = this.randv(); + this.rv = new VectorN(temp); + + // Reset everything to before + this.ta = initial_ta; + + VectorN out = new VectorN(3); + out.x[0] = temp[0]; + out.x[1] = temp[1]; + out.x[2] = temp[2]; + // out.print("sat pos at t"); + return out; + + } + + // propagates from whatever ta is starting from time t=0 relative to the + // starting point + public void propagate(double t0, double tf, Printable pr, boolean print_switch, double steps) { + double[] temp = new double[6]; + // double ta_save = this.ta; + this.steps = steps; + + // Determine step size + double n = this.meanMotion(); + double period = this.period(); + double dt = period / steps; + if ((t0 + dt) > tf) // check to see if we're going past tf + { + dt = tf - t0; + } + + // determine initial E and M + double sqrome2 = Math.sqrt(1.0 - this.e * this.e); + double cta = Math.cos(this.ta); + double sta = Math.sin(this.ta); + double sine0 = (sqrome2 * sta) / (1.0 + this.e * cta); + double cose0 = (this.e + cta) / (1.0 + this.e * cta); + double e0 = Math.atan2(sine0, cose0); + + double ma = e0 - this.e * Math.sin(e0); + + // initialize t + + double t = t0; + + if (print_switch) { + temp = this.randv(); + pr.print(t, temp); + } + + while (t < tf) { + ma = ma + n * dt; + double ea = solveKepler(ma, this.e); + + double sinE = Math.sin(ea); + double cosE = Math.cos(ea); + double den = 1.0 - this.e * cosE; + + double sinv = (sqrome2 * sinE) / den; + double cosv = (cosE - this.e) / den; + + this.ta = Math.atan2(sinv, cosv); + if (this.ta < 0.0) { + this.ta = this.ta + 2.0 * Constants.pi; + } + + t = t + dt; + + temp = this.randv(); + this.rv = new VectorN(temp); + + if (print_switch) { + pr.print(t, temp); + } + + if ((t + dt) > tf) { + dt = tf - t; + } + + } + // Reset everything to before + this.ta = initial_ta; + + } + + public double[] randv(double ta) { + double p = a * (1.0 - e * e); + double cta = Math.cos(ta); + double sta = Math.sin(ta); + double opecta = 1.0 + e * cta; + double sqmuop = Math.sqrt(this.mu / p); + + VectorN xpqw = new VectorN(6); + xpqw.x[0] = p * cta / opecta; + xpqw.x[1] = p * sta / opecta; + xpqw.x[2] = 0.0; + xpqw.x[3] = -sqmuop * sta; + xpqw.x[4] = sqmuop * (e + cta); + xpqw.x[5] = 0.0; + + Matrix cmat = PQW2ECI(); + + VectorN rpqw = new VectorN(xpqw.x[0], xpqw.x[1], xpqw.x[2]); + VectorN vpqw = new VectorN(xpqw.x[3], xpqw.x[4], xpqw.x[5]); + + VectorN rijk = cmat.times(rpqw); + VectorN vijk = cmat.times(vpqw); + + double[] out = new double[6]; + + for (int i = 0; i < 3; i++) { + out[i] = rijk.x[i]; + out[i + 3] = vijk.x[i]; + } + + return out; + } + + public double eccentricAnomaly(double ta) { + double cta = Math.cos(ta); + double e0 = Math.acos((e + cta) / (1.0 + e * cta)); + return e0; + } + + public double meanAnomaly(double t) { + + return 2. * Math.PI * t / period(); + } + + public double t_from_ta() { + + double M = meanAnomaly(); + double P = period(); + + return P * M / 2. / Math.PI; + + } + + public double ta_from_t(double t) { + + double M = meanAnomaly(t); + double ea = solveKepler(M, this.e); + + double sinE = Math.sin(ea); + double cosE = Math.cos(ea); + double den = 1.0 - this.e * cosE; + double sqrome2 = Math.sqrt(1.0 - this.e * this.e); + double sinv = (sqrome2 * sinE) / den; + double cosv = (cosE - this.e) / den; + + double ta = Math.atan2(sinv, cosv); + if (this.ta < 0.0) { + this.ta = this.ta + 2.0 * Constants.pi; + } + + return ta; + } + +} diff --git a/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java new file mode 100644 index 0000000..41a565c --- /dev/null +++ b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java @@ -0,0 +1,47 @@ +package jat.examples.ThreeBodyExample; + +import org.apache.commons.math3.ode.FirstOrderDifferentialEquations; + +import jat.coreNOSA.cm.ThreeBody; +import jat.coreNOSA.math.MatrixVector.data.VectorN; + +public class ThreeBodyAPL extends ThreeBody implements FirstOrderDifferentialEquations{ + + public ThreeBodyAPL(double G, double m1, double m2, double m3) { + super(G, m1, m2, m3); + } + + @Override + public void computeDerivatives(double arg0, double[] arg1, double[] arg2) { + // TODO Auto-generated method stub + + } + + @Override + public int getDimension() { + // returns a dimension of 6 + return 6; + } + public double[] randv() { + double[] randv = new double[18]; + randv[0] = 2.0; + randv[1] = 6.0; + randv[2] = 7.0; + randv[3] = 4.0; + randv[4] = 8.0; + randv[5] = 7.0; + randv[6] = 5.0; + randv[7] = 8.0; + randv[8] = 3.0; + randv[9] = 6.0; + randv[10] = 8.0; + randv[11] = 9.0; + randv[12] = 0.0; + randv[13] = 7.0; + randv[14] = 4.0; + randv[15] = 5.0; + randv[16] = 7.0; + randv[17] = 8.0; + return randv; + } +} diff --git a/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java b/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java new file mode 100644 index 0000000..a77a8db --- /dev/null +++ b/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java @@ -0,0 +1,94 @@ +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 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; + +public class ThreeBodyExample { + + public ThreeBodyExample() { + } + + public static void main(String[] args) { + + ThreeBodyExample x = new ThreeBodyExample(); + + // create a TwoBody orbit using orbit elements + TwoBodyAPL sat = new TwoBodyAPL(7000.0, 0.3, 0.0, 0.0, 0.0, 0.0); + + 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 period of the orbit + double period = sat.period(); + + // set the final time = one orbit period + double tf = period; + + // set the initial time to zero + double t0 = 0.0; + + // 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); + + new FrameView(p).setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); + + System.out.println("end"); + + } +} From 4bdbee0e6b3fb1cf666e18aafa42a7a184b6ba73 Mon Sep 17 00:00:00 2001 From: mxd12001 Date: Fri, 4 Dec 2015 19:41:32 -0500 Subject: [PATCH 4/7] test commit to fix branch merge problem --- src/jat/core/cm/ThreeBodyAPL.java | 274 ------------------------------ 1 file changed, 274 deletions(-) delete mode 100644 src/jat/core/cm/ThreeBodyAPL.java diff --git a/src/jat/core/cm/ThreeBodyAPL.java b/src/jat/core/cm/ThreeBodyAPL.java deleted file mode 100644 index 18f186f..0000000 --- a/src/jat/core/cm/ThreeBodyAPL.java +++ /dev/null @@ -1,274 +0,0 @@ -/* JAT: Java Astrodynamics Toolkit - * - Copyright 2012 Tobias Berthold - - Licensed under the Apache License, Version 2.0 (the "License"); - you may not use this file except in compliance with the License. - You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - - Unless required by applicable law or agreed to in writing, software - distributed under the License is distributed on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - See the License for the specific language governing permissions and - limitations under the License. - */ - -package jat.core.cm; - -import jat.coreNOSA.algorithm.integrators.Printable; -import jat.coreNOSA.cm.Constants; -import jat.coreNOSA.cm.TwoBody; -import jat.coreNOSA.math.MatrixVector.data.Matrix; -import jat.coreNOSA.math.MatrixVector.data.VectorN; - -import java.util.ArrayList; - -import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; -import org.apache.commons.math3.ode.FirstOrderDifferentialEquations; -import org.apache.commons.math3.ode.sampling.StepHandler; -import org.apache.commons.math3.ode.sampling.StepInterpolator; - -public class ThreeBodyAPL extends ThreeBody implements FirstOrderDifferentialEquations { - double initial_ta; - public ArrayList time = new ArrayList(); - public ArrayList xsol = new ArrayList(); - public ArrayList ysol = new ArrayList(); - public ArrayList zsol = new ArrayList(); - - public ThreeBodyAPL(double mu, VectorN r, VectorN v) { - super(mu, r, v); - initial_ta = ta; - } - - public ThreeBodyAPL(double a, double e, double i, double raan, double w, double ta) { - super(a, e, i, raan, w, ta); - } - - @Override - public void computeDerivatives(double t, double[] y, double[] yDot) { - - Vector3D r = new Vector3D(y[0], y[1], y[2]); - double rnorm = r.getNorm(); - double r3 = rnorm * rnorm * rnorm; - double k = -1. * this.mu / r3; - yDot[0] = y[3]; - yDot[1] = y[4]; - yDot[2] = y[5]; - yDot[3] = k * y[0]; - yDot[4] = k * y[1]; - yDot[5] = k * y[2]; - } - - public StepHandler stepHandler = new StepHandler() { - public void init(double t0, double[] y0, double t) { - } - - public void handleStep(StepInterpolator interpolator, boolean isLast) { - double t = interpolator.getCurrentTime(); - double[] y = interpolator.getInterpolatedState(); - // System.out.println(t + " " + y[0] + " " + y[1]+ " " + y[2]); - time.add(t); - xsol.add(y[0]); - ysol.add(y[1]); - zsol.add(y[2]); - } - }; - - @Override - public int getDimension() { - return 6; - } - - public VectorN position(double t) { - - double[] temp = new double[6]; - - // Determine step size - double n = this.meanMotion(); - - // determine initial E and M - double sqrome2 = Math.sqrt(1.0 - this.e * this.e); - double cta = Math.cos(this.ta); - double sta = Math.sin(this.ta); - double sine0 = (sqrome2 * sta) / (1.0 + this.e * cta); - double cose0 = (this.e + cta) / (1.0 + this.e * cta); - double e0 = Math.atan2(sine0, cose0); - - double ma = e0 - this.e * Math.sin(e0); - - ma = ma + n * t; - double ea = solveKepler(ma, this.e); - - double sinE = Math.sin(ea); - double cosE = Math.cos(ea); - double den = 1.0 - this.e * cosE; - - double sinv = (sqrome2 * sinE) / den; - double cosv = (cosE - this.e) / den; - - this.ta = Math.atan2(sinv, cosv); - if (this.ta < 0.0) { - this.ta = this.ta + 2.0 * Constants.pi; - } - - temp = this.randv(); - this.rv = new VectorN(temp); - - // Reset everything to before - this.ta = initial_ta; - - VectorN out = new VectorN(3); - out.x[0] = temp[0]; - out.x[1] = temp[1]; - out.x[2] = temp[2]; - // out.print("sat pos at t"); - return out; - - } - - // propagates from whatever ta is starting from time t=0 relative to the - // starting point - public void propagate(double t0, double tf, Printable pr, boolean print_switch, double steps) { - double[] temp = new double[6]; - // double ta_save = this.ta; - this.steps = steps; - - // Determine step size - double n = this.meanMotion(); - double period = this.period(); - double dt = period / steps; - if ((t0 + dt) > tf) // check to see if we're going past tf - { - dt = tf - t0; - } - - // determine initial E and M - double sqrome2 = Math.sqrt(1.0 - this.e * this.e); - double cta = Math.cos(this.ta); - double sta = Math.sin(this.ta); - double sine0 = (sqrome2 * sta) / (1.0 + this.e * cta); - double cose0 = (this.e + cta) / (1.0 + this.e * cta); - double e0 = Math.atan2(sine0, cose0); - - double ma = e0 - this.e * Math.sin(e0); - - // initialize t - - double t = t0; - - if (print_switch) { - temp = this.randv(); - pr.print(t, temp); - } - - while (t < tf) { - ma = ma + n * dt; - double ea = solveKepler(ma, this.e); - - double sinE = Math.sin(ea); - double cosE = Math.cos(ea); - double den = 1.0 - this.e * cosE; - - double sinv = (sqrome2 * sinE) / den; - double cosv = (cosE - this.e) / den; - - this.ta = Math.atan2(sinv, cosv); - if (this.ta < 0.0) { - this.ta = this.ta + 2.0 * Constants.pi; - } - - t = t + dt; - - temp = this.randv(); - this.rv = new VectorN(temp); - - if (print_switch) { - pr.print(t, temp); - } - - if ((t + dt) > tf) { - dt = tf - t; - } - - } - // Reset everything to before - this.ta = initial_ta; - - } - - public double[] randv(double ta) { - double p = a * (1.0 - e * e); - double cta = Math.cos(ta); - double sta = Math.sin(ta); - double opecta = 1.0 + e * cta; - double sqmuop = Math.sqrt(this.mu / p); - - VectorN xpqw = new VectorN(6); - xpqw.x[0] = p * cta / opecta; - xpqw.x[1] = p * sta / opecta; - xpqw.x[2] = 0.0; - xpqw.x[3] = -sqmuop * sta; - xpqw.x[4] = sqmuop * (e + cta); - xpqw.x[5] = 0.0; - - Matrix cmat = PQW2ECI(); - - VectorN rpqw = new VectorN(xpqw.x[0], xpqw.x[1], xpqw.x[2]); - VectorN vpqw = new VectorN(xpqw.x[3], xpqw.x[4], xpqw.x[5]); - - VectorN rijk = cmat.times(rpqw); - VectorN vijk = cmat.times(vpqw); - - double[] out = new double[6]; - - for (int i = 0; i < 3; i++) { - out[i] = rijk.x[i]; - out[i + 3] = vijk.x[i]; - } - - return out; - } - - public double eccentricAnomaly(double ta) { - double cta = Math.cos(ta); - double e0 = Math.acos((e + cta) / (1.0 + e * cta)); - return e0; - } - - public double meanAnomaly(double t) { - - return 2. * Math.PI * t / period(); - } - - public double t_from_ta() { - - double M = meanAnomaly(); - double P = period(); - - return P * M / 2. / Math.PI; - - } - - public double ta_from_t(double t) { - - double M = meanAnomaly(t); - double ea = solveKepler(M, this.e); - - double sinE = Math.sin(ea); - double cosE = Math.cos(ea); - double den = 1.0 - this.e * cosE; - double sqrome2 = Math.sqrt(1.0 - this.e * this.e); - double sinv = (sqrome2 * sinE) / den; - double cosv = (cosE - this.e) / den; - - double ta = Math.atan2(sinv, cosv); - if (this.ta < 0.0) { - this.ta = this.ta + 2.0 * Constants.pi; - } - - return ta; - } - -} From e889a08a35544c9f88d7ae45db8141c9a38a63de Mon Sep 17 00:00:00 2001 From: jic13003 Date: Sun, 6 Dec 2015 20:46:10 -0500 Subject: [PATCH 5/7] test commit --- src/jat/coreNOSA/cm/ThreeBody.java | 6 +- .../ThreeBodyExample/ThreeBodyAPL.java | 81 +++++++++++++ .../ThreeBodyExample/ThreeBodyExample.java | 114 ++++++++++++++++++ 3 files changed, 198 insertions(+), 3 deletions(-) create mode 100644 src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java create mode 100644 src/jat/examples/ThreeBodyExample/ThreeBodyExample.java diff --git a/src/jat/coreNOSA/cm/ThreeBody.java b/src/jat/coreNOSA/cm/ThreeBody.java index a9eaabd..ef46b20 100644 --- a/src/jat/coreNOSA/cm/ThreeBody.java +++ b/src/jat/coreNOSA/cm/ThreeBody.java @@ -34,9 +34,9 @@ import jat.coreNOSA.math.MatrixVector.data.VectorN; public class ThreeBody implements Derivatives //, Printable { - - private double G; // Gravitational constant - private double m1, m2, m3; // masses + protected double steps = 500.; + protected double G; // Gravitational constant + protected double m1, m2, m3; // masses /** * Method ThreeBody. diff --git a/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java new file mode 100644 index 0000000..f1d5d07 --- /dev/null +++ b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java @@ -0,0 +1,81 @@ +package jat.examples.ThreeBodyExample; + +import java.util.ArrayList; + +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.commons.math3.ode.FirstOrderDifferentialEquations; +import org.apache.commons.math3.ode.sampling.StepHandler; +import org.apache.commons.math3.ode.sampling.StepInterpolator; + +import jat.coreNOSA.algorithm.integrators.Printable; +import jat.coreNOSA.cm.Constants; +import jat.coreNOSA.cm.ThreeBody; +import jat.coreNOSA.math.MatrixVector.data.VectorN; + +public class ThreeBodyAPL extends ThreeBody implements FirstOrderDifferentialEquations{ + public ArrayList time = new ArrayList(); + public ArrayList xsol = new ArrayList(); + public ArrayList ysol = new ArrayList(); + public ArrayList zsol = new ArrayList(); + protected double[] randv[]; + public ThreeBodyAPL(double G, double m1, double m2, double m3) { + super(G, m1, m2, m3); + this.randv(); + } + + @Override + public void computeDerivatives(double t, double[] y, double[] yDot) { + // returns the derivatives of the ThreeBody problem + yDot = this.derivs(t, y); + } + + @Override + public int getDimension() { + // returns a dimension of 6 + return 18; + } + public double[] randv() { + // returns the position and velocity of the three objects + // vector r1 is x[0], x[1], x[2] + // vector v1 is x[3], x[4], x[5] + // vector r2 is x[6], x[7], x[8] + // vector v2 is x[9], x[10], x[11] + // vector r3 is x[12], x[13], x[14] + // vector v3 is x[15], x[16], x[17] + double[] randv = new double[18]; + randv[0] = 2.0; + randv[1] = 6.0; + randv[2] = 7.0; + randv[3] = 4.0; + randv[4] = 8.0; + randv[5] = 7.0; + randv[6] = 5.0; + randv[7] = 8.0; + randv[8] = 3.0; + randv[9] = 6.0; + randv[10] = 8.0; + randv[11] = 9.0; + randv[12] = 0.0; + randv[13] = 7.0; + randv[14] = 4.0; + randv[15] = 5.0; + randv[16] = 7.0; + randv[17] = 8.0; + return randv; + } + + public StepHandler stepHandler = new StepHandler() { + public void init(double t0, double[] y0, double t) { + } + + public void handleStep(StepInterpolator interpolator, boolean isLast) { + double t = interpolator.getCurrentTime(); + double[] y = interpolator.getInterpolatedState(); + System.out.println(t + " " + y[0] + " " + y[1]+ " " + y[2]); + time.add(t); + xsol.add(y[0]); + ysol.add(y[1]); + zsol.add(y[2]); + } + }; +} \ No newline at end of file diff --git a/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java b/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java new file mode 100644 index 0000000..e9243c8 --- /dev/null +++ b/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java @@ -0,0 +1,114 @@ +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; + +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); + + // 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"); + } +} From 51ab91f1d6d16a4b3f53cf830da013925a4162e7 Mon Sep 17 00:00:00 2001 From: mxd12001 Date: Mon, 7 Dec 2015 15:48:35 -0500 Subject: [PATCH 6/7] Possibly fixed the issues with merging the branch --- .../ThreeBodyExample/ThreeBodyAPL.java | 81 ++++++++++++ .../ThreeBodyExample/ThreeBodyExample.java | 115 ++++++++++++++++++ 2 files changed, 196 insertions(+) create mode 100644 src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java create mode 100644 src/jat/examples/ThreeBodyExample/ThreeBodyExample.java diff --git a/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java new file mode 100644 index 0000000..f22a529 --- /dev/null +++ b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java @@ -0,0 +1,81 @@ +package jat.examples.ThreeBodyExample; + +import java.util.ArrayList; + +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.commons.math3.ode.FirstOrderDifferentialEquations; +import org.apache.commons.math3.ode.sampling.StepHandler; +import org.apache.commons.math3.ode.sampling.StepInterpolator; + +import jat.coreNOSA.algorithm.integrators.Printable; +import jat.coreNOSA.cm.Constants; +import jat.coreNOSA.cm.ThreeBody; +import jat.coreNOSA.math.MatrixVector.data.VectorN; + +public class ThreeBodyAPL extends ThreeBody implements FirstOrderDifferentialEquations{ + public ArrayList time = new ArrayList(); + public ArrayList xsol = new ArrayList(); + public ArrayList ysol = new ArrayList(); + public ArrayList zsol = new ArrayList(); + + public ThreeBodyAPL(double G, double m1, double m2, double m3) { + super(G, m1, m2, m3); + } + + @Override + public void computeDerivatives(double t, double[] y, double[] yDot) { + // returns the derivatives of the ThreeBody problem + yDot = this.derivs(t, y); + } + + @Override + public int getDimension() { + // returns a dimension of 6 + return 18; + } + public double[] randv() { + // returns the position and velocity of the three objects + // vector r1 is x[0], x[1], x[2] + // vector v1 is x[3], x[4], x[5] + // vector r2 is x[6], x[7], x[8] + // vector v2 is x[9], x[10], x[11] + // vector r3 is x[12], x[13], x[14] + // vector v3 is x[15], x[16], x[17] + double[] randv = new double[18]; + randv[0] = 2.0; + randv[1] = 6.0; + randv[2] = 7.0; + randv[3] = 4.0; + randv[4] = 8.0; + randv[5] = 7.0; + randv[6] = 5.0; + randv[7] = 8.0; + randv[8] = 3.0; + randv[9] = 6.0; + randv[10] = 8.0; + randv[11] = 9.0; + randv[12] = 0.0; + randv[13] = 7.0; + randv[14] = 4.0; + randv[15] = 5.0; + randv[16] = 7.0; + randv[17] = 8.0; + return randv; + } + + public StepHandler stepHandler = new StepHandler() { + public void init(double t0, double[] y0, double t) { + } + + public void handleStep(StepInterpolator interpolator, boolean isLast) { + double t = interpolator.getCurrentTime(); + double[] y = interpolator.getInterpolatedState(); + System.out.println(t + " " + y[0] + " " + y[1]+ " " + y[2]); + time.add(t); + xsol.add(y[0]); + ysol.add(y[1]); + zsol.add(y[2]); + } + }; + +} \ No newline at end of file diff --git a/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java b/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java new file mode 100644 index 0000000..0b5585f --- /dev/null +++ b/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java @@ -0,0 +1,115 @@ +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; + +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); + + // 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"); + } +} \ No newline at end of file From 9d8b1db59c52bd17a0215aac8f525cbde33b26ab Mon Sep 17 00:00:00 2001 From: mxd12001 Date: Mon, 7 Dec 2015 16:36:45 -0500 Subject: [PATCH 7/7] Comments on the code are more extensive and clear up some of the confusion we've been having with what things are and how they work as well as why we implemented certains classes and functions the way we did Needs to be further updated as well --- .../ThreeBodyExample/ThreeBodyAPL.java | 17 ++++++++++++-- .../ThreeBodyExample/ThreeBodyExample.java | 22 ++++++++++++++++++- 2 files changed, 36 insertions(+), 3 deletions(-) diff --git a/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java index f22a529..aa379c1 100644 --- a/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java +++ b/src/jat/examples/ThreeBodyExample/ThreeBodyAPL.java @@ -18,16 +18,29 @@ public class ThreeBodyAPL extends ThreeBody implements FirstOrderDifferentialEqu public ArrayList ysol = new ArrayList(); public ArrayList zsol = new ArrayList(); + /* + * Why we created the ThreeBodyAPL: + * The APL is responsible for computing the derivatives used in the ThreeBodyExample and interpolating data + * APL stands for Applied Physics Laboratory + */ + + /* + * What we used from existing code + * The information that this class uses from the ThreeBody class are calculations for: + * Masses, center of mass and energy. + * This is usefull for us because all of this is necessary for our modification + */ + public ThreeBodyAPL(double G, double m1, double m2, double m3) { super(G, m1, m2, m3); } - + @Override public void computeDerivatives(double t, double[] y, double[] yDot) { // returns the derivatives of the ThreeBody problem yDot = this.derivs(t, y); } - + @Override public int getDimension() { // returns a dimension of 6 diff --git a/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java b/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java index 0b5585f..cd9a912 100644 --- a/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java +++ b/src/jat/examples/ThreeBodyExample/ThreeBodyExample.java @@ -19,6 +19,13 @@ 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() { @@ -53,8 +60,21 @@ public class ThreeBodyExample { 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);