This notebook performs the PSF deconvolution code SHARPESST (Plowman et al. 2023) on the SPICE dataset analyzed in the manuscript. We aim to demonstrate that the consistent blue-shifted patches across tens of arcsec are real signals instead of a PSF artifact. Link to Figure C.2.
Note: The internal hyperlink only works on GitHub Pages or nbviewer. Do not click when viewing the notebook on GitHub.
In [42]:
import numpy as np
import matplotlib.pyplot as plt
from sunraster.instr.spice import read_spice_l2_fits
import h5py
import sunpy
import sunpy.map
from sharpesst.correct_2d_psf import get_fwd_matrices, correct_spice_raster
from sharpesst.util import bindown, as_dict, get_iris_data, masked_median_filter, get_mask_errs
from sharpesst.fit_spice_lines import get_overall_center, fit_spice_lines as fsl
import astropy
from astropy.visualization import (ImageNormalize, AsinhStretch)
from astropy import constants as const
import astropy.units as u
from astropy.io import fits
from astropy.wcs import WCS
import juanfit
import importlib
importlib.reload(juanfit)
from juanfit import SpectrumFit2D, gaussian, SpectrumFitSingle, SpectrumFitRow
from scipy.signal import find_peaks
from sospice import spice_error
from scipy.optimize import curve_fit
import os
import itertools
from concurrent.futures import ProcessPoolExecutor
import concurrent
import concurrent.futures
import logging
from matplotlib.patches import Rectangle, Circle, Ellipse
from astropy.convolution import interpolate_replace_nans, Gaussian2DKernel
from fancy_colorbar import plot_colorbar
import matplotlib.patheffects as path_effects
import matplotlib.transforms as transforms
from scipy.io import readsav
from copy import deepcopy
from glob import glob
import sys
sys.path.append("/cluster/home/zhuyin/scripts/spice-line-fits/linefit_modules/")
from skew_correction import full_correction
from skew_parameter_search import detrend_dopp
In [2]:
spice_file = "/cluster/home/zhuyin/Solar/SPICE_psf_202210/src/solo_L2_spice-n-ras_20221024T231535_V22_150995398-000.fits"
with fits.open(spice_file) as hdul:
# hdul.info()
spice_NeVIII_data = hdul[3].data[0].copy()
spice_NeVIII_header = hdul[3].header.copy()
spice_NeVIII_data[spice_NeVIII_data < -0] = np.nan
spice_window_mask = np.isfinite(spice_NeVIII_data[:,120:700,:]).any(axis=(1,2))
spice_NeVIII_wvl = (np.arange(1, spice_NeVIII_header["NAXIS3"] + 1) - spice_NeVIII_header["CRPIX3"])*\
spice_NeVIII_header["CDELT3"] + spice_NeVIII_header["CRVAL3"]
spice_NeVIII_wvl = spice_NeVIII_wvl[spice_window_mask]*10 # nm to AA
start_wvl_index, end_wvl_index = np.where(spice_window_mask == True)[0][np.r_[0,-1]]
spice_NeVIII_header_new_naxis3 = np.sum(spice_window_mask)
spice_NeVIII_header_new_crpix3 = (1 + spice_NeVIII_header_new_naxis3)/2.
spice_NeVIII_header_new_crval3 = ((start_wvl_index + end_wvl_index)/2. + 1 - spice_NeVIII_header['CRPIX3'])*spice_NeVIII_header['CDELT3'] + \
spice_NeVIII_header['CRVAL3']
spice_NeVIII_header["NAXIS3"] = spice_NeVIII_header_new_naxis3
spice_NeVIII_header["CRPIX3"] = spice_NeVIII_header_new_crpix3
spice_NeVIII_header["CRVAL3"] = spice_NeVIII_header_new_crval3
spice_NeVIII_data_bg = deepcopy(spice_NeVIII_data)
spice_NeVIII_data_bg[spice_NeVIII_data_bg > np.nanpercentile(spice_NeVIII_data, 15, axis=0)[np.newaxis,:,:]] = np.nan
spice_NeVIII_data_bg = np.nanmedian(spice_NeVIII_data_bg, axis=0)
spice_NeVIII_data_bg_rm = spice_NeVIII_data - spice_NeVIII_data_bg[np.newaxis,:,:] #np.nanmin(spice_NeVIII_data, axis=0)[np.newaxis,:,:]
spice_NeVIII_data_bg_rm = spice_NeVIII_data_bg_rm[spice_window_mask,120:700,:]
# kernel = Gaussian2DKernel(3)
# for ii in range(spice_NeVIII_data_bg_rm.shape[2]):
# spice_NeVIII_data_bg_rm[:, :, ii] = interpolate_replace_nans(spice_NeVIII_data_bg_rm[:, :, ii], kernel)
rebin_facs = [1,2,1]
plot_aspect = spice_NeVIII_header['CDELT2']/spice_NeVIII_header['CDELT1']*rebin_facs[1]/rebin_facs[2]
spice_NeVIII_data_bin = bindown(spice_NeVIII_data_bg_rm,np.round(np.array(spice_NeVIII_data_bg_rm.shape)/rebin_facs).astype(np.int32))
spice_NeVIII_window_start = np.argmax(np.isfinite(spice_NeVIII_data_bin), axis=0)
spice_NeVIII_window_end = spice_NeVIII_data_bin.shape[0] - 1 - np.argmax(np.isfinite(spice_NeVIII_data_bin[::-1,:,:]), axis=0)
spice_NeVIII_data_bin[19:21,200:202,:] = np.nan # mask some hot pixels
# spice_NeVIII_window_start = np.zeros(spice_NeVIII_data_bin.shape[1:]).astype(np.int32)
# spice_NeVIII_window_end = (np.ones(spice_NeVIII_data_bin.shape[1:])*(spice_NeVIII_data_bin.shape[0] - 1)).astype(np.int32)
/cluster/software/stacks/2024-05/python-cuda/3.11.6/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1563: RuntimeWarning: All-NaN slice encountered return function_base._ureduce(a, /scratch/tmp.29624082.zhuyin/ipykernel_1348293/3596438582.py:29: RuntimeWarning: All-NaN slice encountered spice_NeVIII_data_bg = np.nanmedian(spice_NeVIII_data_bg, axis=0)
In [3]:
fig, ax = plt.subplots(figsize=(4,6), layout="constrained")
ax.imshow(np.nansum(spice_NeVIII_data_bin, axis=0), origin="lower", aspect=plot_aspect,
norm=ImageNormalize(vmin=np.nanpercentile(np.nansum(spice_NeVIII_data_bin, axis=0), 0.5),
vmax=np.nanpercentile(np.nansum(spice_NeVIII_data_bin, axis=0), 99.5),
stretch=AsinhStretch()))
rect_box_1 = Rectangle((30, 120), 60, 120, linewidth=1, edgecolor='r', facecolor='none')
ax.add_patch(rect_box_1)
Out[3]:
<matplotlib.patches.Rectangle at 0x14ba6f259910>
In [4]:
spice_NeVIII_data_bin_r1 = spice_NeVIII_data_bin[:,120:240,30:90]
In [5]:
fig, ax = plt.subplots(figsize=(3,3), layout="constrained")
ax.imshow(np.nansum(spice_NeVIII_data_bin_r1, axis=0), origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(np.nansum(spice_NeVIII_data_bin_r1, axis=0), 0.5),
vmax=np.nanpercentile(np.nansum(spice_NeVIII_data_bin_r1, axis=0), 99.5),
stretch=AsinhStretch()), aspect=plot_aspect)
Out[5]:
<matplotlib.image.AxesImage at 0x14ba6012d990>
In [6]:
yl_core_xpo = 1.8
yl_wing_xpo = 1
# Rotation angle of the PSF, both core and wings
psf_yl_angle = -25*np.pi/180
# FWHMs of PSF core. First argument is width along y axis before rotation,
# and is in arcseconds. Second is along lambda axis and is in angstrom.
fwhm_core0_yl = np.array([2.5, 1.11])
# This descriptor for plots should be manually edited to reflect the PSF parameters
gaussian_desc = '2-part Gaussian PSF'
fwhm_wing0_yl = np.array([13.0, 4.5]) # FWHMs of PSF wings in arcseconds and angstroms, respectively
desc_str='; standard wing aspect ratio'
# Fraction of overall PSF amplitude in wings (core weight is 1.0 - wing_weight).
# PSFs have unit peak amplitude, -- PLEASE NOTE: they do not integrate to 1.
wing_weight = 0.28
In [7]:
spice_NeVIII_corr_data_r1, spice_NeVIII_cor_chi2s_r1, metadict_r1 = \
correct_spice_raster(spice_NeVIII_data_bin_r1.transpose([2,1,0]), spice_NeVIII_header, fwhm_core0_yl, fwhm_wing0_yl,
psf_yl_angle, wing_weight, yl_core_xpo=yl_core_xpo,
yl_wing_xpo=yl_wing_xpo,
super_fac=1, chi2_th = 0.5,
psf_thold_core=0.0005, spice_bin_facs = rebin_facs)
spice_NeVIII_corr_data_r1 = spice_NeVIII_corr_data_r1.transpose([2,1,0])
Correcting Ne VIII 770 - Peak; ref. wavelength=768.456459599 Computing PSF Core: 4.980204209210252 % done after 0.4538135528564453 seconds 9.98124609293603 % done after 0.9062261581420898 seconds 14.982287976661805 % done after 1.3571484088897705 seconds 19.98332986038758 % done after 1.8096771240234375 seconds 24.984371744113357 % done after 2.261590003967285 seconds 29.985413627839133 % done after 2.7141690254211426 seconds 34.98645551156491 % done after 3.1659111976623535 seconds 39.98749739529069 % done after 3.617821216583252 seconds 44.988539279016464 % done after 4.070712566375732 seconds 49.98958116274224 % done after 4.522996664047241 seconds 54.990623046468016 % done after 4.975438117980957 seconds 59.99166493019379 % done after 5.428099155426025 seconds 64.99270681391957 % done after 5.8819074630737305 seconds 69.99374869764534 % done after 6.333714962005615 seconds 74.99479058137112 % done after 6.785520553588867 seconds 79.9958324650969 % done after 7.236984014511108 seconds 84.99687434882267 % done after 7.6898064613342285 seconds 89.99791623254845 % done after 8.143182516098022 seconds 94.99895811627422 % done after 8.595542192459106 seconds 100.0 % done after 9.048874378204346 seconds Computing PSF Wings: 4.980204209210252 % done after 1.031674861907959 seconds 9.98124609293603 % done after 2.0659804344177246 seconds 14.982287976661805 % done after 3.1016266345977783 seconds 19.98332986038758 % done after 4.133039474487305 seconds 24.984371744113357 % done after 5.164168357849121 seconds 29.985413627839133 % done after 6.199804067611694 seconds 34.98645551156491 % done after 7.234876394271851 seconds 39.98749739529069 % done after 8.267709255218506 seconds 44.988539279016464 % done after 9.299019813537598 seconds 49.98958116274224 % done after 10.337345838546753 seconds 54.990623046468016 % done after 11.370586156845093 seconds 59.99166493019379 % done after 12.402878284454346 seconds 64.99270681391957 % done after 13.434340953826904 seconds 69.99374869764534 % done after 14.469213724136353 seconds 74.99479058137112 % done after 15.502605438232422 seconds 79.9958324650969 % done after 16.534666061401367 seconds 84.99687434882267 % done after 17.575329542160034 seconds 89.99791623254845 % done after 18.607362270355225 seconds 94.99895811627422 % done after 20.025184154510498 seconds 100.0 % done after 21.053448915481567 seconds Computing Symmetric PSF for output: 4.980204209210252 % done after 0.3137040138244629 seconds 9.98124609293603 % done after 0.6276655197143555 seconds 14.982287976661805 % done after 0.9405100345611572 seconds 19.98332986038758 % done after 1.2505545616149902 seconds 24.984371744113357 % done after 1.5593297481536865 seconds 29.985413627839133 % done after 1.8690383434295654 seconds 34.98645551156491 % done after 2.179258108139038 seconds 39.98749739529069 % done after 2.489119529724121 seconds 44.988539279016464 % done after 2.7992663383483887 seconds 49.98958116274224 % done after 3.1098756790161133 seconds 54.990623046468016 % done after 3.420069694519043 seconds 59.99166493019379 % done after 3.7300686836242676 seconds 64.99270681391957 % done after 4.040565013885498 seconds 69.99374869764534 % done after 4.350973129272461 seconds 74.99479058137112 % done after 4.66081166267395 seconds 79.9958324650969 % done after 4.971075534820557 seconds 84.99687434882267 % done after 5.2817604541778564 seconds 89.99791623254845 % done after 5.594995498657227 seconds 94.99895811627422 % done after 5.905846357345581 seconds 100.0 % done after 6.216505527496338 seconds
2025-04-17 14:04:24 - numba.core.transforms - INFO: finding looplift candidates
/cluster/home/zhuyin/scripts/SHARPESST/sharpesst/util.py:83: NumbaWarning:
Compilation is falling back to object mode WITHOUT looplifting enabled because Function "_masked_medfilt_inner" failed type inference due to: Use of unsupported NumPy function 'numpy.unravel_index' or unsupported use of the function.
File "../../../scripts/SHARPESST/sharpesst/util.py", line 86:
def _masked_medfilt_inner(flatinds,data,footprint_inds,footprint_ind_offset,tparg,dat_pad_shape,data_fppad,mask_fppad,data_filt):
<source elided>
for ind in flatinds:
ijkpad = np.unravel_index(ind,data.shape)+footprint_inds-footprint_ind_offset
^
During: typing of get attribute at /cluster/home/zhuyin/scripts/SHARPESST/sharpesst/util.py (86)
File "../../../scripts/SHARPESST/sharpesst/util.py", line 86:
def _masked_medfilt_inner(flatinds,data,footprint_inds,footprint_ind_offset,tparg,dat_pad_shape,data_fppad,mask_fppad,data_filt):
<source elided>
for ind in flatinds:
ijkpad = np.unravel_index(ind,data.shape)+footprint_inds-footprint_ind_offset
^
@numba.jit(fastmath=True, parallel=True, forceobj=True)
Correcting the PSF: 0.4248075485229492 s; Finished slit index 0 of 60 , chi^2= 0.41305386910886277 0.8257670402526855 s; Finished slit index 1 of 60 , chi^2= 0.40159536301564863 1.219766616821289 s; Finished slit index 2 of 60 , chi^2= 0.3696127490735742 1.5971400737762451 s; Finished slit index 3 of 60 , chi^2= 0.3457252122173611 1.984438419342041 s; Finished slit index 4 of 60 , chi^2= 0.3632225035924998 2.383026361465454 s; Finished slit index 5 of 60 , chi^2= 0.4101936083216352 2.7740252017974854 s; Finished slit index 6 of 60 , chi^2= 0.3240936625944776 3.1689224243164062 s; Finished slit index 7 of 60 , chi^2= 0.3725659650738048 3.5648458003997803 s; Finished slit index 8 of 60 , chi^2= 0.36333128794258035 3.9620156288146973 s; Finished slit index 9 of 60 , chi^2= 0.36270863028290384 4.3638646602630615 s; Finished slit index 10 of 60 , chi^2= 0.3134356534586687 4.772308826446533 s; Finished slit index 11 of 60 , chi^2= 0.36938956486080665 5.177753686904907 s; Finished slit index 12 of 60 , chi^2= 0.3514965649548807 5.581541061401367 s; Finished slit index 13 of 60 , chi^2= 0.41476186108424307 6.554769992828369 s; Finished slit index 14 of 60 , chi^2= 0.26219304584231334 7.953392267227173 s; Finished slit index 15 of 60 , chi^2= 0.3349801046966764 9.196560621261597 s; Finished slit index 16 of 60 , chi^2= 0.3525014105250098 10.308531999588013 s; Finished slit index 17 of 60 , chi^2= 0.367622504781808 11.406569957733154 s; Finished slit index 18 of 60 , chi^2= 0.36935886695806686 12.569491863250732 s; Finished slit index 19 of 60 , chi^2= 0.41143251536729153 13.705390691757202 s; Finished slit index 20 of 60 , chi^2= 0.41729297111671304 14.793606996536255 s; Finished slit index 21 of 60 , chi^2= 0.4278929952700844 16.440433502197266 s; Finished slit index 22 of 60 , chi^2= 0.44952353158021796 18.942987203598022 s; Finished slit index 23 of 60 , chi^2= 0.4880126267270865 21.54291534423828 s; Finished slit index 24 of 60 , chi^2= 0.48023270936810314 22.586934328079224 s; Finished slit index 25 of 60 , chi^2= 0.4589220716019457 24.296475887298584 s; Finished slit index 26 of 60 , chi^2= 0.49024621067041757 26.64468002319336 s; Finished slit index 27 of 60 , chi^2= 0.47612800344834616 28.500867128372192 s; Finished slit index 28 of 60 , chi^2= 0.47073016664530765 30.503663778305054 s; Finished slit index 29 of 60 , chi^2= 0.4751191535039403 32.560588121414185 s; Finished slit index 30 of 60 , chi^2= 0.4825207706176365 36.22824573516846 s; Finished slit index 31 of 60 , chi^2= 0.48690789333729473 38.18460726737976 s; Finished slit index 32 of 60 , chi^2= 0.4695902499238668 39.285648822784424 s; Finished slit index 33 of 60 , chi^2= 0.4924975525427254 40.49050211906433 s; Finished slit index 34 of 60 , chi^2= 0.4755146712470065 41.611255168914795 s; Finished slit index 35 of 60 , chi^2= 0.4795503378722582 42.66443610191345 s; Finished slit index 36 of 60 , chi^2= 0.448530594415842 44.71842288970947 s; Finished slit index 37 of 60 , chi^2= 0.4393934940944445 46.717400312423706 s; Finished slit index 38 of 60 , chi^2= 0.4541722572505919 48.84012961387634 s; Finished slit index 39 of 60 , chi^2= 0.4331087951341688 50.048972845077515 s; Finished slit index 40 of 60 , chi^2= 0.4313734562258759 52.050524950027466 s; Finished slit index 41 of 60 , chi^2= 0.47805088661201806 56.60388422012329 s; Finished slit index 42 of 60 , chi^2= 0.4826330252235942 63.10903453826904 s; Finished slit index 43 of 60 , chi^2= 0.4745617040964823 68.0278902053833 s; Finished slit index 44 of 60 , chi^2= 0.4971314811167409 74.89669513702393 s; Finished slit index 45 of 60 , chi^2= 0.48568853352714914 79.15765714645386 s; Finished slit index 46 of 60 , chi^2= 0.4963835813070636 81.01180291175842 s; Finished slit index 47 of 60 , chi^2= 0.4930947191766559 84.2812750339508 s; Finished slit index 48 of 60 , chi^2= 0.48895048940399477 89.70541286468506 s; Finished slit index 49 of 60 , chi^2= 0.4917735250522421 97.63153576850891 s; Finished slit index 50 of 60 , chi^2= 0.49844262606925593 108.3949384689331 s; Finished slit index 51 of 60 , chi^2= 0.4910862479491528 117.33379340171814 s; Finished slit index 52 of 60 , chi^2= 0.49792775919740007 128.46971487998962 s; Finished slit index 53 of 60 , chi^2= 0.490722857864807 137.6218445301056 s; Finished slit index 54 of 60 , chi^2= 0.48635675721189403 144.43010425567627 s; Finished slit index 55 of 60 , chi^2= 0.48959879968515385 152.12741231918335 s; Finished slit index 56 of 60 , chi^2= 0.4848052171869965 158.8964433670044 s; Finished slit index 57 of 60 , chi^2= 0.48521520092683473 167.46382355690002 s; Finished slit index 58 of 60 , chi^2= 0.4789426634436249 179.80132412910461 s; Finished slit index 59 of 60 , chi^2= 0.4956294252736102
In [8]:
def func_NeVIII_770_sg(x, wvl_770, int_770, wid_770, b):
return gaussian(x, wvl_770, int_770, wid_770) + \
b
In [9]:
def work_fit_single_spec(args):
(NeVIII_770_wvl, NeVIII_770_data_rebin,
NeVIII_770_window_start, NeVIII_window_end,
ii, jj) = args
try:
data_slice = NeVIII_770_data_rebin[NeVIII_770_window_start[ii,jj]:NeVIII_window_end[ii,jj]+1,
ii, jj]
p0 = [
NeVIII_770_wvl[np.nanargmax(data_slice)], # wavelength of max
np.nanmax(data_slice), # intensity
0.8, # width
0 # background
]
# Bounds for parameters
bounds = (
(0, 0, 0, -np.inf),
(np.inf, np.inf, np.inf, np.inf)
)
# Perform curve fitting
popt, pcov = curve_fit(
func_NeVIII_770_sg,
NeVIII_770_wvl[NeVIII_770_window_start[ii,jj]:NeVIII_window_end[ii,jj]+1],
data_slice,
p0=p0,
bounds=bounds,
nan_policy="omit"
)
perr = np.sqrt(np.diag(pcov))
if np.any(~np.isfinite(popt)) or np.any(~np.isfinite(perr)) \
or np.any(np.abs(perr[:3]/popt[:3]) > 0.3):
return (ii, jj, None)
else:
return (ii, jj, popt)
except Exception as e:
# Return None or a special marker for failed fits
return (ii, jj, None)
In [10]:
def fit_parallel(NeVIII_770_wvl, NeVIII_770_data_rebin, NeVIII_770_window_start=spice_NeVIII_window_start,
NeVIII_770_window_end=spice_NeVIII_window_end, ncpu=8):
# Prepare input arguments for each position
input_args = [
(NeVIII_770_wvl, NeVIII_770_data_rebin,
NeVIII_770_window_start, NeVIII_770_window_end,
ii, jj)
for ii in range(NeVIII_770_data_rebin.shape[1])
for jj in range(NeVIII_770_data_rebin.shape[2])
]
# Prepare result arrays
shape = NeVIII_770_data_rebin.shape[1:]
NeVIII_fit_wvl = np.ones(shape)*np.nan
NeVIII_fit_int = np.ones(shape)*np.nan
NeVIII_fit_wid = np.ones(shape)*np.nan
NeVIII_fit_b = np.ones(shape)*np.nan
# Use ProcessPoolExecutor to parallelize fitting
with ProcessPoolExecutor(max_workers=ncpu) as executor:
# Submit all tasks
futures = [executor.submit(work_fit_single_spec, args) for args in input_args]
# Process results as they complete
for future in concurrent.futures.as_completed(futures):
result = future.result()
# Check if fitting was successful
if result[2] is not None:
ii, jj, popt = result
# Store results in respective arrays
NeVIII_fit_wvl[ii, jj] = popt[0]
NeVIII_fit_int[ii, jj] = popt[1]
NeVIII_fit_wid[ii, jj] = popt[2]
NeVIII_fit_b[ii, jj] = popt[3]
return (
NeVIII_fit_wvl,
NeVIII_fit_int,
NeVIII_fit_wid,
NeVIII_fit_b
)
In [11]:
(NeVIII_orig_fit_wvl_r1, NeVIII_orig_fit_int_r1, NeVIII_orig_fit_wid_r1,
NeVIII_orig_fit_b_r1) = fit_parallel(spice_NeVIII_wvl,spice_NeVIII_data_bin_r1)
In [12]:
(NeVIII_corr_fit_wvl_r1, NeVIII_corr_fit_int_r1, NeVIII_corr_fit_wid_r1,
NeVIII_corr_fit_b_r1) = fit_parallel(spice_NeVIII_wvl,spice_NeVIII_corr_data_r1)
In [13]:
fig, ((ax1,ax2,ax3), (ax4,ax5,ax6)) = plt.subplots(2,3,figsize=(8,6), layout="constrained")
ax1.imshow(NeVIII_orig_fit_int_r1, origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_orig_fit_int_r1, 0.1),
vmax=np.nanpercentile(NeVIII_orig_fit_int_r1, 99.5),
stretch=AsinhStretch(0.1)), aspect=plot_aspect, interpolation='none')
ax2.imshow((NeVIII_orig_fit_wvl_r1 - np.nanmedian(NeVIII_orig_fit_wvl_r1, axis=0)[np.newaxis,:])/770*3e5,
norm=ImageNormalize(vmin=-40, vmax=40), cmap="bwr", origin="lower",
aspect=plot_aspect, interpolation='none')
ax3.imshow(NeVIII_orig_fit_wid_r1, norm=ImageNormalize(vmin=0.8, vmax=0.9),
origin="lower", aspect=plot_aspect, interpolation='none')
ax4.imshow(NeVIII_corr_fit_int_r1, origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_corr_fit_int_r1, 0.5),
vmax=np.nanpercentile(NeVIII_corr_fit_int_r1, 99.5),
stretch=AsinhStretch(0.1)), aspect=plot_aspect, interpolation='none')
ax5.imshow((NeVIII_corr_fit_wvl_r1 - np.nanmedian(NeVIII_corr_fit_wvl_r1, axis=0)[np.newaxis,:])/770*3e5,
norm=ImageNormalize(vmin=-40, vmax=40), cmap="bwr", origin="lower",
aspect=plot_aspect, interpolation='none')
ax6.imshow(NeVIII_corr_fit_wid_r1, norm=ImageNormalize(vmin=0.3, vmax=0.6),
origin="lower", aspect=plot_aspect, interpolation='none')
Out[13]:
<matplotlib.image.AxesImage at 0x14ba6dec1d90>
In [37]:
with np.load("/cluster/home/zhuyin/work/spice_psf/spice_skew_20221024/spice_20221024T231535_NeVIII_skew_fit_params.npz", allow_pickle=True) as f:
spice_skew_sav = dict(f)
In [41]:
spice_skew_sav['linefits'].item()[0][0]['amplitudes'].data.shape
Out[41]:
(192, 832, 1, 1)
In [15]:
spice_skew_dat = spice_skew_sav['dat_skew'].transpose(2,1,0)
spice_skew_dat_r1 = spice_skew_dat[:, 360:360+240, 30:90]
spice_skew_dat_r1 = np.nanmean(spice_skew_dat_r1.reshape(50, 120, 2, 60), axis=2)
/scratch/tmp.29624082.zhuyin/ipykernel_1348293/1287099549.py:3: RuntimeWarning: Mean of empty slice spice_skew_dat_r1 = np.nanmean(spice_skew_dat_r1.reshape(50, 120, 2, 60), axis=2)
In [43]:
NeVIII_skew_fit_int_r1_load = spice_skew_sav['linefits'].item()[0][0]['amplitudes'].data[30:90,360:360+240,0,0].T
NeVIII_skew_fit_wvl_r1_load = detrend_dopp(spice_skew_sav['linefits'].item()[0][0]['centers'])[30:90,360:360+240].T
In [16]:
(NeVIII_skew_fit_wvl_r1, NeVIII_skew_fit_int_r1, NeVIII_skew_fit_wid_r1,
NeVIII_skew_fit_b_r1) = fit_parallel(spice_skew_sav['waves'],spice_skew_dat_r1)
In [17]:
fig, ((ax1,ax2,ax3), (ax4,ax5,ax6),
(ax7,ax8,ax9)) = plt.subplots(3,3,figsize=(8,8), layout="constrained")
ax1.imshow(NeVIII_orig_fit_int_r1, origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_orig_fit_int_r1, 0.1),
vmax=np.nanpercentile(NeVIII_orig_fit_int_r1, 99.5),
stretch=AsinhStretch(0.1)), aspect=plot_aspect, interpolation='none')
ax2.imshow((NeVIII_orig_fit_wvl_r1 - np.nanmedian(NeVIII_orig_fit_wvl_r1, axis=0)[np.newaxis,:])/770*3e5,
norm=ImageNormalize(vmin=-40, vmax=40), cmap="bwr", origin="lower",
aspect=plot_aspect, interpolation='none')
ax3.imshow(NeVIII_orig_fit_wid_r1, norm=ImageNormalize(vmin=0.8, vmax=0.9),
origin="lower", aspect=plot_aspect, interpolation='none')
ax4.imshow(NeVIII_corr_fit_int_r1, origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_corr_fit_int_r1, 0.5),
vmax=np.nanpercentile(NeVIII_corr_fit_int_r1, 99.5),
stretch=AsinhStretch(0.1)), aspect=plot_aspect, interpolation='none')
ax5.imshow((NeVIII_corr_fit_wvl_r1 - np.nanmedian(NeVIII_corr_fit_wvl_r1, axis=0)[np.newaxis,:])/770*3e5,
norm=ImageNormalize(vmin=-40, vmax=40), cmap="bwr", origin="lower",
aspect=plot_aspect, interpolation='none')
ax6.imshow(NeVIII_corr_fit_wid_r1, norm=ImageNormalize(vmin=0.3, vmax=0.6),
origin="lower", aspect=plot_aspect, interpolation='none')
ax7.imshow(NeVIII_skew_fit_int_r1, origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_skew_fit_int_r1, 0.5),
vmax=np.nanpercentile(NeVIII_skew_fit_int_r1, 99.5),
stretch=AsinhStretch(0.1)), aspect=plot_aspect, interpolation='none')
ax8.imshow((NeVIII_skew_fit_wvl_r1 - np.nanmedian(NeVIII_skew_fit_wvl_r1, axis=0)[np.newaxis,:])/770*3e5,
norm=ImageNormalize(vmin=-40, vmax=40), cmap="bwr", origin="lower",
aspect=plot_aspect, interpolation='none')
ax9.imshow(NeVIII_skew_fit_wid_r1, norm=ImageNormalize(vmin=0.6, vmax=0.9),
origin="lower", aspect=plot_aspect, interpolation='none')
Out[17]:
<matplotlib.image.AxesImage at 0x14ba8cb25c90>
In [18]:
ms_style_dict = {'text.usetex': True, 'font.family': 'serif', 'axes.linewidth': 1.2,
'xtick.major.width': 1.2, 'xtick.major.size': 4,
'ytick.major.width': 1.2, 'ytick.major.size': 4,
'xtick.minor.width': 1.2, 'xtick.minor.size': 2,
'ytick.minor.width': 1.2, 'ytick.minor.size': 2,
'xtick.direction': 'in', 'ytick.direction': 'in',
'text.latex.preamble': r'\usepackage[T1]{fontenc}'
r'\usepackage{amsmath}' r'\usepackage{siunitx}'
r'\sisetup{detect-all=True}'}
In [45]:
with plt.rc_context(ms_style_dict):
example_x_index = 51
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6),
(ax7,ax8,ax9)) = plt.subplots(3,3,figsize=(7,9), layout="constrained",
gridspec_kw={'width_ratios': [1,1,0.55]})
im1 = ax1.imshow(NeVIII_orig_fit_int_r1, origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_orig_fit_int_r1, 0.1),
vmax=np.nanpercentile(NeVIII_orig_fit_int_r1, 99.9),
stretch=AsinhStretch(0.1)), aspect=plot_aspect, interpolation="none",
cmap="sdoaia171")
clb1, clb_ax1 = plot_colorbar(im1, ax1, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
clb_ax1.set_title(r"Ne\,\textsc{viii} 77.04\,nm $I_{\rm tot}$" + "\n" + \
r"$\mathrm{\left(W\,m^{-2}\,sr^{-1}\right)}$", fontsize=10)
im2 = ax2.imshow((NeVIII_orig_fit_wvl_r1 - np.nanmedian(NeVIII_orig_fit_wvl_r1, axis=0)[np.newaxis,:])/770.428*const.c.to_value(u.km/u.s),
norm=ImageNormalize(vmin=-40, vmax=40), cmap="coolwarm", origin="lower",
aspect=plot_aspect, interpolation="none")
clb2, clb_ax2 = plot_colorbar(im2, ax2, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
clb_ax2.set_title(r"Ne\,\textsc{viii} 77.04\,nm $v_{\rm Doppler}$" + "\n" + \
r"$\mathrm{\left(km\,s^{-1}\right)}$", fontsize=10)
im3 = ax3.imshow(spice_NeVIII_data_bin_r1[:, :, example_x_index].T, origin="lower", aspect="auto",
norm=ImageNormalize(vmin=np.nanpercentile(spice_NeVIII_data_bin_r1[:, :, example_x_index],0.5),
vmax=np.nanpercentile(spice_NeVIII_data_bin_r1[:, :, example_x_index],99.9),
stretch=AsinhStretch(0.3)), interpolation="none", cmap="magma_r")
clb3, clb_ax3 = plot_colorbar(im3, ax3, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
clb_ax3.set_title(r"Ne\,\textsc{viii} $I_{\lambda}$" + "\n" + \
r"$\mathrm{\left(W\,m^{-2}\,sr^{-1}\,nm^{-1}\right)}$", fontsize=10)
im4 = ax4.imshow(NeVIII_corr_fit_int_r1, origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_corr_fit_int_r1, 0.5),
vmax=np.nanpercentile(NeVIII_corr_fit_int_r1, 99.9),
stretch=AsinhStretch(0.1)), aspect=plot_aspect, interpolation="none",
cmap="sdoaia171")
clb4, clb_ax4 = plot_colorbar(im4, ax4, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
im5 = ax5.imshow((NeVIII_corr_fit_wvl_r1 - np.nanmedian(NeVIII_corr_fit_wvl_r1, axis=0)[np.newaxis,:])/770.428*const.c.to_value(u.km/u.s),
norm=ImageNormalize(vmin=-40, vmax=40), cmap="coolwarm", origin="lower",
aspect=plot_aspect, interpolation="none")
clb5, clb_ax5 = plot_colorbar(im5, ax5, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
im6 = ax6.imshow(spice_NeVIII_corr_data_r1[:, :, example_x_index].T, origin="lower", aspect="auto",
norm=ImageNormalize(vmin=np.nanpercentile(spice_NeVIII_corr_data_r1[:, :, example_x_index],0.1),
vmax=np.nanpercentile(spice_NeVIII_corr_data_r1[:, :, example_x_index],99.9),
stretch=AsinhStretch(0.3)), interpolation="none", cmap="magma_r")
clb6, clb_ax6 = plot_colorbar(im6, ax6, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
im7 = ax7.imshow(NeVIII_skew_fit_int_r1, origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_skew_fit_int_r1, 0.5),
vmax=np.nanpercentile(NeVIII_skew_fit_int_r1, 99.9),
stretch=AsinhStretch(0.1)), aspect=plot_aspect, interpolation="none",
cmap="sdoaia171")
# im7 = ax7.imshow(NeVIII_skew_fit_int_r1_load, origin="lower",
# norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_skew_fit_int_r1, 0.5),
# vmax=np.nanpercentile(NeVIII_skew_fit_int_r1, 99.9),
# stretch=AsinhStretch(0.1)), aspect=plot_aspect/2, interpolation="none",
# cmap="sdoaia171")
clb7, clb_ax7 = plot_colorbar(im7, ax7, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
im8 = ax8.imshow((NeVIII_skew_fit_wvl_r1 - np.nanmedian(NeVIII_skew_fit_wvl_r1, axis=0)[np.newaxis,:])/770.428*const.c.to_value(u.km/u.s),
norm=ImageNormalize(vmin=-40, vmax=40), cmap="coolwarm", origin="lower",
aspect=plot_aspect, interpolation="none")
# im8 = ax8.imshow((NeVIII_skew_fit_wvl_r1_load - np.nanmedian(NeVIII_skew_fit_wvl_r1_load))/770.428*const.c.to_value(u.km/u.s),
# norm=ImageNormalize(vmin=-40, vmax=40), cmap="coolwarm", origin="lower",
# aspect=plot_aspect/2, interpolation="none")
clb8, clb_ax8 = plot_colorbar(im8, ax8, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
im9 = ax9.imshow(spice_skew_dat_r1[:, :, example_x_index].T, origin="lower", aspect="auto",
norm=ImageNormalize(vmin=np.nanpercentile(spice_skew_dat_r1[:, :, example_x_index],0.1),
vmax=np.nanpercentile(spice_skew_dat_r1[:, :, example_x_index],99.9),
stretch=AsinhStretch(0.3)), interpolation="none", cmap="magma_r")
clb9, clb_ax9 = plot_colorbar(im9, ax9, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
ax1.tick_params(labelbottom=False, top=True, right=True)
ax2.tick_params(labelleft=False, labelbottom=False, top=True, right=True)
ax3.tick_params(labelleft=False, labelbottom=False, top=True, right=True)
ax4.tick_params(top=True, right=True)
ax5.tick_params(labelleft=False, top=True, right=True)
ax6.tick_params(labelleft=False, top=True, right=True)
ax9.tick_params(labelleft=False, top=True, right=True)
ax1.set_ylabel(r"Detector-$y$ Bin$\times$2 (Pixel)")
ax4.set_ylabel(r"Detector-$y$ Bin$\times$2 (Pixel)")
ax7.set_ylabel(r"Detector-$y$ Bin$\times$2 (Pixel)")
ax7.set_xlabel(r"Raster Step")
ax8.set_xlabel(r"Raster Step")
ax9.set_xlabel(r"Detector-$\lambda$ (Pixel)")
for clb_ax_ in (clb_ax1, clb_ax2, clb_ax3, clb_ax4, clb_ax5, clb_ax6,
clb_ax7, clb_ax8, clb_ax9):
clb_ax_.tick_params(top=True, labeltop=True, labelbottom=False, bottom=False, which="both")
for ax_ in (ax1,ax2,ax4,ax5,ax7,ax8):
ax_.axvline(example_x_index, color="w", lw=2, ls="--", alpha=0.7)
for ax_, text_ in zip((ax1,ax2,ax4,ax5,ax7,ax8),
("a", "b", "d", "e", "g", "h")):
ax_.text(0.03,0.03,r"\textbf{{{:})}}".format(text_), transform=ax_.transAxes,
va="bottom", ha="left", color="w", fontsize=10,
path_effects=[path_effects.Stroke(linewidth=1.5, foreground='black'),path_effects.Normal()])
for ax_, text_ in zip((ax3,ax6,ax9),("c","f","i")):
ax_.text(0.06,0.03,r"\textbf{{{:})}}".format(text_), transform=ax_.transAxes,
va="bottom", ha="left", color="w", fontsize=10,
path_effects=[path_effects.Stroke(linewidth=1.5, foreground='black'),path_effects.Normal()])
ax1.text(0.05,0.95,r"\textbf{Raw}", transform=ax1.transAxes,
va="top", ha="left", color="w", fontsize=12,
path_effects=[path_effects.Stroke(linewidth=1.5, foreground='black'),path_effects.Normal()])
ax4.text(0.05,0.95,r"\textbf{PSF Deconv.}", transform=ax4.transAxes,
va="top", ha="left", color="w", fontsize=12,
path_effects=[path_effects.Stroke(linewidth=1.5, foreground='black'),path_effects.Normal()])
ax7.text(0.05,0.95,r"\textbf{Skew Corr.}", transform=ax7.transAxes,
va="top", ha="left", color="w", fontsize=12,
path_effects=[path_effects.Stroke(linewidth=1.5, foreground='black'),path_effects.Normal()])
target_r1_param = {'xy':(50, 100), 'width': 12, 'height': 12/plot_aspect, 'angle':0,
'ec':"#1B813E", 'fill': False, 'lw':2}
target_r2_param = {'xy':(54, 37), 'width': 18, 'height': 5/plot_aspect, 'angle':60,
'ec':"#F75C2F", 'fill': False, 'lw':2}
# the patches.Ellipse cannot handle the rotation
# in axes with different scales so we need to
# plot it explicitly
r2_t_ = np.linspace(0, 2*np.pi, 360)
r2_x_ = target_r2_param['xy'][0] + \
target_r2_param['width']*np.cos(r2_t_)*np.cos(np.deg2rad(target_r2_param['angle'])) - \
target_r2_param['height']*plot_aspect*np.sin(r2_t_)*np.sin(np.deg2rad(target_r2_param['angle']))
r2_y_ = target_r2_param['xy'][1] + \
target_r2_param['width']*np.cos(r2_t_)*np.sin(np.deg2rad(target_r2_param['angle'])) + \
target_r2_param['height']*plot_aspect*np.sin(r2_t_)*np.cos(np.deg2rad(target_r2_param['angle']))
target_r3_param = {'xy':(49, 5), 'width': 12, 'height': 12/plot_aspect, 'angle':0,
'ec':"#E03C8A", 'fill': False, 'lw':2}
for ax_ in (ax1,ax2,ax4,ax5,ax7,ax8):
ax_lim_ = ax_.axis()
circle_r1_ = Ellipse(**target_r1_param)
circle_r3_ = Ellipse(**target_r3_param)
ax_.add_patch(circle_r1_)
ax_.add_patch(circle_r3_)
ax_.plot(r2_x_, r2_y_, lw=target_r2_param['lw'], color=target_r2_param['ec'])
ax_.axis(ax_lim_)
# fig.savefig("../figs/spice_psf_deconv.pdf", dpi=600, format="pdf", bbox_inches="tight")
plt.show()
In [50]:
with plt.rc_context(ms_style_dict):
example_x_index = 51
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6),
(ax7,ax8,ax9)) = plt.subplots(3,3,figsize=(7,9), layout="constrained",
gridspec_kw={'width_ratios': [1,1,0.55]})
im1 = ax1.imshow(NeVIII_orig_fit_int_r1, origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_orig_fit_int_r1, 0.1),
vmax=np.nanpercentile(NeVIII_orig_fit_int_r1, 99.9),
stretch=AsinhStretch(0.1)), aspect=plot_aspect, interpolation="none",
cmap="sdoaia171")
clb1, clb_ax1 = plot_colorbar(im1, ax1, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
clb_ax1.set_title(r"Ne\,\textsc{viii} 77.04\,nm $I_{\rm tot}$" + "\n" + \
r"$\mathrm{\left(W\,m^{-2}\,sr^{-1}\right)}$", fontsize=10)
im2 = ax2.imshow((NeVIII_orig_fit_wvl_r1 - np.nanmedian(NeVIII_orig_fit_wvl_r1, axis=0)[np.newaxis,:])/770.428*const.c.to_value(u.km/u.s),
norm=ImageNormalize(vmin=-40, vmax=40), cmap="coolwarm", origin="lower",
aspect=plot_aspect, interpolation="none")
clb2, clb_ax2 = plot_colorbar(im2, ax2, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
clb_ax2.set_title(r"Ne\,\textsc{viii} 77.04\,nm $v_{\rm Doppler}$" + "\n" + \
r"$\mathrm{\left(km\,s^{-1}\right)}$", fontsize=10)
im3 = ax3.imshow(spice_NeVIII_data_bin_r1[:, :, example_x_index].T, origin="lower", aspect="auto",
norm=ImageNormalize(vmin=np.nanpercentile(spice_NeVIII_data_bin_r1[:, :, example_x_index],0.5),
vmax=np.nanpercentile(spice_NeVIII_data_bin_r1[:, :, example_x_index],99.9),
stretch=AsinhStretch(0.3)), interpolation="none", cmap="magma_r")
clb3, clb_ax3 = plot_colorbar(im3, ax3, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
clb_ax3.set_title(r"Ne\,\textsc{viii} $I_{\lambda}$" + "\n" + \
r"$\mathrm{\left(W\,m^{-2}\,sr^{-1}\,nm^{-1}\right)}$", fontsize=10)
im4 = ax4.imshow(NeVIII_corr_fit_int_r1, origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_corr_fit_int_r1, 0.5),
vmax=np.nanpercentile(NeVIII_corr_fit_int_r1, 99.9),
stretch=AsinhStretch(0.1)), aspect=plot_aspect, interpolation="none",
cmap="sdoaia171")
clb4, clb_ax4 = plot_colorbar(im4, ax4, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
im5 = ax5.imshow((NeVIII_corr_fit_wvl_r1 - np.nanmedian(NeVIII_corr_fit_wvl_r1, axis=0)[np.newaxis,:])/770.428*const.c.to_value(u.km/u.s),
norm=ImageNormalize(vmin=-40, vmax=40), cmap="coolwarm", origin="lower",
aspect=plot_aspect, interpolation="none")
clb5, clb_ax5 = plot_colorbar(im5, ax5, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
im6 = ax6.imshow(spice_NeVIII_corr_data_r1[:, :, example_x_index].T, origin="lower", aspect="auto",
norm=ImageNormalize(vmin=np.nanpercentile(spice_NeVIII_corr_data_r1[:, :, example_x_index],0.1),
vmax=np.nanpercentile(spice_NeVIII_corr_data_r1[:, :, example_x_index],99.9),
stretch=AsinhStretch(0.3)), interpolation="none", cmap="magma_r")
clb6, clb_ax6 = plot_colorbar(im6, ax6, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
# im7 = ax7.imshow(NeVIII_skew_fit_int_r1, origin="lower",
# norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_skew_fit_int_r1, 0.5),
# vmax=np.nanpercentile(NeVIII_skew_fit_int_r1, 99.9),
# stretch=AsinhStretch(0.1)), aspect=plot_aspect, interpolation="none",
# cmap="sdoaia171")
im7 = ax7.imshow(NeVIII_skew_fit_int_r1_load, origin="lower",
norm=ImageNormalize(vmin=np.nanpercentile(NeVIII_skew_fit_int_r1, 0.5),
vmax=np.nanpercentile(NeVIII_skew_fit_int_r1, 99.9),
stretch=AsinhStretch(0.1)), aspect=plot_aspect/2, interpolation="none",
cmap="sdoaia171")
clb7, clb_ax7 = plot_colorbar(im7, ax7, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
# im8 = ax8.imshow((NeVIII_skew_fit_wvl_r1 - np.nanmedian(NeVIII_skew_fit_wvl_r1, axis=0)[np.newaxis,:])/770.428*const.c.to_value(u.km/u.s),
# norm=ImageNormalize(vmin=-40, vmax=40), cmap="coolwarm", origin="lower",
# aspect=plot_aspect, interpolation="none")
im8 = ax8.imshow((NeVIII_skew_fit_wvl_r1_load - np.nanmedian(NeVIII_skew_fit_wvl_r1_load))/770.428*const.c.to_value(u.km/u.s),
norm=ImageNormalize(vmin=-40, vmax=40), cmap="coolwarm", origin="lower",
aspect=plot_aspect/2, interpolation="none")
clb8, clb_ax8 = plot_colorbar(im8, ax8, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
im9 = ax9.imshow(spice_skew_dat_r1[:, :, example_x_index].T, origin="lower", aspect="auto",
norm=ImageNormalize(vmin=np.nanpercentile(spice_skew_dat_r1[:, :, example_x_index],0.1),
vmax=np.nanpercentile(spice_skew_dat_r1[:, :, example_x_index],99.9),
stretch=AsinhStretch(0.3)), interpolation="none", cmap="magma_r")
clb9, clb_ax9 = plot_colorbar(im9, ax9, bbox_to_anchor=(0.03, 1.02, 0.94, 0.08),
fontsize=10, orientation="horizontal")
ax1.tick_params(labelbottom=False, top=True, right=True)
ax2.tick_params(labelleft=False, labelbottom=False, top=True, right=True)
ax3.tick_params(labelleft=False, labelbottom=False, top=True, right=True)
ax4.tick_params(top=True, right=True)
ax5.tick_params(labelleft=False, top=True, right=True)
ax6.tick_params(labelleft=False, top=True, right=True)
ax9.tick_params(labelleft=False, top=True, right=True)
ax1.set_ylabel(r"Detector-$y$ Bin$\times$2 (Pixel)")
ax4.set_ylabel(r"Detector-$y$ Bin$\times$2 (Pixel)")
ax7.set_ylabel(r"Detector-$y$ (Pixel)")
ax7.set_xlabel(r"Raster Step")
ax8.set_xlabel(r"Raster Step")
ax9.set_xlabel(r"Detector-$\lambda$ (Pixel)")
for clb_ax_ in (clb_ax1, clb_ax2, clb_ax3, clb_ax4, clb_ax5, clb_ax6,
clb_ax7, clb_ax8, clb_ax9):
clb_ax_.tick_params(top=True, labeltop=True, labelbottom=False, bottom=False, which="both")
for ax_ in (ax1,ax2,ax4,ax5,ax7,ax8):
ax_.axvline(example_x_index, color="w", lw=2, ls="--", alpha=0.7)
for ax_, text_ in zip((ax1,ax2,ax4,ax5,ax7,ax8),
("a", "b", "d", "e", "g", "h")):
ax_.text(0.03,0.03,r"\textbf{{{:})}}".format(text_), transform=ax_.transAxes,
va="bottom", ha="left", color="w", fontsize=10,
path_effects=[path_effects.Stroke(linewidth=1.5, foreground='black'),path_effects.Normal()])
for ax_, text_ in zip((ax3,ax6,ax9),("c","f","i")):
ax_.text(0.06,0.03,r"\textbf{{{:})}}".format(text_), transform=ax_.transAxes,
va="bottom", ha="left", color="w", fontsize=10,
path_effects=[path_effects.Stroke(linewidth=1.5, foreground='black'),path_effects.Normal()])
ax1.text(0.05,0.95,r"\textbf{Raw}", transform=ax1.transAxes,
va="top", ha="left", color="w", fontsize=12,
path_effects=[path_effects.Stroke(linewidth=1.5, foreground='black'),path_effects.Normal()])
ax4.text(0.05,0.95,r"\textbf{PSF Deconv.}", transform=ax4.transAxes,
va="top", ha="left", color="w", fontsize=12,
path_effects=[path_effects.Stroke(linewidth=1.5, foreground='black'),path_effects.Normal()])
ax7.text(0.05,0.95,r"\textbf{Skew Corr.}", transform=ax7.transAxes,
va="top", ha="left", color="w", fontsize=12,
path_effects=[path_effects.Stroke(linewidth=1.5, foreground='black'),path_effects.Normal()])
target_r1_param = {'xy':(50, 100), 'width': 12, 'height': 12/plot_aspect, 'angle':0,
'ec':"#1B813E", 'fill': False, 'lw':2}
target_r2_param = {'xy':(54, 37), 'width': 18, 'height': 5/plot_aspect, 'angle':60,
'ec':"#F75C2F", 'fill': False, 'lw':2}
target_r1_param_nobin = {'xy':(50, 100*2), 'width': 12, 'height': 12/plot_aspect*2, 'angle':0,
'ec':"#1B813E", 'fill': False, 'lw':2}
target_r2_param_nobin = {'xy':(54, 37*2), 'width': 18, 'height': 5*2/plot_aspect, 'angle':60,
'ec':"#F75C2F", 'fill': False, 'lw':2}
# the patches.Ellipse cannot handle the rotation
# in axes with different scales so we need to
# plot it explicitly
r2_t_ = np.linspace(0, 2*np.pi, 360)
r2_x_ = target_r2_param['xy'][0] + \
target_r2_param['width']*np.cos(r2_t_)*np.cos(np.deg2rad(target_r2_param['angle'])) - \
target_r2_param['height']*plot_aspect*np.sin(r2_t_)*np.sin(np.deg2rad(target_r2_param['angle']))
r2_y_ = target_r2_param['xy'][1] + \
target_r2_param['width']*np.cos(r2_t_)*np.sin(np.deg2rad(target_r2_param['angle'])) + \
target_r2_param['height']*plot_aspect*np.sin(r2_t_)*np.cos(np.deg2rad(target_r2_param['angle']))
target_r3_param = {'xy':(49, 5), 'width': 12, 'height': 12/plot_aspect, 'angle':0,
'ec':"#E03C8A", 'fill': False, 'lw':2}
target_r3_param_nobin = {'xy':(49, 5*2), 'width': 12, 'height': 12*2/plot_aspect, 'angle':0,
'ec':"#E03C8A", 'fill': False, 'lw':2}
for ax_ in (ax1,ax2,ax4,ax5):
ax_lim_ = ax_.axis()
circle_r1_ = Ellipse(**target_r1_param)
circle_r3_ = Ellipse(**target_r3_param)
ax_.add_patch(circle_r1_)
ax_.add_patch(circle_r3_)
ax_.plot(r2_x_, r2_y_, lw=target_r2_param['lw'], color=target_r2_param['ec'])
ax_.axis(ax_lim_)
for ax_ in (ax7,ax8):
ax_lim_ = ax_.axis()
circle_r1_ = Ellipse(**target_r1_param_nobin)
circle_r3_ = Ellipse(**target_r3_param_nobin)
ax_.add_patch(circle_r1_)
ax_.add_patch(circle_r3_)
ax_.plot(r2_x_, r2_y_*2, lw=target_r2_param['lw'], color=target_r2_param['ec'])
ax_.axis(ax_lim_)
fig.savefig("../figs/spice_psf_deconv.pdf", dpi=600, format="pdf", bbox_inches="tight")
plt.show()
