UniversalETC


Introduction

Sensitivity related calculations including sensitivity limits and data simulation. Developed for the Australian Space Eye but much of this should have general applicability.

Code preamble


In [1]:
from __future__ import division,unicode_literals # Half arsed attempt at Python 2 compatibility

In [2]:
%matplotlib inline

In [3]:
import numpy as np
from scipy.interpolate import RectSphereBivariateSpline, SmoothBivariateSpline, InterpolatedUnivariateSpline
from matplotlib import pyplot as plt
from matplotlib import rcParams
from matplotlib import colors
from astropy.io import fits
from astropy.table import Table
from astropy import units as u
from astropy.constants import h,c
import healpy as hp

In [4]:
rcParams['figure.figsize'] = 12, 9

Instrument efficiency & light gathering power

End-to-end throughput of the instrument from the entrance aperture through to photo-electrons in the image sensor

Entrance pupil area

  • Cubesat format effectively limits us to a maximum aperture diameter of 90 mm
  • To avoid diffraction spike we want to a circular pupil
  • Pupil will be unobscured

In [5]:
pupil_area = (np.pi * (90 * u.mm)**2 / 4).si
pupil_area


Out[5]:
$0.0063617251 \; \mathrm{m^{2}}$

Throughput of optics

  • With appropriate glasses there will be negligible bulk absorption in the visible wavelength range
  • With anti-reflection coatings optimised for a specific band can expect at least 0.995 transmission at each surface throughout each band
  • Baseline design is 5 elements in 5 groups, i.e. 10 vacuum-glass interfaces
  • Last surface accounted for in next section

In [6]:
n_surfaces = 9 # Doesn't include last surface
throughput = 0.995**n_surfaces * u.dimensionless_unscaled
throughput


Out[6]:
$0.95588958 \; \mathrm{}$

Filters


In [7]:
# Load DECam filter data
filters = [Table.read('http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.u'), \
           Table.read('http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.g'), \
           Table.read('http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.r'), \
           Table.read('http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.i'), \
           Table.read('http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.z')]

for f in filters:
    # Set units correctly
    f['Wavelength'] = f['Wavelength'].to(u.um)
    f['Wavelength'].unit = u.um
    f['Transmission'].unit = u.dimensionless_unscaled


Downloading http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.u [Done]
WARNING: W42: None:2:0: W42: No XML namespace specified [astropy.io.votable.tree]
WARNING:astropy:W42: None:2:0: W42: No XML namespace specified
Downloading http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.g [Done]
Downloading http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.r [Done]
Downloading http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.i [Done]
Downloading http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.z [Done]

In [8]:
for f in filters:
    plt.plot(f['Wavelength'], f['Transmission'])
    plt.xlim(0.25,1.05)
    plt.xlabel('Wavelength/$\mu$m')
    plt.ylabel('Transmission')
    plt.title('DECam filter transmission')


Image sensor

Baseline concept uses an e2V CIS115, a 2000 x 1504 pixel back-illuminated monolithic CMOS image sensor with 7 micron pixel pitch. Numerical QE data is hard to come by but there are some plots available from which we can extract approximate date with http://arohatgi.info/WebPlotDigitizer/, the most useful of which appears to be in Soman et al (2014), http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=1891353, which has full wavelength range predicted QE at both +20$^\circ$C and -40$^\circ$C. These data seem quite conservative as Wang et al (2014), http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=1891414, get higher peak QE from both e2v's and their own measurements on a prototype device.


In [11]:
qe_p20C = Table.read('resources/cis115_QE_20C.csv')
qe_m40C = Table.read('resources/cis115_QE_-40C.csv')

qe_p20C['Wavelength'].unit = u.um
qe_m40C['Wavelength'].unit = u.um

qe_p20C['QE'].unit = u.electron / u.photon
qe_m40C['QE'].unit = u.electron / u.photon

In [12]:
plt.plot(qe_p20C['Wavelength'], qe_p20C['QE'], 'r-', label='$+20^\circ$C')
plt.plot(qe_m40C['Wavelength'], qe_m40C['QE'], 'b-', label='$-40^\circ$C')
plt.xlim(0.25,1.05)
plt.ylim(0,1)
plt.xlabel('Wavelength/$\mu$m')
plt.ylabel('QE / e photon$^{-1}$')
plt.legend(loc='best')
plt.title('e2v CIS115 predicted QE (Soman et al)')
plt.savefig('QE.png')


