"""Produces a graph of the visibility & accessible position angles
for a given RA & DEC, and prints out corresponding information,
including the ranges of accessible and inaccessible PAs.
Usage: python visibilityPA.py RA DEC [targetName]
if targetName is specified, then the figure is saved
-Created by David Lafreniere, March 2016
-makes use of (and hacks) several scripts created by Pierre Ferruit
that are part of the JWST Python tools JWSTpylib and JWSTpytools
"""
import datetime
import math
import pkg_resources
from astropy.table import Table
from astropy.time import Time
from bokeh.models import HoverTool, ranges
from bokeh.plotting import figure, ColumnDataSource
import matplotlib.dates as mdates
import numpy as np
from . import ephemeris_old2x as EPH
# from jwst_gtvt.find_tgt_info import get_table
from ..utils import blockPrint, enablePrint
D2R = math.pi / 180. # degrees to radians
R2D = 180. / math.pi # radians to degrees
[docs]
def checkVisPA(ra, dec, targetName=None, ephFileName=None, fig=None):
"""Check the visibility at a range of position angles.
Parameters
----------
ra: str
The RA of the target in hh:mm:ss.s or dd:mm:ss.s or representing
a float
dec: str
The Dec of the target in hh:mm:ss.s or dd:mm:ss.s or
representing a float
targetName: str
The target name
ephFileName: str
The filename of the ephemeris file
fig: bokeh.plotting.figure
The figure to plot to
Returns
-------
paGood : float
The good position angle.
paBad : float
The bad position angle.
gd : matplotlib.dates object
The greogrian date.
fig : bokeh.plotting.figure object
The plotted figure.
"""
if ephFileName is None:
eph_file = 'data/contam_visibility/JWST_ephem_short.txt'
ephFileName = pkg_resources.resource_filename('exoctk', eph_file)
if ra.find(':') > -1: # format is hh:mm:ss.s or dd:mm:ss.s
ra = convert_ddmmss_to_float(ra) * 15. * D2R
dec = convert_ddmmss_to_float(dec) * D2R
else: # format is decimal
ra = float(ra) * D2R
dec = float(dec) * D2R
# load ephemeris
eclFlag = False
eph = EPH.Ephemeris(ephFileName, eclFlag)
# convert dates from MJD to Gregorian calendar dates
mjd = np.array(eph.datelist)
d = mdates.julian2num(mjd + 2400000.5)
gd = mdates.num2date(d)
# loop through dates and determine VIS and PAs (nominal, min, max)
vis = np.empty(mjd.size, dtype=bool)
paNom = np.empty(mjd.size)
paMin = np.empty(mjd.size)
paMax = np.empty(mjd.size)
for i in range(mjd.size):
# is it visible?
vis[i] = eph.in_FOR(mjd[i], ra, dec)
# nominal PA at this date
pa = eph.normal_pa(mjd[i], ra, dec)
# search for minimum PA allowed by roll
pa0 = pa
while eph.is_valid(mjd[i], ra, dec, pa0 - 0.002):
pa0 -= 0.002
# search for maximum PA allowed by roll
pa1 = pa
while eph.is_valid(mjd[i], ra, dec, pa1 + 0.002):
pa1 += 0.002
paNom[i] = (pa * R2D) % 360
paMin[i] = (pa0 * R2D) % 360
paMax[i] = (pa1 * R2D) % 360
# does PA go through 360 deg?
wrap = np.any(np.abs(np.diff(paNom[np.where(vis)[0]])) > 350)
# Determine good and bad PA ranges
# Good PAs
i, = np.where(vis)
pa = np.concatenate((paNom[i], paMin[i], paMax[i]))
if wrap:
pa = np.append(pa, (0., 360.))
pa.sort()
i1, = np.where(np.diff(pa) > 10)
i0 = np.insert(i1 + 1, 0, 0)
i1 = np.append(i1, -1)
paGood = np.dstack((pa[i0], pa[i1])).round(1).reshape(-1, 2).tolist()
# bad PAs (complement of the good PAs)
paBad = []
if paGood[0][0] > 0:
paBad.append([0., paGood[0][0]])
for i in range(1, len(paGood)):
paBad.append([paGood[i - 1][1], paGood[i][0]])
if paGood[-1][1] < 360.:
paBad.append([paGood[-1][1], 360.])
# Make a figure
if fig is None or fig:
tools = 'crosshair, reset, hover, save'
radec = ', '.join([str(ra), str(dec)])
fig = figure(tools=tools, width=800, height=400,
x_axis_type='datetime',
title=targetName or radec)
# Do all figure calculations
iBad, = np.where(vis == False)
paMasked = np.copy(paNom)
paMasked[iBad] = np.nan
gdMasked = np.copy(gd)
i = np.argmax(paNom)
if paNom[i + 1] < 10:
i += 1
paMasked = np.insert(paMasked, i, np.nan)
gdMasked = np.insert(gdMasked, i, gdMasked[i])
i = np.argmax(paMin)
goUp = paMin[i - 2] < paMin[i - 1] # PA going up at wrap point?
# Top part
i0_top = 0 if goUp else i
i1_top = i if goUp else paMin.size - 1
paMaxTmp = np.copy(paMax)
paMaxTmp[np.where(paMin > paMax)[0]] = 360
# Bottom part
i = np.argmin(paMax)
i0_bot = i if goUp else 0
i1_bot = paMin.size - 1 if goUp else i
paMinTmp = np.copy(paMin)
paMinTmp[np.where(paMin > paMax)[0]] = 0
# Convert datetime to a number for Bokeh
gdMaskednum = [datetime.date(2019, 6, 1) + datetime.timedelta(days=n)
for n, d in enumerate(gdMasked)]
color = 'green'
# Draw the curve and error
try:
fig.line(
gdMaskednum,
paMasked,
legend_label='cutoff',
line_color=color)
except AttributeError:
fig.line(gdMaskednum, paMasked, legend='cutoff', line_color=color)
# Top
terr_y = np.concatenate([paMin[i0_top:i1_top + 1],
paMaxTmp[i0_top:i1_top + 1][::-1]])
terr_x = np.concatenate([gdMaskednum[i0_top:i1_top + 1],
gdMaskednum[i0_top:i1_top + 1][::-1]])
fig.patch(terr_x, terr_y, color=color, fill_alpha=0.2, line_alpha=0)
# Bottom
berr_y = np.concatenate([paMinTmp[i0_bot:i1_bot + 1],
paMax[i0_bot:i1_bot + 1][::-1]])
berr_x = np.concatenate([gdMaskednum[i0_bot:i1_bot + 1],
gdMaskednum[i0_bot:i1_bot + 1][::-1]])
fig.patch(berr_x, berr_y, color='red', fill_alpha=0.2, line_alpha=0)
# Plot formatting
fig.xaxis.axis_label = 'Date'
fig.yaxis.axis_label = 'Aperture Position Angle (degrees)'
return paGood, paBad, gd, fig
[docs]
def using_gtvt(ra, dec, instrument, targetName='noName', ephFileName=None, output='bokeh'):
"""Plot the visibility (at a range of position angles) against time.
Parameters
----------
ra : str
The RA of the target (in degrees) hh:mm:ss.s or dd:mm:ss.s or
representing a float
dec : str
The Dec of the target (in degrees) hh:mm:ss.s or dd:mm:ss.s or
representing a float
instrument : str
Name of the instrument. Can either be (case-sensitive):
'NIRISS', 'NIRCam', 'MIRI', 'FGS', or 'NIRSpec'
ephFileName : str
The filename of the ephemeris file.
output : str
Switches on plotting with Bokeh. Parameter value must be
'bokeh'.
Returns
-------
paGood : float
The good position angle.
paBad : float
The bad position angle.
gd : matplotlib.dates object
The gregorian date.
fig : bokeh.plotting.figure object
The plotted figure.
"""
# Getting calculations from GTVT (General Target Visibility Tool)
# blockPrint()
tab = get_table(ra, dec)
# enablePrint()
gd = tab['Date']
paMin = tab[str(instrument) + ' min']
paMax = tab[str(instrument) + ' max']
paNom = tab[str(instrument) + ' nom']
v3min = tab['V3PA min']
v3max = tab['V3PA max']
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~NOTE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Addressing NIRSpec issue*
# *the issue that NIRSpec's angle goes beyond 360 degrees with some targs,
# thus resetting back to 0 degrees, which can make the plot look weird
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
index = np.arange(0, len(paNom), 1)
for idx in index:
try:
a1 = paNom[idx]
b1 = paNom[idx + 1]
if (np.isfinite(a1)) & (np.isfinite(b1)):
delta = np.abs(a1 - b1)
if delta > 250:
gd = np.insert(gd, idx + 1, np.nan)
paMin = np.insert(paMin, idx + 1, np.nan)
paMax = np.insert(paMax, idx + 1, np.nan)
paNom = np.insert(paNom, idx + 1, np.nan)
v3min = np.insert(v3min, idx + 1, np.nan)
v3max = np.insert(v3min, idx + 1, np.nan)
except BaseException:
pass
# Setting up HoverTool parameters & other variables
COLOR = 'green'
TOOLS = 'pan, wheel_zoom, box_zoom, reset, save'
SOURCE = ColumnDataSource(data=dict(pamin=paMin,
panom=paNom,
pamax=paMax,
date=gd))
# Time to plot
if output == 'bokeh':
fig = figure(tools=TOOLS,
width=800,
height=400,
x_axis_type='datetime',
title='{} Visibility with {}'.format(targetName,
instrument))
# Draw the curve and PA min/max circles
try:
nom = fig.line('date', 'panom',
line_color=COLOR,
legend_label='Nominal Aperture PA',
alpha=.5,
source=SOURCE)
except AttributeError:
nom = fig.line('date', 'panom',
line_color=COLOR,
legend='Nominal Aperture PA',
alpha=.5,
source=SOURCE)
fig.circle('date', 'pamin', color=COLOR, size=1, source=SOURCE)
fig.circle('date', 'pamax', color=COLOR, size=1, source=SOURCE)
# Adding HoverTool
fig.add_tools(HoverTool(tooltips=[('Date', '@date{%F}'),
('Maximum Aperture PA', '@pamax'),
('Nominal Aperture PA', '@panom'),
('Minimum Aperture PA', '@pamin')],
formatters={'@date': 'datetime'}))
# Plot formatting
fig.xaxis.axis_label = 'Date'
fig.yaxis.axis_label = 'Aperture Position Angle (degrees)'
fig.y_range = ranges.Range1d(0, 360)
# Making the output table
# Creating new lists w/o the NaN values
v3minnan, v3maxnan, paNomnan, paMinnan, paMaxnan, gdnan = [], [], [], [], [], []
for vmin, vmax, pnom, pmin, pmax, date in zip(
v3min, v3max, paNom, paMin, paMax, gd):
if np.isfinite(pmin):
v3minnan.append(vmin)
v3maxnan.append(vmax)
paNomnan.append(pnom)
paMinnan.append(pmin)
paMaxnan.append(pmax)
gdnan.append(date)
# Adding MJD column
mjdnan = []
for date in gdnan:
t = Time(str(date), format='iso')
mjd = t.mjd
mjdnan.append(mjd)
# Adding lists to a table object
table = Table([v3minnan,
v3maxnan,
paMinnan,
paMaxnan,
paNomnan,
gdnan,
mjdnan],
names=('min_V3_PA',
'max_V3_PA',
'min_Aperture_PA',
'max_Aperture_PA',
'nom_Aperture_PA',
'Gregorian',
'MJD'))
# Getting bad PAs
allPAs = np.arange(0, 360, 1)
badPAs = []
for pa in allPAs:
if (pa not in np.round(paMinnan)) & \
(pa not in np.round(paMaxnan)) & \
(pa not in np.round(paNomnan)):
badPAs.append(pa)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~NOTE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# This addresses a bokeh shading issue that accidentally shades
# accessible PAs (e.g: trappist-1b)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
badPAs = select_badPAs_ge_paNomnan(badPAs, paNomnan)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~NOTE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Grouping the bad PAs into lists within the badPAs list.
# This will make bad PA shading easier in the contamination Bokeh plot
# (sossContamFig.py)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
badPAs = np.sort(badPAs)
if len(badPAs > 0):
grouped_badPAs = [[badPAs[0]]]
for idx in range(1, len(badPAs)):
if ((badPAs[idx - 1] + 1) == badPAs[idx]):
grouped_badPAs[len(grouped_badPAs) - 1].append(badPAs[idx])
elif ((badPAs[idx - 1] + 1) < badPAs[idx]):
grouped_badPAs.append([badPAs[idx]])
# grouped_badPAs = np.asarray(grouped_badPAs)
else: # Accounting for targets with 100% visibility
grouped_badPAs = [] #np.asarray([])
return paMin, paMax, gd, fig, table, grouped_badPAs
[docs]
def select_badPAs_ge_paNomnan(badPAs, paNomnan, threshold=7):
"""Returns the absolute difference between each badPAs and paNomnan
Should be greater than threshold (default=7)
Parameters
----------
badPAs: list
The list of bad position angles
paNomnan: list
The list of nominal PAs
Returns
-------
np.ndarray
The array of PAs
"""
# Reshaping
badPAs_array = np.array(badPAs)[np.newaxis] # (1, len(badPAs))
paNomnan_array = np.array(paNomnan)[np.newaxis].T # (len(paNomnan), 1)
# elementwise absolute difference
diff = np.abs(np.subtract(badPAs_array, paNomnan_array)) # (len(paNomnan), len(badPAs))
# boolean array above threshold
above_thresh = np.all(diff >= threshold, axis=0) # (len(badPAs),)
# index and return those that are above threshold
return badPAs_array[0, above_thresh]