Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DOC: Add ellipticity example #471

Open
wants to merge 16 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file added examples/galaxies/SDSS_ellipticity.npz
Binary file not shown.
126 changes: 126 additions & 0 deletions examples/galaxies/plot_ellipticity.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,126 @@
"""
Galaxy Ellipticity Distributions
================================

This example demonstrate how to sample ellipticity distributions
in SkyPy.

"""


# %%
# 3D Ellipticity Distribution
# ---------------------------
#
# In Ryden 2004 [1]_, the ellipticity is sampled by randomly projecting
# a 3D ellipsoid with principal axes :math:`A > B > C`.
#
# The distribution of the axis ratio :math:`\gamma = C/A` is a truncated
# normal with mean :math:`\mu_\gamma` and standard deviation
# :math:`\sigma_\gamma`.
#
# The distribution of :math:`\epsilon = \log(1 - B/A)` is truncated normal
# with mean :math:`\mu` and standard deviation :math:`\sigma`.
#
#
# The func:`skypy.galaxies.morphology.ryden04_ellipticity()` model samples
# projected 2D axis ratios. Specifically, it samples the axis ratios of the
# 3D ellipsoid according to the description above and
# then randomly projects them using
# :func:`skypy.utils.random.triaxial_axis_ratio()`.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Potentially confusing to mention the triaxial_axis_ratio function. It is hidden from the user who doesn't need to understand the implementation in this context.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure I agree, but I'm open to delete that bit.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To be clear, we still want to explain the model (population of triaxial ellipsoids with some distribution in projection) but the implementation details (ryden04_ellipticity calls another function triaxial_axis_ratio) are at best irrelevant and at worst confusing for the user.

#
#
# Validation plot with SDSS Data
# ------------------------------
#
# Here we validate our function by comparing our ``SkyPy`` simulated galaxy
# ellipticities against observational data from SDSS DR1 in Figure 1 of Ryden 2004.
# You can download the data file
# :download:`SDSS_ellipticity <../../../examples/galaxies/SDSS_ellipticity.npz>`,
# stored as a 2D array: :math:`e_1`, :math:`e_2`.
#

import numpy as np

# Load SDSS data from Fig. 1 in Ryden 2004
data = np.load('SDSS_ellipticity.npz')
e1, e2 = data['sdss']['e1'], data['sdss']['e2']

Ngal = len(e1)
e = np.hypot(e1, e2)
q_amSDSS = np.sqrt((1 - e)/(1 + e))


# %%
#
# In this example, we generate 100 galaxy ellipticity samples using the
# ``SkyPy`` function
# :func:`skypy.galaxies.morphology.ryden04_ellipticity()` and the
# best fit parameters
# given in Ryden 2004 [1]:
# :math:`\mu_\gamma =0.222`, :math:`\sigma_\gamma=0.057`, :math:`\mu =-1.85`
# and :math:`\sigma=0.89`.
#

from skypy.galaxies.morphology import ryden04_ellipticity

# Best fit parameters from Fig. 1 in Ryden 2004
mu_gamma, sigma_gamma, mu, sigma = 0.222, 0.057, -1.85, 0.89

# Binning scheme of Fig. 1
bins = np.linspace(0, 1, 41)
mid = 0.5 * (bins[:-1] + bins[1:])

# Mean and variance of sampling
mean = np.zeros(len(bins)-1)
var = np.zeros(len(bins)-1)

# %%
#

# Create 100 SkyPy realisations
for i in range(100):
# sample ellipticity
ellipticity = ryden04_ellipticity(mu_gamma, sigma_gamma, mu, sigma, size=Ngal)
axis_ratio = (1 - ellipticity) / (1 + ellipticity)
n, _ = np.histogram(axis_ratio, bins=bins)

# update mean and variance
x = n - mean
mean += x/(i+1)
y = n - mean
var += x*y
Lucia-Fonseca marked this conversation as resolved.
Show resolved Hide resolved

# finalise variance and standard deviation
var = var/i
std = np.sqrt(var)


# %%
#
# We now plot the distribution of axis ratio :math:`𝑞_{am}`
# using adaptive moments in the i band, for exponential galaxies in the SDSS DR1
# (solid line). The data points with error bars represent
# the `SkyPy` simulation:
#

import matplotlib.pyplot as plt

plt.hist(q_amSDSS, range=[0, 1], bins=40, histtype='step',
ec='k', lw=0.5, label='SDSS data')
plt.errorbar(mid, mean, yerr=std, fmt='.r', ms=4, capsize=3,
lw=0.5, mew=0.5, label='SkyPy model')

plt.xlabel(r'Axis ratio, $q_{am}$')
plt.ylabel(r'N')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Units? Do we know the area covered by the data?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd say these quantities are dimensionless, right? One is a ratio and the other one is simply a count.
@ntessore ANy idea about the area? According to their paper, DR1 is 2099 deg2.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes it's axis ratio and number of galaxies, the area shouldn't come into it

plt.legend(frameon=False)
plt.show()

# %%
# References
# ----------
#
#
# .. [1] Ryden, Barbara S., 2004, `The Astrophysical Journal, Volume 601, Issue 1, pp. 214-220`_
#
# .. _The Astrophysical Journal, Volume 601, Issue 1, pp. 214-220: https://arxiv.org/abs/astro-ph/0310097