End-to-end efficiency

Combine throughput of optics, filters and image sensor QE (at nominal operating temperature of -40$^\circ$C).


In [13]:
efficiencies = []

for f in filters:
    # Interpolate QE data at same wavelengths as filter transmission data
    i = np.interp(f['Wavelength'], qe_m40C['Wavelength'], qe_m40C['QE']) *  qe_m40C['QE'].unit
    efficiencies.append(throughput * f['Transmission'].quantity * i)

In [14]:
for i, f in enumerate(filters):
    plt.plot(f['Wavelength'], efficiencies[i])
    plt.xlim(0.25,1.05)
    plt.ylim(0,1)
    plt.xlabel('Wavelength/$\mu$m')
    plt.ylabel('Efficiency / e photon$^{-1}$')
    plt.title('End-to-end efficiency')
    plt.savefig('eff.png')


Effective aperture area

Can calculate an effective aperture area as the aperture area multiplied by the end-to-end efficiency


In [15]:
eff_areas = [pupil_area * efficiency for efficiency in efficiencies]

Effective étendue

Similarly can calculate an effective étendue ($A\Omega$) per pixel by multiplying the effective aperture area by the solid angle subtended by each pixel. In the baseline design each pixel is 3 arcseconds across.


In [16]:
eff_etendues = [eff_area * (3 * u.arcsecond)**2 / u.pixel for eff_area in eff_areas]

Sky background

Zodical light

In space the visible wavelength sky background is dominated by zodiacal light. To model the zodiacal light we follow the same prescription as the HST ETC, as described in Giavilisco, Saul & Bohlin (2002), http://www.stsci.edu/hst/wfc3/documents/ISRs/2002/WFC3-2002-12.pdf.

Solar spectrum

In the visible wavelength range the zodiacal light is the result of scattered sunlight. Consequently the starting point is a solar spectrum from Colina, Bohlin & Castelli (1996), http://adslabs.org/adsabs/abs/1996AJ....112..307C/


In [17]:
sun = fits.open('resources/sun_castelli.fits')

In [18]:
sun.info()


Filename: resources/sun_castelli.fits
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU      14   ()              
1                BinTableHDU     16   22999R x 2C   [1E, 1E]   

In [19]:
sun[0].header


Out[19]:
SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                   16 / number of bits per data pixel                  
NAXIS   =                    0 / number of data axes                            
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format defined in Astronomy and
COMMENT   Astrophysics Supplement Series v44/p363, v44/p371, v73/p359, v73/p365.
COMMENT   Contact the NASA Science Office of Standards and Technology for the   
COMMENT   FITS Definition document #100 and other FITS information.             
ORIGIN  = 'STScI-STSDAS/TABLES' / Tables version 1999-03-22                     
FILENAME= 'sun_castelli.fits'  / name of file                                   
HISTORY   Created Thu 15:49:31 16-Nov-95                                        
COMMENT   solar model spectrum calculated by F. Castelli.                       
COMMENT   Absolute flux normalized to a V flux of 184.2 ergs/s/cm^2/A           
COMMENT   For more details see Colina, Bohlin & Castelli 1996 CAL/SCS-008       

In [20]:
sun[1].header


Out[20]:
XTENSION= 'BINTABLE'           / binary table extension                         
BITPIX  =                    8 / 8-bit bytes                                    
NAXIS   =                    2 / 2-dimensional binary table                     
NAXIS1  =                    8 / width of table in bytes                        
NAXIS2  =                22999                                                  
PCOUNT  =                    0 / size of special data area                      
GCOUNT  =                    1 / one data group (required keyword)              
TFIELDS =                    2                                                  
TTYPE1  = 'WAVELENGTH'         / label for field   1                            
TFORM1  = '1E      '           / data format of field: 4-byte REAL              
TUNIT1  = 'ANGSTROMS'          / physical unit of field                         
TTYPE2  = 'FLUX    '           / label for field   2                            
TFORM2  = '1E      '           / data format of field: 4-byte REAL              
TUNIT2  = 'FLAM    '           / physical unit of field                         
TDISP1  = 'G15.7   '           / display format                                 
TDISP2  = 'G15.7   '           / display format                                 

In [21]:
sun_waves = sun[1].data['WAVELENGTH'] * u.Angstrom

In [22]:
sun_flux = sun[1].data['FLUX'] * u.erg * u.second**-1 * u.cm**-2 * u.Angstrom**-1

Normalisation & reddening

To convert the solar spectrum to a zodical light spectrum it is normalised and reddened following the prescription of Leinert et al (1997), http://www.edpsciences.org/10.1051/aas:1998105, with the revised parameters from Aldering (2001), http://www-supernova.lbl.gov/~aldering/zodi.pdf. The result is normalised to the correct surface brightness for the north ecliptic pole (NEP).


In [23]:
# Colina, Bohlin & Castelli solar spectrum normalised to V band flux 
# of 184.2 ergs/s/cm^2/A,
solar_normalisation = 184.2 * u.erg * u.second**-1 * u.cm**-2 * u.Angstrom**-1
# Leinert at al NEP zodical light 1.81e-18 erg/s/cm^2/A/arcsec^2 at 0.5 um, 
# Aldering suggests using 0.01 dex lower.
zl_nep = 1.81e-18 * u.erg * u.second**-1 * u.cm**-2 * u.Angstrom**-1 * \
    u.arcsecond**-2 * 10**(-0.01)

zl_normalisation = zl_nep / solar_normalisation

lambda_c = 0.5 * u.micron

# Aldering reddening parameters
f_blue = 0.9
f_red = 0.48

def solar_to_ZL(waves):
    """
    Calculate factors to convert a Colina et alsolar spectrum to spectrum
    for the zodical light normalised for NEP.  Follows the prescription of 
    Leinert et al (1997) with the revised parameters from Aldering (2001) 
    as used in the HST ETC (Giavalsico, Sahi & Bohlin (2002)).
    """
    rfactor = np.where(waves < lambda_c, \
                       1.0 + f_blue * np.log(waves/lambda_c), \
                       1.0 + f_red * np.log(waves/lambda_c))
    return zl_normalisation * rfactor

In [24]:
zl_waves = sun_waves.to(u.micron)
zl_flux = (sun_flux * solar_to_ZL(sun_waves))
zl_flux = zl_flux.to(u.Watt * u.m**-2 * u.arcsecond**-2 * u.micron**-1)

In [25]:
plt.semilogy(zl_waves, zl_flux)
plt.xlim(0.2,2.5)
plt.ylim(1e-18,3e-17)
plt.xlabel('Wavelength/$\mu$m')
plt.ylabel('Surface brightness / W m$^{-2}$ arcsec$^{-2}$ $\mu$m$^{-1}$')
plt.title('Zodical light at North Ecliptic Pole - energy flux surface brightness units');


Conversion to photon flux surface brightness units

Divide by photon energies to convert to photon based units.


In [26]:
photon_energies = (h * c / sun_waves) / u.photon
zl_photons = zl_flux / photon_energies
zl_photons = zl_photons.to(u.photon * u.second**-1 * u.m**-2 * u.arcsecond**-2 *  u.micron**-1)

In [27]:
plt.semilogy(zl_waves, zl_photons, label='Zodical Light - NEP')
plt.xlim(0.2,2.5)
plt.ylim(2e-0,6e1)
plt.xlabel('Wavelength/$\mu$m')
plt.ylabel('Surface brightness / photons s$^{-1}$ m$^{-2}$ arcsec$^{-2}$ $\mu$m$^{-1}$')
plt.title('Zodical light at North Ecliptic Pole - photon flux surface brightness units');


Detected background

To convert the sky background brightness to electrons per second per pixel in the image sensor we multiply the sky background by the effective étendues then integrate over the wavelength range of each band.


In [28]:
zl_ep = []
for i, f in enumerate(filters):
    # Interpolate etendue data at same wavelengths as zodiacal light data
    eff_etendue =  np.interp(zl_waves, f['Wavelength'], eff_etendues[i]) * \
        eff_etendues[i].unit
    # Multiply zodical light surface brightness by effective etendue
    zl_epw = zl_photons * eff_etendue
    # Integrate over wavelength
    zl_ep.append(np.trapz(zl_epw, x=zl_waves))
    
zl_ep = u.Quantity(zl_ep)
zl_ep


Out[28]:
$[0.029449793,~0.26831068,~0.35355352,~0.25612122,~0.082758802] \; \mathrm{\frac{e^{-}}{pix\,s}}$

Dependence on position

The zodiacal light background calculated so far is for the North Ecliptic Pole. Tables 16 & 17 of Leinert et al (1997) present data showing how the average surface brightness varies with (ecliptic longitude - solar ecliptic longitude) and ecliptic latitude. To allow evaluation of the zodiacal light surface brightness at arbitrary positions on the sky we will interpolate these tabular data and use it to scale the North Ecliptic Pole value.

The zodical light is generally very smooth, effects such as resonant structures, comet tails and asteroid collisions result in variations of less than 1% which are largely confined to low ecliptic latitudes. We will not consider these here.

Leinert et al also present evidence for variations in the degree of reddening with solar elongation however at visible wavelengths this effect will be small unless the telescope is pointed close to the Sun. As positions close to the Sun would not be observed by ASE we will neglect this effect.


In [29]:
# Data from table 17, Leinert et al (1997).
llsun = np.array([0, 5 ,10, 15, 20, 25, 30, 35, 40, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180])
beta = np.array([0, 5, 10, 15, 20, 25, 30, 45, 60, 75])
zl_scale = np.array([[np.NaN, np.NaN, np.NaN, 3140, 1610, 985, 640, 275, 150, 100], \
                     [np.NaN, np.NaN, np.NaN, 2940, 1540, 945, 625, 271, 150, 100], \
                     [np.NaN, np.NaN, 4740, 2470, 1370, 865, 590, 264, 148, 100], \
                     [11500, 6780, 3440, 1860, 1110, 755, 525, 251, 146, 100], \
                     [6400, 4480, 2410, 1410, 910, 635, 454, 237, 141, 99], \
                     [3840, 2830, 1730, 1100, 749, 545, 410, 223, 136, 97], \
                     [2480, 1870, 1220, 845, 615, 467, 365, 207, 131, 95], \
                     [1650, 1270, 910, 680, 510, 397, 320, 193, 125, 93], \
                     [1180, 940, 700, 530, 416, 338, 282, 179, 120, 92], \
                     [910, 730, 555, 442, 356, 292, 250, 166, 116, 90], \
                     [505, 442, 352, 292, 243, 209, 183, 134, 104, 86], \
                     [338, 317, 269, 227, 196, 172, 151, 116, 93, 82], \
                     [259, 251, 225, 193, 166, 147, 132, 104, 86, 79], \
                     [212, 210, 197, 170, 150, 133, 119, 96, 82, 77], \
                     [188, 186, 177, 154, 138, 125, 113, 90, 77, 74], \
                     [179, 178, 166, 147, 134, 122, 110, 90, 77, 73], \
                     [179, 178, 165, 148, 137, 127, 116, 96, 79, 72], \
                     [196, 192, 179, 165, 151, 141, 131, 104, 82, 72], \
                     [230, 212, 195, 178, 163, 148, 134, 105, 83, 72]]).transpose()

# Normalise to a value of 1.0 at the NEP
zl_scale /= 77

# Expand range of angles to cover the full sphere by using symmetry
beta = np.array(np.concatenate((-np.flipud(beta)[:-1], beta))) * u.degree
llsun = np.array(np.concatenate((llsun, 360 - np.flipud(llsun)[1:-1]))) * u.degree
zl_scale = np.concatenate((np.flipud(zl_scale)[:-1],zl_scale))
zl_scale = np.concatenate((zl_scale,np.fliplr(zl_scale)[:,1:-1]),axis=1)

# Convert angles to radians within the required ranges.
beta = beta.to(u.radian).value + np.pi/2
llsun = llsun.to(u.radian).value

The Scipy spherical interpolation functions can't handle the missing values for positions near to the Sun, either the interpolation fails to run or there are ringing artefacts. To overcome this we use an initial cartesian interpolation pass on a subset of the data near to the sun in order to smoothly fill in the mising values, then use spherical interpolation on the full, filled in data set.


