//
// orbit.cpp
//
// This program calculates the orbit of an object around a black hole, or other heavy object,
// using the equations of motion from general relativity.
//
// This program was made for my final project in the course Phys291 at the University of Bergen,
// in the spring of 2020. Instructions for using the program are found at: https://folk.uib.no/hpl001/project/
// I will refer to this webpage in the comments when it is appropriate.
// If you are unsure of the definition of some expression used in the code,
// you will probably find the it there.
//
{
//---------------------------------------------------
//--------------------------------------------------- -->SET NUMBER OF ROTATIONS AND STEP SIZE HERE!<--
//
const int N = 3; // Sets the number of rotations around the black hole that are to be calculated.
const int n = 10000; // Sets the number of steps in the calculation of one rotation.
// The needed number of steps for a good calculation will depend on the parameters
// A very elliptical orbit requires smaller steps
//
//---------------------------------------------------
//---------------------------------------------------
const double pi = 3.14159265359;
const double c = 299792000; // The speed of light in meters/second.
double E2; // The energy -->per mass<-- squared of the orbiting object.
double E; // The energy -->per mass<-- of the orbiting object.
double L2; // The angular momentum -->per mass<-- squared of the orbiting object.
double L; // The angular momentum -->per mass<-- of the orbiting object.
double R; // The Schwarzschild Radius of the source of the gravitational field.
// See the webpage for the definition of the Schwarzschild radius.
double dN = N; // A double-version of N
double dn = n; // A double-version of n
double h = (2*pi)/dn; // The step size in the calculation.
const int steps = n*N; // The total number of steps in the calculation.
const int points = steps + 1;
double t = 0; // The proper time.
double T = 0; // The coordinate time.
double r_apo; // The radius of the apoapsis.
double r_peri; // The radius of the periapsis.
double phi_peri; // The angle of the first periapsis.
double phi_apo; // The angle of the first apoapsis.
double precession; // The precession-angle.
// all the points along the orbit will be put into the double vectors r, x, and y.
// We solve for r(\phi) first, and convert all these points to xy-coordinates afterward.
// ehx and ehy are for plotting the event horizon.
double r[points], x[points], y[points], ehx[n], ehy[n];
cout << " This program calculates orbits around black holes." << endl;
cout << " For more details about this program and the project, visit https://folk.uib.no/hpl001/Project/" << endl;
cout << " There you will also find some values of the input-parameters for different known celestial bodies, " <> r_peri;
cout << " Please enter the apoapsis of the orbit: " << endl;
cin >> r_apo;
cout << " Please enter the Schwarzschild radius:" << endl;
cin >> R;
r_peri = 1000*r_peri; // Converts the inputs to meters
r_apo = 1000*r_apo;
R = 1000*R;
// Calculation of the energy and angular momentum of the orbit.
double numerator = R*(r_peri*r_peri*r_apo - r_apo*r_apo*r_peri);
double denominator = r_peri*r_peri*(1 - R/r_apo) - r_apo*r_apo*(1-R/r_peri);
L2 = numerator/denominator;
L = sqrt(L2);
E2 = (1 - R/r_apo)*(1 + L2/(r_apo*r_apo));
E = sqrt(E2);
// Sets the radius where the calculation starts.
r[0] = 0.5*(r_peri + r_apo);
// sgn keeps track of the sign of dr_dphi.
// The orbit starts out going inward, so sgn = -1 at the beginning.
// See the web-page for more.
int sgn = -1;
// a and b keep track of when the object is between the first apoapsis and
// the first periapsis. At the first apoapsis the program sets b -> 0, and a -> 1,
// and the calculation of the orbital time(s) starts.
// At the first periapsis, a->0, and the time-calculations stop.
int a = 0;
int b = 1;
// s1 and s2 are the first and second derivatives, respectively, of r with respect to the angle phi.
// r_d is used in a test-calculation to ensure that r_{i+1} will not break the formula for dr/dphi.
// See the computation methods-section in the .pdf for details.
double s1, s2, r_d;
// It's time to start the calculation!
// For each loop, a test-iteration is done, to ensure that the calculated r[i+1] lies within the
// possible range of r, where the expression for the derivative actually makes sence.
// If the result is acceptable, the calculation of r[i+1] is done without any change.
// If the the result is UNacceptable, the sign of the first derivative is changed BEFORE the calculation is done.
//
// After the calculation is done, the code checks if there is a turning point along
// the curve between r[i] and r[i+1], using the second order approximation of the curve.
// If there is a turning point, the sign is changed AFTER r[i+1] is calculated,
// and before r[i+2] is calculated.
int i;
for(i = 0; i < steps + 1; i++){
// Calculation of the derivatives.
s1 = sgn*((r[i]*r[i])/L)*sqrt(E2 - (1 - R/r[i])*(L2/(r[i]*r[i]) + 1));
s2 = (2*(E2 - 1)*r[i]*r[i]*r[i])/(L2)
+(3*R*r[i]*r[i])/(2*L2)
- r[i]
+ R/2;
r_d = r[i] + s1*h + 0.5*h*h*s2; // The test-iteration
// Changes the sign of the derivative before the actual calculation if this is needed.
// The condition is that r_d lie outside the range where dr/dphi is properly defined.
if(E2 - (1 - R/r_d)*(1 + (L2)/(r_d*r_d)) < 0){
sgn = -sgn;
s1 = -s1;
// If the orbit hits the first periapsis, a -> 0, so that the time-calculations stop.
if(a == 1){
a = 0;
phi_apo = h*i;
}
// If the orbit hits the first apoapsis, a -> 1, and b ->0, so that the time-calculations start to run.
if(b == 1){
b = 0;
a = 1;
phi_peri = h*i;
}
}
r[i+1] = r[i] + s1*h + 0.5*h*h*s2;
// Calculates the proper time and coordinate time when a == 1.
// a == 1 between the first apoapsis and the first periapsis.
if(a == 1){
t = t + h*(r[i]*r[i])/(c*L);
T = T + h*((r[i]*r[i])/(c*L))*E*(1/(1 - R/r[i]));
}
// These two parts changes the sign before the next calculation if this is needed.
//
if(sgn == 1 and (s1 + s2*h < 0)){
sgn = -1;
if(a == 1){
phi_apo = h*i;
a = 0;
}
}
if(sgn == -1 and (s1 + s2*h > 0)){
sgn = 1;
if( b == 1){
phi_peri = h*i;
b = 0;
a = 1;
}
}
}
// The orbital times calculated in the loop will be half the time of one full orbit.
T = 2*T;
t = 2*t;
precession = 2*(phi_apo - phi_peri - pi);
// Puts all the points in r[i] into xy-coordinates:
// Values are converted back to kilometers
for(i = 0; i < steps + 1; i++){
x[i] = r[i]*sin(h*i - phi_peri);
y[i] = -r[i]*cos(h*i - phi_peri);
}
// Sets up the points in the event horizon.
for(i = 0; i < n; i++){
ehx[i] = R*sin(h*i - phi_peri);
ehy[i] = -R*cos(h*i - phi_peri);
}
TCanvas *can = new TCanvas("can", "can", 0, 0, 400, 400);
TGraph *Orbit = new TGraph(steps + 1, x, y);
TGraph *EventHorizon = new TGraph(n + 1, ehx, ehy);
Orbit->SetTitle("Orbit");
Orbit->Draw();
EventHorizon->Draw("SAME");
cout << " Result of calculation: " << endl;
cout << " Energy per mc^2: " << E << endl;
cout << " Angular momentum per mc: " << L << endl;
cout << " Proper time of orbit: " << t << " Seconds" << endl;
cout << " Coordinate time of orbit: " << T << " Seconds" << endl;
cout << " Time-dilation factor T/t: " << T/t << endl;
cout << " Precession of the orbit: " << precession*(360/(2*pi)) << " Degrees" << endl;
cout << " " << endl;
}