#!/usr/bin/env python
# coding: utf-8

# # Observing Run Preparation Module
# 
# **Lecturer:** Robert Quimby<br>
# **Jupyter Notebook Author:** Shubham Srivastav, Cameron Hummels & Robert Quimby
# 
# This is a Jupyter notebook lesson taken from the GROWTH Summer School 2019.  For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-astro-school-2018-resources.html
# 
# ## Objective
# Demonstrate how to plan observations prior to an observing run.
# 
# ## Key steps
# - Select targets
# - Get visibility and airmass plots
# - Get moon separation angles
# - Calculate exposure times for targets
# 
# ## Required dependencies
# 
# See GROWTH school webpage for detailed instructions on how to install these modules and packages.  Nominally, you should be able to install the python modules with `pip install <module>`.  The external astromatic packages are easiest installed using package managers (e.g., `rpm`, `apt-get`).
# 
# ### Python modules
# * python 3
# * astropy
# * numpy
# * matplotlib
# * astroplan
# * pytz
# 
# ### External packages
# None

# In[1]:


import numpy as np
from astropy import units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord
from astropy.coordinates import EarthLocation
import pytz
get_ipython().run_line_magic('matplotlib', 'inline')
from astroplan import Observer, FixedTarget
from astropy.utils.iers import conf
conf.auto_max_age = None
from astroplan import download_IERS_A 
from astropy.coordinates import get_sun, get_moon, get_body
from astroplan import moon_illumination


# ### Date and Time
# - Dates and times are in UTC
# - Default time is 00:00:00 UTC (verify this)

# In[2]:


date = Time("2018-12-03", format='iso')
print(date)


# ### What is the current UTC?

# In[3]:


now = Time.now()
print(now)
print(now.jd)
print(now.mjd)
print(now.decimalyear)


# ### Exercise
# What time will it be (in UTC) after 1 hour 45 minutes from `now`? Complete the line below to print it out.

# In[4]:


print("In 1 hour and 45 minutes, the time will be {0} UTC".format(now + 1*u.h + 45*u.min))


# ### Using UT1
# - To keep accurate time, the changes in earth's rotation period have to be taken into account.
# - AstroPy does this using a convention called UT1, that is tied to the rotation of earth with respect to the position of distant quasars. IERS - International Earth Rotation and Reference Systems Service keeps continuous tabs on the orientation of the earth and updates the data in the IERS bulletin.
# Update the bulletin:

# In[5]:


download_IERS_A()


# ### Check to see what observatories are available in the database.

# In[6]:


print("Available observatories: \n{0}"
      .format(', '.join(EarthLocation.get_site_names())))


# ### Setting up observatory location

# In[7]:


# Mount Laguna Observatory is not listed in the database, so let's define the location
latitude = 32.842167 * u.deg
longitude = -116.426938 * u.deg
elevation = 1860 * u.m
location = EarthLocation.from_geodetic(longitude, latitude, elevation)
mlo = Observer(location = location, timezone = 'America/Los_Angeles',
                             name = "MLO", description = "MLO 1.0-m telescope")
mlo


# ### Sunset, Sunrise, Midnight

# In[8]:


# ##### just for testing...REMOVE! #####
now = Time('2019-08-04 21:17:59')


# In[9]:


# Calculating the sunset, midnight and sunrise times for our observatory 
# What is astronomical twilight?
sunset_mlo = mlo.sun_set_time(now, which='nearest')
eve_twil_mlo = mlo.twilight_evening_astronomical(now, which='nearest')
midnight_mlo = mlo.midnight(now, which='next')
morn_twil_mlo = mlo.twilight_morning_astronomical(now, which='next')
sunrise_mlo = mlo.sun_rise_time(now, which='next')

print("Sunset at MLO will be at {0.iso} UTC".format(sunset_mlo))
print("Astronomical evening twilight at MLO will be at {0.iso} UTC".format(eve_twil_mlo))
print("Midnight at MLO will be at {0.iso} UTC".format(midnight_mlo))
print("Astronomical morning twilight at MLO will be at {0.iso} UTC".format(morn_twil_mlo))
print("Sunrise at MLO will be at {0.iso} UTC".format(sunrise_mlo))


# ### Exercise 
# Find the effective length of time (in hours) available for optical astronomical observations at MLO tonight

# In[10]:


observing_time = (morn_twil_mlo - eve_twil_mlo).to(u.hour)
print("You can observe for {0:.1f} at MLO tonight".format(observing_time))


