/* ----- sympl_CRTBP.c ----- */ /* This program makes the numerical integration of the Hamiltonian differential equations related to the Circular Restricted Three-Body Problem (hereafter CRTBP). This is made by implementing the symplectic method SBAB3, in such a way that the splitting of the Hamiltonian H = A + B is made so that A is the Keplerian Hamiltonian with all the mass put in the origin of a non-inertial frame, rotating around the z-axis with angular velocity equal to 1. With such a splitting, B / A = O(mu), being mu the ratio between the masses of the primaries. Let us recall that the mass of the third body is assumed to be so negligible that it does not affect the other two bodies, that are moving on circular orbits, by following the Newton laws, so to keep their distance constantly equal to 1 (in the so called adimensional units). Of course, the motion of the third body is ruled by the Hamiltonian H. */ #include #include #include #include #include "basics_cel_mech.c" // This file, including a sort of "basic library" // for Celestial Mechanics, must stay in the // same folder where the present code is located /* Definitions of global variables */ double mu, plus1minmu; /* prototypes of the functions */ double CRTBP_energy(double x[], double p[]); void fix2rotating(double capx[], double capp[], double t, double x[], double p[]); void heliocentric_polar_coord(double x[3], double *rho, double *theta); void main_motion(double x[], double p[], double dt); void perturb_motion(double x[], double p[], double dt); void rotating2fix(double x[], double p[], double t, double capx[], double capp[]); void sbab3_kepler_perturb(double x[], double p[], double c_dt[], double d_dt[]); /* ### */ int main() { int i, j, steps_total_number, steps4sampling; double x[3], p[3]; double initial_time, final_time, sampling_time, time, h; double c[2], d[2], c_dt[2], d_dt[2]; double initial_energy, energy, err, max_err, rho, theta; FILE *ifp, *ofp; /* Input of initial data */ ifp = fopen("invmasse.inp", "r"); fscanf(ifp, "%lf", &mu); fclose(ifp); mu = (1. / mu) / (1. + 1. / mu); plus1minmu = 1. - mu; printf(" Parameters of the model:\n"); printf("mu = %le plus1minmu = %le\n\n", mu, plus1minmu); ifp = fopen("condiniz.inp", "r"); for(i=0; i<3; i++) fscanf(ifp, "%lf", &x[i]); for(i=0; i<3; i++) fscanf(ifp, "%lf", &p[i]); fclose(ifp); ifp = fopen("sympl_CRTBP.inp", "r"); fscanf(ifp, "%lf", &initial_time); fscanf(ifp, "%lf", &final_time); fscanf(ifp, "%lf", &sampling_time); fscanf(ifp, "%lf", &h); fclose(ifp); /* Evaluation of the initial energy of the system */ initial_energy = CRTBP_energy(x,p); printf("The initial energy of the entire system is: %24.16lf\n\n", initial_energy); max_err = 0.; /* Definition of the parameters ruling the numerical integration */ steps_total_number = (final_time - initial_time) / h; steps4sampling = sampling_time / h; /* In particular, the following ones are concerning with the SBAB3 symplectic integration method */ c[0] = 0.5 - sqrt(5.) / 10.; c[1] = sqrt(5.) / 5.; d[0] = 1. / 12.; d[1] = 5. / 12.; for(j=0; j<2; j++) { c_dt[j] = c[j] * h; d_dt[j] = d[j] * h; } // Opening of the output file ofp = fopen("sympl_CRTBP.out", "w"); /* Some more initial messages, preparing the rest of the output on the screen */ printf("--------------------------------------------------\n"); printf(" time Error on Energy\n"); printf("--------------------------------------------------\n"); // Iteration of sbab3_kepler_perturb for steps_total_number times for(j=1; j<=steps_total_number; j++) { sbab3_kepler_perturb(x, p, c_dt, d_dt); if(j % steps4sampling == 0) { /* Writing on screen of the messages about the error on the energy */ time = initial_time + j * h; energy = CRTBP_energy(x,p); err = fabs(energy - initial_energy); printf(" %14.6le %14.6le\n", time, err); if(max_err < err) max_err = err; /* Writing on an output file of the data about position, momenta and heliocentric polar coordinates */ heliocentric_polar_coord(x, &rho, &theta); fprintf(ofp, "%14.6le %14.6le %14.6le %14.6le %14.6le %14.6le %14.6le\n", time, x[0], x[1], p[0], p[1], theta, rho); } } /* Final operations */ fclose(ofp); printf("\n Maximal relative error on the preservation of the energy:\n"); printf("\n %14.6le\n\n", fabs(max_err / initial_energy) ); } /* ### */ double CRTBP_energy(x,p) double x[3], p[3]; { int i; double energy, denom; // printf("Entry CRTBP_energy\n"); // printf("Exit CRTBP_energy\n"); return(energy); } /* ### */ /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIX2ROTATING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PURPOSE: to compute the transformation passing from an inertial frame to the rotating one (with angular velocity equal to 1 and oriented as the versor of the z-axis) having the same origin. input: capx[0],capx[1], capx[2] ==> coordinates in the inertial frame capp[0],capp[1], capp[2] ==> momenta in the inertial frame t ==> time of the transformation (it is essential to individuate the angle of rotation between the two frames) Output: x[0],x[1],x[2] ==> coordinates in the rotating frame p[0],p[1],p[2] ==> momenta in the rotating frame %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ void fix2rotating(capx,capp,t,x,p) double capx[3],capp[3],t,x[3],p[3]; { double cost, sint, xdot[2]; // printf("Entry fix2rotating\n"); // printf("Entry fix2rotating\n"); } /* ### */ void heliocentric_polar_coord(x, rho, theta) double x[3], *rho, *theta; { *rho = sqrt( (x[0]-mu) * (x[0]-mu) + x[1] * x[1]); *theta = atan2(x[1], x[0]-mu); } /* ### */ /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MAIN_MOTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PURPOSE: Given the cartesian coordinates and the related momenta at a fixed time compute the same quantities after a given interval of time. The trajectory is calculated for the main (and integrable) part of a circular restricted three-body planetary system (CRTB), that is the Keplerian orbit due to a field attracting toward the center by a force equal to 1/r^2 (where r is the euclidean norm of the coordinates) in rotating frame with angular velocity equal to 1 and oriented as the z-axis. x[0],x[1],x[2] ==> position of the "massless" third body (in the synodical frame) p[0],p[1],p[2] ==> position of the "massless" third body (in the synodical frame) dt ==> interval of time Output: x[0],x[1],x[2] ==> position of the "massless" third body (after dt time) p[0],p[1],p[2] ==> position of the "massless" third body (after dt time) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ void main_motion(x,p,dt) double x[3], p[3], dt; { double capx[3], capp[3]; // printf("Entry main_motion\n"); // printf("Exit main_motion\n"); } /* ### */ /* ********************************************************************** PERTURB_MOTION ********************************************************************** PURPOSE: Given the cartesian coordinates and the related momenta compute the same quantities after a given interval of time. The trajectory is calculated along the flow induced by the part of the Hamiltonian corresponding to the perturbation due to the smaller of the two primary bodies in the circular restricted three-body problem (CRTB). The perturbing part of the Hamiltonian depends only on the cartesian coordinates. Thus, the cartesian coordinates are constants of motion and the motion of the related momenta is trivial. input: x[0],x[1],x[2] ==> position of the "massless" third body (in the synodical frame) p[0],p[1],p[2] ==> momentum of the "massless" third body (in the synodical frame) dt ==> interval of time mu ==> mass ratio of the smaller primary over the sum of the two primaries (it is a global variable) plus1minmu ==> mass ratio of the larger primary over the sum of the two primaries (it is a global variable) Output: x[0],x[1],x[2] ==> position of the "massless" third body (in the synodical frame, after dt time) p[0],p[1],p[2] ==> momentum of the "massless" third body (in the synodical frame, after dt time) ********************************************************************** */ void perturb_motion(x, p, dt) double x[3], p[3], dt; { int i; double delta_cub; // printf("Entry perturb_motion\n"); // printf("Exit perturb_motion\n"); } /* ### */ /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ROTATING2FIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PURPOSE: computing the transformation passing from a rotating frame (with angular velocity equal to 1 and oriented as the versor of the z-axis) to the inertial one having the same origin. input: x[0],x[1],x[2] ==> coordinates in the rotating frame p[0],p[1],p[2] ==> momenta in the rotating frame t ==> time of the transformation (it is essential to individuate the angle of rotation between the two frames) Output: capx[0],capx[1], capx[2] ==> coordinates in the inertial frame capp[0],capp[1], capp[2] ==> momenta in the inertial frame %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ void rotating2fix(x, p, t, capx, capp) double x[3], p[3], t, capx[3], capp[3]; { double cost, sint, xdot[2]; // printf("Entry rotating2fix\n"); // printf("Exit rotating2fix\n"); } /* ### */ /* sbab3_kepler_perturb: it implements one cycle of the SBAB3 method. */ void sbab3_kepler_perturb(x, p, c_dt, d_dt) double x[3]; double p[3]; double c_dt[2], d_dt[2]; { perturb_motion(x, p, d_dt[0]); main_motion(x, p, c_dt[0]); perturb_motion(x, p, d_dt[1]); main_motion(x, p, c_dt[1]); perturb_motion(x, p, d_dt[1]); main_motion(x, p, c_dt[0]); perturb_motion(x, p, d_dt[0]); }