Skip to content

Commit

Permalink
Merge f5dd9d7 into 6a43f21
Browse files Browse the repository at this point in the history
  • Loading branch information
jlustigy committed Apr 12, 2019
2 parents 6a43f21 + f5dd9d7 commit 4b3f2de
Show file tree
Hide file tree
Showing 18 changed files with 462 additions and 223 deletions.
105 changes: 26 additions & 79 deletions README.md
@@ -1,91 +1,38 @@
# coronagraph
<p align="center">
<img width = "450" src="https://github.com/jlustigy/coronagraph/blob/luvoir_drm/docs/_static/logo6.png?raw=true"/>
<br>
A Python noise model for directly imaging exoplanets with coronagraph equipped telescopes.
</p>

A Python noise model for directly imaging exoplanets with a space based coronagraph.
<p align="center">
<a href="https://travis-ci.org/jlustigy/coronagraph"><img src="http://img.shields.io/travis/jlustigy/coronagraph/master.svg?style=flat"/></a>
<a href="http://coronagraph.readthedocs.io/en/latest/"><img src="https://readthedocs.org/projects/coronagraph/badge/?version=latest"/></a>
<a href="https://github.com/jlustigy/coronagraph/blob/master/LICENSE"><img src="http://img.shields.io/badge/license-MIT-blue.svg?style=flat"/></a>
<a href="https://coveralls.io/github/jlustigy/coronagraph?branch=master"><img src="https://coveralls.io/repos/github/jlustigy/coronagraph/badge.svg?branch=master&style=flat"/></a>
</p>