# ### Local Sidereal Time (LST)

# In[11]:


#What is the LST now at MLO?
#What would the LST be at MLO at local midnight?
lst_now = mlo.local_sidereal_time(now)
lst_mid = mlo.local_sidereal_time(midnight_mlo)
print("LST at MLO now is {0:.2f}".format(lst_now))
print("LST at MLO at local midnight will be {0:.2f}".format(lst_mid))


# ### Choosing targets for observations
# Targets can be defined by name or coordinates.
# 

# In[12]:


# using coordinates
coords = SkyCoord('18h53m35.097s +33d01m44.8831s', frame='icrs') # coordinates of the Ring Nebula (M57) 
m57 = FixedTarget(name = 'M57', coord=coords)
m57.ra.hms


# In[13]:


# by name
target = FixedTarget.from_name('m57') # Messier 57
target.coord


# Check to see if target is "up" at evening twilight (assume "up" means more than 30 degrees above the horizon). Also check if target is available at midnight and morning twilight.

# In[14]:


# check if the target is up
print(mlo.target_is_up(eve_twil_mlo, m57, horizon=30*u.deg))
print(mlo.target_is_up(midnight_mlo, m57, horizon=30*u.deg))
print(mlo.target_is_up(morn_twil_mlo, m57, horizon=30*u.deg))


# In[15]:


# Altitude and Azimuth of target at evening twilight
aa = mlo.altaz(eve_twil_mlo, m57)
aa.alt.degree, aa.az.degree


# Determine the time at which the target rises

# In[16]:


m57rise = mlo.target_rise_time(now, m57, which = 'next', horizon=0*u.deg)
print(m57rise.iso)  #default format is JD


# ### Dealing with moving targets

# In[17]:


get_body('jupiter', now)


# In[18]:


# get moon position at midnight 
get_moon(midnight_mlo)


# In[19]:


# How bright is the moon at midnight?
moon_illumination(midnight_mlo)


# In[20]:


# We can turn solar system objects into 'pseudo-fixed' targets to plan observations
saturn_midnight = FixedTarget(name = 'Saturn', coord = get_body('saturn', midnight_mlo))
saturn_midnight.coord


# ### Airmass
# - Ideally, targets should be observed when they have the least airmass. Airmass ranges from 1 (zenith) to ~38 at the horizon.
# - Airmass is 2.0 at alt=30, 2.9 at alt=20 and 3.9 at alt=15 degrees
# - As a general rule of thumb, try observing targets when airmass < 2
# - Let us find the airmass of M57 at midnight at MLO

# In[21]:


#Is the target up at MLO at midnight?
mlo.target_is_up(midnight_mlo, target)


# In[22]:


#lets check the alt and az of the target at midnight
target_altaz = mlo.altaz(midnight_mlo, target)
target_altaz.altaz


# That is high over head; ideal for observing.

# In[23]:


#Find the airmass
target_altaz.secz


# ### Plots to help planning

# Now we can visualize what we have done so far using some plots

# In[24]:


import matplotlib.pyplot as plt
from astroplan.plots import plot_sky, plot_airmass


# In[25]:


#position of target at midnight
plot_sky(target, mlo, midnight_mlo);


# Now let us see how the target moves over the course of the night

# In[26]:


t_start = eve_twil_mlo
t_end = morn_twil_mlo
t_observe = t_start + (t_end - t_start) * np.linspace(0.0, 1.0, 20)
plot_sky(target, mlo, t_observe);


# Now let's plot the airmass as a function of time

# In[27]:


plot_airmass(target, mlo, t_observe)
plt.grid();


# The airmass is above 2 for the better part of the night, making M57 a good summer target from MLO.
# Note that the default airmass limit is 3 in astroplan, corresponding to ~19 degrees elevation.

# ### Finder Charts

# In[28]:


from astroplan.plots import plot_finder_image
from astroquery.skyview import SkyView


# Load an image of the field in which the target lies.

# In[29]:


# field of view corresponding to the MLO 1.0-m telesocpe
fov = 14*u.arcmin

# plot the image
plot_finder_image(target, fov_radius=fov);  


# Now let's define an array of targets to work with

# In[30]:


target_names = ['vega', 'polaris', 'm1', 'm42', 'm55']
targets = [FixedTarget.from_name(target) for target in target_names]
targets


