Skip to content

Stellar dynamics

mlpoppyns.simulator.stellar_dynamics.coordinate_conversions

Conversions between different coordinate systems and related issues.

Authors:

    Vanessa Graber (graber@ice.csic.es)
    Michele Ronchi (ronchi@ice.csic.es)

check_radial_coordinate(r)

Check that the distance from the origin is not negative.

Parameters:

Name Type Description Default
r ndarray

Distance from the origin in units of length.

 required
Source code in mlpoppyns/simulator/stellar_dynamics/coordinate_conversions.py 
20
21
22
23
24
25
26
27
28
def check_radial_coordinate(r: np.ndarray) -> None:
    """
    Check that the distance from the origin is not negative.

    Args:
        r (np.ndarray): Distance from the origin in units of length.
    """
    if np.any(r < 0):
        raise ValueError("One of the radial coordinates is out of range.")

convert_cylindrical_to_all_coordinates(dict_pop_dyn)

Converts a population's positions and velocities from galactocentric cylindrical (polar) coordinates to ICRS, and galactic coordinate frames.

The resulting sky coordinates, distances, proper motions and line of sight velocity are added to the input dictionary.

Parameters:

Name Type Description Default
dict_pop_dyn dict

A dictionary containing the initial polar coordinates and velocities of a population of neutron stars.

 required

Returns:

Type Description
dict

The input dictionary dict_pop_dyn with the following keys added:

  • 'ra' (np.ndarray): Right Ascension [deg].
  • 'dec' (np.ndarray): Declination [deg].
  • 'dist' (np.ndarray): Heliocentric distance in ICRS frame [kpc].
  • 'pm_ra' (np.ndarray): Proper motion in RA [mas/yr].
  • 'pm_dec' (np.ndarray): Proper motion in Dec [mas/yr].
  • 'l_gal' (np.ndarray): Galactic longitude [deg].
  • 'b_gal' (np.ndarray): Galactic latitude [deg].
  • 'pm_l' (np.ndarray): Proper motion in Galactic longitude [mas/yr].
  • 'pm_b' (np.ndarray): Proper motion in Galactic latitude [mas/yr].
Source code in mlpoppyns/simulator/stellar_dynamics/coordinate_conversions.py 
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
def convert_cylindrical_to_all_coordinates(dict_pop_dyn: dict) -> dict:
    """
    Converts a population's positions and velocities from galactocentric cylindrical
    (polar) coordinates to ICRS, and galactic coordinate frames.

    The resulting sky coordinates, distances, proper motions and line of sight velocity
    are added to the input dictionary.

    Args:
        dict_pop_dyn (dict): A dictionary containing the initial polar coordinates and velocities
            of a population of neutron stars.

    Returns:
        (dict): The input dictionary `dict_pop_dyn` with the following keys added:

            - 'ra' (np.ndarray): Right Ascension [deg].
            - 'dec' (np.ndarray): Declination [deg].
            - 'dist' (np.ndarray): Heliocentric distance in ICRS frame [kpc].
            - 'pm_ra' (np.ndarray): Proper motion in RA [mas/yr].
            - 'pm_dec' (np.ndarray): Proper motion in Dec [mas/yr].
            - 'l_gal' (np.ndarray): Galactic longitude [deg].
            - 'b_gal' (np.ndarray): Galactic latitude [deg].
            - 'pm_l' (np.ndarray): Proper motion in Galactic longitude [mas/yr].
            - 'pm_b' (np.ndarray): Proper motion in Galactic latitude [mas/yr].
    """

    r = dict_pop_dyn["r"]
    phi = dict_pop_dyn["phi"]
    z = dict_pop_dyn["z"]
    v_r = dict_pop_dyn["v_r"]
    v_phi = dict_pop_dyn["v_phi"]
    v_z = dict_pop_dyn["v_z"]

    # Convert from polar coordinates to Cartesian coordinates.
    x, y = polar_to_cartesian(r, phi)

    # Convert velocity components from galactocentric cylindrical coordinates
    # to galactocentric Cartesian coordinates.
    (
        v_x,
        v_y,
        v_z,
    ) = speed_cylindrical_to_cartesian(v_r, v_phi, v_z, phi)

    # Convert galactocentric coordinates and velocities into ICRS and galactic frame.
    (
        ra,
        dec,
        sun_dist_icrs,
        pm_ra,
        pm_dec,
        v_ls_icrs,
    ) = galactocentric_to_icrs(x, y, z, v_x, v_y, v_z)

    (
        l_gal,
        b_gal,
        sun_dist_gal,
        pm_l,
        pm_b,
        v_ls_gal,
    ) = galactocentric_to_galactic(x, y, z, v_x, v_y, v_z)

    dict_pop_dyn["ra"] = ra
    dict_pop_dyn["dec"] = dec
    dict_pop_dyn["dist"] = sun_dist_icrs
    dict_pop_dyn["pm_ra"] = pm_ra
    dict_pop_dyn["pm_dec"] = pm_dec
    dict_pop_dyn["l"] = l_gal
    dict_pop_dyn["b"] = b_gal
    dict_pop_dyn["v_ls"] = v_ls_icrs

    return dict_pop_dyn

galactocentric_to_galactic(x, y, z, v_x, v_y, v_z)

Calculating the galactic longitude l, galactic latitude b, distance, pm_l, pm_b proper motion components and line of sight velocity v_ls from the galactocentric coordinates x, y, z and velocity v_x, v_y, v_z. The galactocentric coordinates refer to the galactocentric reference frame defined as a right-handed reference frame with the Sun located at the coordinate point (x = 0 kpc, y = 8.3 kpc, z = 0.02 kpc). We use the astropy.coordinates package that allows automatic conversions between coordinate systems. l = 0 deg, b = 0 deg corresponds to the location of the Galactic center, l increase anticlockwise (in the direction of galactic rotation as seen from the Sun) and ranges in the interval [-180, 180] deg while b is in the range [-90, 90] deg.

Parameters:

Name Type Description Default
x ndarray

x coordinate in [kpc] in the galactocentric reference frame.

 required
y ndarray

y coordinate in [kpc] in the galactocentric reference frame.

 required
z ndarray

z coordinate in [kpc] in the galactocentric reference frame.

 required
v_x ndarray

x component of the velocity in [km/s] in galactocentric reference frame.

 required
v_y ndarray

y component of the velocity in [km/s] in galactocentric reference frame.

 required
v_z ndarray

z component of the velocity in [km/s] in galactocentric reference frame.

 required

Returns:

Type Description
Tuple[ndarray, ndarray, ndarray, ndarray, ndarray, ndarray]

Galactic longitude l, Galactic latitude b in [deg], distance in [kpc], pm_l, pm_b proper motion components in [mas/yr] and line of sight velocity v_ls in [km/s].
Source code in mlpoppyns/simulator/stellar_dynamics/coordinate_conversions.py 
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
def galactocentric_to_galactic(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray,
    v_x: np.ndarray,
    v_y: np.ndarray,
    v_z: np.ndarray,
) -> Tuple[
    np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray
]:
    """
    Calculating the galactic longitude l, galactic latitude b, distance, pm_l, pm_b proper motion components and line
    of sight velocity v_ls from the galactocentric coordinates x, y, z and velocity v_x, v_y, v_z.
    The galactocentric coordinates refer to the galactocentric reference frame defined as a right-handed reference
    frame with the Sun located at the coordinate point (x = 0 kpc, y = 8.3 kpc, z = 0.02 kpc).
    We use the astropy.coordinates package that allows automatic conversions between coordinate systems. l = 0 deg,
    b = 0 deg corresponds to the location of the Galactic center, l increase anticlockwise (in the direction of
    galactic rotation as seen from the Sun) and ranges in the interval [-180, 180] deg while b is in the range
    [-90, 90] deg.

    Args:
        x (np.ndarray): x coordinate in [kpc] in the galactocentric reference frame.
        y (np.ndarray): y coordinate in [kpc] in the galactocentric reference frame.
        z (np.ndarray): z coordinate in [kpc] in the galactocentric reference frame.
        v_x (np.ndarray): x component of the velocity in [km/s] in galactocentric reference frame.
        v_y (np.ndarray): y component of the velocity in [km/s] in galactocentric reference frame.
        v_z (np.ndarray): z component of the velocity in [km/s] in galactocentric reference frame.

    Returns:
        (Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]): Galactic longitude l,
            Galactic latitude b in [deg], distance in [kpc], pm_l, pm_b proper motion components in [mas/yr]
            and line of sight velocity v_ls in [km/s].
    """

    # The galactocentric reference frame used for the input is a right-handed
    # reference frame with the Sun located at the coordinate point (x = 0 kpc,
    # y = 8.3 kpc, z = 0.02 kpc).
    # The module astropy.coordinates.Galactocentric deals with galactocentric
    # coordinates but it is defined with the x, y, axes rotated by 90 degrees
    # clockwise with respect to the galactocentric reference frame used in the simulation.
    # In this new frame the position of the Sun is (x = -8.3 kpc, y = 0 kpc, z = 0.02
    # kpc). We therefore need to convert the galactocentric coordinates we used
    # into the galactocentric frame defined in astropy. To do that we
    # apply the transformation (x_gal -> y_astropy, y_gal -> -x_astropy, z_gal -> z_astropy).

    # Set the astropy galactocentric frame with the parameter
    # values from astropy version 4.0.
    _ = galactocentric_frame_defaults.set("v4.0")

    # Create an object containing the coordinates using the class
    # coordinates.Galactocentric from astropy.
    c = coord.Galactocentric(
        x=-y * u.kpc,
        y=x * u.kpc,
        z=z * u.kpc,
        v_x=-v_y * u.km / u.s,
        v_y=v_x * u.km / u.s,
        v_z=-v_z * u.km / u.s,
        z_sun=cfg["z_sun"] * u.kpc,
        galcen_distance=cfg["R_sun"] * u.kpc,
    )

    # Transform from galactocentric to galactic frame.
    galactic = c.transform_to(coord.Galactic())

    # Make the galactic longitude ranging in [-180, 180] deg with the galactic center at (0, 0) deg.
    l_gal = galactic.l.value
    l_gal[(l_gal > 180.0) & (l_gal <= 360.0)] = (
        l_gal[(l_gal > 180.0) & (l_gal <= 360.0)] - 360.0
    )
    b_gal = galactic.b.value
    distance = galactic.distance.value

    pm_l_cosb = galactic.pm_l_cosb.value
    pm_b = galactic.pm_b.value
    radial_velocity = galactic.radial_velocity.value

    return l_gal, b_gal, distance, pm_l_cosb, pm_b, radial_velocity

galactocentric_to_icrs(x, y, z, v_x, v_y, v_z)

Calculating the ICRS (International Celestial Reference Frame) coordinates RA, DEC, distance and proper velocities v_RA, v_DEC, v_ls from galactocentric spatial coordinates and velocities x, y, z, v_x, v_y and v_z. This galactocentric coordinates refers to the galactocentric reference frame used in the simulation defined as a right-handed reference frame with the Sun located at the coordinate point (x = 0 kpc, y = 8.5 kpc, z = 0.02 kpc). We use the astropy.coordinates package that allows automatic conversions between coordinate systems.

Parameters:

Name Type Description Default
x ndarray

x coordinate in [kpc] in galactocentric reference frame.

 required
y ndarray

y coordinate in [kpc] in galactocentric reference frame.

 required
z ndarray

z coordinate in [kpc] in galactocentric reference frame.

 required
v_x ndarray

x velocity component in [km/s] in galactocentric reference frame.

 required
v_y ndarray

y velocity component in [km/s] in galactocentric reference frame.

 required
v_z ndarray

z velocity component in [km/s] in galactocentric reference frame.

 required

Returns:

Type Description
Tuple[ndarray, ndarray, ndarray, ndarray, ndarray, ndarray]

RA, DEC coordinates in [deg], distance from the ICRS origin in [kpc], proper motion pm_RA, pm_DEC components in [mas/yr] in the ICRS reference frame and the line of sight velocity in [km/s].
Source code in mlpoppyns/simulator/stellar_dynamics/coordinate_conversions.py 
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
def galactocentric_to_icrs(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray,
    v_x: np.ndarray,
    v_y: np.ndarray,
    v_z: np.ndarray,
) -> Tuple[
    np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray
]:
    """
    Calculating the ICRS (International Celestial Reference Frame) coordinates RA, DEC, distance and proper velocities
    v_RA, v_DEC, v_ls from galactocentric spatial coordinates and velocities x, y, z, v_x, v_y and v_z. This
    galactocentric coordinates refers to the galactocentric reference frame used in the simulation defined as a
    right-handed reference frame with the Sun located at the coordinate point (x = 0 kpc, y = 8.5 kpc, z = 0.02 kpc).
    We use the astropy.coordinates package that allows automatic conversions between coordinate systems.

    Args:
        x (np.ndarray): x coordinate in [kpc] in galactocentric reference frame.
        y (np.ndarray): y coordinate in [kpc] in galactocentric reference frame.
        z (np.ndarray): z coordinate in [kpc] in galactocentric reference frame.
        v_x (np.ndarray): x velocity component in [km/s] in galactocentric reference frame.
        v_y (np.ndarray): y velocity component in [km/s] in galactocentric reference frame.
        v_z (np.ndarray): z velocity component in [km/s] in galactocentric reference frame.

    Returns:
        (Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]): RA, DEC coordinates in [deg],
            distance from the ICRS origin in [kpc], proper motion pm_RA, pm_DEC components
            in [mas/yr] in the ICRS reference frame and the line of sight velocity in [km/s].
    """

    # Set the astropy galactocentric frame with the parameter
    # values from astropy version 4.0.
    _ = galactocentric_frame_defaults.set("v4.0")

    # The galactocentric reference frame used in the simulation is a right-handed
    # reference frame with the Sun located at the coordinate point (x = 0 kpc,
    # y = 8.5 kpc, z = 0.02 kpc).
    # The module astropy.coordinates.Galactocentric deals with galactocentric
    # coordinates but it is defined with the x, y, axes rotated of 90 degrees
    # clockwise respect to the galactocentric reference frame used in the simulation.
    # In this new frame the position of the Sun is (x = -8.5 kpc, y = 0 kpc, z = 0.02
    # kpc). We therefore need to convert the galactocentric coordinates we used in
    # the simulation into the galactocentric frame defined in astropy. To do that we
    # apply the transformation (x -> y_gal, y -> -x_gal, z -> z_gal).

    x_gal = -y
    y_gal = x
    z_gal = z

    v_x_gal = -v_y
    v_y_gal = v_x
    v_z_gal = v_z

    # Create an object containing the coordinates using the class
    # coordinates.Galactocentric from astropy.
    gc_coord = coord.Galactocentric(
        x=x_gal * u.kpc,
        y=y_gal * u.kpc,
        z=z_gal * u.kpc,
        v_x=v_x_gal * (u.km / u.s),
        v_y=v_y_gal * (u.km / u.s),
        v_z=v_z_gal * (u.km / u.s),
        z_sun=cfg["z_sun"] * u.kpc,
        galcen_distance=cfg["R_sun"] * u.kpc,
    )

    # Transform from galactocentric to ICRS frame.
    icrs_coord = gc_coord.transform_to(coord.ICRS())

    # Determine RA and DEC in [deg] in the ranges [0, 360] and [-90, 90],
    # respectively, and drop the units to the astropy objects by taking only the values.
    ra = icrs_coord.ra.degree
    dec = icrs_coord.dec.degree
    sun_dist = icrs_coord.distance.value
    pm_ra = icrs_coord.pm_ra_cosdec.value
    pm_dec = icrs_coord.pm_dec.value
    v_ls = icrs_coord.radial_velocity.value

    return ra, dec, sun_dist, pm_ra, pm_dec, v_ls

polar_to_cartesian(r, phi)

Calculating the Cartesian x and y coordinates from plane polar r and phi.

Parameters:

Name Type Description Default
r ndarray

Radial component (magnitude of the vector) in plane polar coordinates, r>0.

 required
phi ndarray

Angular coordinate, [0, 2*pi].

 required

Returns:

Type Description
Tuple[ndarray, ndarray]

