From 68218746d03d8e3e5e68dbe232e13aa40442366f Mon Sep 17 00:00:00 2001 From: watem Date: Mon, 21 May 2018 18:44:01 -0400 Subject: [PATCH 1/4] Add files via upload --- BodyMaths.java | 187 ++++++++++++++++ ElectrostaticMaths.java | 30 +++ Electrostatics.java | 103 +++++++++ Main.java | 32 +++ MainPanel.java | 478 ++++++++++++++++++++++++++++++++++++++++ NBody.java | 324 +++++++++++++++++++++++++++ NBodyGraph.java | 45 ++++ Object.java | 75 +++++++ SolarBody.java | 236 ++++++++++++++++++++ 9 files changed, 1510 insertions(+) create mode 100644 BodyMaths.java create mode 100644 ElectrostaticMaths.java create mode 100644 Electrostatics.java create mode 100644 Main.java create mode 100644 MainPanel.java create mode 100644 NBody.java create mode 100644 NBodyGraph.java create mode 100644 Object.java create mode 100644 SolarBody.java diff --git a/BodyMaths.java b/BodyMaths.java new file mode 100644 index 0000000..8db887f --- /dev/null +++ b/BodyMaths.java @@ -0,0 +1,187 @@ +/** +* This logic class contains all of the equations for three body problems involving gravity or electrostatics. +* This class also contains the formulae for gravitational Lagrange points. +* @author Matthew Williams, Yulia Kosharych +* @version 2018-05-16 +**/ +public class BodyMaths { + + + // import constants + public static final double G = SolarBody.G; + + + /** + * This method calculates the acceleration that one body causes on another due to gravity + * @param obj[][] double these are the positions of the object that the acceleration is calculated on + * @param body[][] double these are the positions of the body that causes the gravitional acceleration on the object + * @param massBody double this is the mass of the body that causes the acceleration + * @param i int this is the current timestep + * @param axis int this is the axis that is currently being calculated (0=x, y=1) + * @return acceleration double this is the value of acceleration along the selected axis + * @version 2018-04-25 + **/ + public static double gravityAccel(double obj[][], double body[][], double massBody, int i, int axis){ + double acceleration; + acceleration = -G*massBody*(obj[axis][i]-body[axis][i])/Math.pow(radius(obj, body, i), 3); + return acceleration; + }//dualBodyAccel() + + + + /** + * This method calculates the magnitude of a two dimensional vector + * @param x double this is the x component of the vector + * @param y double this is the y component of the vector + * @return norm double this is the magnitude of the vector + * @version 2018-04-25 + **/ + public static double norm(double x, double y) { + double norm = Math.sqrt(x*x+y*y); + return norm; + }//norm() + + /** + * This method calculates the distance between two bodies at a specific timestep + * @param body1[][] double these are the positions of first body + * @param body2[][] double these are the positions of second body + * @param i int this is the current timestep + * @return norm this is the magnitude of the distance between the two bodies + * @version 2018-04-25 + **/ + public static double radius(double body1[][], double body2[][], int i) { + double x = body1[0][i] - body2[0][i]; + double y = body1[1][i] - body2[1][i]; + return norm(x, y); + }//radius() + + /** + * This method updates the velcity and positions for a timestep based on the acceleration + * @param obj[][][] double these are the position, velocity, acceleration vectors for a body + * @param dt double this is the timestep length + * @param i int this is the current timestep + * @version 2018-05-16 + **/ + public static void updateBody(double obj[][][], double dt, int i){ + // velocity = vI + a*dt + obj[1][0][0]+=obj[2][0][i]*dt; + obj[1][1][0]+=obj[2][1][i]*dt; + // position = xI + v*dt + obj[0][0][i]=obj[0][0][i-1]+obj[1][0][0]*dt; + obj[0][1][i]=obj[0][1][i-1]+obj[1][1][0]*dt; + }//updateBody() + + + + // Lagrange points + /** + * These methods give the distances to the legrange points. m1>m2 + * + * @param distance double this is the distance between the major bodies + * @param m1 double this is the mass of the body that is much larger than m2 + * @param m2 double this is the mass of the body that is much smaller than m1 + * + * @return r double for L1 this is the distance from the smaller object towards the larger object + * for L2 this is the distance from the smaller object away from the larger object + * for L3 this is the distance from the larger object away from the smaller object + * @version 2018-05-08 + **/ + public static double L1(double distance, double m1, double m2){ + double r = Math.abs(distance)*Math.pow(m2/(3*m1), 1./3.); + return r; + }//L1() + public static double L2(double distance, double m1, double m2){ + double r = distance*Math.pow(m2/(3*m1), 1./3.); + return r; + }//L2() + public static double L3(double distance, double m1, double m2){ + double r = 7./12.*distance*m2/m1; + return r; + }//L3() + + + + /** + * This method calculates the required orbital velocity for a circular orbit + * @param massBig double this is the mass of the driving body + * @param dis double this is the starting distance between the two bodies + * @return velocity double this is velocity required for a circular orbit + * @version 2018-05-16 + **/ + public static double circleVelocityG(double massBig, double dis){ + double velocity = Math.sqrt(G*massBig/dis); + return velocity; + }//circleVelocityG() + + /** + * This method calculates the the angularVelocity of a body in circular gravitational motion + * @param massBig double this is the mass of the driving body + * @param dis double this is the distance between the two bodies + * @return omega double this is the angular velocity of the body + * @version 2018-05-16 + **/ + public static double angularVelocity(double massBig, double dis){ + double omega = circleVelocityG(massBig, dis)/dis; + return omega; + }//angularVelocity() + + /** + * This method calculates the orbital period of two bodies in circular gravitational motion + * @param massBig double this is the mass of the driving body + * @param dis double this is the distance between the two bodies + * @return time double this is the orbital period in hours + * @version 2018-05-16 + **/ + public static double orbitalPeriod(double massBig, double dis){ + double time = 2*Math.PI/angularVelocity(massBig, dis); + return time; + }//orbitalPeriod() + + + /** + * This method converts the positions from cartesian to polar coordinates in order to make the plot relative to one body + * @param body[][][][] double these are the position, velocity, acceleration vectors for a body + * @param inertiaNum int this is the body number to use for the reference + * @version 2018-05-18 + **/ + public static void inertialReference(double body[][][][], int inertiaNum){ + int bodies = body.length; + int imax = body[0][0][0].length; + double bodyTheta[][] = new double[bodies][imax]; + double radial[][] = new double[bodies][imax]; + for (int i=0;i magPosFin ) { + + return posFinal; + + }else { + + return posInitial; + } + +} +public static double compareVelocity(double[] velInitial, double[] velFinal) { + + double magVelIn = getMag(velInitial); + double magVelFin = getMag(velFinal); + + + if (magVelIn < magVelFin ) { + + return magVelFin; + + }else { + + return magVelIn; + } + +} + +} diff --git a/NBody.java b/NBody.java new file mode 100644 index 0000000..fdd22ce --- /dev/null +++ b/NBody.java @@ -0,0 +1,324 @@ +import java.io.*; +import java.awt.*; +import javax.swing.*; +import java.util.Scanner; + +/** +* This is the driving class for the N body solution for the problem using the Euler computational analysis method. +* not to be confused with the Euler general solution for 3 body problems involving two fixed bodies. +* @author Matthew Williams, Yulia Kosharych +* @version 2018-05-16 +**/ +public class NBody{ + //presets + public static SolarBody body0 = SolarBody.Sol; + public static SolarBody body1 = SolarBody.Earth; + // public static SolarBody body1 = SolarBody.Jupiter; + // public static SolarBody body2 = new SolarBody("satellite", 100, SolarBody.au/7, SolarBody.vJupiter, 0.001, 3*Math.PI/2); + public static SolarBody body2 = SolarBody.PlanetX; + // public static SolarBody body2 = SolarBody.Horseshoe; + // public static SolarBody body2 = SolarBody.Luna; + + // public static double dis = SolarBody.solDisJupiter+Math.pow(10, 5); + // public static SolarBody body2 = new SolarBody("Tadpole orbit", 10, dis, BodyMaths.circleVelocityG(body0.mass, dis), 0, -Math.PI/6); + // public static SolarBody body2 = SolarBody.IrregularTadpole; + + // public static SolarBody body3 = new SolarBody("probe", -1, -1, SolarBody.vEarth, 0, 0); + public static SolarBody body3 = SolarBody.Jupiter; + + + // if set to true, it shows the path the body without exerting any gravity on other bodies + public static boolean setPath0 = false; + public static boolean setPath1 = false; + public static boolean setPath2 = false; + public static boolean setPath3 = false; + // useful to compare the affects of gravity on another object + + + public static final double dt = 1.; //delta t in hours + // public static final double tmax = 8765.82*5*3.54; //number of hours to run + public static final double tmax =0.1*BodyMaths.orbitalPeriod(body0.mass, body1.solDis); //number of "years" to run (in hours) + public static final int imax =(int) (tmax/dt); //number of time steps + // number of bodies to use + public static final int bodies = 3; + + public static final double G = SolarBody.G; + + + // toggle outputs + + public static final double timeStart = 6*BodyMaths.orbitalPeriod(body0.mass, body1.solDis); + public static final double timeEnd = 7*BodyMaths.orbitalPeriod(body0.mass, body1.solDis);;; + public static final boolean partial = false; + + + public static final boolean graphing = false; //plots locations + public static final boolean print = false; //prints all steps in calculations. not recommended. + public static final boolean verbose = true; //provides details as to what the program is doing + + public static final boolean inertialPlot = true; + public static final int inertiaNum = 1; //which body the inertial frame is in reference to + + + + + public static void main(String[] args) { + final int maxJavaHeap = 55855815; //max number of points that java can handle + int nPoints=bodies*3*2*imax+bodies*5; + // System.out.println(nPoints); + if (nPoints>maxJavaHeap) { + System.out.println("WARNING. Datapoints exceed maximum acceptable number of points within java (trying to use "+nPoints+" with a maximum of "+maxJavaHeap+" points)\nConsider either lowering accuracy or lowering the max time of simulation"); + // System.exit(0); + } + // body2.setName("Tadpole"); + // body2.setDisAU(); + // body2.addDis(); + // body2.addVel(); + + + + + + if (setPath0) { + body0.setPath(); + } + if (setPath1) { + body1.setPath(); + } + if (setPath2) { + body2.setPath(); + } + if (setPath3) { + body3.setPath(); + } + + Scanner kb = new Scanner(System.in); + double[] mass = new double[bodies]; + double[][][][] body= new double[bodies][3][2][imax];//[body][pos][x/y][timestep] + String bodyName[] = new String[bodies]; + double vInit[] = new double[bodies]; + double radius[] = new double[bodies]; + double theta[] = new double[bodies]; + double distance; + + + for (int i =0;i2) { + bodyName[2]=body2.name; + } + if (bodyName.length>3) { + bodyName[3]=body3.name; + } + for (int i = 4;i2) { + mass[2]=body2.mass; + } + if (mass.length>3) { + mass[3]=body3.mass; + } + for (int i=4;i2) { + theta[2]=body2.theta; + } + if (theta.length>3) { + theta[3]=body3.theta; + } + for (int i=4;i2){ + while (body2.solDis<0) { + System.out.println("distance from centre for "+body2.name); + body2.solDis = kb.nextDouble(); + } + body[2][0][0][0]=body2.solDis*Math.cos(theta[2]); + body[2][0][1][0]=body2.solDis*Math.sin(theta[2]); + } + if (body.length>3){ + while (body3.solDis<0) { + System.out.println("distance from centre for "+body3.name); + body3.solDis = kb.nextDouble(); + } + body[3][0][0][0]=body3.solDis*Math.cos(theta[3]); + body[3][0][1][0]=body3.solDis*Math.sin(theta[3]); + } + for (int i = 4;i2){ + // while (body2.vInitial<0) { + // System.out.println("initial velocity for "+body2.name); + // body2.vInitial = kb.nextDouble(); + // } + body[2][1][0][0]=-body2.vInitial*Math.sin(theta[2]); + body[2][1][1][0]=body2.vInitial*Math.cos(theta[2]); + } + if (body.length>3){ + // while (body3.vInitial<0) { + // System.out.println("initial velocity for "+body3.name); + // body3.vInitial = kb.nextDouble(); + // } + body[3][1][0][0]=-body3.vInitial*Math.sin(theta[3]); + body[3][1][1][0]=body3.vInitial*Math.cos(theta[3]); + } + for (int i = 4;i Date: Mon, 21 May 2018 18:47:09 -0400 Subject: [PATCH 3/4] Delete ThreeBodyInstantaneousSimulationVersion --- ThreeBodyInstantaneousSimulationVersion | 329 ------------------------ 1 file changed, 329 deletions(-) delete mode 100644 ThreeBodyInstantaneousSimulationVersion diff --git a/ThreeBodyInstantaneousSimulationVersion b/ThreeBodyInstantaneousSimulationVersion deleted file mode 100644 index 22ac7a1..0000000 --- a/ThreeBodyInstantaneousSimulationVersion +++ /dev/null @@ -1,329 +0,0 @@ - import java.util.Scanner; - import java.util.Arrays; - -// TODO: javadocs - - - public class ThreeBody { - - //Final variables declaration - - public static final double G = 8.64960768*Math.pow(10, -13);//gravitational constant (km^3/kg/h^2) - public static final double dt = 0.1; //delta t in hours - public static double tmax ; //number of hours to run - public static double time = 0; //initial time - // public static final double tmax; //number of hours to run {inputs} - public static final int imax =(int) (tmax/dt); //number of time steps - // public static final double imax; //number of time steps {inputs} - public static final boolean graphing = true; - public static final double au = 1.49597870700*Math.pow(10,8);//initial distance of earth from the sun in km - // public static final double radEarth=; //km - // public static final double radSun=; //km - - public static final double massEarth=5.972*Math.pow(10, 24); //kg - public static final double massSun=1.989*Math.pow(10, 30); //kg - public static double massSat; //kg - public static final double velocityEarth = 1.08*Math.pow(10, 5);// initial velocity of Earth - public static Scanner kb = new Scanner(System.in); - - - - public static void main(String[] args) { - - - System.out.print("Please enter the duration of the simulation: "); - tmax = kb.nextDouble(); - - System.out.print("Please enter the mass of the satellite: "); - massSat = kb.nextDouble(); - - double satX, satY, velInitialX, velInitialY; - - // earthX[0][0]=1.495978*Math.pow(10,16); - // earthX[1][0]=0; - // earthY[1][0]=108000; - // earthY[0][0]=0; - // System.out.println("please enter the mass of the satellite: "); - // massSat = kb.nextDouble(); - System.out.print("Please enter the initial x position of the satellite: "); - satX = kb.nextDouble(); - System.out.print("Please enter the initial y position of the satellite: "); - satY = kb.nextDouble(); - System.out.print("Please enter the initial x velocity of the satellite: "); - velInitialX = kb.nextDouble(); - System.out.print("Please enter the initial y velocity of the satellite: "); - velInitialY = kb.nextDouble(); - - - - - - - double [][][] object = new double [3][3][2]; - - // object[a] : a = 0(Earth), a=1(Sun), a=2(Sat) -// object[][b] : b = 0(position), b=1(velocity), b=2(acceleration) -//object[][][c] : c = 0(x position), c=1(y position) - - //Earth - object[0][0][0] = au; - object[0][0][1] = 0; - - object[0][1][0] = 0; - object[0][1][1] = velocityEarth; - - object[0][2][0] = 0; - object[0][2][1] = 0; - - //Sun - object[1][0][0] = 0; - object[1][0][1] = 0; - - object[1][1][0] = 0; - object[1][1][1] = 0; - - object[1][2][0] = 0; - object[1][2][1] = 0; - - //Sat - object[2][0][0] = satX; - object[2][0][1] = satY; - - object[2][0][0] = velInitialX; - object[2][0][1] = velInitialY; - - object[2][2][0] = 0; - object[2][2][1] = 0; - - - - - System.out.println(Arrays.toString(accelSunAnyObj(massSun,object[0][0],object[1][0]))); - - double [] maxAccel, minPosSun, minPosEarth; - maxAccel = object[2][2]; - - minPosSun = differencePos(object[2][0],object[1][0]); - minPosEarth = differencePos(object[2][0],object[0][0]); - - while (time magPosFin ) { - - return posFinal; - - }else { - - return posInitial; - } - -} - - - - public static double compareVelocity(double[] velInitial, double[] velFinal) { - - double magVelIn = getMag(velInitial); - double magVelFin = getMag(velFinal); - - - if (magVelIn < magVelFin ) { - - return magVelFin; - - }else { - - return magVelIn; - } - -} - - -//TODO: graphing method() -// public static void graph(double[][] satX,double[][] satY,double[][] sunX,double[][] sunY,double[][] earthX,double[][] earthY) - - - - -}//end class ThreeBody From a1c145c36f65321227cae262bcf68e3279337699 Mon Sep 17 00:00:00 2001 From: watem Date: Mon, 21 May 2018 22:25:10 -0400 Subject: [PATCH 4/4] Add files via upload --- ElectrostaticMaths.java | 5 ++-- MainPanel.java | 57 ++++++++++++++++++++++++++++++++++++----- 2 files changed, 53 insertions(+), 9 deletions(-) diff --git a/ElectrostaticMaths.java b/ElectrostaticMaths.java index d0a753f..2e03891 100644 --- a/ElectrostaticMaths.java +++ b/ElectrostaticMaths.java @@ -4,7 +4,7 @@ * @author Matthew Williams, Yulia Kosharych * @version 2018-05-16 **/ -public class BodyMaths { +public class ElectrostaticMaths { // import constants public static final double Ke = Electrostatics.Ke; @@ -23,8 +23,7 @@ public class BodyMaths { **/ public static double electrostaticAccel(double obj[][], double body[][], double qObj, double qBody, double massObj, int i, int axis){ double acceleration; - acceleration = Ke*qObj*qBody*(obj[axis][i]-body[axis][i])/Math.pow(Bodymaths.radius(obj, body, i), 3)/massObj; + acceleration = Ke*qObj*qBody*(obj[axis][i]-body[axis][i])/Math.pow(BodyMaths.radius(obj, body, i), 3)/massObj; return acceleration; }//electrostaticAccel() } - diff --git a/MainPanel.java b/MainPanel.java index 5ad2a02..61b8a86 100644 --- a/MainPanel.java +++ b/MainPanel.java @@ -41,6 +41,9 @@ public class MainPanel extends JPanel{ public static final double au = SolarBody.au;//initial distance of earth from the sun in km + // public static final double massEarth=SolarBody.massEarth; //kg + // public static final double velocityEarth = SolarBody.vEarth;// initial velocity of Earth + public static final double massEarth=SolarBody.massEarth; //kg public static final double velocityEarth = SolarBody.vEarth;// initial velocity of Earth @@ -48,9 +51,9 @@ public class MainPanel extends JPanel{ public static final double velocityJupiter = SolarBody.vJupiter;//velocity jupiter - public static final double massL5 = SolarBody.SE_L5; //kg + public static final double massL5 = SolarBody.SE_L5.mass; //kg - public static final double massSun=SolarBody.massSun; //kg + public static final double massSun=SolarBody.massSol; //kg public static double massSat; //kg public static Scanner kb = new Scanner(System.in); @@ -69,6 +72,7 @@ public class MainPanel extends JPanel{ public void initSunEarthLuna() { scaling = 350; + //Sat Conditions massSat = SolarBody.massLuna; satX = SolarBody.solDisLuna; @@ -84,17 +88,36 @@ public void initSunEarthLuna() velInitialPlanetY = velocityEarth; } + public void initSunJupiterCallisto() + { + scaling = 75; + //Sat Conditions + massSat = SolarBody.massCallisto; + satX = SolarBody.solDisCallisto+9*SolarBody.jupDisCallisto; + satY = 0; + velInitialX = 0; + velInitialY = SolarBody.vJupiter+BodyMaths.circleVelocityG(massEarth, 10*SolarBody.jupDisCallisto); + + //Jupiter Condition + massPlanet = massJupiter; + planetX = SolarBody.solDisJupiter; + planetY = 0; + velInitialPlanetX = 0; + velInitialPlanetY = velocityJupiter; + + } public void initSunEarthL5() { scaling = 350; //Sat Conditions + double theta = -Math.PI/3; massSat = massL5; - satX = au*Math.cos(-Math.PI/3); - satY = au*Math.sin(-Math.PI/3); - velInitialX = -velocityEarth*Math.sin(theta[0]); - velInitialY = velocityEarth*Math.cos(theta[0]); + satX = au*Math.cos(theta); + satY = au*Math.sin(theta); + velInitialX = -velocityEarth*Math.sin(theta); + velInitialY = velocityEarth*Math.cos(theta); //Planet Condition massPlanet = massEarth; @@ -131,6 +154,26 @@ public void initSunJupiterProbe() } + public void initHorseshoe() + { + scaling = 75; + //Sat Conditions + double theta = SolarBody.Horseshoe.theta; + massSat = SolarBody.Horseshoe.mass; + satX = SolarBody.Horseshoe.solDis*Math.cos(theta); + satY = SolarBody.Horseshoe.solDis*Math.sin(theta); + velInitialX = -SolarBody.Horseshoe.vInitial*Math.sin(theta); + velInitialY = SolarBody.Horseshoe.vInitial*Math.cos(theta); + + //Planet Condition + massPlanet = massJupiter; + planetX = SolarBody.solDisJupiter; + planetY = 0; + velInitialPlanetX = 0; + velInitialPlanetY =velocityJupiter; + + } + public void Algorithm() { @@ -182,7 +225,9 @@ public MainPanel() if (sim == 3) initSunJupiterProbe(); + if (sim == 4) initSunJupiterCallisto(); + if (sim == 5) initHorseshoe(); Algorithm();