Quickstart

Here we take a look at a quick use case for lbparticles, namely finding a particle’s position and velocity some time after initialization.

from lbparticles import Precomputer, Particle, LogPotential, PotentialWrapper
import numpy as np

G= 0.00449987 # in pc^3 / (Msun Myr^2)
iHaveACopy=False # we don't yet have a copy of the precomputer, which we will need.

The first thing we’ll need after the imports above is the precomputer object.

The easiest thing to do is just to download it (this is what the next cell does). If you don’t run this cell the code later in the notebook will initialize a Precomputer object for you – this will take O(10 minutes) to run.

import urllib.request
lbprefn, headers = urllib.request.urlretrieve('https://www.dropbox.com/scl/fi/yo7ud56pxap10facsqk6k/big_10_1000_alpha2p2_lbpre.pickle?rlkey=dix9e1a5yn7gqjlw2sdtmfbnx&dl=1')
iHaveACopy=True

if iHaveACopy:
    #lbpre = Precomputer.load('../../big_10_1000_alpha2p2_lbpre.pickle') # adjust path to where you saved the pickle file.
    lbpre = Precomputer.load(lbprefn)
else:
    lbpre = Precomputer() # Takes a bit of time to run.
    lbpre.save() # Save a copy
    lbprefn = lbpre.identifier + '_lbpre.pickle' # remember the filename so we can load it in if this cell is called again
iHaveACopy=True

Now we initialize a particle’s position and velocity:

# the Sun's location in x,y,z cartesian coordinates (in parsecs) relative to the Galaxy centre
xcart = [8100.0, 0.0, 90.0] 
# similar to the Sun's velocity in vx, vy, vz (given the position xcart) in units of pc/Myr.
vcart = [-11.1, 12.24 + 220.0, 7.25] 

Now we need to set up the potential in which the particle will orbit. This is split into a vertical component and a radial component. For the vertical component, we assume that the vertical oscillation frequency is

$\nu_0 = \sqrt{4\pi G \cdot 0.2}$

a value about right for the self-gravity of the stellar disk in the Solar neighborhood.

We then assume that we can describe the radial variation of this vertical frequency with a powerlaw (though any function is supported).

$\nu = \nu_0 (r/8100\ \mathrm{pc})^{-\alpha/2}$

The powerlaw $\alpha$ describes the variation of the midplane density with radius, hence the factor of 2 in the exponent. This is a reasonable description of the variation near the Solar circle for $\alpha\approx 2.2$.

The radial component of the potential is described with a LogPotential, which takes as its argument the circular velocity. This potential produces a flat rotation curve.

# Set up the particle's vertical oscillation frequency as a function of radius
# vertical oscillation frequency at r=8100 pc.
nu0 = np.sqrt( 4.0*np.pi * G * 0.2) 
# powerlaw slope of the midplane density with radius, so that nu = nu0 (r/r0)^(-alpha/2)
alpha = 2.2 
# The potential, including the radial component (a logarithmic potential with a circular velocity of 220 pc/Myr) and a vertical component as above.
psir = PotentialWrapper(LogPotential(220.0), nur=lambda r:nu0*(r/8100.0)**(-alpha/2.0)) 

Finally we make a choice for how many terms to use in the expansions, and initialize the particle.

# number of terms used in the series to find the tangential position of the particle
ordershape = 10 
# number of terms used in the series to find the relationship between the particle's phase in its radial oscillation and the current time.
ordertime = 5 

# initialize the particle
part = Particle( xcart, vcart, psir, lbpre, ordershape=ordershape, ordertime=ordertime)

With the particle initialized, we can find its position and velocity at any time:

# find the particle's position and velocity 100 Myr later.
X,V = part.xvabs(100) 

print(X, V)

[-7899.86418169  4605.27704349  -121.26918128] [-104.14461122 -177.4118112     0.24529073]

# We get
# [-7899.86418167  4605.27704353  -121.26918128] [-104.14461122 -177.4118112     0.24529073]