Skip to content

Commit

Permalink
radius inflation, wrong interpolation
Browse files Browse the repository at this point in the history
  • Loading branch information
Tiger Lu committed Feb 19, 2024
1 parent 838ba6a commit 5a42e09
Show file tree
Hide file tree
Showing 2 changed files with 38 additions and 32 deletions.
64 changes: 35 additions & 29 deletions examples/hatp11/problem.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,14 +14,17 @@ double obl(struct reb_vec3d v1, struct reb_vec3d v2){
return acos(reb_vec3d_dot(v1,v2) / (sqrt(reb_vec3d_length_squared(v1)) * sqrt(reb_vec3d_length_squared(v2))));
}

char title[100] = "test_bigr";
//char title_stats[100] = "zlk1212_stats";
char title_remove[100] = "rm -v test_bigr";
//char title[100] = "test_bigr";
char title_stats[100] = "zlk218_stats";
//char title_remove[100] = "rm -v test_bigr";
int ind;
int printed_stats=1;
double planet_a;
double rfac;
double angle;
double planet_k2;

double tmax = 1e6*2*M_PI;
double tmax = 1e9*2*M_PI;

// eccentricity functions
double h3(double e){
Expand Down Expand Up @@ -57,13 +60,14 @@ int main(int argc, char* argv[]){
// Yee et al 2018
struct reb_particle planet = {0};
double planet_m = reb_random_uniform(sim, 0.0736 - 0.0047, 0.0736 + 0.0047) * 9.55e-4;
double planet_r = 10. * 4.2588e-5;//reb_random_uniform(sim, 4.36 - 0.06, 4.36 + 0.06) * 4.2588e-5;
planet_a = 0.5;//reb_random_uniform(sim, 0.5, 1.5);
rfac = reb_random_uniform(sim, 1.5, 2.5);
double planet_r = rfac * 4.36 * 4.2588e-5;//reb_random_uniform(sim, 4.36 - 0.06, 4.36 + 0.06) * 4.2588e-5;
planet_a = reb_random_uniform(sim, 0.3, 0.7);
double planet_e = 0.01;//reb_random_uniform(sim, 0.01, 0.1);
double planet_Omega = (117.1 - 180.) * (M_PI / 180.); //reb_random_uniform(sim, 0., 2 * M_PI);
double planet_omega = 0.;//reb_random_uniform(sim, 0.0, 2 * M_PI);
double planet_f = 0;//reb_random_uniform(sim, 0.0, 2 * M_PI);
double planet_inc = 30. * M_PI/180.;//reb_random_uniform(sim, 6. * M_PI/180., 39. * M_PI/180.);
double planet_inc = reb_random_uniform(sim, 6. * M_PI/180., 39. * M_PI/180.);
//double planet_inc = 106. * M_PI/180.;

// HAT-P-11c - treated as a point particle
Expand All @@ -81,8 +85,8 @@ int main(int argc, char* argv[]){
// Reset until over Kozai inclination
struct reb_orbit o1 = reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]);
struct reb_orbit o2 = reb_orbit_from_particle(sim->G, sim->particles[2], sim->particles[0]);
double angle = acos(reb_vec3d_dot(o1.hvec, o2.hvec) / (sqrt(reb_vec3d_length_squared(o1.hvec)) * sqrt(reb_vec3d_length_squared(o2.hvec))));
printf("%f\n", angle * 180./M_PI);
angle = acos(reb_vec3d_dot(o1.hvec, o2.hvec) / (sqrt(reb_vec3d_length_squared(o1.hvec)) * sqrt(reb_vec3d_length_squared(o2.hvec))));
//printf("%f\n", angle * 180./M_PI);

// Initial conditions
// Setup constants
Expand All @@ -108,10 +112,10 @@ int main(int argc, char* argv[]){
const double solar_spin_period = 25. * 2. * M_PI / 365.;
const double solar_spin = (2 * M_PI) / solar_spin_period;
rebx_set_param_vec3d(rebx, &sim->particles[0].ap, "Omega", (struct reb_vec3d){.z=solar_spin}); // Omega_x = Omega_y = 0 by default
rebx_set_param_double(rebx, &sim->particles[0].ap, "tau", 1e-8);
rebx_set_param_double(rebx, &sim->particles[0].ap, "tau", 1e-10);

// Planet
double planet_k2 = 0.4;//reb_random_uniform(sim, 0.4, 0.8);
planet_k2 = reb_random_uniform(sim, 0.1, 0.8);
rebx_set_param_double(rebx, &sim->particles[1].ap, "k2", planet_k2);
rebx_set_param_double(rebx, &sim->particles[1].ap, "I", 0.25 * planet_m * planet_r * planet_r);

Expand All @@ -121,7 +125,7 @@ int main(int argc, char* argv[]){
const double phi_p = 0. * M_PI / 180;
struct reb_vec3d Omega_sv = reb_tools_spherical_to_xyz(spin_p, theta_p, phi_p);
rebx_set_param_vec3d(rebx, &sim->particles[1].ap, "Omega", Omega_sv);
rebx_set_param_double(rebx, &sim->particles[1].ap, "tau", 1e-4);
rebx_set_param_double(rebx, &sim->particles[1].ap, "tau", 1e-7);


// add GR precession:
Expand All @@ -144,24 +148,18 @@ int main(int argc, char* argv[]){
//system("rm -v test.txt"); // delete previous output file
//system(title_remove);

FILE* of = fopen(title, "w");
fprintf(of, "#Seed: %d,%e,%e,%e,%e,%e,%e,%e,%e\n", index, planet_m, planet_r, planet_a, planet_e, planet_omega, planet_inc, planet_f, planet_k2);
fprintf(of, "t,a1,i1,e1,p_ob,a2,i2,e2,pert_ob,mi\n");
//FILE* of = fopen(title, "w");
//fprintf(of, "#Seed: %d,%e,%e,%e,%e,%e,%e,%e,%e\n", index, planet_m, planet_r, planet_a, planet_e, planet_omega, planet_inc, planet_f, planet_k2);
//fprintf(of, "t,a1,i1,e1,p_ob,a2,i2,e2,pert_ob,mi\n");
//fprintf(of, "t,a1,i1,e1,p_ob,mag_p,theta_p,phi_p\n");
//"t,ssx,ssy,ssz,mag1,theta1,phi1,a1,e1,nx1,ny1,nz1,nOm1,pom1,a2,e2,i2,Om2,pom2,nx2,ny2,nz2,p_ob,pert_ob,mi\n");
fclose(of);
//fclose(of);

reb_simulation_integrate(sim, tmax);
/*
for (unsigned int i = 0; i < sim->N; i++){
struct reb_particle* p = &sim->particles[i];
printf("%d %e %e %e %e %e %e\n", i, p->x,p->y,p->z,p->vx,p->vy,p->vz);
if (i == 0 || i == 1){
struct reb_vec3d* Omega = rebx_get_param(rebx, p->ap, "Omega");
printf("%e %e %e\n", Omega->x, Omega->y, Omega->z);
}
}
*/
struct reb_orbit orb = reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]);
FILE* sf = fopen(title_stats, "a");
fprintf(sf, "%d,%f,%f,%f,%f,%f,1\n", ind, orb.e, sim->particles[1].r / planet_r, planet_a, planet_k2, angle*180./M_PI);
fclose(sf);
rebx_free(rebx);
reb_simulation_free(sim);
}
Expand Down Expand Up @@ -193,7 +191,7 @@ void heartbeat(struct reb_simulation* sim){
double rdot = prefactor * (ms / mp) * pow((re / ob.a), 5.);
*/
// Output spin and orbital information to file

/*
if(reb_simulation_output_check(sim, 10. * 2 * M_PI)){ // outputs every 100 years
struct rebx_extras* const rebx = sim->extras;
Expand Down Expand Up @@ -252,9 +250,17 @@ void heartbeat(struct reb_simulation* sim){
fclose(of);
}
*/

//if(reb_simulation_output_check(sim, 20.*M_PI)){ // outputs to the screen
// reb_simulation_output_timing(sim, tmax);
//}

if(reb_simulation_output_check(sim, 20.*M_PI)){ // outputs to the screen
reb_simulation_output_timing(sim, tmax);
struct reb_orbit orb = reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]);
if (orb.a < 0.05){
FILE* sf = fopen(title_stats, "a");
fprintf(sf, "%d,%f,%f,%f,%f,%f,0\n", ind, orb.e, rfac, planet_a, planet_k2, angle*180./M_PI);
fclose(sf);
exit(1);
}
}
6 changes: 3 additions & 3 deletions src/tides_spin.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ double rebx_curlyL(double r){
double rrj = log(r/4.676e-4);

double exponent = A * rrj * rrj + B * rrj + C;
double lsun = 1.12234e-11;
double lsun = 2.84293e-13; // 10^27 erg/s to rebound units

return pow(10., exponent) * lsun;
}
Expand Down Expand Up @@ -166,7 +166,7 @@ double rebx_calculate_radius_inflation(struct reb_particle* source, struct reb_p
double denominator = 0.5 * G * ms * ms / (Rs * Rs) + alpha * ms * Rs * magp * magp;

// Curly L
double curlyL = -1. * rebx_curlyL(Rs);
double curlyL = -1.0 * rebx_curlyL(Rs);

// Tidal force
struct reb_vec3d v3_1 = reb_vec3d_mul(tidal_force, mu_ij);
Expand Down Expand Up @@ -280,7 +280,7 @@ static void rebx_spin_derivatives(struct reb_ode* const ode, double* const yDot,
unsigned int Nspins = 0;
const int N_real = sim->N - sim->N_var;
for (int i=0; i<N_real; i++){
struct reb_particle* pi = &sim->particles[i]; // target particle
struct reb_particle* pi = &sim->particles[i]; // source particle
const double* k2 = rebx_get_param(rebx, pi->ap, "k2"); // This is slow
const double* tau = rebx_get_param(rebx, pi->ap, "tau");
const double* I = rebx_get_param(rebx, pi->ap, "I");
Expand Down

0 comments on commit 5a42e09

Please sign in to comment.