[![Build Status](http://img.shields.io/travis/jlustigy/coronagraph/master.svg?style=flat)](https://travis-ci.org/jlustigy/coronagraph)
[![Docs Status](https://readthedocs.org/projects/coronagraph/badge/?version=latest)](http://coronagraph.readthedocs.io/en/latest/)
[![MIT licensed](http://img.shields.io/badge/license-MIT-blue.svg?style=flat)](https://github.com/jlustigy/coronagraph/blob/master/LICENSE)
[![Coverage Status](https://coveralls.io/repos/github/jlustigy/coronagraph/badge.svg?branch=master&style=flat)](https://coveralls.io/github/jlustigy/coronagraph?branch=master)
Documentation
=============

If you use this model in your own research please cite [Robinson et al (2016)](http://adsabs.harvard.edu/cgi-bin/bib_query?arXiv:1507.00777) and include the following acknowledgement: "This work made use of the Python coronagraph noise model, developed by J. Lustig-Yaeger and available at https://github.com/jlustigy/coronagraph/".

## Install
Please refer to the [code documentation](http://coronagraph.readthedocs.io/en/latest/) on Read the Docs.

* Clone this github repository:
```bash
git clone git@github.com:jlustigy/coronagraph.git
```
* (optional) Add to python path to use `coronagraph` in any dirctory. This is ideal if you want to be able to access the core coronagraph functions from different project directories. If you're on a mac, then open up the `.bash_profile` file and add the following lines:
```bash
# create python path if it isn't already an environment variable
export PYTHONPATH=$PATH
# add this repo to the python path
export PYTHONPATH=${PYTHONPATH}:/Users/Name/Folder/coronagraph/
```
If you use the `.cshrc` file it might look something like this:
```shell
# For .cshrc (I know, terrible...)
setenv PYTHONPATH ${PYTHONPATH}:/Users/Name/Folder/coronagraph/
```
Installation
============

## Examples
The simplest option is to install using `pip`:

#### Running `luvoir_demo.py`
You can start to get a feel for how the coronagraph noise model works by examining and running `luvoir_demo.py`. From the command line type:
```bash
python luvoir_demo.py
pip install coronagraph
```
You'll get an output file that looks like this:

<img src="https://github.com/jlustigy/coronagraph/blob/master/scripts/plots/luvoir_demo_.png" width="60%" height="60%" align="middle" />

Now you can create your own scripts like this to play with all the tunable parameters!

#### Simulate observation with the Integral Field Spectrograph (IFS)

```python
# Import coronagraph package
import coronagraph as cg

# Initialize Telescope, Planet, and Star objects
telescope = cg.Telescope()
planet = cg.Planet()
star = cg.Star()

# Read-in wavelength, reflectance model
model = np.loadtxt('planets/earth_quadrature_radiance_refl.dat', skiprows=8)
lam = model[:,0] # wavelength (microns)
refl = np.pi * model[:,3] # geometric albedo
solhr = model[:,2] # solar flux
Users may also clone the repository on Github:

# Specify telescope integration time in hours
integration_time = 10.0

# Observe!
lam, dlam, Cratio, spec, sig, SNR = \
cg.generate_observation(lam, refl, solhr, integration_time, telescope, planet, star)

```
<img src="https://github.com/jlustigy/coronagraph/blob/master/scripts/plots/earth_quad_R70.png" width="100%" height="100%" align="middle" />

#### Simulate observation with the Imaging camera

```python
# Set telescope to 'Imaging' mode
telescope.mode = 'Imaging'

# Load Filter Wheel for obsevation (the default filters are the Johnson-Counsins UBVRI filters)
landsat = cg.filters.landsat()
jc = cg.filters.johnson_cousins2()

# Add Filter Wheel to Telescope
telescope.filter_wheel = jc

# Observe!
lam, spec, sig = cg.generate_observation(lam, refl, integration_time, telescope, planet, star)
```bash
git clone https://github.com/jlustigy/coronagraph.git
cd coronagraph
python setup.py
```
<img src="https://github.com/jlustigy/coronagraph/blob/master/scripts/plots/earth_quad_jc.png" width="100%" height="100%" align="middle" />

## Notes

* Check out the [SVO Filter Profile Service](http://svo2.cab.inta-csic.es/svo/theory/fps3/index.php?id=2MASS/2MASS.J&&mode=browse&gname=2MASS&gname2=2MASS#filter)
Note
====
If you use this model in your own research please cite [Robinson et al (2016)](http://adsabs.harvard.edu/cgi-bin/bib_query?arXiv:1507.00777) and include the following acknowledgement: "This work made use of the Python coronagraph noise model, developed by J. Lustig-Yaeger and available at https://github.com/jlustigy/coronagraph/".
3 changes: 2 additions & 1 deletion coronagraph/Noise.py
Expand Up @@ -7,7 +7,7 @@ class Output(object):

def __init__(self, lam=None, dlam=None, A=None, q=None, Cratio=None,
cp=None, csp=None, cz=None, cez=None, cD=None, cR=None,
cth=None, DtSNR=None):
cth=None, cc=None, DtSNR=None):

self.lam = lam
self.dlam = dlam
Expand All @@ -21,4 +21,5 @@ def __init__(self, lam=None, dlam=None, A=None, q=None, Cratio=None,
self.cD = cD
self.cR = cR
self.cth = cth
self.cc = cc
self.DtSNR = DtSNR
113 changes: 82 additions & 31 deletions coronagraph/count_rates.py
Expand Up @@ -133,6 +133,8 @@ def run_count_rates(self, Ahr, lamhr, solhr):
Dark current photon count rate [photons/s]
cR : array
Read noise photon count rate [photons/s]
cc : array
Clock induced charge photon count rate [photons/s]
cb : array
Total background photon noise count rate [photons/s]
DtSNR : array
Expand Down Expand Up @@ -165,7 +167,7 @@ def run_count_rates(self, Ahr, lamhr, solhr):
assert False, "telescope.aperture is invalid"

# Call count_rates
lam, dlam, A, q, Cratio, cp, csp, cz, cez, cD, cR, cth, DtSNR = \
lam, dlam, A, q, Cratio, cp, csp, cz, cez, cD, cR, cth, cc, DtSNR = \
count_rates(Ahr, lamhr, solhr,
alpha = self.planet.alpha,
Phi = self.planet.Phi,
Expand All @@ -191,14 +193,20 @@ def run_count_rates(self, Ahr, lamhr, solhr):
De = self.telescope.darkcurrent,
DNHpix = self.telescope.DNHpix,
Re = self.telescope.readnoise,
Rc = self.telescope.Rc,
Dtmax = self.telescope.Dtmax,
X = self.telescope.X,
qe = self.telescope.qe,
MzV = self.planet.MzV,
MezV = self.planet.MezV,
A_collect = self.telescope.A_collect,
diam_circumscribed = self.telescope.diam_circumscribed,
diam_inscribed = self.telescope.diam_inscribed,
lam = self.telescope.lam,
dlam = self.telescope.dlam,
Tput_lam = self.telescope.Tput_lam,
qe_lam = self.telescope.qe_lam,
lammin_lenslet = self.telescope.lammin_lenslet,
NIR = self.NIR,
GROUND = self.GROUND,
THERMAL = self.THERMAL,
Expand All @@ -208,19 +216,20 @@ def run_count_rates(self, Ahr, lamhr, solhr):
)

# Save output arrays
self.lam = lam
self.dlam = dlam
self.A = A
self.Cratio = Cratio
self.cp = cp
self.csp =csp
self.cz = cz
self.cez = cez
self.cD = cD
self.cR = cR
self.cth = cth
self.cb = cz + cez + csp + cD + cR + cth
self.DtSNR = DtSNR
self.lam = lam
self.dlam = dlam
self.A = A
self.Cratio = Cratio
self.cp = cp
self.csp = csp
self.cz = cz
self.cez = cez
self.cD = cD
self.cR = cR
self.cth = cth
self.cc = cc
self.cb = cz + cez + csp + cD + cR + cth + cc
self.DtSNR = DtSNR

# Flip the switch
self._computed = True
Expand Down Expand Up @@ -476,14 +485,20 @@ def count_rates(Ahr, lamhr, solhr,
De = 1e-4,
DNHpix = 3.0,
Re = 0.1,
Rc = 0.0,
Dtmax = 1.0,
X = 1.5,
qe = 0.9,
MzV = 23.0,
MezV = 22.0,
A_collect = None,
diam_circumscribed = None,
diam_inscribed = None,
lam = None,
dlam = None,
Tput_lam = None,
qe_lam = None,
lammin_lenslet = None,
wantsnr=10.0, FIX_OWA = False, COMPUTE_LAM = False,
SILENT = False, NIR = False, THERMAL = False, GROUND = False,
vod=False, set_fpa=None, CIRC = True, roll_maneuver = True):
Expand Down Expand Up @@ -547,6 +562,8 @@ def count_rates(Ahr, lamhr, solhr,
Number of horizontal/spatial pixels for dispersed spectrum
Re : float, optional
Read noise counts per pixel
Rc : float, optional
Clock induced charge [counts/pixel/photon]
Dtmax : float, optional
Detector maximum exposure time [hours]
X : float, optional
Expand All @@ -559,10 +576,24 @@ def count_rates(Ahr, lamhr, solhr,
V-band exozodiacal light surface brightness [mag/arcsec**2]
A_collect : float, optional
Mirror collecting area (m**2) (uses :math:`\pi(D/2)^2` by default)
diam_circumscribed : float, optional
Circumscribed telescope diameter [m] used for IWA and OWA (uses `diam`
if `None` provided)
diam_inscribed : float, optional
Inscribed telescope diameter [m] used for lenslet calculations
(uses `diam` if `None` provided)
lam : array-like, optional
Wavelength grid for spectrograph [microns] (uses ``lammin``, ``lammax``,
and ``resolution`` to determine if ``None`` provided)
dlam : array-like, optional
Wavelength grid `widths` for spectrograph [microns] (uses ``lammin``, ``lammax``,
and ``resolution`` to determine if ``None`` provided)
Tput_lam : tuple of arrays
Wavelength-dependent throughput e.g. ``(wls, tputs)``
qe_lam : tuple of arrays
Wavelength-dependent throughput e.g. ``(wls, qe)``
lammin_lenslet : float, optional
Minimum wavelength to use for lenslet calculation (default is ``lammin``)
wantsnr : float, optional
Desired signal-to-noise ratio in each pixel
FIX_OWA : bool, optional
Expand Down Expand Up @@ -605,32 +636,42 @@ def count_rates(Ahr, lamhr, solhr,
Cratio : ndarray
Planet-star contrast ratio
cp : ndarray
Planetary photon count rate on detector
Planetary photon count rate on detector [1/s]
csp : ndarray
Speckle photon count rate on detector
Speckle photon count rate on detector [1/s]
cz : ndarray
Zodiacal photon count rate on detector
Zodiacal photon count rate on detector [1/s]
cez : ndarray
Exozodiacal photon count rate on detector
Exozodiacal photon count rate on detector [1/s]
cD : ndarray
Dark current photon count rate on detector
Dark current photon count rate on detector [1/s]
cR : ndarray
Read noise photon count rate on detector
Read noise photon count rate on detector [1/s]
cth : ndarray
Instrument thermal photon count rate on detector
Instrument thermal photon count rate on detector [1/s]
cc : ndarray
Clock induced charge photon count rate [1/s]
DtSNR : ndarray
Exposure time required to get desired S/N (wantsnr) [hours]
"""

convolution_function = downbin_spec
#convolution_function = degrade_spec

# Define a diameter for IWA (inscribed area) and collecting area
diam_inscribed = diam
# Define a diameter for IWA (circumscribed),
# collecting area, and lenslet (inscribed)
if diam_inscribed is None:
# Defaults to diam
diam_inscribed = diam
if A_collect is None:
# Defaults to diam
diam_collect = diam
else:
# Calculated from provided collecting area
diam_collect = 2. * (A_collect / np.pi)**0.5
if diam_circumscribed is None:
# Defaults to diam
diam_circumscribed = diam

# Configure for different telescope observing modes
if mode == 'Imaging':
Expand Down Expand Up @@ -658,8 +699,10 @@ def count_rates(Ahr, lamhr, solhr,
fpa = set_fpa * f_airy(X)

# Set wavelength grid
# GENERALIZE THIS:
if COMPUTE_LAM:
lam, dlam = construct_lam(lammin, lammax, Res)
if (lam is None) or (dlam is None):
lam, dlam = construct_lam(lammin, lammax, Res)
elif IMAGE:
pass
else:
Expand All @@ -675,11 +718,12 @@ def count_rates(Ahr, lamhr, solhr,
Re = set_read_noise(lam, Re, NIR=NIR)

# Set Angular size of lenslet
theta = set_lenslet(lam, lammin, diam_inscribed, X, NIR=True)
if lammin_lenslet is None: lammin_lenslet = lammin
theta = set_lenslet(lam, lammin_lenslet, diam_inscribed, X, NIR=True)

# Set throughput (for inner and outer working angle cutoffs)
sep = r/d*np.sin(alpha*np.pi/180.)*np.pi/180./3600. # separation in radians
T = set_throughput(lam, Tput, diam_inscribed, sep, IWA, OWA, lammin, FIX_OWA=FIX_OWA, SILENT=SILENT)
T = set_throughput(lam, Tput, diam_circumscribed, sep, IWA, OWA, lammin, FIX_OWA=FIX_OWA, SILENT=SILENT)

# Apply wavelength-dependent throuput, if needed
if Tput_lam is not None:
Expand Down Expand Up @@ -708,7 +752,7 @@ def count_rates(Ahr, lamhr, solhr,

# Degrade albedo and stellar spectrum
if COMPUTE_LAM:
A = convolution_function(Ahr,lamhr,lam,dlam=dlam)
A = convolution_function(Ahr, lamhr, lam, dlam=dlam)
Fs = convolution_function(solhr, lamhr, lam, dlam=dlam)
elif IMAGE:
# Convolve with filter response
Expand All @@ -727,7 +771,7 @@ def count_rates(Ahr, lamhr, solhr,
cp = cplan(q, fpa, T, lam, dlam, Fp, diam_collect) # planet count rate
cz = czodi(q, X, T, lam, dlam, diam_collect, MzV) # solar system zodi count rate
cez = cezodi(q, X, T, lam, dlam, diam_collect, r, \
Fstar(lam,Teff,Rs,1.,AU=True), Nez, MezV) # exo-zodi count rate
Fstar(lam, Teff, Rs,1. , AU=True), Nez, MezV) # exo-zodi count rate
csp = cspeck(q, T, C, lam, dlam, Fstar(lam,Teff,Rs,d), diam_collect) # speckle count rate
cD = cdark(De, X, lam, diam_collect, theta, DNHpix, IMAGE=IMAGE) # dark current count rate
cR = cread(Re, X, lam, diam_collect, theta, DNHpix, Dtmax, IMAGE=IMAGE) # readnoise count rate
Expand Down Expand Up @@ -759,8 +803,15 @@ def count_rates(Ahr, lamhr, solhr,
plt.show()
'''

# Clock induced charge photon count rate
# Calculate photon count rate in the scene (everything except read noise)
cscene = cp + cz + cez + csp + cD + cth
# Calculate the clock induced charge photon count rate
cc = ccic(Rc, cscene, X, lam, diam_collect, theta, DNHpix, Dtmax,
IMAGE=IMAGE, CIRC=CIRC)

# Calculate total background counts
cb = (cz + cez + csp + cD + cR + cth)
cb = (cz + cez + csp + cD + cR + cth + cc)

# Use telescope roll maneuver
if roll_maneuver:
Expand All @@ -774,7 +825,7 @@ def count_rates(Ahr, lamhr, solhr,
cnoise = cp + roll_factor*cb

# Calculate total counts
ctot = cp + cz + cez + csp + cD + cR + cth
ctot = cp + cz + cez + csp + cD + cR + cth + cc

'''
Giada: where does the factor of 2 come from [above]?
Expand All @@ -796,4 +847,4 @@ def count_rates(Ahr, lamhr, solhr,
# Exposure time to SNR
DtSNR = exptime_element(lam, cp, cnoise, wantsnr)

return lam, dlam, A, q, Cratio, cp, csp, cz, cez, cD, cR, cth, DtSNR
return lam, dlam, A, q, Cratio, cp, csp, cz, cez, cD, cR, cth, cc, DtSNR

0 comments on commit 4b3f2de

Please sign in to comment.