Commit 5c02300a authored by alexandre beelen's avatar alexandre beelen
Browse files

Added Gaussian ratio

parent 6c3f4a11
Pipeline #2719 failed with stage
in 0 seconds
......@@ -5,6 +5,7 @@ Basic Usage of powspec
"""
import matplotlib.pyplot as plt
# sphinx_gallery_thumbnail_path = '_static/demo.png'
import astropy.units as u
from astropy.convolution import Gaussian2DKernel
......@@ -29,9 +30,7 @@ res = 1 * u.arcmin
alphas = [-1, -2, -3]
images = []
for alpha in alphas:
images.append(
gen_pkfield(npix=1024, alpha=alpha, fknee=0.1 / u.arcmin, res=res) * u.MJy
)
images.append(gen_pkfield(npix=1024, alpha=alpha, fknee=0.1 / u.arcmin, res=res) * u.MJy)
# %%
# Compute P(k)
......@@ -71,12 +70,17 @@ plt.show()
# Create a fake catalog of sources
n_pix = 512
n_sources = 128*5
n_sources = 128 * 5
positions = np.random.uniform(0, n_pix, size=(2, n_sources))
fluxes = np.random.uniform(1, 10, n_sources)
sigma = 10 # pixels
sigma = 10 # pixels
images = [
gen_psffield(positions, fluxes, n_pix, kernel=Gaussian2DKernel(sigma)) * u.Jy / u.beam,
gen_psffield(positions, fluxes, n_pix, kernel=Gaussian2DKernel(sigma, x_size=n_pix // 2)) * u.Jy / u.beam,
]
labels = ["", "x_size"]
image = gen_psffield(positions, fluxes, n_pix, kernel=Gaussian2DKernel(sigma), factor=4) * u.Jy / u.beam
# %%
# Compute P(k)
......@@ -84,22 +88,48 @@ image = gen_psffield(positions, fluxes, n_pix, kernel=Gaussian2DKernel(sigma), f
#
# Compute power spectra of each images
#
powspec, k = power_spectral_density(image, res=res, range='tight', bins='auto')
powspecs = []
for image in images:
powspec, k = power_spectral_density(image, res=res, range="tight", bins=n_pix // 2)
powspecs.append(powspec)
# powspec, k = power_spectral_density(image, res=res, range='tight', bins='auto')
k_mid = np.mean(u.Quantity([k[1:], k[:-1]]), axis=0)
# %%
# Plots
# -----
fig, axes = plt.subplots(ncols=2)
axes[0].imshow(image.value, origin='lower')
axes[1].loglog(k_mid, powspec)
def gaussian_pk(k, sigma):
return np.exp(- (np.pi * k)**2 * (2 * sigma**2) * 2 )
return np.exp(-((np.pi * k) ** 2) * (2 * sigma ** 2) * 2)
axes[1].plot(k_mid, gaussian_pk(k_mid, sigma * res) * powspec.max())
axes[1].set_ylim(np.min(powspec), np.max(powspec))
fig = plt.figure()
gs = fig.add_gridspec(ncols=2, nrows=len(images))
# ax_pk = fig.add_subplot(gs[:, 0])
ax_pk = fig.add_subplot(gs[: len(images) // 2, 0])
ax_pk_ratio = fig.add_subplot(gs[len(images) // 2 :, 0])
for i, (image, powspec, label) in enumerate(zip(images, powspecs, labels)):
ax_pk.loglog(k_mid.to(u.arcmin ** -1), powspec.to(u.Jy ** 2 / u.beam ** 2 * u.arcmin ** 2), label=label)
ax_pk_ratio.semilogx(
k_mid.to(u.arcmin ** -1),
powspec.to(u.Jy ** 2 / u.beam ** 2 * u.arcmin ** 2) / gaussian_pk(k_mid, sigma * res),
label=label,
)
ax = fig.add_subplot(gs[i, 1])
ax.imshow(image.value, origin="lower")
ax_pk.plot(
k_mid.to(u.arcmin ** -1),
gaussian_pk(k_mid, sigma * res) * (u.Jy ** 2 / u.beam ** 2 * u.arcmin ** 2),
linestyle="--",
label="analytical",
)
ax_pk.legend()
ax_pk.set_ylim(1e-15, np.max(fluxes * 2 * np.pi * sigma ** 2) ** 2)
ratio = np.sum((2 * np.pi * sigma ** 2 * fluxes * u.Jy / u.beam) ** 2) * res ** 2 / n_pix ** 2
ax_pk_ratio.axhline(ratio, linestyle="--", label="analytical")
ax_pk_ratio.legend()
ax_pk_ratio.set_ylim(0, 3 * ratio)
fig.show()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment