DJA Imaging Data Products

imaging demo jwst

(This page is auto-generated from the Jupyter notebook image-data-products.ipynb.)

Here we summarize the files available for imaging datasets, for example in the v7 data release.

%matplotlib inline
import os
import warnings
warnings.filterwarnings('ignore')

import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage as nd

import astropy.io.fits as pyfits
import astropy.units as u

import sep

import grizli
from grizli import utils
print(f'grizli version: {grizli.__version__}')

BASE_URL = 'https://s3.amazonaws.com/grizli-v2/JwstMosaics/v7/'
grizli version: 1.10.dev3+g341a999

File extensions

Generally, for a given root and filter combination, the following files are available:

  • {root}-{filter}_drc_sci.fits.gz: Science image
  • {root}-{filter}_drc_wht.fits.gz: Inverse variance weight image (sky + readnoise)
  • {root}-{filter}_drc_exp.fits.gz: Exposure-time map
  • {root}-{filter}_wcs.csv: Table summarizing individual exposures that contribute to the mosaic

Notes

  1. All mosaics are created with the legacy drizzlepac.adrizzle.do_driz drizzle implementation that works interchangeably with JWST and HST.
    • grizli generates WCS for each exposure that follow the SIP-WCS convention and that match the newer gwcs JWST wcs at the level of 1e-4 pixels or better
  2. All NIRCam LW, NIRISS (and HST) mosaics are created with 40 mas pixels
  3. Most fields have 20 mas pixels for the NIRCam SW images that exacly subsample the LW grid 2x2.
    • The very large primer-cosmos and primer-uds SW mosaics have 40 mas pixels
  4. All sci mosaics have intensity units of 10 nJy / pix, corresponding to an AB magnitude zeropoint 28.9. This has the slightly desirable property that the image pixel values are not too different from unity.
    • These are not the same as the surface brightness units of the JWST pipeline!
  5. The exp exposure time images have units of seconds rounded to the nearest integer
    • Subsampled to 4x4 of the parent mosaic to keep the file sizes small
    • The exp images are created directly from the footprints of the constituent exposures and don’t account for masked pixels within an exposure.
# Example
root = 'smacs0723-grizli-v7.0'

filter = 'f444w-clear'
# Open the files directly from the web

img = {}

print('# File shape')

for ext in ['sci','wht','exp']:
    _file = f'{root}-{filter}_drc_{ext}.fits.gz'
    img[ext] = pyfits.open(os.path.join(BASE_URL, _file))
    print(f'{_file} : {img[ext][0].data.shape}')
# File shape
smacs0723-grizli-v7.0-f444w-clear_drc_sci.fits.gz : (12000, 12000)
smacs0723-grizli-v7.0-f444w-clear_drc_wht.fits.gz : (12000, 12000)
smacs0723-grizli-v7.0-f444w-clear_drc_exp.fits.gz : (3000, 3000)
# Files have a single PrimaryHDU
img['sci'].info()
Filename: /Users/gbrammer/.astropy/cache/download/url/b0380671ce11dec1c5653485f66f705c/contents
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      93   (12000, 12000)   float32   

Primary sci header

img['sci'][0].header
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -32 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                12000                                                  
NAXIS2  =                12000                                                  
WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =               4591.5 / Pixel coordinate of reference point            
CRPIX2  =               6515.5 / Pixel coordinate of reference point            
CD1_1   = -1.1111111111111E-05 / Coordinate transformation matrix element       
CD2_2   =  1.1111111111111E-05 / Coordinate transformation matrix element       
CDELT1  =                  1.0 / [deg] Coordinate increment at reference point  
CDELT2  =                  1.0 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'RA---TAN'           / Right ascension, gnomonic projection           
CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection               
CRVAL1  =            110.83403 / [deg] Coordinate value at reference point      
CRVAL2  =            -73.45429 / [deg] Coordinate value at reference point      
LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
LATPOLE =            -73.45429 / [deg] Native latitude of celestial pole        
MJDREF  =                  0.0 / [d] MJD of fiducial time                       
DATE-OBS= '2022-06-07'         / ISO-8601 time of observation                   
MJD-OBS =              59737.0 / [d] MJD of observation                         
RADESYS = 'ICRS'               / Equatorial coordinate system                   
CD1_2   =                  0.0                                                  
CD2_1   =                  0.0                                                  
DRIZKERN= 'square  '           / Drizzle kernel                                 
DRIZPIXF=                 0.75 / Drizzle pixfrac                                
EXPTIME =    15074.42400000001                                                  
NDRIZIM =                   18                                                  
PIXFRAC =                 0.75                                                  
KERNEL  = 'square  '                                                            
OKBITS  =                    4 / FLT bits treated as valid                      
PHOTSCAL=    1.001401962747847 / Scale factor applied                           
GRIZLIV = '1.8.16.dev12+g86ad0c1' / Grizli code version                         
WHTTYPE = 'jwst    '           / Exposure weighting strategy                    
RNPERC  =                   99 / VAR_RNOISE clip percentile for JWST            
TELESCOP= 'JWST    '                                                            
FILTER  = 'F444W   '                                                            
PUPIL   = 'CLEAR   '                                                            
DETECTOR= 'NRCALONG'                                                            
INSTRUME= 'NIRCAM  '                                                            
PHOTFLAM= 1.54184756289340E-22                                                  
PHOTPLAM=    44036.71097714713                                                  
PHOTFNU =                1E-08                                                  
EXPSTART=    59737.22120032604                                                  
EXPEND  =    59737.23101751157                                                  
TIME-OBS= '05:18:31.708'                                                        
UPDA_CTX= 'jwst_0995.pmap'                                                      
CRDS_CTX= 'jwst_1041.pmap'                                                      
R_DISTOR= 'jwst_nircam_distortion_0141.asdf'                                    
R_PHOTOM= 'jwst_nircam_photom_0111.fits'                                        
R_FLAT  = 'jwst_nircam_flat_0574.fits'                                          
PHOTMJSR=   0.3925000131130219                                                  
PIXAR_SR=             9.31E-14                                                  
FLT00001= 'jw02736001001_02105_00001_nrcalong_rate.fits'                        
WHT00001=      14335.236328125 / Median weight of exposure 1                    
FLT00002= 'jw02736001001_02105_00001_nrcblong_rate.fits'                        
WHT00002=     14818.2607421875 / Median weight of exposure 2                    
FLT00003= 'jw02736001001_02105_00002_nrcalong_rate.fits'                        
WHT00003=      14245.970703125 / Median weight of exposure 3                    
FLT00004= 'jw02736001001_02105_00002_nrcblong_rate.fits'                        
WHT00004=     14709.4404296875 / Median weight of exposure 4                    
FLT00005= 'jw02736001001_02105_00003_nrcalong_rate.fits'                        
WHT00005=     14335.7177734375 / Median weight of exposure 5                    
FLT00006= 'jw02736001001_02105_00003_nrcblong_rate.fits'                        
WHT00006=     14799.9052734375 / Median weight of exposure 6                    
FLT00007= 'jw02736001001_02105_00004_nrcalong_rate.fits'                        
WHT00007=       14333.94921875 / Median weight of exposure 7                    
FLT00008= 'jw02736001001_02105_00004_nrcblong_rate.fits'                        
WHT00008=       14851.24609375 / Median weight of exposure 8                    
FLT00009= 'jw02736001001_02105_00005_nrcalong_rate.fits'                        
WHT00009=     14321.4130859375 / Median weight of exposure 9                    
FLT00010= 'jw02736001001_02105_00005_nrcblong_rate.fits'                        
WHT00010=      14816.029296875 / Median weight of exposure 10                   
FLT00011= 'jw02736001001_02105_00006_nrcalong_rate.fits'                        
WHT00011=     14384.9931640625 / Median weight of exposure 11                   
FLT00012= 'jw02736001001_02105_00006_nrcblong_rate.fits'                        
WHT00012=      14822.169921875 / Median weight of exposure 12                   
FLT00013= 'jw02736001001_02105_00007_nrcalong_rate.fits'                        
WHT00013=       14322.65234375 / Median weight of exposure 13                   
FLT00014= 'jw02736001001_02105_00007_nrcblong_rate.fits'                        
WHT00014=     14808.3759765625 / Median weight of exposure 14                   
FLT00015= 'jw02736001001_02105_00008_nrcalong_rate.fits'                        
WHT00015=     14282.8134765625 / Median weight of exposure 15                   
FLT00016= 'jw02736001001_02105_00008_nrcblong_rate.fits'                        
WHT00016=     14759.8623046875 / Median weight of exposure 16                   
FLT00017= 'jw02736001001_02105_00009_nrcalong_rate.fits'                        
WHT00017=     14353.7353515625 / Median weight of exposure 17                   
FLT00018= 'jw02736001001_02105_00009_nrcblong_rate.fits'                        
WHT00018=     14783.2705078125 / Median weight of exposure 18                   
OPHOTFNU= 9.30775449348276E-08 / Original PHOTFNU before scaling                
BUNIT   = '10.0*nanoJansky'                                                     
# Images have units of 10 nJy / pix
for k in ('FILTER','PHOTFNU','PHOTPLAM','BUNIT'):
    print(f"{k:>8}: {img['sci'][0].header[k]}")
  FILTER: F444W
 PHOTFNU: 1e-08
PHOTPLAM: 44036.71097714713
   BUNIT: 10.0*nanoJansky

Primary exp header

img['exp'][0].header
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                   32 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                 3000                                                  
NAXIS2  =                 3000                                                  
WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =             1147.875 / Pixel coordinate of reference point            
CRPIX2  =             1628.875 / Pixel coordinate of reference point            
CD1_1   = -4.4444444444444E-05 / Coordinate transformation matrix element       
CD2_2   =  4.4444444444444E-05 / Coordinate transformation matrix element       
CDELT1  =                  1.0 / [deg] Coordinate increment at reference point  
CDELT2  =                  1.0 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'RA---TAN'           / Right ascension, gnomonic projection           
CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection               
CRVAL1  =            110.83403 / [deg] Coordinate value at reference point      
CRVAL2  =            -73.45429 / [deg] Coordinate value at reference point      
LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
LATPOLE =            -73.45429 / [deg] Native latitude of celestial pole        
MJDREF  =                  0.0 / [d] MJD of fiducial time                       
DATE-OBS= '2022-06-07'         / ISO-8601 time of observation                   
MJD-OBS =              59737.0 / [d] MJD of observation                         
RADESYS = 'ICRS'               / Equatorial coordinate system                   
CD1_2   =                  0.0                                                  
CD2_1   =                  0.0                                                  
DRIZKERN= 'square  '           / Drizzle kernel                                 
DRIZPIXF=                 0.75 / Drizzle pixfrac                                
EXPTIME =    15074.42400000001                                                  
NDRIZIM =                   18                                                  
PIXFRAC =                 0.75                                                  
KERNEL  = 'square  '                                                            
OKBITS  =                    4 / FLT bits treated as valid                      
PHOTSCAL=    1.001401962747847 / Scale factor applied                           
GRIZLIV = '1.8.16.dev12+g86ad0c1' / Grizli code version                         
WHTTYPE = 'jwst    '           / Exposure weighting strategy                    
RNPERC  =                   99 / VAR_RNOISE clip percentile for JWST            
TELESCOP= 'JWST    '                                                            
FILTER  = 'F444W   '                                                            
PUPIL   = 'CLEAR   '                                                            
DETECTOR= 'NRCALONG'                                                            
INSTRUME= 'NIRCAM  '                                                            
PHOTFLAM= 1.54184756289340E-22                                                  
PHOTPLAM=    44036.71097714713                                                  
PHOTFNU =                1E-08                                                  
EXPSTART=    59737.22120032604                                                  
EXPEND  =    59737.23101751157                                                  
TIME-OBS= '05:18:31.708'                                                        
UPDA_CTX= 'jwst_0995.pmap'                                                      
CRDS_CTX= 'jwst_1041.pmap'                                                      
R_DISTOR= 'jwst_nircam_distortion_0141.asdf'                                    
R_PHOTOM= 'jwst_nircam_photom_0111.fits'                                        
R_FLAT  = 'jwst_nircam_flat_0574.fits'                                          
PHOTMJSR=   0.3925000131130219                                                  
PIXAR_SR=             9.31E-14                                                  
FLT00001= 'jw02736001001_02105_00001_nrcalong_rate.fits'                        
WHT00001=      14335.236328125 / Median weight of exposure 1                    
FLT00002= 'jw02736001001_02105_00001_nrcblong_rate.fits'                        
WHT00002=     14818.2607421875 / Median weight of exposure 2                    
FLT00003= 'jw02736001001_02105_00002_nrcalong_rate.fits'                        
WHT00003=      14245.970703125 / Median weight of exposure 3                    
FLT00004= 'jw02736001001_02105_00002_nrcblong_rate.fits'                        
WHT00004=     14709.4404296875 / Median weight of exposure 4                    
FLT00005= 'jw02736001001_02105_00003_nrcalong_rate.fits'                        
WHT00005=     14335.7177734375 / Median weight of exposure 5                    
FLT00006= 'jw02736001001_02105_00003_nrcblong_rate.fits'                        
WHT00006=     14799.9052734375 / Median weight of exposure 6                    
FLT00007= 'jw02736001001_02105_00004_nrcalong_rate.fits'                        
WHT00007=       14333.94921875 / Median weight of exposure 7                    
FLT00008= 'jw02736001001_02105_00004_nrcblong_rate.fits'                        
WHT00008=       14851.24609375 / Median weight of exposure 8                    
FLT00009= 'jw02736001001_02105_00005_nrcalong_rate.fits'                        
WHT00009=     14321.4130859375 / Median weight of exposure 9                    
FLT00010= 'jw02736001001_02105_00005_nrcblong_rate.fits'                        
WHT00010=      14816.029296875 / Median weight of exposure 10                   
FLT00011= 'jw02736001001_02105_00006_nrcalong_rate.fits'                        
WHT00011=     14384.9931640625 / Median weight of exposure 11                   
FLT00012= 'jw02736001001_02105_00006_nrcblong_rate.fits'                        
WHT00012=      14822.169921875 / Median weight of exposure 12                   
FLT00013= 'jw02736001001_02105_00007_nrcalong_rate.fits'                        
WHT00013=       14322.65234375 / Median weight of exposure 13                   
FLT00014= 'jw02736001001_02105_00007_nrcblong_rate.fits'                        
WHT00014=     14808.3759765625 / Median weight of exposure 14                   
FLT00015= 'jw02736001001_02105_00008_nrcalong_rate.fits'                        
WHT00015=     14282.8134765625 / Median weight of exposure 15                   
FLT00016= 'jw02736001001_02105_00008_nrcblong_rate.fits'                        
WHT00016=     14759.8623046875 / Median weight of exposure 16                   
FLT00017= 'jw02736001001_02105_00009_nrcalong_rate.fits'                        
WHT00017=     14353.7353515625 / Median weight of exposure 17                   
FLT00018= 'jw02736001001_02105_00009_nrcblong_rate.fits'                        
WHT00018=     14783.2705078125 / Median weight of exposure 18                   
OPHOTFNU= 9.30775449348276E-08 / Original PHOTFNU before scaling                
BUNIT   = 'second  '                                                            
SAMPLE  =                    4 / Sampling factor                                
NXORIG  =                12000                                                  
NYORIG  =                12000                                                  
MOSPSCL =  0.03999999999999958 / Mosaic pixel scale arcsec                      
ORIGPSCL=   0.0629361212228063 / Original detector pixel scale arcsec           
DNTOEPS =    58.62417855098175 / Inverse flux conversion back to e per second   
BSCALE  =                    1                                                  
BZERO   =           2147483648                                                  

WCS log

The wcs.csv files contain the full SIP header of each exposure that contributes to the mosaic, along with some epoch information.

_file = f'{root}-{filter}_wcs.csv'
wcs = utils.read_catalog(os.path.join(BASE_URL, _file))
print(wcs.colnames)
['file', 'ext', 'exptime', 'wcsaxes', 'crpix1', 'crpix2', 'cd1_1', 'cd1_2', 'cd2_1', 'cd2_2', 'cdelt1', 'cdelt2', 'cunit1', 'cunit2', 'ctype1', 'ctype2', 'crval1', 'crval2', 'lonpole', 'latpole', 'wcsname', 'mjdref', 'date-beg', 'mjd-beg', 'date-avg', 'mjd-avg', 'date-end', 'mjd-end', 'xposure', 'telapse', 'obsgeo-x', 'obsgeo-y', 'obsgeo-z', 'radesys', 'velosys', 'a_order', 'a_0_2', 'a_0_3', 'a_0_4', 'a_0_5', 'a_1_1', 'a_1_2', 'a_1_3', 'a_1_4', 'a_2_0', 'a_2_1', 'a_2_2', 'a_2_3', 'a_3_0', 'a_3_1', 'a_3_2', 'a_4_0', 'a_4_1', 'a_5_0', 'b_order', 'b_0_2', 'b_0_3', 'b_0_4', 'b_0_5', 'b_1_1', 'b_1_2', 'b_1_3', 'b_1_4', 'b_2_0', 'b_2_1', 'b_2_2', 'b_2_3', 'b_3_0', 'b_3_1', 'b_3_2', 'b_4_0', 'b_4_1', 'b_5_0', 'naxis', 'naxis1', 'naxis2', 'sipcrpx1', 'sipcrpx2']
# First few lines
wcs['file','ext','exptime','mjd-avg','date-avg','crpix1','crpix2','crval1','crval2'][:4]
GTable length=4
fileextexptimemjd-avgdate-avgcrpix1crpix2crval1crval2
str44int64float64float64str23float64float64float64float64
jw02736001001_02105_00001_nrcalong_rate.fits1837.46859737.2261089192022-06-07T05:25:35.8111024.6471024.66110.68612332528-73.481395298838
jw02736001001_02105_00001_nrcblong_rate.fits1837.46859737.2261052152022-06-07T05:25:35.4911024.4891024.662110.8276309259-73.453986741388
jw02736001001_02105_00002_nrcalong_rate.fits1837.46859737.2367955742022-06-07T05:40:59.1381024.6471024.66110.68291892635-73.47999884627
jw02736001001_02105_00002_nrcblong_rate.fits1837.46859737.236791872022-06-07T05:40:58.8181024.4891024.662110.82443032455-73.452590822835

Compare the WHT and EXP images


exts = ['sci','wht','exp']

fig, axes = plt.subplots(2,len(exts),figsize=(3*len(exts),6))

for j, ext in enumerate(exts):
    msk = img[ext][0].data != 0
    wmax = np.nanpercentile(img[ext][0].data[msk], 95)
    for i in [0,1]:
        axes[i][j].imshow(img[ext][0].data, vmin=0, vmax=wmax, origin='lower', cmap='magma')
        axes[i][j].grid()
        if i == 0:
            axes[i][j].set_title(ext)
        
xy = 2600, 6100, 256

for j, p in enumerate([0,0,1]):
    axes[0][j].set_xlim(*(xy[0] + np.array([-1,1])*xy[2])/4**p)
    axes[0][j].set_ylim(*(xy[1] + np.array([-1,1])*xy[2])/4**p)
    
fig.tight_layout(pad=1)

png

Make a full variance image including the Poisson component from the sources themselves

The mosaics are created by weighting each input exposure by a factor like 1/wht = VAR_RNOISE + median(VAR_POISSON) from the JWST exposure files. The first term incorporates pixel-to-pixel variations resulting from pixels where some fraction of the reads may have been masked as saturated or affected by cosmic rays. The second term effectively provides the noise from the sky background, but without including the Poisson term for individual sources.

Certain applications like photometry or morphology fitting with galfit may require a variance / sigma image that includes the poisson term from the individual sources. This can be generated from the exp maps as shown below.

# Grow the exposure map to the original frame
full_exp = np.zeros(img['sci'][0].data.shape, dtype=int)
full_exp[2::4,2::4] += img['exp'][0].data*1
full_exp = nd.maximum_filter(full_exp, 4)

img['Full exp'] = pyfits.HDUList([pyfits.PrimaryHDU(data=full_exp)])
# Show the full exposure map
exts = ['wht','Full exp','exp']

fig, axes = plt.subplots(2,len(exts),figsize=(3*len(exts),6))

for j, ext in enumerate(exts):
    msk = img[ext][0].data != 0
    wmax = np.nanpercentile(img[ext][0].data[msk], 95)
    for i in [0,1]:
        axes[i][j].imshow(img[ext][0].data, vmin=0, vmax=wmax, origin='lower', cmap='magma')
        axes[i][j].grid()
        if i == 0:
            axes[i][j].set_title(ext)
        
xy = 2600, 6100, 256

for j, p in enumerate([0,0,1]):
    axes[0][j].set_xlim(*(xy[0] + np.array([-1,1])*xy[2])/4**p)
    axes[0][j].set_ylim(*(xy[1] + np.array([-1,1])*xy[2])/4**p)
    
fig.tight_layout(pad=1)

png

Effective “gain”

To convert the Poisson variance associated with the sci image, write down the multiplicative factors that had been applied to the original count-rate data in the pipeline rate files. PHOTMJSR is the original flux calibration to units of “MJy/sr” and PHOTSCL is any (small) additional photometric term that was included by grizli. The PHOTFNU / OPHOTFNU term accounts for any final scaling of the output mosaic and the ratio of the original and mosaic pixel areas.

# Scale factors
phot_scale = 1 / (PHOTMJSR * PHOTSCAL) * PHOTFNU / OPHOTFNU

# Effective gain e-/DN, including exposure time
effective_gain = phot_scale * exposure_time_map

# Variance in electrons = counts in electrons
var_poisson_elec = sci * effective_gain

# Variance in mosaic DN
var_poisson_dn = var_poisson_elec / effective_gain**2 = sci / effective_gain
header = img['exp'][0].header

# Multiplicative factors that have been applied since the original count-rate images
phot_scale = 1.

for k in ['PHOTMJSR','PHOTSCAL']:
    print(f'{k} {header[k]:.3f}')
    phot_scale /= header[k]

# Unit and pixel area scale factors
if 'OPHOTFNU' in header:
    phot_scale *= header['PHOTFNU'] / header['OPHOTFNU']

# "effective_gain" = electrons per DN of the mosaic
effective_gain = (phot_scale * full_exp)

# Poisson variance in mosaic DN
var_poisson_dn = np.maximum(img['sci'][0].data, 0) / effective_gain

# Original variance from the `wht` image = RNOISE + BACKGROUND
var_wht = 1/img['wht'][0].data

# New total variance
var_total = var_wht + var_poisson_dn
full_wht = 1 / var_total

# Null weights
full_wht[var_total <= 0] = 0

img['Full wht'] = pyfits.HDUList([pyfits.PrimaryHDU(data=full_wht, header=img['wht'][0].header)])
PHOTMJSR 0.393
PHOTSCAL 1.001
# Compare science and weight arrays, where you can now 
# "see" the sources in the Full weight array

exts = ['sci','Full wht','wht']

fig, axes = plt.subplots(2,len(exts),figsize=(3*len(exts),6))

for j, ext in enumerate(exts):
    msk = img[ext][0].data != 0
    wmax = np.nanpercentile(img[ext][0].data[msk], 95)
    for i in [0,1]:
        axes[i][j].imshow(img[ext][0].data, vmin=0, vmax=wmax, origin='lower', cmap='magma')
        axes[i][j].grid()
        if i == 0:
            axes[i][j].set_title(ext)
        
xy = 2600, 6100, 256

for j, p in enumerate([0,0,0]):
    axes[0][j].set_xlim(*(xy[0] + np.array([-1,1])*xy[2])/4**p)
    axes[0][j].set_ylim(*(xy[1] + np.array([-1,1])*xy[2])/4**p)
    
fig.tight_layout(pad=1)

png