From 03fe0c5a395d5879bd16776fc4fe38734c8e6611 Mon Sep 17 00:00:00 2001 From: mxd12001 Date: Fri, 4 Dec 2015 19:38:52 -0500 Subject: [PATCH] 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"); + + } +}