In [30]:
# For initial cartesian interpolation want the hole in the middle of the data,
# i.e. want to remap longitudes onto -180 to 180, not 0 to 360 degrees.
# Only want the region of closely spaced data point near to the Sun.
beta_c = beta[3:-3]
nl = len(llsun)
llsun_c = np.concatenate((llsun[nl//2+9:]-2*np.pi,llsun[:nl//2-8]))
zl_scale_c = np.concatenate((zl_scale[3:-3,nl//2+9:],zl_scale[3:-3,:nl//2-8]),axis=1)

# Convert everthing to 1D arrays (lists) of x, y and z coordinates.
llsuns, betas = np.meshgrid(llsun_c, beta_c)
llsuns_c = llsuns.ravel()
betas_c = betas.ravel()
zl_scale_cflat = zl_scale_c.ravel()

# Indices of the non-NaN points
good = np.where(np.isfinite(zl_scale_cflat))

# 2D cartesian interpolation function
zl_scale_c_interp = SmoothBivariateSpline(betas_c[good], llsuns_c[good], zl_scale_cflat[good])

# Calculate interpolated values
zl_scale_fill = zl_scale_c_interp(beta_c, llsun_c)

# Remap the interpolated values back to original ranges of coordinates
zl_patch = np.zeros(zl_scale.shape)
zl_patch[3:16,0:10] = zl_scale_fill[:,9:]
zl_patch[3:16,-9:] = zl_scale_fill[:,:9]

# Fill the hole in the original data with values from the cartesian interpolation
zl_patched = np.where(np.isfinite(zl_scale),zl_scale,zl_patch)

# Spherical interpolation function from the full, filled data set
zl_scale_interp = RectSphereBivariateSpline(beta, llsun, zl_patched, \
                                            pole_continuity=(False,False), pole_values=(1.0, 1.0), \
                                            pole_exact=True, pole_flat=False)

In [31]:
# Check the process with some plots
n_rows = 4
n_cols = 2

plt.subplot(n_rows,n_cols,1)
plt.imshow(zl_scale, norm=colors.LogNorm(vmin=0.7,vmax=300), \
           interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Leinert et al data, normalised, full sphere')

plt.subplot(n_rows,n_cols,3)
plt.imshow(zl_scale_c, norm=colors.LogNorm(vmin=0.7,vmax=300), interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Centred and clipped')

plt.subplot(n_rows,n_cols,5)
plt.imshow(zl_scale_fill, norm=colors.LogNorm(vmin=0.7,vmax=300), interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Centred, clipped, cartesian interpolated')

plt.subplot(n_rows,n_cols,7)
plt.imshow(zl_patch, norm=colors.LogNorm(vmin=0.7,vmax=300), interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Cartesian interpolation remapped to original coords')

plt.subplot(n_rows,n_cols,2)
plt.imshow(zl_patched, norm=colors.LogNorm(vmin=0.7,vmax=300), interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Leinert et al data, filled with cartesian interpolation')

plt.subplot(n_rows,n_cols,4)
plt.imshow(zl_scale_interp(beta, llsun), norm=colors.LogNorm(vmin=0.7,vmax=300), \
           interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Spherical interpolation evaluated at original points')

plt.subplot(n_rows,n_cols,6)
plt.imshow((zl_scale_interp(beta, llsun) - zl_scale)/zl_scale, interpolation='none',)
plt.colorbar()
plt.title('Relative deviation of interpolation from original')

plt.subplot(n_rows,n_cols,8)
plt.imshow(zl_scale_interp(np.linspace(0,np.pi,180),np.linspace(0,2*np.pi,360)), \
           origin='lower', extent=[0,360,-90,90], norm=colors.LogNorm(vmin=0.7,vmax=300), \
           interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Spherical interpolation')
plt.xlabel('$\lambda - \lambda_{\odot}$')
plt.ylabel('$\\beta$')

plt.tight_layout()
plt.gcf().set_size_inches(12,12)


/home/ajh/Documents/virtualenvs/sci-latest-python3.4.1/lib/python3.4/site-packages/matplotlib/colors.py:994: RuntimeWarning: invalid value encountered in less_equal
  mask |= resdat <= 0

We will precalculate an all-sky map of the zodical light surface brightness distribution using the HEALPix pixelation (http://healpix.sourceforge.net/). We can use the healpy tools to extract projections of this map for signal to noise calculations.

Want to calculate these in galactic coordinates for later combination with dust maps


In [31]:
# Calculate all sky map
# *WARNING* This takes a long time, generally will just want to load it from zl_map.fits instead,
# see cell two below.
nside = 2048
zl_map=np.empty(hp.nside2npix(nside))
for i in range(zl_map.size):
    theta, phi = hp.pix2ang(nside, i)
    zl_map[i] = zl_scale_interp(theta, phi)


---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-31-1a97ce598b86> in <module>()
      3 zl_map=np.empty(hp.nside2npix(nside))
      4 for i in range(zl_map.size):
----> 5     theta, phi = hp.pix2ang(nside, i)
      6     zl_map[i] = zl_scale_interp(theta, phi)

/home/ajh/Documents/virtualenvs/sci-latest-python3.4.1/lib/python3.4/site-packages/healpy/pixelfunc.py in pix2ang(nside, ipix, nest)
    421     (array([ 2.30052398,  0.84106867,  0.41113786,  0.2044802 ]), array([ 5.49778714,  5.89048623,  5.89048623,  5.89048623]))
    422     """
--> 423     check_nside(nside)
    424     if nest:
    425         return pixlib._pix2ang_nest(nside, ipix)

/home/ajh/Documents/virtualenvs/sci-latest-python3.4.1/lib/python3.4/site-packages/healpy/pixelfunc.py in check_nside(nside)
   1014 def check_nside(nside):
   1015     """Raises exception is nside is not valid"""
-> 1016     if not np.all(isnsideok(nside)):
   1017         raise ValueError("%s is not a valid nside parameter (must be a power of 2, less than 2**30)" % str(nside))
   1018 

/home/ajh/Documents/virtualenvs/sci-latest-python3.4.1/lib/python3.4/site-packages/numpy/core/fromnumeric.py in all(a, axis, out, keepdims)
   1916 
   1917     """
-> 1918     arr = asanyarray(a)
   1919 
   1920     try:

KeyboardInterrupt: 

In [32]:
# Write the zodiacal light map to file (it takes a long time to calculate, don't want to have to redo it all the time).
hp.write_map('resources/zl_map.fits', zl_map, fits_IDL=False, coord='E', column_names=('ZL_rel_SB',))


WARNING: AstropyDeprecationWarning: The new_table function is deprecated and may be removed in a future version.
        Use :meth:`BinTableHDU.from_columns` for new BINARY tables or :meth:`TableHDU.from_columns` for new ASCII tables instead. [healpy.fitsfunc]
WARNING:astropy:AstropyDeprecationWarning: The new_table function is deprecated and may be removed in a future version.
        Use :meth:`BinTableHDU.from_columns` for new BINARY tables or :meth:`TableHDU.from_columns` for new ASCII tables instead.
WARNING: AstropyDeprecationWarning: The use of header.update() to add new keywords to a header is deprecated.  Instead, use either header.set() or simply `header[keyword] = value` or `header[keyword] = (value, comment)`.  header.set() is only necessary to use if you also want to use the before/after keyword arguments. [astropy.io.fits.header]
WARNING:astropy:AstropyDeprecationWarning: The use of header.update() to add new keywords to a header is deprecated.  Instead, use either header.set() or simply `header[keyword] = value` or `header[keyword] = (value, comment)`.  header.set() is only necessary to use if you also want to use the before/after keyword arguments.

In [34]:
# Read zodical light map from file
zl_map = hp.read_map('resources/zl_map.fits')


NSIDE = 2048
ORDERING = RING in fits file

In [35]:
hp.mollview(zl_map, coord=('E', 'C'), cmap='cubehelix_r', min=0.7, max=10, norm='log', \
            title='Zodiacal light surface brightness distribution, March equinox', \
            unit='Surface brightness relative to NEP')
hp.graticule(coord='C', local=True)

hp.projscatter(0, 90, lonlat=True, marker='x', coord='E')
hp.projtext(0, 90, 'NEP', lonlat=True, coord='E')
hp.projscatter(0, -90, lonlat=True, marker='x', coord='E')
hp.projtext(0, -90, 'SEP', lonlat=True, coord='E')

plt.gcf().set_size_inches(12,7.5)
plt.savefig('ZL_March.png')


0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.

Seasonal variations

The main effect is simply due to the movement of the Sun around the ecliptic equator over the course of the year. As a result the distribution of zodiacal light on the sky rotates about an axis through the ecliptic poles.

In addition the zodiacal light exhibits seasonal variations due to the orbit of the Earth taking it above and below the symmetry plane of the zodiacal dust. At high ecliptic latitudes this manifests as an approximately 10% sinusoidal variation in surface brightness. For now we will ignore this.

Otherwise the zodiacal light is essentially time invariant. Resonant structures, comet tails and asteroid collisions result in variations of less than 1% and are largely confined to low ecliptic latitude.


In [36]:
# Very wasteful of RAM but for now just implement this using rotated copies of the ZL map

def map_rotate_polar(input_map, angle):
    '''
    Function to return a copy of a healpy map rotated about its polar axis
    '''
    npix = input_map.size
    nside = hp.npix2nside(npix)
    # Get coordinates of each pixel in original map
    indices = np.arange(npix)
    theta,phi = hp.pix2ang(nside, indices)
    # New empty map for rotated copy
    output_map = np.empty_like(input_map)
    # Rotator object to transform pixel coordinates
    rotator = hp.Rotator(deg=True, rot=[angle,0])
    # For each pixel in original map calculate new coordinates and copy
    # into appropriate pixel in new map.
    output_map[hp.ang2pix(nside, *rotator(theta, phi))] = input_map
    return output_map

In [37]:
zl_map_90 = map_rotate_polar(zl_map, 90)
zl_map_180 = map_rotate_polar(zl_map, 180)
zl_map_270 = map_rotate_polar(zl_map, 270)

In [38]:
hp.mollview(zl_map, coord=('E', 'C'), cmap='cubehelix_r', min=0.7, max=10, \
            norm='log', sub=(2,2,1), cbar=False, \
            title='Zodiacal light surface brightness distribution, March equinox', \
            unit='Surface brightness relative to NEP')
hp.graticule(coord='C', local=True)
hp.projscatter(0, 90, lonlat=True, marker='x', coord='E')
hp.projtext(0, 90, 'NEP', lonlat=True, coord='E')
hp.projscatter(0, -90, lonlat=True, marker='x', coord='E')
hp.projtext(0, -90, 'SEP', lonlat=True, coord='E')

hp.mollview(zl_map_270, coord=('E', 'C'), cmap='cubehelix_r', min=0.7, max=10, \
            norm='log', sub=(2,2,2), cbar=False, \
            title='Zodiacal light surface brightness distribution, June solstice', \
            unit='Surface brightness relative to NEP')
hp.graticule(coord='C', local=True)
hp.projscatter(0, 90, lonlat=True, marker='x', coord='E')
hp.projtext(0, 90, 'NEP', lonlat=True, coord='E')
hp.projscatter(0, -90, lonlat=True, marker='x', coord='E')
hp.projtext(0, -90, 'SEP', lonlat=True, coord='E')

hp.mollview(zl_map_180, coord=('E', 'C'), cmap='cubehelix_r', min=0.7, max=10, \
            norm='log', sub=(2,2,3), cbar=False, \
            title='Zodiacal light surface brightness distribution, September equinox', \
            unit='Surface brightness relative to NEP')
hp.graticule(coord='C', local=True)
hp.projscatter(0, 90, lonlat=True, marker='x', coord='E')
hp.projtext(0, 90, 'NEP', lonlat=True, coord='E')
hp.projscatter(0, -90, lonlat=True, marker='x', coord='E')
hp.projtext(0, -90, 'SEP', lonlat=True, coord='E')

hp.mollview(zl_map_90, coord=('E', 'C'), cmap='cubehelix_r', min=0.7, max=10, \
            norm='log', sub=(2,2,4), cbar=False, \
            title='Zodiacal light surface brightness distribution, December solstice', \
            unit='Surface brightness relative to NEP')
hp.graticule(coord='C', local=True)
hp.projscatter(0, 90, lonlat=True, marker='x', coord='E')
hp.projtext(0, 90, 'NEP', lonlat=True, coord='E')
hp.projscatter(0, -90, lonlat=True, marker='x', coord='E')
hp.projtext(0, -90, 'SEP', lonlat=True, coord='E')

plt.gcf().set_size_inches(12,7.5)
plt.savefig('ZL_seasonal.png')


0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.