Mixing of Particles' Indices

[ruaSulaiman]

Hello,

I'm using Rebound to simulate 1000 semi-active particles, and I've disabled collision detection by setting it to None. I'm wondering if I can assume that the particles at specific indices in the orbits file (that is generated using the function reb_output_orbits) don't mix with each other. This is considering:

1. The simulation runs for 10 days, with each day cut into 24-hour segments. I need to restart the simulation from the archive file at the start of each new day.
2. The initial conditions are generated using the Python version of Rebound and stored in a file named "initial_conditions.txt."

This is because I would like to track the final semi-major axis of each particle.
I would appreciate your help. Here is my code:

#########################################
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include "rebound.h"

double rand_range(double min, double max);
void heartbeat(struct reb_simulation* r);
double E0;
double ang_m0[5];
double ang_m[5];
double com_rel[6];
double L_x;
double L_y;
double L_z;
double com_tot_vmag;
double com_vmag;
struct reb_particle com0;
const int N_planetesimals = 1000;
void calculate_com(struct reb_simulation* const r, struct reb_particle com0, double com_rel[]);
void calculate_angular_momentum(struct reb_simulation* const r, double L_tot[]);

int main(int argc, char* argv[]){

struct reb_simulation* r = reb_create_simulation();

reb_simulationarchive_automate_step(r, "archive.bin", 1e5);
// Simulation Setup
r->integrator = REB_INTEGRATOR_BS;
r->ri_bs.eps_rel = 1e-10;
r->ri_bs.eps_abs = 1e-10;
r->heartbeat    = heartbeat;

r->testparticle_type = 1;

// Collisions
r->collision = REB_COLLISION_NONE;

int i;
int j;
int particle_N = N_planetesimals + 2;
double** mat=malloc(particle_N*sizeof(double*)); 
FILE *file;
file=fopen("initial_conditions.txt", "r");
for(i=0;i<particle_N;++i)
mat[i]=malloc(7*sizeof(double));

for(i = 0; i < particle_N; i++){
	for(j = 0; j < 7; j++){
		if (!fscanf(file, "%lf", &mat[i][j])) 
		break;
}
}
fclose(file);

for (int i=0;i<particle_N;i++){
	struct reb_particle p = {0};
	p.x = mat[i][1];	 	p.y = mat[i][2];		p.z  = mat[i][3];
	p.vx = mat[i][4]; 		p.vy = mat[i][5];	 	p.vz = mat[i][6];
p.m  = mat[i][0];		
reb_add(r, p); 
}
r->N_active = 2;

srand(12);
double a_scat_planet = 10;
double a_disk = 20;
r->dt = 6.283*pow(a_scat_planet,1.5)/50; //2*pi()a^(1.5)

reb_move_to_com(r);
E0 = reb_tools_energy(r);
system("rm -rf energy.txt");

calculate_angular_momentum(r, ang_m0);
system("rm -rf angular_m.txt");

com0 = reb_get_com(r);
system("rm -rf com.txt");

FILE* f5 = fopen("conservation.txt","a");
reb_integrator_synchronize(r);
fprintf(f5,"%e %e %e %e %e %e %e %e %e %e %e\n",E0, com0.x, com0.y, com0.z, com0.vx, com0.vy, com0.vz, ang_m0[0], ang_m0[1], ang_m0[2], ang_m0[3]);
fclose(f5);

reb_integrate(r,2*M_PI*5*1e7);
reb_free_simulation(r);

}

void heartbeat(struct reb_simulation* r){
if (reb_output_check(r, 10000)){ // every ~ 1600 years
//relative energy error

    FILE* f = fopen("energy.txt","a");
    reb_integrator_synchronize(r);
    double E = reb_tools_energy(r);
    fprintf(f,"%e %e %d\n",r->t, fabs((E-E0)/E0), r->N);
    fclose(f);
    
    FILE* f2 = fopen("angular_m.txt","a");
    reb_integrator_synchronize(r);
    calculate_angular_momentum(r, ang_m);
    fprintf(f2,"t = %e, L_ini_me = %e, L = %e, L_mag_rel = %e, L_x_rel = %e, L_y_rel = %e, L_z_rel = %e \n", r->t, ang_m0[3], ang_m[3], fabs((ang_m0[3] - ang_m[3])/ang_m0[3]), fabs((ang_m0[0] - ang_m[0])/ang_m0[0]), fabs((ang_m0[1] - ang_m[1])/ang_m0[1]),fabs((ang_m0[2] - ang_m[2])/ang_m0[2]));
    fclose(f2);
    
    FILE* f3 = fopen("com.txt","a");
    reb_integrator_synchronize(r);
    calculate_com(r, com0, com_rel);
    fprintf(f3,"t = %e, com_vmag_rel = %e, com0_vmag = %e, com_vmag = %e\n", r->t, com_rel[3], com_rel[4], com_rel[5]);
    fclose(f3);

    //get orbital elements
    struct reb_particle p = r->particles[1];
    struct reb_particle star = r->particles[0];
    struct reb_orbit o = reb_tools_particle_to_orbit(r->G,p,star);
    
    printf("a_p=%f,e_p=%f,N=%d\n",o.a,o.e,r->N);
}
if (reb_output_check(r, 1e5)){

    reb_integrator_synchronize(r);
    reb_output_orbits(r,"orbits.txt");
    
    FILE *file;
    file = fopen("orbits.txt","a");
    fprintf(file,"Next \n");
    fclose(file);
}

}

void calculate_angular_momentum(struct reb_simulation* const r, double L_tot[]) {
//reb_move_to_com(r);
L_x = 0;
L_y = 0;
L_z = 0;
struct reb_particle p;
for (int i = 0; i < r->N; i++) {
p = r->particles[i];
L_x += p.m*(p.yp.vz - p.zp.vy);
L_y += -p.m*(p.xp.vz - p.zp.vx);
L_z += p.m*(p.xp.vy - p.yp.vx);
}
double L_t = sqrt(L_zL_z + L_xL_x + L_y*L_y);
L_tot[0] = L_x;
L_tot[1] = L_y;
L_tot[2] = L_z;
L_tot[3] = L_t;
}

void calculate_com(struct reb_simulation* const r, struct reb_particle com0, double com_rel[]) {

struct reb_particle com = reb_get_com(r);
com_rel[0] = fabs((com0.x - com.x)/com0.x);
com_rel[1] = fabs((com0.y - com.y)/com0.y);
com_rel[2] = fabs((com0.z - com.z)/com0.z);
com_tot_vmag = sqrt(com0.vx*com0.vx + com0.vy*com0.vy + com0.vz*com0.vz);
com_vmag =  sqrt(com.vx*com.vx + com.vy*com.vy + com.vz*com.vz);
com_rel[3] = fabs((com_tot_vmag - com_vmag)/com_tot_vmag);
com_rel[4] = com_tot_vmag;
com_rel[5] = com_vmag;

}

[hannorein]

Yes. The order will be fixed in your case.

[hannorein]

I am closing this for now, but if you have further questions, feel free to reopen it.