Hi all,
I have resumed and old project about tides forces. I coded the equations for a reboundx module but I have a doubt in the last part, after compute the tides forces.
The point is that the final particle acceleration is given by:
Where F^{T} are the tides forces and F^{GR} and F^{R} are not considered in this case. So, please ignore them.
As one can notice, the main problem is that every planet acceleration depends on the forces of the others planets in the problem. So, when I coded the module (following the examples and the guide)
I have something like that:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "reboundx.h"
static void rebx_calculate_tides_forces(struct reb_simulation* const sim, const double K2s, const double K2sDt, const double Ss_z, const double Rs, const int source_index){
struct reb_particle* const particles = sim->particles;
struct reb_particle* const source = &particles[source_index];
const double G = sim->G;
const int _N_real = sim->N - sim->N_var;
for (int i=0;i<_N_real;i++){
const double* K2p = rebx_get_param_check(&particles[i], "K2p", REBX_TYPE_DOUBLE);
if(K2p == NULL) continue; // only particles with K2p set feel tidal forces
const double* K2pDt = rebx_get_param_check(&particles[i], "K2pDt", REBX_TYPE_DOUBLE);
if(K2pDt == NULL) continue; //
const double* Sp_z = rebx_get_param_check(&particles[i], "Sp_z", REBX_TYPE_DOUBLE);
if(Sp_z == NULL) continue; //
const double* Rp = rebx_get_param_check(&particles[i], "Rp", REBX_TYPE_DOUBLE);
if(Rp == NULL) continue; //
const struct reb_particle p = particles[i];
Then, I computed the forces for every i-planet inside the for loop. I don't copy here all the equations involved to avoid the mesh. At the end, we have:
double F_x = (F_tr+(Ptop+Ptos)*(v_dot_er/rps)) * er_x + Ptop*(Sp_cross_er_x - (1./rps_2)*triple_product_x) + Ptos*(Ss_cross_er_x-(1./rps_2)*triple_product_x);
double F_y = (F_tr+(Ptop+Ptos)*(v_dot_er/rps)) * er_y + Ptop*(Sp_cross_er_y - (1./rps_2)*triple_product_y) + Ptos*(Ss_cross_er_y-(1./rps_2)*triple_product_y);
double F_z = (F_tr+(Ptop+Ptos)*(v_dot_er/rps)) * er_z + Ptop*(Sp_cross_er_z - (1./rps_2)*triple_product_z) + Ptos*(Ss_cross_er_z-(1./rps_2)*triple_product_z);
// Compute the accelaration
F_x *= reduced_mass;
F_y *= reduced_mass;
F_z *= reduced_mass;
particles[i].ax += F_x;
particles[i].ay += F_y;
particles[i].az += F_z;
}
}
void rebx_tides_forces(struct reb_simulation* const sim, struct rebx_effect* const tides_forces){
double* K2s = rebx_get_param_check(tides_forces, "K2s", REBX_TYPE_DOUBLE);
if (K2s == NULL){
reb_error(sim, "Need to set potential love number of the star (K2s) in tides_forces effect.\n");
}
double* K2sDt = rebx_get_param_check(tides_forces, "K2sDt", REBX_TYPE_DOUBLE);
if (K2sDt == NULL){
reb_error(sim, "Need to set potential love number times constant time lag of the star (K2sDt) in tides_forces effect.\n");
}
double* Ss_z = rebx_get_param_check(tides_forces, "Ss_z", REBX_TYPE_DOUBLE);
if (Ss_z == NULL){
reb_error(sim, "Need to set the Spin in Z direction of the star (Ss_z) in tides_forces effect.\n");
}
double* Rs = rebx_get_param_check(tides_forces, "Rs", REBX_TYPE_DOUBLE);
if (Rs == NULL){
reb_error(sim, "Need to set the radio of the star (Rs) in tides_forces effect.\n");
}
const int N_real = sim->N - sim->N_var;
struct reb_particle* const particles = sim->particles;
int source_found=0;
for (int i=0; i<N_real; i++){
if (rebx_get_param_check(&particles[i], "Star", REBX_TYPE_INT) != NULL){
source_found = 1;
rebx_calculate_tides_forces(sim, *K2s, *K2sDt, *Ss_z, *Rs, i);
}
}
if (!source_found){
rebx_calculate_tides_forces(sim, *K2s, *K2sDt, *Ss_z, *Rs, 0); // default source to index 0 if "radiation_source" not found on any particle
}
}
However, with this formulation, I computed for i-planet the tidal forces and its acceleration, but without taken into account the forces corresponding to others planets. So, basically I'm computing only the
first term in the right hand of the equation showed at the beginning, and it is missing the second term.
Any suggestion about how to include this missing part?
Sorry in advance if it's not very clear the problem I have. I will be glad to re-explain if necessary.
Thank you very much =)
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent