"""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
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
from ..pkgdata import resource_filename
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 = 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]