modify mass C equation about reboundx

[dtamayo]

Hi Mileva, That code is taking the differential equation dm/dt = m/tau, which has the solution m(t) = m(0)e^(t/tau), but solving it approximately as delta m = m/tau * delta t. But you could implement it exactly, or however you want. The advantage of this is that this is machine independent, the exponential function will depend on your hardware/compiler etc. Dan

…

On Aug 28, 2017, at 4:11 PM, MilevaE ***@***.***> wrote: Hi all, I'm using the external module modify_mass.c in REBOUNDx. What I want is to implement another kind of mass loss instead of the e-folding. Thus, I went to the src code in order to learn how implement this change. However, I don't understand very well the following part of the code: if (tau_mass != NULL){ p->m += p->m*dt/(*tau_mass); I know that is equivalent in python to: M0=initial mass M = M0*np.e**(times/tau_mass) However, I don't understand where is implemented the exponential in the original C code. What is dt in that equation? I know that my problem here is that I'm bad C code person....yet =) Thanks in advance for the help! Best, Mileva — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#21>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ADQAlxNP8vSFOOo9DXZzzRmkdcliJNyRks5scx7cgaJpZM4PFAQD>.

[MilevaE]

Hi Dan,

thanks you very much for the explanation!, now it's clear =)

Cheers,

Mileva

[MilevaE]

Hi again,

sorry for re-open when was already closed =(

As I mentioned, I'm learning about how to implement my own effects. Just for learn, I try to implement a linear mass loss/growth, based in the e-folding code already available. I followed all the steps, which I think that are ok. Everything compile correctly. However, when I call from my python code, the integration never starts =(

This is the C code for the effect (copy from modify_mass.c and change a couple of things):

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "rebound.h"
#include "reboundx.h"

void rebx_modify_mass_linear(struct reb_simulation* const sim, struct rebx_effect* const effect, const double dt, enum rebx_timing timing){
    const int _N_real = sim->N - sim->N_var;
	for(int i=0; i<_N_real; i++){
		struct reb_particle* const p = &sim->particles[i];
        const double* const M_dot = rebx_get_param_check(p, "M_dot", REBX_TYPE_DOUBLE);
        const double* const initial_mass = rebx_get_param_check(p, "initial_mass", REBX_TYPE_DOUBLE);
        if (M_dot != NULL){
		    p->m += *initial_mass + (*M_dot*dt);
        }
	}
    reb_move_to_com(sim);
} 

I have defined M_dot in units of M_sun/year, i.e., the rate of mass loss. And Mo the initial mass of
the star in M_sun. So, what I think is doing is:

m = Mo + M_dot*dt

Thus, in principle, just a linear (and simple) equation.

From python, I call doing:

Mo = 3.04
rebx = reboundx.Extras(sim)
mass_effect_linear = rebx.add("modify_mass_linear")

ps[0].params["M_dot"] = -1.e-5
ps[0].params["initial_mass"] = Mo

Nout = 1e5
mass = np.zeros(Nout)
years=1e5
times = np.linspace(0., years, Nout)
time_step = np.zeros((0,Nout))
semimayor = np.zeros((2,Nout))
ecc = np.zeros((2,Nout))
sim.integrator = "ias15"
#sim.integrator = "whfast"


for i, timex in enumerate(times):
    sim.integrate(timex)
    mass[i] = sim.particles[0].m
    semimayor[0][i] = sim.particles[1].a
    semimayor[1][i] = sim.particles[2].a
    ecc[0][i] = sim.particles[1].e
    ecc[1][i] = sim.particles[2].e
    time_step = np.append(time_step,sim.dt)
    print i

So, any advice what I'm doing wrong and why the integration never starts.. :/

Thanks a lot for the help!

Cheers,

Mileva

[dtamayo]

Hi Mileva,

what do you mean it doesn't start? What error do you get in the notebook and/or in the terminal from which you launched the notebook?

Also, I know it's just a test, but I think your equation updating the masses will always update to the same value assuming a constant timestep. I think you want p->m += M_dotdt;

[MilevaE]

Hi Dan,

yes, with your advice works perfect. Thanks!

When I said "the integration never starts" I mean, that, as you can see in the code:

for i, timex in enumerate(times):
    sim.integrate(timex)
    mass[i] = sim.particles[0].m
    semimayor[0][i] = sim.particles[1].a
    semimayor[1][i] = sim.particles[2].a
    ecc[0][i] = sim.particles[1].e
    ecc[1][i] = sim.particles[2].e
    time_step = np.append(time_step,sim.dt)
    print i

I inserted a "print i" at the end, in order to follow what is going on. However, in this case, never start to count, doesn't print anything. With your correction yes, it prints =)

Thanks!

Best,

Mileva