# Which of these targets are up now?

# In[31]:


mlo.target_is_up(now, targets)


# Which of these targets are up at local midnight?

# In[32]:


mlo.target_is_up(midnight_mlo, targets)


# ### Exercise
# Find out the times at which the targets rise to an elevation of 10 degrees. Use target_rise_time. 

# In[33]:


for target in targets:
    print(mlo.target_rise_time(now, target, which = 'next', horizon = 10*u.deg).iso)    


# How high is Vega above the horizion now?

# In[34]:


mlo.altaz(now, targets[0])


# Now let's plot the elevation of Vega to see how it varies over the night

# In[35]:


times = (t_start - 0.5 * u.h) + (t_end - t_start + 1 * u.h) * np.linspace(0.0, 1.0, 40)
elevations = mlo.altaz(times, targets[0]).alt
ax = plt.gca()
ax.plot_date(times.plot_date, elevations.deg)
ax.set(xlabel = 'Time UTC [MM-DD HH]', ylabel = 'Altitude [deg]')
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.grid()


# ### Exercise
# Plot the altitude as a function of time for tonight for each of the targets in a single plot

# In[36]:


ax = plt.gca()
times = t_start + (t_end - t_start) * np.linspace(0.0, 1.0, 20)
for target in targets: 
    elevation = mlo.altaz(times, target).alt
    ax.plot_date(times.plot_date, elevation, label=target.name)
    
ax.set(xlabel = 'Time UTC [MM-DD HH]', ylabel = 'Altitude [deg]')
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.legend()
plt.grid()


# ### Exercise
# Plot sky positions for each target using plot_sky for tonight at MLO in a single plot.

# In[37]:


times = (t_start - 0.5 * u.h) + (t_end - t_start + 1 * u.h) * np.linspace(0.0, 1.0, 20)
for target in targets:
    plot_sky(target, mlo, times)
plt.legend(loc=[1.1,0]);


# ### Exercise
# Plot airmass vs time for each target in targets for tonight at MLO.

# In[38]:


for target in targets:
    plot_airmass(target, mlo, t_observe)
plt.ylim(4,0.5)    
plt.legend()
plt.grid()


# ### Observational Constraints
# You can set specific constraints that define when a target is "observable"
# - twilight level (e.g., "Civil")
# - airmass
# - altitude limits

# In[39]:


from astroplan import (AltitudeConstraint, AirmassConstraint,
                       AtNightConstraint, MoonSeparationConstraint)
constraints = [AltitudeConstraint(15*u.deg, 84*u.deg),
               AirmassConstraint(3), AtNightConstraint.twilight_civil(), MoonSeparationConstraint(min = 10 * u.deg)]
t_range = Time([t_start - 0.5 * u.hour, t_end + 0.5 * u.hour])


# In[40]:


from astroplan import is_observable, is_always_observable, months_observable

# Are targets ever observable in the time range?
ever_observable = is_observable(constraints, mlo, targets, time_range=t_range)
print(ever_observable)

# Are targets always observable in the time range?
always_observable = is_always_observable(constraints, mlo, targets, time_range=t_range)
print(always_observable)


# The functions is_observable and ever_observable return boolean arrays. Let's print their output in tabular form.

# In[41]:


from astropy.table import Table
observability_table = Table()
observability_table['targets'] = [target.name for target in targets]
observability_table['ever_observable'] = ever_observable
observability_table['always_observable'] = always_observable
print(observability_table)


# Or we could do this directly using the observability_table function 

# In[42]:


from astroplan import observability_table
table = observability_table(constraints, mlo, targets, time_range = t_range)
print(table)


# In[43]:


# During what months are the targets ever observable?
months_observable(constraints, mlo, targets)


# ### Exercise
# - Create a list of your favourite targets and store it in a text file with 3 columns - name, RA and Dec. Or you could    use 'targetlists.txt' which already contains a list of targets. 
# - Read the text file, and store the targets as FixedTarget objects.
# - Get observability tables for all the targets for different moon separation angles (10, 20, 30... degrees)
# - Plot airmass and sky position as a function of time for tonight for all your targets.

# In[44]:


from astropy.io import ascii
table = ascii.read('data/targetlist.txt')
targets = [FixedTarget(coord=SkyCoord(ra=ra*u.deg, dec=dec*u.deg), name=name) for name, ra, dec in table]
#The recipe for the remaining part of the exercise is in the previous solved exercises