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 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