POSYDON-code / POSYDON

POSYDON is a next-generation single and binary-star population synthesis code incorporating full stellar structure and evolution modeling with the use of MESA.
BSD 3-Clause "New" or "Revised" License
25 stars 19 forks source link

Adding option to set q_max to a user specified value #302

Open astroJeff opened 1 month ago

astroJeff commented 1 month ago

This update allows for a user to change the default initial minimum and maximum mass ratio of the two stars in the binary. Note that both the mass ratio limits and the secondary mass limits are simultaneously satisfied. This has been tested and is ready to go.

Defaults are: q_min = 0.05 q_max = 1.0 secondary_mass_min=0.35 secondary_mass_max=120

An isolated test of the function:

In [22]: M1 = generate_primary_masses(number_of_binaries=1000)

In [23]: M2 = generate_secondary_masses(M1)

In [24]: q = M2/M1

In [25]: print(np.min(q), np.max(q), np.min(M2), np.max(M2))
0.055159202746417196 0.9978797254217177 0.4031951085983805 92.65663979838298

In [26]: M2 = generate_secondary_masses(M1, q_min=0.01, q_max=0.5)

In [27]: q = M2/M1

In [28]: print(np.min(q), np.max(q), np.min(M2), np.max(M2))
0.011247410652580144 0.49886272923081393 0.3504552319194787 54.75864549435622

In [29]: 
mkruckow commented 1 month ago

It is still failing. Here an example:

In [8]: M1 = generate_primary_masses(number_of_binaries=1000)

In [9]: M2 = generate_secondary_masses(M1, secondary_mass_max=2.0)

In [10]: q = M2/M1

In [11]: print(np.min(q), np.max(q), np.min(M2), np.max(M2))
0.01747108587586304 0.2718665208879266 0.3753374686092538 5.4458254240205015

When you check all the 24 possible situations you'll see that there are still two cases, where the user input is OK, (because q_min<q_max and m2_min<m2_max; the user can't check more because m1 is not set at this moment, hence can't be considered), while the resulting range isn't OK (because an actual max value is smaller than a min value).

size order sampling range user
m2_min/m1 < q_min < m2_max/m1 < q_max [q_min,m2_max/m1] OK OK
m2_min/m1 < q_min < q_max < m2_max/m1 [q_min,q_max] OK OK
m2_min/m1 < m2_max/m1 < q_min < q_max [m2_max/m1,q_min] fail OK
m2_min/m1 < m2_max/m1 < q_max < q_min [m2_max/m1,q_min] fail mistake
m2_min/m1 < q_max < q_min < m2_max/m1 [q_max,q_min] fail mistake
m2_min/m1 < q_max < m2_max/m1 < q_min [q_max,q_min] fail mistake
q_min < m2_min/m1 < m2_max/m1 < q_max [m2_min/m1,m2_max/m1] OK OK
q_min < m2_min/m1 < q_max < m2_max/m1 [m2_min/m1,q_max] OK OK
q_min < m2_max/m1 < m2_min/m1 < q_max [m2_max/m1,m2_min/m1] fail mistake
q_min < m2_max/m1 < q_max < m2_min/m1 [m2_max/m1,m2_min/m1] fail mistake
q_min < q_max < m2_min/m1 < m2_max/m1 [q_max,m2_min/m1] fail OK
q_min < q_max < m2_max/m1 < m2_min/m1 [q_max,m2_min/m1] fail mistake
m2_max/m1 < m2_min/m1 < q_min < q_max [m2_max/m1,q_min] fail mistake
m2_max/m1 < m2_min/m1 < q_max < q_min [m2_max/m1,q_min] fail mistake
m2_max/m1 < q_min < m2_min/m1 < q_max [m2_max/m1,m2_min/m1] fail mistake
m2_max/m1 < q_min < q_max < m2_min/m1 [m2_max/m1,m2_min/m1] fail mistake
m2_max/m1 < q_max < m2_min/m1 < q_min [m2_max/m1,q_min] fail mistake
m2_max/m1 < q_max < q_min < m2_min/m1 [m2_max/m1,m2_min/m1] fail mistake
q_max < m2_min/m1 < q_min < m2_max/m1 [q_max,q_min] fail mistake
q_max < m2_min/m1 < m2_max/m1 < q_min [q_max,q_min] fail mistake
q_max < q_min < m2_min/m1 < m2_max/m1 [q_max,m2_min/m1] fail mistake
q_max < q_min < m2_max/m1 < m2_min/m1 [q_max,m2_min/m1] fail mistake
q_max < m2_max/m1 < m2_min/m1 < q_min [q_max,q_min] fail mistake
q_max < m2_max/m1 < q_min < m2_min/m1 [q_max,m2_min/m1] fail mistake

Depending on the drawn m1, different cases are in action for m2 at the same time, e.g. m1=10: 0.35/10<0.05<2/10<1.0 vs. m1=50: 0.35/50<2/50<0.05<1.0.

mkruckow commented 1 month ago

Please don't forget to add the new limit to the population_params_default.ini, otherwise a user won't know where the limits should be put and how they have to be called.

mkruckow commented 1 month ago

The more I'm thinking about introducing mass ratio limits, it feels that it gets less intuitive for a user that he needs to care about both. I guess, several users will assume, that either secondary-mass limits or mass-ratio limits will be applied, but not both at the same time. May explaining more, why we currently have a mass ratio limit: only because our HMS-HMS gird has this limit and we don't go to detached there when outside the range, because we don't support currently and HMS-HMS_RLO grid, which we could get from our HMS-HMS data. But best would be to add a lower mass ratio slice to the HMS-HMS gird.