******** Examples ******** For the examples below the following imports have already been executed:: import astropy.units as u import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np from gala.potential import scf .. _coeff-particle: Computing expansion coefficients from particle positions -------------------------------------------------------- To compute expansion coefficients for a distribution of particles or discrete samples from a density distribution, use ~gala.potential.scf.compute_coeffs_discrete. In this example, we will generate particle positions from a Plummer density profile, compute the expansion coefficients assuming spherical symmetry, then re-compute the expansion coefficients and variances (Weinberg 1996; [W96]_) allowing for non-spherical terms (e.g., :math:l,m>0). We'll start by generating samples from a Plummer sphere (see Section 3 of [HMV11]_ for more details). To do this, we will use inverse transform sampling by inverting the cumulative mass function (in this case, the mass enclosed): .. math:: \rho(r) &= \frac{M}{\frac{4}{3}\pi a^3} \, \left(1 + \frac{r^2}{a^2}\right)^{-5/2} m(0. We'll then plot the magnitude of the coefficients as a function of :math:n (but we'll ignore the sine terms, :math:T_{nlm} for this example):: mass = np.ones(n_samples) / n_samples S,T = scf.compute_coeffs_discrete(xyz, mass=mass, nmax=16, lmax=0, r_s=1.) plt.semilogy(np.abs(S[:,0,0]), marker=None, lw=2) plt.xlabel("$n$") plt.ylabel("$S_{n00}$") plt.tight_layout() .. plot:: :align: center :context: close-figs from gala.potential import scf mass = np.ones(n_samples) / n_samples S,T = scf.compute_coeffs_discrete(xyz, mass=mass, nmax=20, lmax=0, r_s=1.) plt.figure(figsize=(6,4)) plt.semilogy(np.abs(S[:,0,0]), marker=None, lw=2) plt.xlabel("$n$") plt.ylabel("$S_{n00}$") plt.tight_layout() In addition to computing the coefficient values, we can also compute the variances of the coefficients. Here we will relax the assumption about spherical symmetry by setting :math:l_{\rm max}=4. By computing the variance of each coefficient, we can estimate the signal-to-noise ratio of each expansion term and use this to help decide when to truncate the expansion (see [W96]_ for the methodology and reasoning behind this):: S, T, Cov = scf.compute_coeffs_discrete( xyz, mass=mass, r_s=1., nmax=10, lmax=4, skip_m=True, compute_var=True ) signal_to_noise = np.sqrt(S**2 / Cov[0, 0]) for l in range(S.shape[1]): plt.semilogy(signal_to_noise[:,l,0], marker=None, lw=2, alpha=0.5, label='l={}'.format(l)) plt.axhline(1., linestyle='dashed') plt.xlabel("$n$") plt.ylabel("$S/N$") plt.legend() .. plot:: :align: center :context: close-figs S, T, Cov = scf.compute_coeffs_discrete( xyz, mass=mass, r_s=1., nmax=10, lmax=4, skip_m=True, compute_var=True ) signal_to_noise = np.sqrt(S**2 / Cov[0, 0]) plt.figure(figsize=(6,4)) for l in range(S.shape[1]): plt.semilogy(signal_to_noise[:,l,0], marker=None, lw=2, alpha=0.5, label='l={}'.format(l)) plt.axhline(1., linestyle='dashed') plt.xlabel("$n$") plt.ylabel("$S/N$") plt.legend() plt.tight_layout() The horizontal line in the plot above is for a signal-to-noise ratio of 1 -- any coefficients with a SNR near or below this line are suspect and likely just adding noise to the expansion. Note that all of the SNR values for :math:l > 0 hover around 1 -- this is a good indication that we only need the :math:l=0 terms to accurately represent the density distribution of the particles. .. _coeff-analytic: Computing expansion coefficients for an analytic density -------------------------------------------------------- To compute expansion coefficients for an analytic density profile, use ~gala.potential.scf.compute_coeffs. In this example, we will write a function to evaluate an oblate density distribution and compute the expansion coefficients. We'll use a flattened Hernquist profile as our density profile: .. math:: \rho(s) &= \frac{M \, a}{2\pi} \, \frac{1}{s (s+a)^3} s^2 &= x^2 + y^2 + \frac{z^2}{q^2} In code:: def hernquist_density(r, M, a): return M*a / (2*np.pi) / (r*(r+a)**3) def flattened_hernquist_density(x, y, z, M, a, q): s = np.sqrt(x**2 + y**2 + (z/q)**2) return hernquist_density(s, M, a) The function to evaluate the density must take at least 3 arguments: the cartesian coordinates x, y, z. We'll again set :math:M=a=1 and we'll use a flattening :math:q=0.8. Let's visualize this by plotting isodensity contours in the :math:x-:math:z plane: .. plot:: :align: center :context: reset import astropy.units as u import matplotlib.pyplot as plt import matplotlib as mpl from matplotlib import ticker import numpy as np from gala.potential import scf def hernquist_density(r, M, a): return M*a / (2*np.pi) / (r*(r+a)**3) def flattened_hernquist_density(x, y, z, M, a, q): s = np.sqrt(x**2 + y**2 + (z/q)**2) return hernquist_density(s, M, a) M = 1. a = 1. q = 0.8 x,z = np.meshgrid(np.linspace(-10., 10., 128), np.linspace(-10., 10., 128)) y = np.zeros_like(x) dens = flattened_hernquist_density(x, y, z, M, a, q) plt.figure(figsize=(6,6)) plt.contourf(x, z, dens, cmap='magma', levels=np.logspace(np.log10(dens.min()), np.log10(dens.max()), 32), locator=ticker.LogLocator()) plt.title("Isodensity") plt.xlabel("$x$", fontsize=22) plt.ylabel("$z$", fontsize=22) plt.tight_layout() To compute the expansion coefficients, we pass the flattened_hernquist_density() function in to ~gala.potential.scf.compute_coeffs. Because this is an axisymmetric density, we will ignore terms with :math:m>0 by setting skip_m=True:: M = 1. a = 1. q = 0.8 coeff = scf.compute_coeffs(flattened_hernquist_density, nmax=8, lmax=8, M=M, r_s=a, args=(M,a,q), skip_m=True) (S,Serr),(T,Terr) = coeff Computing the coefficients involves a numerical integration that uses scipy.integrate.quad, which simultaneously estimates the error in the computed integral. ~gala.potential.scf.compute_coeffs returns the coefficient arrays and these error estimates. Now that we have the coefficients in hand, we can visualize their magnitudes:: plt.figure(figsize=(6,4)) plt.semilogy(np.abs(S[:,0,0]), marker=None, lw=2) plt.xlabel("$n$") plt.ylabel("$S_{n00}$") .. plot:: :align: center :context: close-figs nmax = 8 lmax = 8 coeff = scf.compute_coeffs(flattened_hernquist_density, nmax=nmax, lmax=lmax, M=M, r_s=a, args=(M,a,q), skip_m=True) (S,Serr),(T,Terr) = coeff plt.figure(figsize=(6,4)) plt.semilogy(np.abs(S[:,0,0]), marker=None, lw=2) plt.xlabel("$n$") plt.ylabel("$S_{n00}$") plt.tight_layout() Because we ignored any :math:m terms, the coefficients are computed in a 2D grid in :math:n,l: we can visualize their magnitude by coloring points on such a grid:: nl_grid = np.mgrid[0:lmax+1, 0:nmax+1] plt.figure(figsize=(5,4)) plt.scatter(nl_grid[0].ravel(), nl_grid[1].ravel(), c=np.abs(S[:,:,0].ravel()), norm=mpl.colors.LogNorm(), cmap='viridis', s=80) plt.xlabel('$n$') plt.ylabel('$l$') plt.colorbar() .. plot:: :align: center :context: close-figs nl_grid = np.mgrid[0:lmax+1, 0:nmax+1] plt.figure(figsize=(5,4)) plt.scatter(nl_grid[0].ravel(), nl_grid[1].ravel(), c=np.abs(S[:,:,0].ravel()), norm=mpl.colors.LogNorm(), cmap='viridis', s=80) plt.xlabel('$n$') plt.ylabel('$l$') plt.colorbar() plt.tight_layout() .. _potential-class: Using ~gala.potential.scf.SCFPotential to evaluate the density, potential, gradient ------------------------------------------------------------------------------------- In this example we'll continue where the :ref:previous example  left off: we now have computed expansion coefficients for a given density function and we would like to evaluate the gradient of the gravitational potential at various locations. We will use gala to integrate an orbit in the expansion potential. From the previous example, we have a set of cosine and sine coefficients (S and T) for an SCF representation of a flattened (oblate) Hernquist density profile. First, we'll create an ~gala.potential.scf.SCFPotential object using these coefficients:: potential = scf.SCFPotential(Snlm=S, Tnlm=T, m=M, r_s=a) # M=a=1 Let's compare how our expansion density to the true density by recreating the above isodensity contour figure with SCF density contours overlaid:: x,z = np.meshgrid(np.linspace(-10., 10., 128), np.linspace(-10., 10., 128)) y = np.zeros_like(x) true_dens = flattened_hernquist_density(x, y, z, M, a, q) # we need an array of positions with shape (3,n_samples) for SCFPotential xyz = np.vstack((x.ravel(),y.ravel(),z.ravel())) scf_dens = potential.density(xyz).value # log-spaced contour levels levels = np.logspace(np.log10(true_dens.min()), np.log10(true_dens.max()), 16) plt.figure(figsize=(6,6)) plt.contourf(x, z, true_dens, cmap='magma', levels=levels, locator=ticker.LogLocator()) plt.contour(x, z, scf_dens.reshape(x.shape), colors='w', levels=levels, locator=ticker.LogLocator()) plt.title("Isodensity") plt.xlabel("$x$", fontsize=22) plt.ylabel("$z$", fontsize=22) .. plot:: :align: center :context: close-figs potential = scf.SCFPotential(Snlm=S, Tnlm=T, m=M, r_s=a) # M=a=1 # we need an array of positions with shape (3,n_samples) for SCFPotential xyz = np.vstack((x.ravel(),y.ravel(),z.ravel())) scf_dens = potential.density(xyz).value # log-spaced contour levels true_dens = flattened_hernquist_density(x, y, z, M, a, q) levels = np.logspace(np.log10(true_dens.min()), np.log10(true_dens.max()), 16) plt.figure(figsize=(6,6)) plt.contourf(x, z, true_dens, cmap='magma', levels=levels, locator=ticker.LogLocator()) plt.contour(x, z, scf_dens.reshape(x.shape), colors='w', levels=levels, locator=ticker.LogLocator()) plt.title("Isodensity") plt.xlabel("$x$", fontsize=22) plt.ylabel("$z$", fontsize=22) plt.tight_layout() By eye, the SCF representation looks pretty good. Let's now create a plot of equipotential contours using the ~gala.potential.scf.SCFPotential instance:: scf_pot = np.abs(potential.energy(xyz)) scf_pot = scf_pot.value # get numerical value from ~astropy.units.Quantity # log-spaced contour levels levels = np.logspace(np.log10(scf_pot.min()), np.log10(scf_pot.max()), 16) plt.figure(figsize=(6,6)) plt.contour(x, z, scf_pot.reshape(x.shape), cmap='inferno_r', levels=levels, locator=ticker.LogLocator()) plt.title("Equipotential") plt.xlabel("$x$", fontsize=22) plt.ylabel("$z$", fontsize=22) .. plot:: :align: center :context: close-figs scf_pot = np.abs(potential.energy(xyz)) scf_pot = scf_pot.value # get numerical value from Astropy Quantity # log-spaced contour levels levels = np.logspace(np.log10(scf_pot.min()), np.log10(scf_pot.max()), 16) plt.figure(figsize=(6,6)) plt.contour(x, z, scf_pot.reshape(x.shape), cmap='inferno_r', levels=levels, locator=ticker.LogLocator()) plt.title("Equipotential") plt.xlabel("$x$", fontsize=22) plt.ylabel("$z$", fontsize=22) plt.tight_layout() (the above is actually provided as a convenience method of any ~gala.potential.PotentialBase subclass -- see ~gala.potential.PotentialBase.plot_contours). Now let's integrate an orbit in this potential. We'll use the orbit integration framework from gala.integrate and the convenience method ~gala.potential.scf.SCFPotential.integrate_orbit` to do this:: import gala.dynamics as gd # when using dimensionless units, we don't need to specify units for the # initial conditions w0 = gd.PhaseSpacePosition(pos=[1.,0,0.25], vel=[0.,0.3,0.]) # by default this uses Leapfrog integration orbit = potential.integrate_orbit(w0, dt=0.1, n_steps=10000) fig = orbit_l.plot(marker=',', linestyle='none', alpha=0.5) .. plot:: :align: center :context: close-figs import gala.dynamics as gd # when using dimensionless units, we don't need to specify units for the # initial conditions w0 = gd.PhaseSpacePosition(pos=[1.,0,0.25], vel=[0.,0.3,0.]) # by default this uses Leapfrog integration orbit = potential.integrate_orbit(w0, dt=0.1, n_steps=10000) fig = orbit.plot(marker=',', linestyle='none', alpha=0.5) References ---------- .. [W96] http://dx.doi.org/10.1086/177902 .. [HMV11] http://www.artcompsci.org/kali/vol/plummer/volume11.pdf