x and y coordinates in a Cartesian system.
Source code in mlpoppyns/simulator/stellar_dynamics/coordinate_conversions.py 
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
def polar_to_cartesian(
    r: np.ndarray, phi: np.ndarray
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Calculating the Cartesian x and y coordinates from plane polar r and phi.

    Args:
        r (np.ndarray): Radial component (magnitude of the vector) in plane polar coordinates, r>0.
        phi (np.ndarray): Angular coordinate, [0, 2*pi].

    Returns:
        (Tuple[np.ndarray, np.ndarray]): x and y coordinates in a Cartesian system.
    """

    x = r * np.cos(phi)
    y = r * np.sin(phi)

    return x, y

speed_cylindrical_to_cartesian(v_r, v_phi, v_z, phi)

Calculating the galactocentric Cartesian v_x, v_y and v_z velocity components from cylindrical galactocentric components v_r, v_phi and v_z.

Parameters:

Name Type Description Default
v_r ndarray

Radial velocity component in a cylindrical galactocentric frame.

 required
v_phi ndarray

Azimuthal velocity component in a cylindrical galactocentric frame.

 required
v_z ndarray

z velocity component in a cylindrical galactocentric frame.

 required
phi ndarray

Azimuthal angle [0, 2*pi] in cylindrical coordinates.

 required

Returns:

Type Description
Tuple[ndarray, ndarray, ndarray]

v_x, v_y and v_z velocity components in a Cartesian galactocentric frame.
Source code in mlpoppyns/simulator/stellar_dynamics/coordinate_conversions.py 
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
def speed_cylindrical_to_cartesian(
    v_r: np.ndarray, v_phi: np.ndarray, v_z: np.ndarray, phi: np.ndarray
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Calculating the galactocentric Cartesian v_x, v_y and v_z velocity components
    from cylindrical galactocentric components v_r, v_phi and v_z.

    Args:
        v_r (np.ndarray): Radial velocity component in a cylindrical galactocentric frame.
        v_phi (np.ndarray): Azimuthal velocity component in a cylindrical galactocentric frame.
        v_z (np.ndarray): z velocity component in a cylindrical galactocentric frame.
        phi (np.ndarray): Azimuthal angle [0, 2*pi] in cylindrical coordinates.

    Returns:
        (Tuple[np.ndarray, np.ndarray, np.ndarray]): v_x, v_y and v_z velocity components in a Cartesian
            galactocentric frame.
    """

    v_x = v_r * np.cos(phi) - v_phi * np.sin(phi)
    v_y = v_r * np.sin(phi) + v_phi * np.cos(phi)

    return v_x, v_y, v_z

spherical_to_cartesian(r, theta, psi)

Calculating the Cartesian x, y and z coordinates from spherical coordinates r, theta and psi.

Parameters:

Name Type Description Default
r ndarray

Radial component (magnitude of the vector) in spherical coordinates, r>0.

 required
theta ndarray

Polar angle, [0, pi].

 required
psi ndarray

Azimuthal angle, [0, 2*pi].

 required

Returns:

Type Description
Tuple[ndarray, ndarray, ndarray]

x, y and z coordinates in a Cartesian system.
Source code in mlpoppyns/simulator/stellar_dynamics/coordinate_conversions.py 
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def spherical_to_cartesian(
    r: np.ndarray, theta: np.ndarray, psi: np.ndarray
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Calculating the Cartesian x, y and z coordinates from spherical coordinates
    r, theta and psi.

    Args:
        r (np.ndarray): Radial component (magnitude of the vector) in spherical coordinates, r>0.
        theta (np.ndarray): Polar angle, [0, pi].
        psi (np.ndarray): Azimuthal angle, [0, 2*pi].

    Returns:
        (Tuple[np.ndarray, np.ndarray, np.ndarray]): x, y and z coordinates in a Cartesian system.
    """

    x = r * np.sin(theta) * np.cos(psi)
    y = r * np.sin(theta) * np.sin(psi)
    z = r * np.cos(theta)

    return x, y, z

mlpoppyns.simulator.stellar_dynamics.dynamical_evolution

Dynamical evolution of the neutron stars in the galactic potential.

We solve the system of dynamical differential equations in cylindrical coordinates, using a galactocentric reference frame. Here we are using the scipy.integrate.odeint package which uses the method 'LSODA' (Adams/BDF method with automatic stiffness detection and switching) from the Fortran library ODEPACK.

We improve performance with Numba, which allows just-in-time (JIT) compilation.

Authors:

Vanessa Graber (graber@ice.csic.es)
Michele Ronchi (ronchi@ice.csic.es)

check_angular_momentum_energy_conservation(dict_pop_initial_dyn, dict_pop_final_dyn, logger)

This function computes the total energy and angular momentum (L_z) of a stellar population before and after a dynamical evolution. It reports the percentage variation of each quantity to assess conservation.

Parameters:

Name Type Description Default
dict_pop_initial_dyn dict

Dictionary containing the initial dynamical properties of the population.

 required
dict_pop_final_dyn dict

Dictionary containing the final dynamical properties of the population.

 required
logger Logger

Logger instance used to output conservation diagnostics (energy and angular momentum variations).

 required
Source code in mlpoppyns/simulator/stellar_dynamics/dynamical_evolution.py 
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
def check_angular_momentum_energy_conservation(
    dict_pop_initial_dyn: dict,
    dict_pop_final_dyn: dict,
    logger: logging.Logger,
) -> None:
    """
    This function computes the total energy and angular momentum (L_z) of a
    stellar population before and after a dynamical evolution. It reports the
    percentage variation of each quantity to assess conservation.

    Args:
        dict_pop_initial_dyn (dict): Dictionary containing the initial dynamical
            properties of the population.
        dict_pop_final_dyn (dict): Dictionary containing the final dynamical
            properties of the population.
        logger (logging.Logger): Logger instance used to output conservation
            diagnostics (energy and angular momentum variations).
    """

    r_initial = dict_pop_initial_dyn["r"]
    z_initial = dict_pop_initial_dyn["z"]
    v_r_initial = dict_pop_initial_dyn["v_r"]
    v_phi_initial = dict_pop_initial_dyn["v_phi"]
    v_z_initial = dict_pop_initial_dyn["v_z"]

    # Compute the magnitude of the initial velocity vector for each star.
    v_initial = (
        np.sqrt(v_r_initial**2 + v_phi_initial**2 + v_z_initial**2)
        * const.KPC_TO_KM
        / const.YR_TO_S
    )

    # Compute the total initial energy of the system.
    total_energy_initial = gm.galactic_model.total_energy(
        v_initial, r_initial, z_initial
    )

    # Compute the initial z-component of the total angular momentum of the system.
    L_z_initial = gm.galactic_model.total_angular_momentum_z(
        v_phi_initial * const.KPC_TO_KM / const.YR_TO_S, r_initial
    )

    r_final = dict_pop_final_dyn["r"]
    z_final = dict_pop_final_dyn["z"]
    v_r_final = dict_pop_final_dyn["v_r"]
    v_phi_final = dict_pop_final_dyn["v_phi"]
    v_z_final = dict_pop_final_dyn["v_z"]

    # Compute the magnitude of the final velocity vector for each star.
    v_final = np.sqrt(v_r_final**2 + v_phi_final**2 + v_z_final**2)

    # Compute the total energy of the system after the dynamical evolution.
    total_energy_final = gm.galactic_model.total_energy(
        v_final, r_final, z_final
    )

    # Compute the percentage variation in total energy during the simulation
    # with respect to the initial total energy.
    delta_energy_percentage = (
        (total_energy_final - total_energy_initial)
        / total_energy_initial
        * 100.0
    )

    logger.info(
        f"Percentage variation of total energy of the system: {delta_energy_percentage} %"
    )

    # Compute the final z-component of the total angular momentum of the system.
    L_z_final = gm.galactic_model.total_angular_momentum_z(
        v_phi_final, r_final
    )

    # Compute the percentage variation in total energy during the simulation
    # with respect to the initial total energy.
    delta_Lz_percentage = (L_z_final - L_z_initial) / L_z_initial * 100.0

    logger.info(
        f"Percentage variation of z-component of total angular momentum of the system: {delta_Lz_percentage} %"
    )

dynamical_eq_system(t, initial_cond, galactic_model)

System of dynamical equations to solve to determine the orbits of the neutron stars in the galactic potential. The differential equation are written in cylindrical galactocentric coordinates (r, phi, z).

Parameters:

Name Type Description Default
t float

Unused time variable, required for the integration below.

 required
initial_cond ndarray

Array of 6 components defining the initial conditions in cylindrical coordinates (r0, phi0, z0, v_r0, omega0, v_z0) with the following units ([kpc], [rad], [kpc], [kpc/yr], [rad/yr], [kpc/yr]).

 required
galactic_model GalaxyModelBase

A galactic model to calculate the needed potential.

 required

Returns:

Type Description
ndarray

Array of 6 values of the first order and second order derivatives at each time step.
Source code in mlpoppyns/simulator/stellar_dynamics/dynamical_evolution.py 
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
@jit(
    float64[:](
        float64,
        float64[:],
        gm.galactic_model._numba_type_.class_type.instance_type,
    )
)
def dynamical_eq_system(
    t: float, initial_cond: np.ndarray, galactic_model: gm.GalaxyModelBase
) -> np.ndarray:
    """
    System of dynamical equations to solve to determine the orbits of the neutron stars in the galactic potential.
    The differential equation are written in cylindrical galactocentric coordinates (r, phi, z).

    Args:
        t (float): Unused time variable, required for the integration below.
        initial_cond (np.ndarray): Array of 6 components defining the initial conditions in cylindrical coordinates
            (r0, phi0, z0, v_r0, omega0, v_z0) with the following units ([kpc], [rad], [kpc], [kpc/yr], [rad/yr], [kpc/yr]).
        galactic_model (gm.GalaxyModelBase): A galactic model to calculate the needed potential.

    Returns:
         (np.ndarray): Array of 6 values of the first order and second order derivatives at each time step.

    """

    r = initial_cond[0]
    z = initial_cond[2]

    gradient_mw_pot = galactic_model.cylind_coord_gradient_mw_potential(r, z)

    # First derivatives.
    dr_dt = initial_cond[3]
    dphi_dt = initial_cond[4]
    dz_dt = initial_cond[5]

    # Second derivatives.
    d2r_dt2 = r * dphi_dt**2 - gradient_mw_pot[0]
    d2phi_dt2 = -2 * dr_dt * dphi_dt / r - gradient_mw_pot[1]
    d2z_dt2 = -gradient_mw_pot[2]

    derivatives = np.array(
        [dr_dt, dphi_dt, dz_dt, d2r_dt2, d2phi_dt2, d2z_dt2]
    )

    return derivatives

dynamical_evolution(initial_cond, t_age)

Performing the dynamical evolution of the neutron star population for a given galactic potential, starting from a set of initial conditions.

Parameters:

Name Type Description Default
initial_cond ndarray

Array of 6 components defining the initial conditions in cylindrical coordinates (r0, phi0, z0, v_r0, omega0, v_z0) with the following units ([kpc], [rad], [kpc], [kpc/yr], [rad/yr], [kpc/yr]).

 required
t_age ndarray

Array of neutron star ages in [yr].

 required

Returns:

Type Description
Tuple[ndarray, dict]

Tuple consisting of a two-dimensional array of shape (NS_number, 6) defining the neutron stars' final positions r [kpc], phi [rad], z [kpc] and velocities in [kpc/yr] in cylindrical coordinates and a dictionary containing the time evolution of these quantities for each neutron star (if the option to save the time evolution is enabled).
Source code in mlpoppyns/simulator/stellar_dynamics/dynamical_evolution.py 
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
def dynamical_evolution(
    initial_cond: np.ndarray, t_age: np.ndarray
) -> Tuple[np.ndarray, dict]:
    """
    Performing the dynamical evolution of the neutron star population for a given
    galactic potential, starting from a set of initial conditions.

    Args:
        initial_cond (np.ndarray): Array of 6 components defining the initial
            conditions in cylindrical coordinates (r0, phi0, z0, v_r0, omega0, v_z0)
            with the following units ([kpc], [rad], [kpc], [kpc/yr], [rad/yr], [kpc/yr]).
        t_age (np.ndarray): Array of neutron star ages in [yr].

    Returns:
        (Tuple[np.ndarray, dict]): Tuple consisting of a two-dimensional array of shape (NS_number, 6)
            defining the neutron stars' final positions r [kpc], phi [rad], z [kpc] and velocities
            in [kpc/yr] in cylindrical coordinates and a dictionary containing the time evolution of
            these quantities for each neutron star (if the option to save the time evolution is enabled).
    """

    # Save the number of simulated neutron stars, which is flexible depending on whether
    # they are simulated all at once or one by one.
    n = len(t_age)

    # Initialization of a dictionary that will contain the evolution in time of
    # positions and velocities.
    evolution_dictionary = {}

    # Initialize the arrays that will contain the final positions
    # and velocities of the neutron stars.
    r_final = np.zeros(n)
    phi_final = np.zeros(n)
    z_final = np.zeros(n)
    v_r_final = np.zeros(n)
    v_phi_final = np.zeros(n)
    v_z_final = np.zeros(n)

    for i in range(n):
        # Linear time grid in years over which the dynamical evolution is performed;
        # each star's position and velocity is evolved for a time equal to its age.
        time_grid = np.append(
            np.arange(0.0, t_age[i], cfg["dyn_time_step"]),
            t_age[i],
        )

        # Save the odeint output which is a two-dimensional array of
        # shape (len(time_grid), 6).
        # We set tfirst=True to unify the structure of the input ODEs in order to be able
        # to compare different scipy functions to solve the ODEs.
        evol_output = np.array(
            odeint(
                dynamical_eq_system,
                y0=initial_cond[i],
                t=time_grid,
                args=(gm.galactic_model,),
                tfirst=True,
            )
        )

        if cfg["save_dyn_evolution"]:
            v_r_evol = evol_output[:, 3] * const.KPC_TO_KM / const.YR_TO_S
            v_phi_evol = (
                evol_output[:, 0]
                * evol_output[:, 4]
                * const.KPC_TO_KM
                / const.YR_TO_S
            )
            v_z_evol = evol_output[:, 5] * const.KPC_TO_KM / const.YR_TO_S

            # Save the evolution output of the i-th neutron star in a dictionary.
            evolution = {
                i: {
                    "t": time_grid.tolist(),
                    "r(t)": evol_output[:, 0].tolist(),
                    "phi(t)": evol_output[:, 1].tolist(),
                    "z(t)": evol_output[:, 2].tolist(),
                    "v_r(t)": v_r_evol.tolist(),
                    "v_phi(t)": v_phi_evol.tolist(),
                    "v_z(t)": v_z_evol.tolist(),
                }
            }
            # Update the dictionary containing the evolution information of all the neutron stars.
            evolution_dictionary = {**evolution_dictionary, **evolution}

        # Save the final position and velocity.
        # Note: We save directly the v_phi velocity component and not the angular velocity omega.
        r_final[i] = evol_output[-1, 0]
        phi_final[i] = evol_output[-1, 1]
        z_final[i] = evol_output[-1, 2]
        v_r_final[i] = evol_output[-1, 3]
        omega_final = evol_output[-1, 4]
        v_phi_final[i] = omega_final * r_final[i]
        v_z_final[i] = evol_output[-1, 5]

    final_population = np.array(
        [r_final, phi_final, z_final, v_r_final, v_phi_final, v_z_final]
    ).T

    return final_population, evolution_dictionary

evolve_population_dyn(dict_pop_initial_dyn, output_path)

Evolve the dynamical properties of a neutron star population over time based on initial conditions.

Parameters:

Name Type Description Default
dict_pop_initial_dyn dict

Dictionary containing initial magneto-rotational properties of the population.

 required
output_path Path

The path where the evolution data will be saved if enabled in the configuration.

 required

Returns:

Type Description
dict

A dictionary containing the properties of the evolved neutron star population.
Source code in mlpoppyns/simulator/stellar_dynamics/dynamical_evolution.py 
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
def evolve_population_dyn(
    dict_pop_initial_dyn: dict,
    output_path: pathlib.Path,
) -> dict:
    """
    Evolve the dynamical properties of a neutron star population over time based on initial conditions.

    Args:
        dict_pop_initial_dyn (dict): Dictionary containing initial magneto-rotational properties of the population.
        output_path (pathlib.Path): The path where the evolution data will be saved if enabled in the configuration.

    Returns:
        (dict): A dictionary containing the properties of the evolved neutron star population.
    """
    age = dict_pop_initial_dyn["age"]
    r_initial = dict_pop_initial_dyn["r"]
    phi_initial = dict_pop_initial_dyn["phi"]
    z_initial = dict_pop_initial_dyn["z"]
    v_r_initial = dict_pop_initial_dyn["v_r"]
    v_phi_initial = dict_pop_initial_dyn["v_phi"]
    v_z_initial = dict_pop_initial_dyn["v_z"]

    omega_initial = v_phi_initial / r_initial

    # Define the initial conditions for the dynamical evolution.
    initial_cond = np.array(
        [
            r_initial,
            phi_initial,
            z_initial,
            v_r_initial,
            omega_initial,
            v_z_initial,
        ]
    ).T

    # Determine the evolved positions and velocities.
    final_population, dyn_evol_dict = dynamical_evolution(
        initial_cond,
        age,
    )

    r_final = final_population[:, 0]
    phi_final = final_population[:, 1]
    z_final = final_population[:, 2]
    v_r_final = final_population[:, 3]
    v_phi_final = final_population[:, 4]
    v_z_final = final_population[:, 5]

    # Convert velocities from [kpc/yr] into [km/s].
    v_r_final = v_r_final * const.KPC_TO_KM / const.YR_TO_S
    v_phi_final = v_phi_final * const.KPC_TO_KM / const.YR_TO_S
    v_z_final = v_z_final * const.KPC_TO_KM / const.YR_TO_S

    if cfg["save_dyn_evolution"]:
        # Save dictionary containing evolution information to output path in a .json file.
        dyn_evolution_dump_path = pathlib.Path().joinpath(
            output_path, "dyn_evolution.json"
        )

        with open(dyn_evolution_dump_path, "wb") as f:
            f.write(
                orjson.dumps(
                    dict(dyn_evol_dict),
                    option=orjson.OPT_SERIALIZE_NUMPY
                    | orjson.OPT_NON_STR_KEYS
                    | orjson.OPT_SORT_KEYS,
                )
            )

    dictionary_final_pop_dyn = {
        "age": age,
        "r": r_final,
        "phi": phi_final,
        "z": z_final,
        "v_r": v_r_final,
        "v_phi": v_phi_final,
        "v_z": v_z_final,
    }

    return dictionary_final_pop_dyn

initialize_population_dyn(age)

Initialize the dynamical properties of a neutron star population.

Parameters:

Name Type Description Default
age ndarray

An array of ages in [yr] for the neutron stars in the population to be initialized.

 required

Returns:

Type Description
dict

A dictionary containing the initialized dynamical properties of the neutron star population.
Source code in mlpoppyns/simulator/stellar_dynamics/dynamical_evolution.py 
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def initialize_population_dyn(age: np.ndarray) -> dict:
    """
    Initialize the dynamical properties of a neutron star population.

    Args:
        age (np.ndarray): An array of ages in [yr] for the neutron stars in the population to be initialized.

    Returns:
        (dict): A dictionary containing the initialized dynamical properties of the neutron star population.
    """

    # Initialize neutron star population properties.
    pop_initial = ipop.InitialNeutronStarPopulation(NS_number=len(age))

    # Generating initial positions.
    (
        r_initial,
        phi_initial,
        z_initial,
    ) = pop_initial.position(t_age=age)

    # Generating initial velocities by summing the kick
    # velocities at birth and the orbital velocities.
    (
        vk_r,
        vk_phi,
        vk_z,
    ) = pop_initial.kick_velocity()

    v_orb = pop_initial.orbital_velocity(r_initial, z_initial)

    v_r_initial = vk_r
    v_phi_initial = vk_phi + v_orb
    v_z_initial = vk_z

    dictionary_initial_pop_dyn = {
        "age": age,
        "r": r_initial,
        "phi": phi_initial,
        "z": z_initial,
        "v_r": v_r_initial,
        "v_phi": v_phi_initial,
        "v_z": v_z_initial,
        "v_orb": v_orb,
    }

    return dictionary_initial_pop_dyn

mlpoppyns.simulator.stellar_dynamics.galactic_model

Model for the Milky Way gravitational potential.

We consider two different models:

1) gmFK06: A galactic structure as in Faucher-Giguère & Kaspi (2006). Their model consists of three components: a disk-halo, a bulge and a nucleus. The parameters of the model are taken from Table B1 in Kuijken & Gilmore (1989).

2) gmM19: The galaxy model from Marchetti et al. (2019). This is a four-component galactic potential model consisting of a Hernquist bulge and nucleus (Hernquist 1990), a Miyamoto-Nagai disk (Miyamoto & Nagai 1975) and a Navarro-Frenk-White halo (Navarro et al. 1996). The parameters of the model are taken from Table 1 in Marchetti et al. (2019) and are chosen to fit the enclosed mass profile of the Milky Way (Bovy 2015).

To improve performance when evolving the neutron stars' position in the galactic potential (see dynamical_evolution.py), we add Numba's jit decorator to all functions.

Authors:

    Vanessa Graber (graber@ice.csic.es)
    Michele Ronchi (ronchi@ice.csic.es)

GalaxyModelBase

Base class for any galaxy model to ensure that a common interface between all of them is respected. The following abstract methods must be implemented or an error will be raised.

Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
class GalaxyModelBase:
    """
    Base class for any galaxy model to ensure that a common interface between all of
    them is respected. The following abstract methods must be implemented or an error
    will be raised.
    """

    @abc.abstractmethod
    def MW_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
        raise NotImplementedError("Please implement the method MW_potential.")

    @abc.abstractmethod
    def cylind_coord_gradient_mw_potential(
        self, r: float, z: float
    ) -> np.ndarray:
        raise NotImplementedError(
            "Please implement the method cylind_coord_gradient_mw_potential."
        )

    def total_energy(
        self, v: np.ndarray, r: np.ndarray, z: np.ndarray
    ) -> float:
        """
        Value of the total energy of the system, sum of the total kinetic energy and
        the total gravitational potential energy. We assume here that all the stars
        have unit mass.

        Args:
            v (np.ndarray): Array of magnitudes of the speed of the stars in [km/s].
            r (np.ndarray): Array of distances from the galactic axis in [kpc].
            z (np.ndarray): Array of distances from the galactic disk in [kpc].

        Returns:
            (float): Value of the total energy of the system in [erg].
        """
        # Convert speeds into [cm/s].
        v = v * const.KM_TO_CM

        tot_kin_energy = 0.5 * np.sum(v**2.0)
        tot_pot_energy = np.sum(self.MW_potential(r, z))

        tot_energy = tot_kin_energy + tot_pot_energy

        return tot_energy

    def total_angular_momentum_z(
        self,
        v_phi: np.ndarray,
        r: np.ndarray,
    ) -> float:
        """
        Value of the z-component of the total angular momentum of the system, which is a
        conserved quantity in axisymmetric potentials. We assume here that all the stars
        have unit mass.

        Args:
            v_phi (np.ndarray): Array of magnitudes of the orbital speed of the stars in [km/s].
            r (np.ndarray): Array of distances from the galactic axis in [kpc].

        Returns:
            (float): Value of the total energy of the system in [erg].
        """
        # Convert speeds into [cm/s].
        v_phi = v_phi * const.KM_TO_CM
        # Convert distances into [cm].
        r = r * const.KPC_TO_CM

        L_z = float(np.sum(r * v_phi))

        return L_z

total_angular_momentum_z(v_phi, r)

Value of the z-component of the total angular momentum of the system, which is a conserved quantity in axisymmetric potentials. We assume here that all the stars have unit mass.

Parameters:

Name Type Description Default
v_phi ndarray

Array of magnitudes of the orbital speed of the stars in [km/s].

 required
r ndarray

Array of distances from the galactic axis in [kpc].

 required

Returns:

Type Description
float

Value of the total energy of the system in [erg].
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
def total_angular_momentum_z(
    self,
    v_phi: np.ndarray,
    r: np.ndarray,
) -> float:
    """
    Value of the z-component of the total angular momentum of the system, which is a
    conserved quantity in axisymmetric potentials. We assume here that all the stars
    have unit mass.

    Args:
        v_phi (np.ndarray): Array of magnitudes of the orbital speed of the stars in [km/s].
        r (np.ndarray): Array of distances from the galactic axis in [kpc].

    Returns:
        (float): Value of the total energy of the system in [erg].
    """
    # Convert speeds into [cm/s].
    v_phi = v_phi * const.KM_TO_CM
    # Convert distances into [cm].
    r = r * const.KPC_TO_CM

    L_z = float(np.sum(r * v_phi))

    return L_z

total_energy(v, r, z)

Value of the total energy of the system, sum of the total kinetic energy and the total gravitational potential energy. We assume here that all the stars have unit mass.

Parameters:

Name Type Description Default
v ndarray

Array of magnitudes of the speed of the stars in [km/s].

 required
r ndarray

Array of distances from the galactic axis in [kpc].

 required
z ndarray

Array of distances from the galactic disk in [kpc].

 required

Returns:

Type Description
float

Value of the total energy of the system in [erg].
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def total_energy(
    self, v: np.ndarray, r: np.ndarray, z: np.ndarray
) -> float:
    """
    Value of the total energy of the system, sum of the total kinetic energy and
    the total gravitational potential energy. We assume here that all the stars
    have unit mass.

    Args:
        v (np.ndarray): Array of magnitudes of the speed of the stars in [km/s].
        r (np.ndarray): Array of distances from the galactic axis in [kpc].
        z (np.ndarray): Array of distances from the galactic disk in [kpc].

    Returns:
        (float): Value of the total energy of the system in [erg].
    """
    # Convert speeds into [cm/s].
    v = v * const.KM_TO_CM

    tot_kin_energy = 0.5 * np.sum(v**2.0)
    tot_pot_energy = np.sum(self.MW_potential(r, z))

    tot_energy = tot_kin_energy + tot_pot_energy

    return tot_energy

GalaxyModelFK06

Bases: GalaxyModelBase

Galaxy model from Faucher-Giguère & Kaspi (2006). This model consists of a disk-halo component, a bulge component, and a nucleus component. The parameters of the model are taken from Table 1 in Kuijken & Gilmore (1989) (in Faucher-Giguère & Kaspi (2006) the nucleus and bulge are erroneously inverted).

Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
@jitclass(galaxyModelFK06_spec)
class GalaxyModelFK06(GalaxyModelBase):
    """
    Galaxy model from Faucher-Giguère & Kaspi (2006). This model consists of a
    disk-halo component, a bulge component, and a nucleus component. The
    parameters of the model are taken from Table 1 in Kuijken & Gilmore (1989)
    (in Faucher-Giguère & Kaspi (2006) the nucleus and bulge are erroneously
    inverted).
    """

    def __init__(self) -> None:
        # Parameter values from Table B1 in Kuijken & Gilmore (1989).
        self.a_d = 2.4  # Scale length of the disk in [kpc].
        self.h = np.array(
            [0.325, 0.090, 0.125]
        )  # Array of disk components' scale heights in [kpc].
        self.beta = np.array(
            [0.4, 0.5, 0.1]
        )  # Array of weights for the disk components.
        self.M_dh = 1.45e11 * const.M_SUN  # Disk+halo mass in [g].
        self.b_dh = 5.5  # Core radius of the halo component in [kpc].
        self.M_b = 1.0e10 * const.M_SUN  # Bulge mass in [g].
        self.b_b = 1.5  # Core radius of the bulge component in [kpc].
        self.M_n = 9.3e9 * const.M_SUN  # Nucleus mass in [g].
        self.b_n = 0.25  # Core radius of the nucleus component in [kpc].

    def shape_parameter(self, z: float) -> Tuple[float, float]:
        """
        Shape parameter for the disk-halo potential from Carlberg & Innamen
        (1987) and its derivative with respect to the height z from the galactic
        disk. First term in the denominator of eq. (13) in Faucher-Giguère &
        Kaspi (2006).

        Args:
            z (float): Height from the galactic disk in [kpc].

        Returns:
            (Tuple[float, float]): Value of the shape parameter and its derivative with
                respect to z.
        """

        a_d = self.a_d
        beta = self.beta
        h = self.h

        K = (
            a_d
            + beta[0] * np.sqrt(z**2.0 + h[0] ** 2.0)
            + beta[1] * np.sqrt(z**2.0 + h[1] ** 2.0)
            + beta[2] * np.sqrt(z**2.0 + h[2] ** 2.0)
        )

        dK_dz = (
            beta[0] * z / np.sqrt(z**2.0 + h[0] ** 2.0)
            + beta[1] * z / np.sqrt(z**2.0 + h[1] ** 2.0)
            + beta[2] * z / np.sqrt(z**2.0 + h[2] ** 2.0)
        )

        return K, dK_dz

    def dh_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
        """
        A cylindrical disk-halo component gravitational potential as defined in eq. (14) in
        Faucher-Giguère & Kaspi (2006).

        Args:
            r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (np.ndarray): Height from the galactic disk in [kpc].

        Returns:
            (np.ndarray): Value of the disk-halo potential in [erg/g].
        """

        K, _ = self.shape_parameter(z)

        M_dh = self.M_dh
        b_dh = self.b_dh

        pot_dh = (
            -const.G
            * M_dh
            / (np.sqrt(K**2.0 + b_dh**2.0 + r**2.0) * const.KPC_TO_CM)
        )

        return pot_dh

    def b_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
        """
        A spherical bulge component gravitational potential as defined in eq. (15) in
        Faucher-Giguère & Kaspi (2006). Note that we assume this potential to be
        spherical following Carlberg & Innanen (1987), while Faucher-Giguère & Kaspi take r
        to be the cylindrical coordinate, not the spherical one.

        Args:
            r (np.ndarray): Distance in the galactic disk from the galactic centre in [kpc].
            z (np.ndarray): Height from the galactic disk in [kpc].

        Returns:
            (np.ndarray): Value of the bulge potential in [erg/g].
        """
        M_b = self.M_b
        b_b = self.b_b

        pot_b = (
            -const.G
            * M_b
            / (np.sqrt(b_b**2 + r**2 + z**2) * const.KPC_TO_CM)
        )

        return pot_b

    def n_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
        """
        A spherical nucleus component gravitational potential as defined in eq. (15) in
        Faucher-Giguère & Kaspi (2006). Note that we assume this potential to be
        spherical following Carlberg & Innanen (1987), while Faucher-Giguère & Kaspi take r
        to be the cylindrical coordinate, not the spherical one.

        Args:
            r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (np.ndarray): Height from the galactic disk in [kpc].

        Returns:
            (np.ndarray): Value of the nucleus potential.
        """
        M_n = self.M_n
        b_n = self.b_n

        pot_n = (
            -const.G
            * M_n
            / (np.sqrt(b_n**2.0 + r**2.0 + z**2.0) * const.KPC_TO_CM)
        )

        return pot_n

    def MW_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
        """
        Total Milky Way gravitational potential defined in eq. (13) in
        Faucher-Giguère & Kaspi (2006).

        Args:
            r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (np.ndarray): Height from the galactic disk in [kpc].

        Returns:
            (np.ndarray): Value of the Galactic potential in [erg].
        """

        MW_pot = (
            self.dh_potential(r, z)
            + self.b_potential(r, z)
            + self.n_potential(r, z)
        )

        return MW_pot

    def r_z_derivatives_dh_potential(
        self, r: float, z: float
    ) -> Tuple[float, float]:
        """
        Derivatives with respect to r and z of the disk-halo component gravitational
        potential defined in eq. (14) in Faucher-Giguère & Kaspi (2006).

        Args:
            r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (float): Height from the galactic disk in [kpc].

        Returns:
            (Tuple[float, float]): Derivatives with respect to r and z of the disk-halo
                potential.
        """

        M_dh = self.M_dh
        b_dh = self.b_dh

        K, dK_dz = self.shape_parameter(z)
        dpot_dh_dr = (
            const.G_KPC_YR
            * M_dh
            * r
            * (K**2.0 + b_dh**2.0 + r**2.0) ** (-3.0 / 2.0)
        )
        dpot_dh_dz = (
            const.G_KPC_YR
            * M_dh
            * (K**2.0 + b_dh**2.0 + r**2.0) ** (-3.0 / 2.0)
            * K
            * dK_dz
        )

        return dpot_dh_dr, dpot_dh_dz

    def r_z_derivatives_b_potential(
        self, r: float, z: float
    ) -> Tuple[float, float]:
        """
        Derivatives with respect to r and z of the bulge component gravitational potential
        defined in eq. (15) in Faucher-Giguère & Kaspi (2006).

        Args:
            r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (float): Height from the galactic disk in [kpc].

        Returns:
            (Tuple[float, float]): Derivatives with respect to r and z of the bulge potential.
        """

        M_b = self.M_b
        b_b = self.b_b

        dpot_b_dr = (
            const.G_KPC_YR
            * M_b
            * r
            * (b_b**2.0 + r**2.0 + z**2.0) ** (-3.0 / 2.0)
        )

        dpot_b_dz = (
            const.G_KPC_YR
            * M_b
            * z
            * (b_b**2.0 + r**2.0 + z**2.0) ** (-3.0 / 2.0)
        )

        return dpot_b_dr, dpot_b_dz

    def r_z_derivatives_n_potential(
        self, r: float, z: float
    ) -> Tuple[float, float]:
        """
        Derivatives with respect to r and z of the nucleus component gravitational potential
        defined in eq. (15) in Faucher-Giguère & Kaspi (2006).

        Args:
            r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (float): Height from the galactic disk in [kpc].

        Returns:
            (Tuple[float, float]): Derivatives with respect to r and z of the nucleus potential.
        """

        M_n = self.M_n
        b_n = self.b_n

        dpot_n_dr = (
            const.G_KPC_YR
            * M_n
            * r
            * (b_n**2.0 + r**2.0 + z**2.0) ** (-3.0 / 2.0)
        )

        dpot_n_dz = (
            const.G_KPC_YR
            * M_n
            * z
            * (b_n**2.0 + r**2.0 + z**2.0) ** (-3.0 / 2.0)
        )

        return dpot_n_dr, dpot_n_dz

    def cylind_coord_gradient_mw_potential(
        self, r: float, z: float
    ) -> np.ndarray:
        """
        Gradient in cylindrical coordinates of the Milky Way gravitational potential
        defined in eq. (13) in Faucher-Giguère & Kaspi (2006).

        Args:
            r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (float): Height from the galactic disk in [kpc].

        Returns:
            (np.ndarray): Gradient of the galactic potential in cylindrical
                coordinates.
        """

        dpot_dh_dr, dpot_dh_dz = self.r_z_derivatives_dh_potential(r, z)
        dpot_b_dr, dpot_b_dz = self.r_z_derivatives_b_potential(r, z)
        dpot_n_dr, dpot_n_dz = self.r_z_derivatives_n_potential(r, z)

        dpot_mw_dr = dpot_dh_dr + dpot_b_dr + dpot_n_dr
        dpot_mw_dphi = 0.0
        dpot_mw_dz = dpot_dh_dz + dpot_b_dz + dpot_n_dz

        pot_mw_gradient = np.array([dpot_mw_dr, dpot_mw_dphi, dpot_mw_dz])

        return pot_mw_gradient

MW_potential(r, z)

Total Milky Way gravitational potential defined in eq. (13) in Faucher-Giguère & Kaspi (2006).

Parameters:

Name Type Description Default
r ndarray

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z ndarray

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
ndarray

Value of the Galactic potential in [erg].
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
def MW_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
    """
    Total Milky Way gravitational potential defined in eq. (13) in
    Faucher-Giguère & Kaspi (2006).

    Args:
        r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (np.ndarray): Height from the galactic disk in [kpc].

    Returns:
        (np.ndarray): Value of the Galactic potential in [erg].
    """

    MW_pot = (
        self.dh_potential(r, z)
        + self.b_potential(r, z)
        + self.n_potential(r, z)
    )

    return MW_pot

b_potential(r, z)

A spherical bulge component gravitational potential as defined in eq. (15) in Faucher-Giguère & Kaspi (2006). Note that we assume this potential to be spherical following Carlberg & Innanen (1987), while Faucher-Giguère & Kaspi take r to be the cylindrical coordinate, not the spherical one.

Parameters:

Name Type Description Default
r ndarray

Distance in the galactic disk from the galactic centre in [kpc].

 required
z ndarray

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
ndarray

Value of the bulge potential in [erg/g].
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
def b_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
    """
    A spherical bulge component gravitational potential as defined in eq. (15) in
    Faucher-Giguère & Kaspi (2006). Note that we assume this potential to be
    spherical following Carlberg & Innanen (1987), while Faucher-Giguère & Kaspi take r
    to be the cylindrical coordinate, not the spherical one.

    Args:
        r (np.ndarray): Distance in the galactic disk from the galactic centre in [kpc].
        z (np.ndarray): Height from the galactic disk in [kpc].

    Returns:
        (np.ndarray): Value of the bulge potential in [erg/g].
    """
    M_b = self.M_b
    b_b = self.b_b

    pot_b = (
        -const.G
        * M_b
        / (np.sqrt(b_b**2 + r**2 + z**2) * const.KPC_TO_CM)
    )

    return pot_b

cylind_coord_gradient_mw_potential(r, z)

Gradient in cylindrical coordinates of the Milky Way gravitational potential defined in eq. (13) in Faucher-Giguère & Kaspi (2006).

Parameters:

Name Type Description Default
r float

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z float

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
ndarray

Gradient of the galactic potential in cylindrical coordinates.
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
def cylind_coord_gradient_mw_potential(
    self, r: float, z: float
) -> np.ndarray:
    """
    Gradient in cylindrical coordinates of the Milky Way gravitational potential
    defined in eq. (13) in Faucher-Giguère & Kaspi (2006).

    Args:
        r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (float): Height from the galactic disk in [kpc].

    Returns:
        (np.ndarray): Gradient of the galactic potential in cylindrical
            coordinates.
    """

    dpot_dh_dr, dpot_dh_dz = self.r_z_derivatives_dh_potential(r, z)
    dpot_b_dr, dpot_b_dz = self.r_z_derivatives_b_potential(r, z)
    dpot_n_dr, dpot_n_dz = self.r_z_derivatives_n_potential(r, z)

    dpot_mw_dr = dpot_dh_dr + dpot_b_dr + dpot_n_dr
    dpot_mw_dphi = 0.0
    dpot_mw_dz = dpot_dh_dz + dpot_b_dz + dpot_n_dz

    pot_mw_gradient = np.array([dpot_mw_dr, dpot_mw_dphi, dpot_mw_dz])

    return pot_mw_gradient

dh_potential(r, z)

A cylindrical disk-halo component gravitational potential as defined in eq. (14) in Faucher-Giguère & Kaspi (2006).

Parameters:

Name Type Description Default
r ndarray

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z ndarray

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
ndarray

Value of the disk-halo potential in [erg/g].
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
def dh_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
    """
    A cylindrical disk-halo component gravitational potential as defined in eq. (14) in
    Faucher-Giguère & Kaspi (2006).

    Args:
        r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (np.ndarray): Height from the galactic disk in [kpc].

    Returns:
        (np.ndarray): Value of the disk-halo potential in [erg/g].
    """

    K, _ = self.shape_parameter(z)

    M_dh = self.M_dh
    b_dh = self.b_dh

    pot_dh = (
        -const.G
        * M_dh
        / (np.sqrt(K**2.0 + b_dh**2.0 + r**2.0) * const.KPC_TO_CM)
    )

    return pot_dh

n_potential(r, z)

A spherical nucleus component gravitational potential as defined in eq. (15) in Faucher-Giguère & Kaspi (2006). Note that we assume this potential to be spherical following Carlberg & Innanen (1987), while Faucher-Giguère & Kaspi take r to be the cylindrical coordinate, not the spherical one.

Parameters:

Name Type Description Default
r ndarray

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z ndarray

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
ndarray

Value of the nucleus potential.
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
def n_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
    """
    A spherical nucleus component gravitational potential as defined in eq. (15) in
    Faucher-Giguère & Kaspi (2006). Note that we assume this potential to be
    spherical following Carlberg & Innanen (1987), while Faucher-Giguère & Kaspi take r
    to be the cylindrical coordinate, not the spherical one.

    Args:
        r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (np.ndarray): Height from the galactic disk in [kpc].

    Returns:
        (np.ndarray): Value of the nucleus potential.
    """
    M_n = self.M_n
    b_n = self.b_n

    pot_n = (
        -const.G
        * M_n
        / (np.sqrt(b_n**2.0 + r**2.0 + z**2.0) * const.KPC_TO_CM)
    )

    return pot_n

r_z_derivatives_b_potential(r, z)

Derivatives with respect to r and z of the bulge component gravitational potential defined in eq. (15) in Faucher-Giguère & Kaspi (2006).

Parameters:

Name Type Description Default
r float

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z float

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
Tuple[float, float]

Derivatives with respect to r and z of the bulge potential.
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
def r_z_derivatives_b_potential(
    self, r: float, z: float
) -> Tuple[float, float]:
    """
    Derivatives with respect to r and z of the bulge component gravitational potential
    defined in eq. (15) in Faucher-Giguère & Kaspi (2006).

    Args:
        r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (float): Height from the galactic disk in [kpc].

    Returns:
        (Tuple[float, float]): Derivatives with respect to r and z of the bulge potential.
    """

    M_b = self.M_b
    b_b = self.b_b

    dpot_b_dr = (
        const.G_KPC_YR
        * M_b
        * r
        * (b_b**2.0 + r**2.0 + z**2.0) ** (-3.0 / 2.0)
    )

    dpot_b_dz = (
        const.G_KPC_YR
        * M_b
        * z
        * (b_b**2.0 + r**2.0 + z**2.0) ** (-3.0 / 2.0)
    )

    return dpot_b_dr, dpot_b_dz

r_z_derivatives_dh_potential(r, z)

Derivatives with respect to r and z of the disk-halo component gravitational potential defined in eq. (14) in Faucher-Giguère & Kaspi (2006).

Parameters:

Name Type Description Default
r float

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z float

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
Tuple[float, float]

Derivatives with respect to r and z of the disk-halo potential.
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
def r_z_derivatives_dh_potential(
    self, r: float, z: float
) -> Tuple[float, float]:
    """
    Derivatives with respect to r and z of the disk-halo component gravitational
    potential defined in eq. (14) in Faucher-Giguère & Kaspi (2006).

    Args:
        r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (float): Height from the galactic disk in [kpc].

    Returns:
        (Tuple[float, float]): Derivatives with respect to r and z of the disk-halo
            potential.
    """

    M_dh = self.M_dh
    b_dh = self.b_dh

    K, dK_dz = self.shape_parameter(z)
    dpot_dh_dr = (
        const.G_KPC_YR
        * M_dh
        * r
        * (K**2.0 + b_dh**2.0 + r**2.0) ** (-3.0 / 2.0)
    )
    dpot_dh_dz = (
        const.G_KPC_YR
        * M_dh
        * (K**2.0 + b_dh**2.0 + r**2.0) ** (-3.0 / 2.0)
        * K
        * dK_dz
    )

    return dpot_dh_dr, dpot_dh_dz

r_z_derivatives_n_potential(r, z)

Derivatives with respect to r and z of the nucleus component gravitational potential defined in eq. (15) in Faucher-Giguère & Kaspi (2006).

Parameters:

Name Type Description Default
r float

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z float

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
Tuple[float, float]

Derivatives with respect to r and z of the nucleus potential.
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
def r_z_derivatives_n_potential(
    self, r: float, z: float
) -> Tuple[float, float]:
    """
    Derivatives with respect to r and z of the nucleus component gravitational potential
    defined in eq. (15) in Faucher-Giguère & Kaspi (2006).

    Args:
        r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (float): Height from the galactic disk in [kpc].

    Returns:
        (Tuple[float, float]): Derivatives with respect to r and z of the nucleus potential.
    """

    M_n = self.M_n
    b_n = self.b_n

    dpot_n_dr = (
        const.G_KPC_YR
        * M_n
        * r
        * (b_n**2.0 + r**2.0 + z**2.0) ** (-3.0 / 2.0)
    )

    dpot_n_dz = (
        const.G_KPC_YR
        * M_n
        * z
        * (b_n**2.0 + r**2.0 + z**2.0) ** (-3.0 / 2.0)
    )

    return dpot_n_dr, dpot_n_dz

shape_parameter(z)

Shape parameter for the disk-halo potential from Carlberg & Innamen (1987) and its derivative with respect to the height z from the galactic disk. First term in the denominator of eq. (13) in Faucher-Giguère & Kaspi (2006).

Parameters:

Name Type Description Default
z float

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
Tuple[float, float]

Value of the shape parameter and its derivative with respect to z.
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
def shape_parameter(self, z: float) -> Tuple[float, float]:
    """
    Shape parameter for the disk-halo potential from Carlberg & Innamen
    (1987) and its derivative with respect to the height z from the galactic
    disk. First term in the denominator of eq. (13) in Faucher-Giguère &
    Kaspi (2006).

    Args:
        z (float): Height from the galactic disk in [kpc].

    Returns:
        (Tuple[float, float]): Value of the shape parameter and its derivative with
            respect to z.
    """

    a_d = self.a_d
    beta = self.beta
    h = self.h

    K = (
        a_d
        + beta[0] * np.sqrt(z**2.0 + h[0] ** 2.0)
        + beta[1] * np.sqrt(z**2.0 + h[1] ** 2.0)
        + beta[2] * np.sqrt(z**2.0 + h[2] ** 2.0)
    )

    dK_dz = (
        beta[0] * z / np.sqrt(z**2.0 + h[0] ** 2.0)
        + beta[1] * z / np.sqrt(z**2.0 + h[1] ** 2.0)
        + beta[2] * z / np.sqrt(z**2.0 + h[2] ** 2.0)
    )

    return K, dK_dz

GalaxyModelM19

Bases: GalaxyModelBase

Galaxy model from Marchetti et al. (2019). This is a four-component galactic potential model consisting of a Hernquist bulge and nucleus (Hernquist 1990), a Miyamoto-Nagai disk (Miyamoto & Nagai 1975) and a Navarro-Frenk-White halo (Navarro et al. 1996). The parameters of the model are taken from Table 1 in Marchetti et al. (2019) and are chosen to fit the enclosed mass profile of the Milky Way (Bovy 2015).

Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
@jitclass(galaxyModelM19_spec)
class GalaxyModelM19(GalaxyModelBase):
    """
    Galaxy model from Marchetti et al. (2019). This is a four-component galactic
    potential model consisting of a Hernquist bulge and nucleus (Hernquist 1990),
    a Miyamoto-Nagai disk (Miyamoto & Nagai 1975) and a Navarro-Frenk-White halo
    (Navarro et al. 1996). The parameters of the model are taken from Table 1 in
    Marchetti et al. (2019) and are chosen to fit the enclosed mass profile of
    the Milky Way (Bovy 2015).
    """

    def __init__(self) -> None:
        # Parameters of the model, values from Table 1 in Marchetti et al. (2019).
        self.a_d = 3.0  # Scale length of the disk in [kpc].
        self.b_d = 0.28  # Scale height for the disk.
        self.M_d = 6.8e10 * const.M_SUN  # Disk+halo mass in [g].
        self.M_b = 5.0e9 * const.M_SUN  # Bulge mass in [g].
        self.r_b = 1.0  # Core radius of the bulge component in [kpc].
        self.M_n = 1.71e9 * const.M_SUN  # Nucleus mass in [g].
        self.r_n = 0.07  # Core radius of the nucleus component in [kpc].
        self.M_h = 5.4e11 * const.M_SUN  # Halo mass in [g].
        self.r_h = 15.62  # Core radius of the halo component in [kpc].

    def shape_parameter(self, z: float) -> Tuple[float, float]:
        """
        Shape parameter for the disk potential from Marchetti et al. (2019) and
        its derivative with respect to the height z from the galactic disk.
        Second term in the denominator of eq. (8) in Marchetti et al. (2019).

        Args:
            z (float): Height from the galactic disk in [kpc].

        Returns:
            (Tuple[float, float]): Value of the shape parameter and its derivative with
                respect to z.
        """

        a_d = self.a_d
        b_d = self.b_d

        K = a_d + np.sqrt(z**2.0 + b_d**2.0)

        dK_dz = z / np.sqrt(z**2.0 + b_d**2.0)

        return K, dK_dz

    def d_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
        """
        A cylindrical Miyamoto-Nagai disk component gravitational potential as defined in eq. (8) in
        Marchetti et al. (2019).

        Args:
            r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (np.ndarray): Height from the galactic disk in [kpc].

        Returns:
            (np.ndarray): Value of the disk-halo potential in [erg/g].
        """

        K, _ = self.shape_parameter(z)

        M_d = self.M_d

        pot_d = (
            -const.G * M_d / (np.sqrt(K**2.0 + r**2.0) * const.KPC_TO_CM)
        )

        return pot_d

    def b_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
        """
        A spherical Hernquist bulge component gravitational potential as defined in eq. (7) in
        Marchetti et al. (2019).

        Args:
            r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (np.ndarray): Height from the galactic disk in [kpc].

        Returns:
            (np.ndarray): Value of the bulge potential in [erg/g].
        """
        M_b = self.M_b
        r_b = self.r_b

        pot_b = (
            -const.G
            * M_b
            / ((r_b + np.sqrt(r**2.0 + z**2.0)) * const.KPC_TO_CM)
        )

        return pot_b

    def n_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
        """
        A spherical Hernquist nucleus component gravitational potential as defined in eq. (7) in
        Marchetti et al. (2019).

        Args:
            r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (np.ndarray): Height from the galactic disk in [kpc].

        Returns:
            (np.ndarray): Value of the bulge potential in [erg/g].
        """
        M_n = self.M_n
        r_n = self.r_n

        pot_n = (
            -const.G
            * M_n
            / ((r_n + np.sqrt(r**2.0 + z**2.0)) * const.KPC_TO_CM)
        )

        return pot_n

    def h_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
        """
        A spherical Navarro-Frenk-White halo component gravitational potential as defined in
        eq. (9) in Marchetti et al. (2019).

        Args:
            r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (np.ndarray): Height from the galactic disk in [kpc].

        Returns:
            (np.ndarray): Value of the bulge potential in [erg/g].
        """
        M_h = self.M_h
        r_h = self.r_h

        pot_h = (
            -const.G
            * M_h
            / (np.sqrt(r**2.0 + z**2.0) * const.KPC_TO_CM)
            * np.log(1.0 + np.sqrt(r**2.0 + z**2.0) / r_h)
        )

        return pot_h

    def MW_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
        """
        Total Milky Way gravitational potential in Marchetti et al. (2019).

        Args:
            r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (np.ndarray): Height from the galactic disk in [kpc].

        Returns:
            (np.ndarray): Value of the Galactic potential in [erg].
        """

        MW_pot = (
            self.d_potential(r, z)
            + self.b_potential(r, z)
            + self.n_potential(r, z)
            + self.h_potential(r, z)
        )

        return MW_pot

    def r_z_derivatives_d_potential(
        self, r: float, z: float
    ) -> Tuple[float, float]:
        """
        Derivatives with respect to r and z of the disk component gravitational
        potential defined in eq. (8) in Marchetti et al. (2019).

        Args:
            r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (float): Height from the galactic disk in [kpc].

        Returns:
            (Tuple[float, float]): Derivatives with respect to r and z of the disk potential.
        """

        M_d = self.M_d

        K, dK_dz = self.shape_parameter(z)
        dpot_d_dr = (
            const.G_KPC_YR * M_d * r * (r**2.0 + K**2.0) ** (-3.0 / 2.0)
        )
        dpot_d_dz = (
            const.G_KPC_YR
            * M_d
            * (r**2.0 + K**2.0) ** (-3.0 / 2.0)
            * K
            * dK_dz
        )

        return dpot_d_dr, dpot_d_dz

    def r_z_derivatives_b_potential(
        self, r: float, z: float
    ) -> Tuple[float, float]:
        """
        Derivatives with respect to r and z of the bulge component gravitational potential
        defined in eq. (7) in Marchetti et al. (2019).

        Args:
            r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (float): Height from the galactic disk in [kpc].
        Returns:
            (Tuple[float, float]): Derivatives with respect to r and z of the bulge potential.
        """

        M_b = self.M_b
        r_b = self.r_b

        dpot_b_dr = (
            const.G_KPC_YR
            * M_b
            * r
            * ((r_b + np.sqrt(r**2.0 + z**2.0)) ** (-2.0))
            * ((r**2.0 + z**2.0) ** (-1.0 / 2.0))
        )
        dpot_b_dz = (
            const.G_KPC_YR
            * M_b
            * z
            * ((r_b + np.sqrt(r**2.0 + z**2.0)) ** (-2.0))
            * ((r**2.0 + z**2.0) ** (-1.0 / 2.0))
        )

        return dpot_b_dr, dpot_b_dz

    def r_z_derivatives_n_potential(
        self, r: float, z: float
    ) -> Tuple[float, float]:
        """
        Derivatives with respect to r and z of the nucleus component gravitational potential
        defined in eq. (7) in Marchetti et al. (2019).

        Args:
            r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (float): Height from the galactic disk in [kpc].
        Returns:
            (Tuple[float, float]): Derivatives with respect to r and z of the nucleus potential.
        """

        M_n = self.M_n
        r_n = self.r_n

        dpot_n_dr = (
            const.G_KPC_YR
            * M_n
            * r
            * ((r_n + np.sqrt(r**2.0 + z**2.0)) ** (-2.0))
            * ((r**2.0 + z**2.0) ** (-1.0 / 2.0))
        )
        dpot_n_dz = (
            const.G_KPC_YR
            * M_n
            * z
            * ((r_n + np.sqrt(r**2.0 + z**2.0)) ** (-2.0))
            * ((r**2.0 + z**2.0) ** (-1.0 / 2.0))
        )

        return dpot_n_dr, dpot_n_dz

    def r_z_derivatives_h_potential(
        self, r: float, z: float
    ) -> Tuple[float, float]:
        """
        Derivatives with respect to r and z of the halo component gravitational potential
        defined in eq. (9) in Marchetti et al. (2019).

        Args:
            r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (float): Height from the galactic disk in [kpc].

        Returns:
            (Tuple[float, float]): Derivatives with respect to r and z of the halo potential.
        """

        M_h = self.M_h
        r_h = self.r_h

        dpot_h_dr = (const.G_KPC_YR * M_h * r / (r**2.0 + z**2.0)) * (
            np.log(1.0 + np.sqrt(r**2.0 + z**2.0) / r_h)
            / (np.sqrt(r**2.0 + z**2.0))
            - 1.0 / (r_h + np.sqrt(r**2.0 + z**2.0))
        )

        dpot_h_dz = (const.G_KPC_YR * M_h * z / (r**2.0 + z**2.0)) * (
            np.log(1.0 + np.sqrt(r**2.0 + z**2.0) / r_h)
            / (np.sqrt(r**2.0 + z**2.0))
            - 1.0 / (r_h + np.sqrt(r**2.0 + z**2.0))
        )

        return dpot_h_dr, dpot_h_dz

    def cylind_coord_gradient_mw_potential(
        self, r: float, z: float
    ) -> np.ndarray:
        """
        Gradient in cylindrical coordinates of the Milky Way gravitational potential for
        the components defined in eq. (7,8,9) in Marchetti et al. (2019).

        Args:
            r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
            z (float): Height from the galactic disk in [kpc].

        Returns:
            (np.ndarray): Gradient of the galactic potential in cylindrical
                coordinates.
        """

        dpot_d_dr, dpot_d_dz = self.r_z_derivatives_d_potential(r, z)
        dpot_b_dr, dpot_b_dz = self.r_z_derivatives_b_potential(r, z)
        dpot_n_dr, dpot_n_dz = self.r_z_derivatives_n_potential(r, z)
        dpot_h_dr, dpot_h_dz = self.r_z_derivatives_h_potential(r, z)

        dpot_mw_dr = dpot_d_dr + dpot_b_dr + dpot_n_dr + dpot_h_dr
        dpot_mw_dphi = 0.0
        dpot_mw_dz = dpot_d_dz + dpot_b_dz + dpot_n_dz + dpot_h_dz

        pot_mw_gradient = np.array([dpot_mw_dr, dpot_mw_dphi, dpot_mw_dz])

        return pot_mw_gradient

MW_potential(r, z)

Total Milky Way gravitational potential in Marchetti et al. (2019).

Parameters:

Name Type Description Default
r ndarray

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z ndarray

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
ndarray

Value of the Galactic potential in [erg].
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
def MW_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
    """
    Total Milky Way gravitational potential in Marchetti et al. (2019).

    Args:
        r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (np.ndarray): Height from the galactic disk in [kpc].

    Returns:
        (np.ndarray): Value of the Galactic potential in [erg].
    """

    MW_pot = (
        self.d_potential(r, z)
        + self.b_potential(r, z)
        + self.n_potential(r, z)
        + self.h_potential(r, z)
    )

    return MW_pot

b_potential(r, z)

A spherical Hernquist bulge component gravitational potential as defined in eq. (7) in Marchetti et al. (2019).

Parameters:

Name Type Description Default
r ndarray

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z ndarray

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
ndarray

Value of the bulge potential in [erg/g].
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
def b_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
    """
    A spherical Hernquist bulge component gravitational potential as defined in eq. (7) in
    Marchetti et al. (2019).

    Args:
        r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (np.ndarray): Height from the galactic disk in [kpc].

    Returns:
        (np.ndarray): Value of the bulge potential in [erg/g].
    """
    M_b = self.M_b
    r_b = self.r_b

    pot_b = (
        -const.G
        * M_b
        / ((r_b + np.sqrt(r**2.0 + z**2.0)) * const.KPC_TO_CM)
    )

    return pot_b

cylind_coord_gradient_mw_potential(r, z)

Gradient in cylindrical coordinates of the Milky Way gravitational potential for the components defined in eq. (7,8,9) in Marchetti et al. (2019).

Parameters:

Name Type Description Default
r float

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z float

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
ndarray

Gradient of the galactic potential in cylindrical coordinates.
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
def cylind_coord_gradient_mw_potential(
    self, r: float, z: float
) -> np.ndarray:
    """
    Gradient in cylindrical coordinates of the Milky Way gravitational potential for
    the components defined in eq. (7,8,9) in Marchetti et al. (2019).

    Args:
        r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (float): Height from the galactic disk in [kpc].

    Returns:
        (np.ndarray): Gradient of the galactic potential in cylindrical
            coordinates.
    """

    dpot_d_dr, dpot_d_dz = self.r_z_derivatives_d_potential(r, z)
    dpot_b_dr, dpot_b_dz = self.r_z_derivatives_b_potential(r, z)
    dpot_n_dr, dpot_n_dz = self.r_z_derivatives_n_potential(r, z)
    dpot_h_dr, dpot_h_dz = self.r_z_derivatives_h_potential(r, z)

    dpot_mw_dr = dpot_d_dr + dpot_b_dr + dpot_n_dr + dpot_h_dr
    dpot_mw_dphi = 0.0
    dpot_mw_dz = dpot_d_dz + dpot_b_dz + dpot_n_dz + dpot_h_dz

    pot_mw_gradient = np.array([dpot_mw_dr, dpot_mw_dphi, dpot_mw_dz])

    return pot_mw_gradient

d_potential(r, z)

A cylindrical Miyamoto-Nagai disk component gravitational potential as defined in eq. (8) in Marchetti et al. (2019).

Parameters:

Name Type Description Default
r ndarray

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z ndarray

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
ndarray

Value of the disk-halo potential in [erg/g].
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
def d_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
    """
    A cylindrical Miyamoto-Nagai disk component gravitational potential as defined in eq. (8) in
    Marchetti et al. (2019).

    Args:
        r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (np.ndarray): Height from the galactic disk in [kpc].

    Returns:
        (np.ndarray): Value of the disk-halo potential in [erg/g].
    """

    K, _ = self.shape_parameter(z)

    M_d = self.M_d

    pot_d = (
        -const.G * M_d / (np.sqrt(K**2.0 + r**2.0) * const.KPC_TO_CM)
    )

    return pot_d

h_potential(r, z)

A spherical Navarro-Frenk-White halo component gravitational potential as defined in eq. (9) in Marchetti et al. (2019).

Parameters:

Name Type Description Default
r ndarray

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z ndarray

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
ndarray

Value of the bulge potential in [erg/g].
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
def h_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
    """
    A spherical Navarro-Frenk-White halo component gravitational potential as defined in
    eq. (9) in Marchetti et al. (2019).

    Args:
        r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (np.ndarray): Height from the galactic disk in [kpc].

    Returns:
        (np.ndarray): Value of the bulge potential in [erg/g].
    """
    M_h = self.M_h
    r_h = self.r_h

    pot_h = (
        -const.G
        * M_h
        / (np.sqrt(r**2.0 + z**2.0) * const.KPC_TO_CM)
        * np.log(1.0 + np.sqrt(r**2.0 + z**2.0) / r_h)
    )

    return pot_h

n_potential(r, z)

A spherical Hernquist nucleus component gravitational potential as defined in eq. (7) in Marchetti et al. (2019).

Parameters:

Name Type Description Default
r ndarray

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z ndarray

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
ndarray

Value of the bulge potential in [erg/g].
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
def n_potential(self, r: np.ndarray, z: np.ndarray) -> np.ndarray:
    """
    A spherical Hernquist nucleus component gravitational potential as defined in eq. (7) in
    Marchetti et al. (2019).

    Args:
        r (np.ndarray): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (np.ndarray): Height from the galactic disk in [kpc].

    Returns:
        (np.ndarray): Value of the bulge potential in [erg/g].
    """
    M_n = self.M_n
    r_n = self.r_n

    pot_n = (
        -const.G
        * M_n
        / ((r_n + np.sqrt(r**2.0 + z**2.0)) * const.KPC_TO_CM)
    )

    return pot_n

r_z_derivatives_b_potential(r, z)

Derivatives with respect to r and z of the bulge component gravitational potential defined in eq. (7) in Marchetti et al. (2019).

Parameters:

Name Type Description Default
r float

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z float

Height from the galactic disk in [kpc].

 required

Returns: (Tuple[float, float]): Derivatives with respect to r and z of the bulge potential.

Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
def r_z_derivatives_b_potential(
    self, r: float, z: float
) -> Tuple[float, float]:
    """
    Derivatives with respect to r and z of the bulge component gravitational potential
    defined in eq. (7) in Marchetti et al. (2019).

    Args:
        r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (float): Height from the galactic disk in [kpc].
    Returns:
        (Tuple[float, float]): Derivatives with respect to r and z of the bulge potential.
    """

    M_b = self.M_b
    r_b = self.r_b

    dpot_b_dr = (
        const.G_KPC_YR
        * M_b
        * r
        * ((r_b + np.sqrt(r**2.0 + z**2.0)) ** (-2.0))
        * ((r**2.0 + z**2.0) ** (-1.0 / 2.0))
    )
    dpot_b_dz = (
        const.G_KPC_YR
        * M_b
        * z
        * ((r_b + np.sqrt(r**2.0 + z**2.0)) ** (-2.0))
        * ((r**2.0 + z**2.0) ** (-1.0 / 2.0))
    )

    return dpot_b_dr, dpot_b_dz

r_z_derivatives_d_potential(r, z)

Derivatives with respect to r and z of the disk component gravitational potential defined in eq. (8) in Marchetti et al. (2019).

Parameters:

Name Type Description Default
r float

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z float

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
Tuple[float, float]

Derivatives with respect to r and z of the disk potential.
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
def r_z_derivatives_d_potential(
    self, r: float, z: float
) -> Tuple[float, float]:
    """
    Derivatives with respect to r and z of the disk component gravitational
    potential defined in eq. (8) in Marchetti et al. (2019).

    Args:
        r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (float): Height from the galactic disk in [kpc].

    Returns:
        (Tuple[float, float]): Derivatives with respect to r and z of the disk potential.
    """

    M_d = self.M_d

    K, dK_dz = self.shape_parameter(z)
    dpot_d_dr = (
        const.G_KPC_YR * M_d * r * (r**2.0 + K**2.0) ** (-3.0 / 2.0)
    )
    dpot_d_dz = (
        const.G_KPC_YR
        * M_d
        * (r**2.0 + K**2.0) ** (-3.0 / 2.0)
        * K
        * dK_dz
    )

    return dpot_d_dr, dpot_d_dz

r_z_derivatives_h_potential(r, z)

Derivatives with respect to r and z of the halo component gravitational potential defined in eq. (9) in Marchetti et al. (2019).

Parameters:

Name Type Description Default
r float

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z float

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
Tuple[float, float]

Derivatives with respect to r and z of the halo potential.
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
def r_z_derivatives_h_potential(
    self, r: float, z: float
) -> Tuple[float, float]:
    """
    Derivatives with respect to r and z of the halo component gravitational potential
    defined in eq. (9) in Marchetti et al. (2019).

    Args:
        r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (float): Height from the galactic disk in [kpc].

    Returns:
        (Tuple[float, float]): Derivatives with respect to r and z of the halo potential.
    """

    M_h = self.M_h
    r_h = self.r_h

    dpot_h_dr = (const.G_KPC_YR * M_h * r / (r**2.0 + z**2.0)) * (
        np.log(1.0 + np.sqrt(r**2.0 + z**2.0) / r_h)
        / (np.sqrt(r**2.0 + z**2.0))
        - 1.0 / (r_h + np.sqrt(r**2.0 + z**2.0))
    )

    dpot_h_dz = (const.G_KPC_YR * M_h * z / (r**2.0 + z**2.0)) * (
        np.log(1.0 + np.sqrt(r**2.0 + z**2.0) / r_h)
        / (np.sqrt(r**2.0 + z**2.0))
        - 1.0 / (r_h + np.sqrt(r**2.0 + z**2.0))
    )

    return dpot_h_dr, dpot_h_dz

r_z_derivatives_n_potential(r, z)

Derivatives with respect to r and z of the nucleus component gravitational potential defined in eq. (7) in Marchetti et al. (2019).

Parameters:

Name Type Description Default
r float

Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.

 required
z float

Height from the galactic disk in [kpc].

 required

Returns: (Tuple[float, float]): Derivatives with respect to r and z of the nucleus potential.

Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
def r_z_derivatives_n_potential(
    self, r: float, z: float
) -> Tuple[float, float]:
    """
    Derivatives with respect to r and z of the nucleus component gravitational potential
    defined in eq. (7) in Marchetti et al. (2019).

    Args:
        r (float): Distance from the galactic rotation axis in [kpc], i.e., cylindrical r coordinate.
        z (float): Height from the galactic disk in [kpc].
    Returns:
        (Tuple[float, float]): Derivatives with respect to r and z of the nucleus potential.
    """

    M_n = self.M_n
    r_n = self.r_n

    dpot_n_dr = (
        const.G_KPC_YR
        * M_n
        * r
        * ((r_n + np.sqrt(r**2.0 + z**2.0)) ** (-2.0))
        * ((r**2.0 + z**2.0) ** (-1.0 / 2.0))
    )
    dpot_n_dz = (
        const.G_KPC_YR
        * M_n
        * z
        * ((r_n + np.sqrt(r**2.0 + z**2.0)) ** (-2.0))
        * ((r**2.0 + z**2.0) ** (-1.0 / 2.0))
    )

    return dpot_n_dr, dpot_n_dz

shape_parameter(z)

Shape parameter for the disk potential from Marchetti et al. (2019) and its derivative with respect to the height z from the galactic disk. Second term in the denominator of eq. (8) in Marchetti et al. (2019).

Parameters:

Name Type Description Default
z float

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
Tuple[float, float]

Value of the shape parameter and its derivative with respect to z.
Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
def shape_parameter(self, z: float) -> Tuple[float, float]:
    """
    Shape parameter for the disk potential from Marchetti et al. (2019) and
    its derivative with respect to the height z from the galactic disk.
    Second term in the denominator of eq. (8) in Marchetti et al. (2019).

    Args:
        z (float): Height from the galactic disk in [kpc].

    Returns:
        (Tuple[float, float]): Value of the shape parameter and its derivative with
            respect to z.
    """

    a_d = self.a_d
    b_d = self.b_d

    K = a_d + np.sqrt(z**2.0 + b_d**2.0)

    dK_dz = z / np.sqrt(z**2.0 + b_d**2.0)

    return K, dK_dz

initialize_galactic_model()

Initializing the galactic model employed in the simulation. We have implemented two versions, i.e., the galactic structure of Faucher-Giguère & Kaspi (2006) (with parameters from Kuijken & Gilmore (1989)) and Galaxy model from Marchetti et al. (2019). The galactic_model variable is made available on a global level.

Source code in mlpoppyns/simulator/stellar_dynamics/galactic_model.py 
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
def initialize_galactic_model() -> None:
    """
    Initializing the galactic model employed in the simulation. We have implemented two
    versions, i.e., the galactic structure of Faucher-Giguère & Kaspi (2006) (with
    parameters from Kuijken & Gilmore (1989)) and Galaxy model from Marchetti et al.
    (2019). The galactic_model variable is made available on a global level.
    """
    global galactic_model

    if cfg["galactic_model"] == "gmM19":
        galactic_model = GalaxyModelM19()
    elif cfg["galactic_model"] == "gmFK06":
        galactic_model = GalaxyModelFK06()
    else:
        raise ValueError(
            "The galactic model does not exist. Choose between gmFK06 or gmM19."
        )

mlpoppyns.simulator.stellar_dynamics.initial_position

Initial galactocentric position for the stellar population.

We follow Faucher-Giguère & Kaspi (2006) and choose a galactocentric coordinate system, where the Galactic center is located at the origin. In terms of Galactic latitude l and longitude b, the x-,y-, and z-axes are parallel to (l, b) = (90, 0), (180, 0) and (0, 90), respectively, forming a right-handed Cartesian frame. This implies that the Sun is positioned at (x=0, y=8.5 kpc). Moreover, we define r = (x^2 + y^2)^0.5 as the distance from the Galactic center in the Galactic plane and phi = arctan(y/x). Here, the angle phi is the same as theta in Faucher-Giguère & Kaspi (2006). We reserve the variable theta for the polar angle in a spherical coordinate system.

Authors:

Vanessa Graber (graber@ice.csic.es)
Michele Ronchi (ronchi@ice.csic.es)

calculate_noise_for_coordinates(r, NS_number)

Calculating noise for the angular and radial coordinate to smear out the distribution and avoid artificial features near the Galactic center; see Sec. 3.2.1 in Faucher-Giguère & Kaspi (2006) for details.

Parameters:

Name Type Description Default
r ndarray

Array of distances from the Galactic center in [kpc].

 required
NS_number int

Total number of neutron stars created in the simulation.

 required

Returns:

Type Description
Tuple[ndarray, ndarray]

Array of noise for the galactocentric coordinates phi [rad], r [kpc].
Source code in mlpoppyns/simulator/stellar_dynamics/initial_position.py 
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
def calculate_noise_for_coordinates(
    r: np.ndarray, NS_number: int
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Calculating noise for the angular and radial coordinate to smear out the
    distribution and avoid artificial features near the Galactic center;
    see Sec. 3.2.1 in Faucher-Giguère & Kaspi (2006) for details.

    Args:
        r (np.ndarray): Array of distances from the Galactic center in [kpc].
        NS_number (int): Total number of neutron stars created in the simulation.

    Returns:
        (Tuple[np.ndarray, np.ndarray]): Array of noise for the galactocentric coordinates
            phi [rad], r [kpc].

    """

    phi_corr = np.random.uniform(0, 2 * np.pi, NS_number) * np.exp(-0.35 * r)
    r_corr = np.random.normal(0, 0.07 * r, NS_number)

    return phi_corr, r_corr

calculate_r_phi_electron_density(t_age)

Calculating the position at birth of each random neutron star in a cylindrical reference frame according to the Galactic electron density distribution ymw16 (see Yao et al. 2017).

Using the notebook ns_distribution_ne_model.ipynb, we create a 2D NumPy array containing the electron density distribution in polar coordinates (r, phi). This 2D array is saved in the file YMW16_density_model.npy in mlpoppyns/simulator/stellar_dynamics, and it is used to sample the neutron star positions in the Galaxy.

Parameters:

Name Type Description Default
t_age ndarray

Array of neutron star ages in [yr].

 required

Returns:

Type Description
Tuple[ndarray, ndarray]

Polar r and phi coordinates in [kpc] and [rad], respectively, for each generated neutron star.
Source code in mlpoppyns/simulator/stellar_dynamics/initial_position.py 
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
def calculate_r_phi_electron_density(
    t_age: np.ndarray,
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Calculating the position at birth of each random neutron star in a cylindrical reference frame according to the
    Galactic electron density distribution ymw16 (see Yao et al. 2017).

    Using the notebook ns_distribution_ne_model.ipynb, we create a 2D NumPy array containing the electron density
    distribution in polar coordinates (r, phi). This 2D array is saved in the file YMW16_density_model.npy in
    mlpoppyns/simulator/stellar_dynamics, and it is used to sample the neutron star positions in the Galaxy.

    Args:
        t_age (np.ndarray): Array of neutron star ages in [yr].

    Returns:
        (Tuple[np.ndarray, np.ndarray]): Polar r and phi coordinates in [kpc] and [rad], respectively, for each
            generated neutron star.
    """

    # Load the neutron star density model table.
    # The model table has been generated through the Jupyter notebook located in
    # tutorials/analysis_notebooks/ns_distribution_ne_model.ipynb.
    # It contains an 2D array of density rho in cylindrical coordinates (r, phi).
    # The density in the table is already multiplied by the galactocentric distance r
    # to take into account the element of area correction.

    file = pathlib.Path().joinpath(
        cfg["path_to_software"],
        "mlpoppyns/simulator/stellar_dynamics/YMW16_density_model.npy",
    )
    NS_density_model = np.load(file)

    # Define the grid of coordinates.
    r_grid = np.linspace(0.0, cfg["r_extent"], NS_density_model.shape[0])
    phi_grid = np.linspace(0.0, 2.0 * np.pi, NS_density_model.shape[1])

    # Drawing a random distance from the Galactic center in [kpc] and a random
    # azimuthal angle in [rad] according to the 2d density model.
    r_rand, phi_rand = rs.random_from_pdf_2d(
        r_grid, phi_grid, NS_density_model, len(t_age)
    )

    # Propagating the azimuthal coordinate of each object backwards in time
    # (according to its age) to account for the rotation of the Galactic arms;
    # we assume that the density structure itself remains rigid.
    phi_rand = spiral_arm_time_evol(phi_rand, t_age)

    return r_rand, phi_rand

calculate_r_phi_spiral_model(t_age, spiral_model)

Calculating the position at birth of each random neutron star in a cylindrical reference frame for a given radial and spiral arm model.

Parameters:

Name Type Description Default
t_age ndarray

Array of neutron star ages in [yr].

 required
spiral_model SpiralModelBase

A class specifying the spiral arm structure model.

 required

Returns:

Type Description
Tuple[ndarray, ndarray]

Polar r and phi coordinates in [kpc] and [rad], respectively, for each generated neutron star.
Source code in mlpoppyns/simulator/stellar_dynamics/initial_position.py 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
def calculate_r_phi_spiral_model(
    t_age: np.ndarray, spiral_model: sm.SpiralModelBase
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Calculating the position at birth of each random neutron star in a cylindrical reference frame for a given radial
    and spiral arm model.

    Args:
        t_age (np.ndarray): Array of neutron star ages in [yr].
        spiral_model (sm.SpiralModelBase): A class specifying the spiral arm structure model.

    Returns:
        (Tuple[np.ndarray, np.ndarray]): Polar r and phi coordinates in [kpc] and [rad], respectively, for each
            generated neutron star.
    """

    # Initializing the radial density prescription.
    radial_model = cfg["radial_model"]
    if radial_model == "rmYK04":
        pdf_radial = pdf_radial_density_YK04
    elif radial_model == "rmVV21":
        pdf_radial = pdf_radial_density_VV21
    else:
        raise ValueError(
            "The radial density model pdf does not exist. Choose between rmYK04 or rmVV21."
        )

    # Randomly associate a spiral arm to each neutron star.
    arm_index_rand = spiral_model.generate_arm_index(
        cfg["arm_number"], len(t_age)
    )
    # Count the number of stars in the Local arm.
    NS_local = len(arm_index_rand[arm_index_rand == 5])

    # Drawing a random distance from the Galactic center in [kpc] for
    # each neutron star according to the radial stellar density.
    r_grid = np.logspace(
        np.log10(0.0001), np.log10(cfg["r_extent"]), cfg["resolution"]
    )

    r_pdf_rand = np.zeros(len(t_age))
    r_pdf_rand[arm_index_rand != 5] = rs.random_from_pdf(
        r_grid, pdf_radial, len(t_age) - NS_local
    )

    # Since the local arm has a shorter extent in r coordinate compared to the other spiral arms,
    # we need to draw the r position of stars in the local arm in the range [local_r_min, local_r_max].
    if NS_local != 0:
        r_grid_local = np.logspace(
            np.log10(spiral_model.local_r_min),
            np.log10(spiral_model.local_r_max),
            cfg["resolution"],
        )

        r_pdf_rand[arm_index_rand == 5] = rs.random_from_pdf(
            r_grid_local, pdf_radial, NS_local
        )

    # Evaluate the angular phi coordinate for each neutron star and
    # add noise to both galactocentric coordinates.
    phi = spiral_model.calculate_phi(r_pdf_rand, arm_index_rand)
    phi_rand, r_rand = smear_initial_coordinates(r_pdf_rand, phi, len(t_age))

    # Propagating the azimuthal coordinate of each object backwards in time
    # (according to its age) to account for the rotation of the Galactic arms;
    # we assume that the arm structure itself remains rigid.
    phi_rand = spiral_arm_time_evol(phi_rand, t_age)

    return r_rand, phi_rand

calculate_z(NS_number)

Drawing a random distance from the Galactic plane in [kpc] for each neutron star according to the probability density function for the height, and randomly distribute the stars above and below the Galactic plane, since we expect that this distribution to be symmetric.

Parameters:

Name Type Description Default
NS_number int

Number of neutron stars for which to draw a distance from the Galactic plane.

 required

Returns:

Type Description
ndarray

Array of distances z from the Galactic plane in [kpc] for the simulated neutron star.
Source code in mlpoppyns/simulator/stellar_dynamics/initial_position.py 
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
def calculate_z(NS_number: int) -> np.ndarray:
    """
    Drawing a random distance from the Galactic plane in [kpc] for each neutron
    star according to the probability density function for the height, and randomly distribute the stars above and below
    the Galactic plane, since we expect that this distribution to be symmetric.

    Args:
        NS_number (int): Number of neutron stars for which to draw a distance from the Galactic plane.

    Returns:
        (np.ndarray): Array of distances z from the Galactic plane in [kpc] for the simulated neutron star.
    """
    z_grid = np.logspace(
        np.log10(0.0001), np.log10(cfg["z_extent"]), cfg["resolution"]
    )
    z_pdf_rand = rs.random_from_pdf(z_grid, pdf_initial_height, NS_number)

    # Randomly distribute the stars above and below the Galactic plane.
    z_rand = random_scatter_about_plane(z_pdf_rand, NS_number)

    return z_rand

pdf_initial_height(z)

Probability density function for the height from the Galactic equatorial plane according to eq. (2) in Gullon et al. (2014).

Parameters:

Name Type Description Default
z ndarray

Distance from the Galactic plane in [kpc].

 required

Returns:

Type Description
ndarray

Distribution of stars per kpc in z direction.
Source code in mlpoppyns/simulator/stellar_dynamics/initial_position.py 
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
def pdf_initial_height(z: np.ndarray) -> np.ndarray:
    """
    Probability density function for the height from the Galactic equatorial plane
    according to eq. (2) in Gullon et al. (2014).

    Args:
        z (np.ndarray): Distance from the Galactic plane in [kpc].

    Returns:
        (np.ndarray): Distribution of stars per kpc in z direction.

    """

    # We use an exponential distribution as given by Wainscoat et al. (1992)
    # and choose a mean scale height characteristic for a young distribution as
    # obtained by Gullon et al. (2014).

    h_c = cfg["h_c"]
    pdf_z = 1.0 / h_c * np.exp(-z / h_c)

    return pdf_z

pdf_radial_density_VV21(r)

The Milky Way's radial density of supernova remnants including both core-collapse and thermonuclear supernovae. The model distribution is the exponential model from eq. (9) from the work of Verberne & Vink (2021).

Parameters:

Name Type Description Default
r ndarray

Distance from the Galactic center in [kpc].

 required

Returns:

Type Description
ndarray

SNR radial density in [1/kpc].
Source code in mlpoppyns/simulator/stellar_dynamics/initial_position.py 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def pdf_radial_density_VV21(r: np.ndarray) -> np.ndarray:
    """
    The Milky Way's radial density of supernova remnants including both core-collapse
    and thermonuclear supernovae. The model distribution is the exponential model from eq. (9)
    from the work of Verberne & Vink (2021).

    Args:
        r (np.ndarray): Distance from the Galactic center in [kpc].

    Returns:
        (np.ndarray): SNR radial density in [1/kpc].

    """

    # Check range of input.
    coco.check_radial_coordinate(r)

    # Here we keep R_sun = 8.0 kpc for consistency with the results
    # of Verberne & Vink (2021).
    rsun = 8.0  # Sun's distance from the Galactic center in [kpc].
    b = 2.46  # +0.39 -0.33

    # SNR surface density following eq. (9) of Verberne & Vink (2021).
    rho = np.exp(-b * (r - rsun) / rsun)

    # Multiply the stellar surface density with the area element in polar coordinates.
    pdf_r = 2 * np.pi * r * rho

    return pdf_r

pdf_radial_density_YK04(r)

The Milky Way's stellar radial density in the Galactic plane according to eq. (15) of Yusifov & Küçük (2004).

Parameters:

Name Type Description Default
r ndarray

Distance from the Galactic center in [kpc].

 required

Returns:

Type Description
ndarray

Stellar radial density in [1/kpc].
Source code in mlpoppyns/simulator/stellar_dynamics/initial_position.py 
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
def pdf_radial_density_YK04(r: np.ndarray) -> np.ndarray:
    """
    The Milky Way's stellar radial density in the Galactic plane according
    to eq. (15) of Yusifov & Küçük (2004).

    Args:
        r (np.ndarray): Distance from the Galactic center in [kpc].

    Returns:
        (np.ndarray): Stellar radial density in [1/kpc].

    """

    # check range of input
    coco.check_radial_coordinate(r)

    # Here we keep R_sun = 8.5 kpc for consistency with the results
    # of Yusifov & Küçük (2004)
    rsun = 8.5  # Sun's distance from the Galactic center in [kpc].
    A = 37.6  # +- 1.90 [1/kpc^2]
    a = 1.64  # +-0.11
    b = 4.01  # +-0.24
    r1 = 0.55  # +- 0.10 [kpc]

    # Stellar surface density following eq. (15) of Yusifov & Küçük (2004).
    rho = (
        A
        * ((r + r1) / (rsun + r1)) ** a
        * np.exp(-b * (r - rsun) / (rsun + r1))
    )

    # Multiply the stellar surface density with the area element in polar coordinates.
    pdf_r = 2 * np.pi * r * rho

    return pdf_r

random_scatter_about_plane(z, NS_number)

Randomly distribute positive height values within z about the Galactic plane located at z=0.

Parameters:

Name Type Description Default
z ndarray

Array of heights in [kpc] with positive values.

 required
NS_number int

Total number of neutron stars created in the simulation.

 required

Returns:

Type Description
ndarray

Array of heights in [kpc] randomly scattered above or below 0.
Source code in mlpoppyns/simulator/stellar_dynamics/initial_position.py 
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
def random_scatter_about_plane(z: np.ndarray, NS_number: int) -> np.ndarray:
    """
    Randomly distribute positive height values within z about the Galactic plane
    located at z=0.

    Args:
        z (np.ndarray): Array of heights in [kpc] with positive values.
        NS_number (int): Total number of neutron stars created in the simulation.

    Returns:
        (np.ndarray): Array of heights in [kpc] randomly scattered above or below 0.

    """

    # Check that z has the length of the number of neutron stars simulated.
    if len(z) != NS_number:
        raise ValueError("Input array has the wrong length")

    # For each neutron star determine if it is above (False) or below (True) the Galactic plane.
    if_below = np.random.choice([False, True], size=NS_number)
    z[if_below] = -z[if_below]

    return z

smear_initial_coordinates(r, phi, NS_number)

Smear the initial radial and angular coordinates in the galactocentric frame by adding noise.

Parameters:

Name Type Description Default
r ndarray

Distances from the Galactic center in [kpc].

 required
phi ndarray

Azimuthal coordinate of the stars on the spiral arms [rad].

 required
NS_number int

Total number of neutron stars created in the simulation.

 required

Returns:

Type Description
Tuple[ndarray, ndarray]

Galactocentric coordinates phi [rad], r [kpc] with noise applied.
Source code in mlpoppyns/simulator/stellar_dynamics/initial_position.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def smear_initial_coordinates(
    r: np.ndarray, phi: np.ndarray, NS_number: int
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Smear the initial radial and angular coordinates in the galactocentric frame by adding noise.

    Args:
        r (np.ndarray): Distances from the Galactic center in [kpc].
        phi (np.ndarray): Azimuthal coordinate of the stars on the spiral arms [rad].
        NS_number (int): Total number of neutron stars created in the simulation.

    Returns:
        (Tuple[np.ndarray, np.ndarray]): Galactocentric coordinates phi [rad], r [kpc] with
            noise applied.

    """

    phi_corr, r_corr = calculate_noise_for_coordinates(r, NS_number)

    phi = phi + phi_corr
    r = r + r_corr

    return phi, r

spiral_arm_time_evol(phi0, t)

Evolving the spiral arm positions backward for a given age t. We assume that the Galactic spiral structure rotates rigidly in clockwise direction with a period of 250 Myr (see 'A guided map to the spiral arms in the Galactic disk of the Milky Way' by Vallée (2017)).

Parameters:

Name Type Description Default
phi0 ndarray

Current angular positions in [rad] for the chosen spiral pattern.

 required
t ndarray

Times in [yr] to propagate backward.

 required

Returns:

Type Description
ndarray

Angular positions in [rad] for the spiral pattern as they were t years ago.
Source code in mlpoppyns/simulator/stellar_dynamics/initial_position.py 
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
def spiral_arm_time_evol(phi0: np.ndarray, t: np.ndarray) -> np.ndarray:
    """
    Evolving the spiral arm positions backward for a given age t. We assume that the
    Galactic spiral structure rotates rigidly in clockwise direction with a
    period of 250 Myr (see 'A guided map to the spiral arms in the Galactic disk
    of the Milky Way' by Vallée (2017)).

    Args:
        phi0 (np.ndarray): Current angular positions in [rad] for the chosen
            spiral pattern.
        t (np.ndarray): Times in [yr] to propagate backward.

    Returns:
        (np.ndarray): Angular positions in [rad] for the spiral pattern as they were
            t years ago.

    """

    # Evaluate the angular velocity of rotation of the spiral pattern;
    # T is the period of rotation in [yr].
    T = 2.5e8
    omega_spiral_arms = 2.0 * np.pi / T

    # Find the values of the angles phi t years ago.
    phi_t = phi0 + omega_spiral_arms * t

    return phi_t

mlpoppyns.simulator.stellar_dynamics.initial_velocity

Initial velocity distribution for the stellar population.

The velocity is composed of two contributions, the neutron stars' proper motion caused by kicks during the supernova as well as the motion of the galaxy itself. For the former, we follow Gullon et al. (2014).

Authors:

Vanessa Graber (graber@ice.csic.es)
Michele Ronchi (ronchi@ice.csic.es)

circular_velocity(r, z)

Circular velocity in [kpc/yr] for a circular orbit at a distance r from the galactic center and at a height z from the galactic plane. This velocity is evaluated by assuming equilibrium between the gravitational acceleration in the r direction due to the galactic potential and the centrifugal acceleration due to rotation.

Parameters:

Name Type Description Default
r float

Distance in the galactic disk from the galactic center in [kpc].

 required
z float

Height from the galactic disk in [kpc].

 required

Returns:

Type Description
float

Value of the circular velocity in [kpc/yr].
Source code in mlpoppyns/simulator/stellar_dynamics/initial_velocity.py 
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
def circular_velocity(r: float, z: float) -> float:
    """
    Circular velocity in [kpc/yr] for a circular orbit at a distance r from the
    galactic center and at a height z from the galactic plane. This velocity is
    evaluated by assuming equilibrium between the gravitational acceleration in the
    r direction due to the galactic potential and the centrifugal acceleration due
    to rotation.

    Args:
        r (float): Distance in the galactic disk from the galactic center in [kpc].
        z (float): Height from the galactic disk in [kpc].

    Returns:
        (float): Value of the circular velocity in [kpc/yr].
    """

    pot_mw_gradient = gm.galactic_model.cylind_coord_gradient_mw_potential(
        r, z
    )
    v_circular = np.sqrt(r * pot_mw_gradient[0])

    return v_circular

kick_velocity_double_maxwell(NS_number)

Draw random kick velocity values from a double Maxwell distribution as suggested in Igoshev (2020).

Parameters:

Name Type Description Default
NS_number int

Total number of neutron stars created in the simulation.

 required

Returns:

Type Description
ndarray

Kick velocities in [kpc/yr].
Source code in mlpoppyns/simulator/stellar_dynamics/initial_velocity.py 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
def kick_velocity_double_maxwell(NS_number: int) -> np.ndarray:
    """
    Draw random kick velocity values from a double Maxwell distribution as suggested in Igoshev (2020).

    Args:
        NS_number (int): Total number of neutron stars created in the simulation.

    Returns:
        (np.ndarray): Kick velocities in [kpc/yr].
    """

    # Drawing a random magnitude of the birth kick velocity in [km/s] for each
    # neutron star according to the underlying velocity probability density
    # function.
    vk_grid = np.linspace(0.0, cfg["vk_extent"], cfg["resolution"])
    vk_rand = rs.random_from_pdf(
        vk_grid, pdf_kick_velocity_double_maxwell, NS_number
    )
    # Convert from [km/s] to [kpc/yr].
    vk_rand = vk_rand * const.YR_TO_S / const.KPC_TO_KM

    return vk_rand

kick_velocity_exp(NS_number)

Draw random kick velocity values from an exponential distribution as suggested in Faucher-Giguère amd Kaspi (2006).

Parameters:

Name Type Description Default
NS_number int

Total number of neutron stars created in the simulation.

 required

Returns:

Type Description
ndarray

Kick velocities in [kpc/yr].
Source code in mlpoppyns/simulator/stellar_dynamics/initial_velocity.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
def kick_velocity_exp(NS_number: int) -> np.ndarray:
    """
    Draw random kick velocity values from an exponential distribution as suggested in Faucher-Giguère amd Kaspi (2006).

    Args:
        NS_number (int): Total number of neutron stars created in the simulation.

    Returns:
        (np.ndarray): Kick velocities in [kpc/yr].
    """

    # Drawing a random magnitude of the birth kick velocity in [km/s] for each
    # neutron star according to the underlying velocity probability density
    # function.
    vk_grid = np.linspace(0.0, cfg["vk_extent"], cfg["resolution"])
    vk_rand = rs.random_from_pdf(vk_grid, pdf_kick_velocity_exp, NS_number)
    # Convert from [km/s] to [kpc/yr].
    vk_rand = vk_rand * const.YR_TO_S / const.KPC_TO_KM

    return vk_rand

kick_velocity_lognormal(mean, sigma, NS_number)

Draw random kick velocity values from a double Maxwell distribution as suggested in Disberg and Mandel (2025).

Parameters:

Name Type Description Default
mean float

Mean of the log-normal distribution.

 required
sigma float

Standard deviation of the log-normal distribution, in [s].

 required
NS_number int

Total number of neutron stars created in the simulation.

 required

Returns:

Type Description
ndarray

Kick velocities in [kpc/yr].
Source code in mlpoppyns/simulator/stellar_dynamics/initial_velocity.py 
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
def kick_velocity_lognormal(
    mean: float, sigma: float, NS_number: int
) -> np.ndarray:
    """
    Draw random kick velocity values from a double Maxwell distribution as suggested in Disberg and Mandel (2025).

    Args:
        mean (float): Mean of the log-normal distribution.
        sigma (float): Standard deviation of the log-normal distribution, in [s].
        NS_number (int): Total number of neutron stars created in the simulation.

    Returns:
        (np.ndarray): Kick velocities in [kpc/yr].
    """

    vk_rand = np.e ** np.random.normal(mean, sigma, NS_number)

    # Convert from [km/s] to [kpc/yr].
    vk_rand = vk_rand * const.YR_TO_S / const.KPC_TO_KM

    return vk_rand

kick_velocity_maxwell(NS_number)

Draw random kick velocity values from a Maxwell distribution as suggested in Hobbs et al. (2005).

Parameters:

Name Type Description Default
NS_number int

Total number of neutron stars created in the simulation.

 required

Returns:

Type Description
ndarray

Kick velocities in [kpc/yr].
Source code in mlpoppyns/simulator/stellar_dynamics/initial_velocity.py 
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
def kick_velocity_maxwell(NS_number: int) -> np.ndarray:
    """
    Draw random kick velocity values from a Maxwell distribution as suggested in Hobbs et al. (2005).

    Args:
        NS_number (int): Total number of neutron stars created in the simulation.

    Returns:
        (np.ndarray): Kick velocities in [kpc/yr].
    """

    # Drawing a random magnitude of the birth kick velocity in [km/s] for each
    # neutron star according to the underlying velocity probability density
    # function.
    vk_grid = np.linspace(0.0, cfg["vk_extent"], cfg["resolution"])
    vk_rand = rs.random_from_pdf(vk_grid, pdf_kick_velocity_maxwell, NS_number)
    # Convert from [km/s] to [kpc/yr].
    vk_rand = vk_rand * const.YR_TO_S / const.KPC_TO_KM

    return vk_rand

pdf_kick_velocity_double_maxwell(v)

Double Maxwell probability density function for the neutron stars' initial kick velocity magnitude following eq. (5) and section 4.2 in Igoshev (2020).

Parameters:

Name Type Description Default
v ndarray

Initial kick velocity magnitude in [km/s].

 required

Returns:

Type Description
ndarray

Stellar kick velocity distribution in [1/(km/s)].
Source code in mlpoppyns/simulator/stellar_dynamics/initial_velocity.py 
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
def pdf_kick_velocity_double_maxwell(v: np.ndarray) -> np.ndarray:
    """
    Double Maxwell probability density function for the neutron stars' initial kick
    velocity magnitude following eq. (5) and section 4.2 in Igoshev (2020).

    Args:
        v (np.ndarray): Initial kick velocity magnitude in [km/s].

    Returns:
        (np.ndarray): Stellar kick velocity distribution in [1/(km/s)].
    """

    # Define the dispersions of the two Maxwellian components.
    sigma_1 = cfg["sigma_k_comp1"]
    sigma_2 = cfg["sigma_k_comp2"]
    # Define the fractional contribution of the first Maxwellian.
    w = cfg["kick_weight_comp1"]

    if (w < 0) or (w > 1):
        raise ValueError(
            "The relative weight parameter of the km_double_maxwell kick velocity model must be in "
            "the range 0 and 1."
        )

    pdf_maxwell_1 = (
        np.sqrt(2 / np.pi)
        * v**2
        / (sigma_1**3)
        * np.exp(-(v**2) / (2 * sigma_1**2))
    )

    pdf_maxwell_2 = (
        np.sqrt(2 / np.pi)
        * v**2
        / (sigma_2**3)
        * np.exp(-(v**2) / (2 * sigma_2**2))
    )

    pdf_vk = w * pdf_maxwell_1 + (1.0 - w) * pdf_maxwell_2

    return pdf_vk

pdf_kick_velocity_exp(v)

Spherically symmetric exponential probability density function for the neutron stars' initial kick velocity magnitude following eq. (5) in Ofek (2009).

Parameters:

Name Type Description Default
v ndarray

Initial kick velocity magnitude in [km/s].

 required

Returns:

Type Description
ndarray

Stellar kick velocity distribution in [1/(km/s)].
Source code in mlpoppyns/simulator/stellar_dynamics/initial_velocity.py 
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def pdf_kick_velocity_exp(v: np.ndarray) -> np.ndarray:
    """
    Spherically symmetric exponential probability density function for the neutron stars' initial
    kick velocity magnitude following eq. (5) in Ofek (2009).

    Args:
        v (np.ndarray): Initial kick velocity magnitude in [km/s].

    Returns:
        (np.ndarray): Stellar kick velocity distribution in [1/(km/s)].
    """
    vk_mean = cfg["vk_c"]
    pdf_vk = v / vk_mean**2 * np.exp(-v / vk_mean)

    return pdf_vk

pdf_kick_velocity_maxwell(v)

Maxwell probability density function for the neutron stars' initial kick velocity magnitude following section 6.2 in Hobbs et al. (2005).

Parameters:

Name Type Description Default
v ndarray

Initial kick velocity magnitude in [km/s].

 required

Returns:

Type Description
ndarray

Stellar kick velocity distribution in [1/(km/s)].
Source code in mlpoppyns/simulator/stellar_dynamics/initial_velocity.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
def pdf_kick_velocity_maxwell(v: np.ndarray) -> np.ndarray:
    """
    Maxwell probability density function for the neutron stars' initial kick
    velocity magnitude following section 6.2 in Hobbs et al. (2005).

    Args:
        v (np.ndarray): Initial kick velocity magnitude in [km/s].

    Returns:
        (np.ndarray): Stellar kick velocity distribution in [1/(km/s)].
    """
    sigma = cfg["sigma_k"]
    pdf_vk = (
        np.sqrt(2 / np.pi)
        * v**2
        / (sigma**3)
        * np.exp(-(v**2) / (2 * sigma**2))
    )

    return pdf_vk

mlpoppyns.simulator.stellar_dynamics.load_dynamical_database

Module to load a dynamically evolved population database.

Authors:

Vanessa Graber (graber @ ice.csic.es)
Michele Ronchi (ronchi @ ice.csic.es)
Alberto Garcia-Garcia (garciagarcia @ ice.csic.es)
Celsa Pardo Araujo (pardo @ ice.csic.es)

load_database_dyn(dyn_path, n_batchsize, idx_remove, logger)

Load a batch of a dynamically evolved population from a .csv file, convert coordinates and return a dictionary with the dynamical information of the selected stars.

Parameters:

Name Type Description Default
dyn_path Path

Path to the directory containing the dynamically evolved population data.

 required
n_batchsize int

Batch size of stars to select when loading the data.

 required
idx_remove list

List of indices to remove from the dynamical database.

 required
logger Logger

Logger to use.

 required

Returns:

Type Description
dict

A dictionary containing the data of the selected dynamical population chunk.
Source code in mlpoppyns/simulator/stellar_dynamics/load_dynamical_database.py 
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
def load_database_dyn(
    dyn_path: pathlib.Path,
    n_batchsize: int,
    idx_remove: list,
    logger: logging.Logger,
) -> dict:
    """
    Load a batch of a dynamically evolved population from a .csv file, convert coordinates and return a dictionary
    with the dynamical information of the selected stars.

    Args:
        dyn_path (pathlib.Path): Path to the directory containing the dynamically evolved population data.
        n_batchsize (int): Batch size of stars to select when loading the data.
        idx_remove (list): List of indices to remove from the dynamical database.
        logger (logging.Logger): Logger to use.

    Returns:
        (dict): A dictionary containing the data of the selected dynamical population chunk.
    """

    # Check if the parsed dynamically simulated population directory exists.
    dyn_path = pathlib.Path(dyn_path)
    dyn_config_path = dyn_path / "configuration.json"
    dyn_data_path = dyn_path / "final_pop_dyn.csv"

    if not dyn_data_path.exists():
        logger.error(f"File {dyn_data_path} not found...")
        sys.exit()

    with open(dyn_config_path, "r") as f:
        config_dyn = json.load(f)

    # Load the batch of the file containing the dynamically evolved population parameters.
    df_dyn = mes.select(
        dyn_data_path,
        n_batchsize,
        config_dyn["NS_number"],
        idx_remove,
    )

    # Remove the second line of the header and convert the dataframe to a dictionary.
    df_dyn.columns = df_dyn.columns.get_level_values(0)
    dyn_dict = {col: df_dyn[col].to_numpy() for col in df_dyn.columns}

    # Add the dynamical properties in all coordinates and the index associated to each star and remove unnecessary keys.
    dictionary_dyn_database_chunk = (
        coco.convert_cylindrical_to_all_coordinates(dyn_dict)
    )
    dictionary_dyn_database_chunk["idx"] = df_dyn.index.values
    dictionary_dyn_database_chunk.pop("r", None)
    dictionary_dyn_database_chunk.pop("phi", None)
    dictionary_dyn_database_chunk.pop("z", None)
    dictionary_dyn_database_chunk.pop("v_r", None)
    dictionary_dyn_database_chunk.pop("v_phi", None)
    dictionary_dyn_database_chunk.pop("v_z", None)

    # Update the parameters in the simulation configuration file with the ones of the dynamical database.
    cfg["t_age_max"] = config_dyn["t_age_max"]
    cfg["NS_number"] = config_dyn["NS_number"]
    cfg["kick_model"] = config_dyn["kick_model"]
    if config_dyn["kick_model"] == "km_maxwell":
        cfg["sigma_k"] = config_dyn["sigma_k"]
    elif config_dyn["kick_model"] == "km_exp":
        cfg["vk_c"] = config_dyn["vk_c"]
    elif config_dyn["kick_model"] == "km_double_maxwell":
        cfg["sigma_k_comp1"] = config_dyn["sigma_k_comp1"]
        cfg["sigma_k_comp2"] = config_dyn["sigma_k_comp2"]
        cfg["kick_weight_comp1"] = config_dyn["kick_weight_comp1"]
    elif config_dyn["kick_model"] == "km_log-normal":
        cfg["vk_ln_mean"] = config_dyn["vk_ln_mean"]
        cfg["vk_ln_sigma"] = config_dyn["vk_ln_sigma"]
    cfg["h_c"] = config_dyn["h_c"]

    # Add the path of the dynamical database in the configuration file.
    cfg["dyn_database_path"] = str(dyn_path)

    return dictionary_dyn_database_chunk

mlpoppyns.simulator.stellar_dynamics.spiral_model

Model for the Milky Way spiral arms.

We consider two different models:

1) saFK06: A galactic spiral structure according to eq. (12) of Faucher-Giguère & Kaspi (2006) (see also Wainscoat et al. 1992). Their model consists of four arms plus a Local arm from Wainscoat et al. (1992).

2) saYMW17: A galactic spiral structure according to eq. (12) of Faucher-Giguère & Kaspi (2006) but with parameters re-adapted from Yao et al. (2017). Their model consists of four arms plus a Local arm from Hou et al. (2014).

Authors:

    Michele Ronchi (ronchi@ice.csic.es)

SpiralModelBase

Base class for any spiral model to ensure that a common interface between all of them is respected. The following abstract methods must be implemented or an error will be raised.

Source code in mlpoppyns/simulator/stellar_dynamics/spiral_model.py 
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
class SpiralModelBase:
    """
    Base class for any spiral model to ensure that a common interface between all of
    them is respected. The following abstract methods must be implemented or an error
    will be raised.
    """

    def __init__(self):
        self.arm_param = dict
        self.arm_probability = []
        self.local_r_min = float
        self.local_r_max = float
        super().__init__()

    def generate_arm_index(
        self, arm_number: int, NS_number: int
    ) -> np.ndarray:
        """
        Generate an array of indices related to the spiral arm associated to each star
        according to a probability evaluated taking into account the radial distribution of stars
        and the limited extension of the Local arm with respect to the other arms.

        Args:
            arm_number (int): Number of arms to simulate.
            NS_number (int): Number of stars to simulate.

        Returns:
            (np.ndarray): Array of random indices for the spiral arm associated to each star.

        """

        # Initialize the array of arm indices.
        arm_index_rand = np.zeros(NS_number)

        # If the Local arm is excluded, uniformly distribute the stars on the other arms.
        if arm_number == 4:
            arm_index_rand = np.random.randint(1, arm_number + 1, NS_number)

        # If the Local arm is included distribute the stars on the arms according to the probability
        # specified in arm_probability.
        elif arm_number == 5:
            choices = [1, 2, 3, 4, 5]
            arm_index_rand = np.array(
                random.choices(
                    choices, weights=self.arm_probability, k=NS_number
                )
            )

        return arm_index_rand

    def check_arm_index(self, arm_index: np.ndarray) -> None:
        """
        Check that the index for the spiral galaxy arms is not <1 or >5.

        Args:
            arm_index (np.ndarray): Index for the respective spiral arms.
        """
        if np.any(arm_index < 1) or np.any(arm_index > 5):
            raise ValueError("One of arm indices is out of range.")

    def calculate_phi(
        self, r: np.ndarray, arm_index: np.ndarray
    ) -> np.ndarray:
        """
        Calculating the angular coordinates of a neutron star for a given distance
        from the galactic center incorporating the Milky Way's arm structure according
        to eq. (12) of Faucher-Giguère & Kaspi (2006) (see also Wainscoat et al. 1992).

        Args:
            r (np.ndarray): Distances from the galactic center in [kpc].
            arm_index (np.ndarray): indices for the respective spiral arms,
                0 < arm_index < 6.

        Returns:
            (np.ndarray): Galactocentric phi coordinates in [rad].

        """

        # Check range of input.
        coco.check_radial_coordinate(r)
        self.check_arm_index(arm_index)

        arm_param_vect = np.vstack(
            np.vectorize(self.arm_param.get, otypes=[np.ndarray])(arm_index)
        )

        phi = (
            arm_param_vect[:, 0] * np.log(r / arm_param_vect[:, 1])
            + arm_param_vect[:, 2]
        )

        return phi

calculate_phi(r, arm_index)

Calculating the angular coordinates of a neutron star for a given distance from the galactic center incorporating the Milky Way's arm structure according to eq. (12) of Faucher-Giguère & Kaspi (2006) (see also Wainscoat et al. 1992).

Parameters:

Name Type Description Default
r ndarray

Distances from the galactic center in [kpc].

 required
arm_index ndarray

indices for the respective spiral arms, 0 < arm_index < 6.

 required

Returns:

Type Description
ndarray

Galactocentric phi coordinates in [rad].
Source code in mlpoppyns/simulator/stellar_dynamics/spiral_model.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
def calculate_phi(
    self, r: np.ndarray, arm_index: np.ndarray
) -> np.ndarray:
    """
    Calculating the angular coordinates of a neutron star for a given distance
    from the galactic center incorporating the Milky Way's arm structure according
    to eq. (12) of Faucher-Giguère & Kaspi (2006) (see also Wainscoat et al. 1992).

    Args:
        r (np.ndarray): Distances from the galactic center in [kpc].
        arm_index (np.ndarray): indices for the respective spiral arms,
            0 < arm_index < 6.

    Returns:
        (np.ndarray): Galactocentric phi coordinates in [rad].

    """

    # Check range of input.
    coco.check_radial_coordinate(r)
    self.check_arm_index(arm_index)

    arm_param_vect = np.vstack(
        np.vectorize(self.arm_param.get, otypes=[np.ndarray])(arm_index)
    )

    phi = (
        arm_param_vect[:, 0] * np.log(r / arm_param_vect[:, 1])
        + arm_param_vect[:, 2]
    )

    return phi

check_arm_index(arm_index)

Check that the index for the spiral galaxy arms is not <1 or >5.

Parameters:

Name Type Description Default
arm_index ndarray

Index for the respective spiral arms.

 required
Source code in mlpoppyns/simulator/stellar_dynamics/spiral_model.py 
 98
 99
100
101
102
103
104
105
106
def check_arm_index(self, arm_index: np.ndarray) -> None:
    """
    Check that the index for the spiral galaxy arms is not <1 or >5.

    Args:
        arm_index (np.ndarray): Index for the respective spiral arms.
    """
    if np.any(arm_index < 1) or np.any(arm_index > 5):
        raise ValueError("One of arm indices is out of range.")

generate_arm_index(arm_number, NS_number)

Generate an array of indices related to the spiral arm associated to each star according to a probability evaluated taking into account the radial distribution of stars and the limited extension of the Local arm with respect to the other arms.

Parameters:

Name Type Description Default
arm_number int

Number of arms to simulate.

 required
NS_number int

Number of stars to simulate.

 required

Returns:

Type Description
ndarray

Array of random indices for the spiral arm associated to each star.
Source code in mlpoppyns/simulator/stellar_dynamics/spiral_model.py 
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def generate_arm_index(
    self, arm_number: int, NS_number: int
) -> np.ndarray:
    """
    Generate an array of indices related to the spiral arm associated to each star
    according to a probability evaluated taking into account the radial distribution of stars
    and the limited extension of the Local arm with respect to the other arms.

    Args:
        arm_number (int): Number of arms to simulate.
        NS_number (int): Number of stars to simulate.

    Returns:
        (np.ndarray): Array of random indices for the spiral arm associated to each star.

    """

    # Initialize the array of arm indices.
    arm_index_rand = np.zeros(NS_number)

    # If the Local arm is excluded, uniformly distribute the stars on the other arms.
    if arm_number == 4:
        arm_index_rand = np.random.randint(1, arm_number + 1, NS_number)

    # If the Local arm is included distribute the stars on the arms according to the probability
    # specified in arm_probability.
    elif arm_number == 5:
        choices = [1, 2, 3, 4, 5]
        arm_index_rand = np.array(
            random.choices(
                choices, weights=self.arm_probability, k=NS_number
            )
        )

    return arm_index_rand

SpiralModelFK06

Bases: SpiralModelBase

Spiral structure of Faucher-Giguère & Kaspi (2006) with the addition of the Local arm from Wainscoat et al. (1992). Spiral arm parameters from table 2 in Faucher-Giguère & Kaspi (2006) and table 1 in Wainscoat et al. (1992) assuming a Sun galactocentric distance R_sun = 8.5 kpc. The Local arm has a radial extension ~ 0.55 rad in the range [1.13, 1.68] rad (see Wainscoat et al. 1992).

Source code in mlpoppyns/simulator/stellar_dynamics/spiral_model.py 
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
class SpiralModelFK06(SpiralModelBase):
    """
    Spiral structure of Faucher-Giguère & Kaspi (2006) with the addition of the Local arm
    from Wainscoat et al. (1992).
    Spiral arm parameters from table 2 in Faucher-Giguère & Kaspi (2006) and table 1 in
    Wainscoat et al. (1992) assuming a Sun galactocentric distance R_sun = 8.5 kpc.
    The Local arm has a radial extension ~ 0.55 rad in the range [1.13, 1.68] rad
    (see Wainscoat et al. 1992).
    """

    def __init__(self):
        # Parameters of the model, values from Table 2 in Faucher-Giguère & Kaspi (2006) and
        # table 1 in Wainscoat et al. (1992). Respectively winding constant k [rad], the inner
        # radius r0 [kpc] and the inner angle phi0 [rad]. The phi0 values in Wainscoat et al.
        # (1992) are increased by pi/2 and reported in the range [0, 2pi].
        super().__init__()
        self.arm_param = {
            1: np.array([4.25, 3.48, 1.57]),  # Norma
            2: np.array([4.25, 3.48, 4.71]),  # Carina-Sagittarius
            3: np.array([4.89, 4.90, 4.09]),  # Perseus
            4: np.array([4.89, 4.90, 0.95]),  # Crux-Scutum
            5: np.array([4.57, 8.10, 1.13]),  # Local
        }
        # Probability of a star to be located in each of the spiral arms.
        # This is evaluated considering the stellar density and the length of the spiral arms
        # in the range of galactocentric radii where the local arm extends.
        self.arm_probability = [0.24615, 0.24615, 0.24615, 0.24615, 0.0154]
        # Galactocentric distances where the Local arm starts and ends [kpc].
        # These values are obtained by evaluating the r coordinates from the phi coordinates
        # specified in Wainscoat et al. (2014).
        self.local_r_min = 8.10
        self.local_r_max = 9.14

SpiralModelYMW17

Bases: SpiralModelBase

Spiral structure of Yao et al. (2017), see also Hou et al. (2014). Spiral arm parameters from table 1 in Yao et al. (2017) assuming a Sun galactocentric distance R_sun = 8.3 kpc. The Local arm has a radial extension ~ 1.05 rad in the range [0.87, 1.92] rad (see Hou et al. 2014).

Source code in mlpoppyns/simulator/stellar_dynamics/spiral_model.py 
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
class SpiralModelYMW17(SpiralModelBase):
    """
    Spiral structure of Yao et al. (2017), see also Hou et al. (2014).
    Spiral arm parameters from table 1 in Yao et al. (2017) assuming a Sun
    galactocentric distance R_sun = 8.3 kpc. The Local arm has a radial extension
    ~ 1.05 rad in the range [0.87, 1.92] rad (see Hou et al. 2014).
    """

    def __init__(self):
        # Parameters of the model, values from Table 1 in Yao et al. (2017) re-adapted
        # to match the same logarithmic functional form used in Faucher-Giguère & Kaspi (2006).
        # Respectively winding constant k [rad], the inner radius r0 [kpc] and the inner angle
        # phi0 [rad].
        super().__init__()
        self.arm_param = {
            1: np.array([4.95, 3.35, 0.77]),  # Norma
            2: np.array([5.46, 3.56, 3.82]),  # Carina-Sagittarius
            3: np.array([5.77, 3.71, 2.09]),  # Perseus
            4: np.array([5.37, 3.67, 5.76]),  # Crux-Scutum
            5: np.array([20.67, 8.21, 0.96]),  # Local
        }
        # Probability of a star to be located in each of the spiral arms.
        # This is evaluated considering the stellar density and the length of the spiral arms
        # in the range of galactocentric radii where the local arm extends.
        self.arm_probability = [0.2455, 0.2455, 0.2455, 0.2455, 0.018]
        # Galactocentric distances where the Local arm starts and ends[kpc].
        # These values are obtained by evaluating the r coordinates from the phi coordinates
        # specified in Hou et al. (2014).
        self.local_r_min = 8.17
        self.local_r_max = 8.6

initialize_spiral_model()

Initializing the spiral model employed in the simulation. We have implemented two versions, i.e., the spiral structure of Faucher-Giguère & Kaspi (2006) (with the addition of the Local arm from Wainscoat et al. (1992)) and the spiral model from Yao et al. (2017). The spiral_model variable is made available on a global level.

Source code in mlpoppyns/simulator/stellar_dynamics/spiral_model.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def initialize_spiral_model() -> None:
    """
    Initializing the spiral model employed in the simulation. We have implemented two
    versions, i.e., the spiral structure of Faucher-Giguère & Kaspi (2006) (with
    the addition of the Local arm from Wainscoat et al. (1992)) and the spiral model from
    Yao et al. (2017). The spiral_model variable is made available on a global level.
    """
    global spiral_model

    if cfg["spiral_arms"] == "saFK06":
        spiral_model = SpiralModelFK06()
    elif cfg["spiral_arms"] == "saYMW17":
        spiral_model = SpiralModelYMW17()
    else:
        raise ValueError(
            "The spiral arm model does not exist. Choose between saFK06 or saYMW17."
        )