Skip to content
Snippets Groups Projects
Commit 18c11076 authored by Julian Rüstig's avatar Julian Rüstig
Browse files

multi: a bit more compare

parent 7e7b329a
No related branches found
No related tags found
1 merge request!50Draft: Mosaic imaging
import matplotlib.pyplot as plt
from charm_lensing.utils import load_fits
import numpy as np
from numpy.typing import ArrayLike
from astropy import units as u
from astropy.coordinates import SkyCoord
from charm_lensing.utils import load_fits
from jubik0.jwst.reconstruction_grid import Grid
from jubik0.jwst.rotation_and_shift import build_nufft_rotation_and_shift
from jax import config
print("Using cpu")
config.update('jax_enable_x64', True)
config.update('jax_platform_name', 'cpu')
def prepare_input_comparison(
sky_input: ArrayLike,
input_grid: Grid,
reconstruction_grid: Grid):
interpolation_points = input_grid.wcs.index_from_wl(
reconstruction_grid.wl_coords())[0]
nufft = build_nufft_rotation_and_shift(
1, 1, input_grid.shape, interpolation_points.shape[1:])
interpolated_sky = nufft(sky_input.T, interpolation_points)
return interpolated_sky / input_grid.dvol * reconstruction_grid.dvol
# Prepare input image
incell, incell_unit = 0.03, u.arcsec
inbright = 0.004
sky_input = load_fits('./small/M51ha.fits') * inbright
center = SkyCoord(0*u.rad, 0*u.rad)
fov = [incell*sky_input.shape[i]*incell_unit for i in range(2)]
input_grid = Grid(center, sky_input.shape, fov)
fov = 25 * u.arcsec
npix = 128
npix = 512
# center_pp = SkyCoord(fov/npix, fov/npix)
reconstruction_grid = Grid(center, (npix,)*2, (fov,)*2)
sky_input = prepare_input_comparison(
sky_input, input_grid, reconstruction_grid)
distance = fov/npix
dvol = distance**2
conversion = dvol.to(u.rad**2).value / dvol.value
sky_casa = load_fits(f'./output/combined_{npix}/sky_casa_{npix}.fits')
sky_reso = load_fits(
f'./output/combined_{npix}/sky_reso_{npix}.fits') * conversion
sky_casa = sky_casa[0, 0]
sky_reso = sky_reso[0, 0]
sky_input /= sky_input.max()
sky_casa /= np.nanmax(sky_casa)
sky_reso /= sky_reso.max()
names = [['input', 'resolve', 'input - resolve'],
['input', 'casa', 'input - casa']]
images = [[sky_input, sky_reso, (sky_input-sky_reso)],
[sky_input, sky_casa, (sky_input-sky_casa)]]
cmaps = ['viridis', 'viridis', 'RdBu_r']
vmax = [None, None, +0.5]
vmin = [None, None, -0.5]
fig, axes = plt.subplots(1, 2)
ax, ay = axes
ax.imshow(sky_reso, origin='lower', vmax=10)
ay.imshow(sky_casa, origin='lower', vmax=10)
fig, axes = plt.subplots(2, 3, sharex=True, sharey=True)
for axe, image, name in zip(axes, images, names):
for ax, imag, nam, cmap, vma, vmi in zip(axe, image, name, cmaps, vmax, vmin):
im = ax.imshow(imag, origin='lower', vmax=vma, vmin=vmi, cmap=cmap)
plt.colorbar(im, ax=ax)
ax.set_title(nam)
plt.show()
from astropy.coordinates import SkyCoord
from jubik0 import build_astropy_wcs
from functools import reduce
import resolve as rve
......
[sky]
freq mode = single
polarization=I
space npix x = 128
space npix y = 128
space npix x = 512
space npix y = 512
space fov x = 25as
space fov y = 25as
# space fov x = 30.72as
......
Subproject commit 1fe138b081ba8e89252dbc90a73d818a826f575a
Subproject commit f0bac8b42da119d3febac5b809fe1eac2502dd66
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment