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

# # Observing Run Preparation Module
# 
# **Lecturer:** Robert Quimby<br>
# **Live Session:** Robert Quimby, Varun Bhalerao<br>
# **Jupyter Notebook Author:** Shubham Srivastav, Cameron Hummels & Robert Quimby
# 
# This is a Jupyter notebook lesson taken from the GROWTH Summer School 2020.  For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-astro-school-2020-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[ ]:


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[ ]:


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


# ### What is the current UTC?

# In[ ]:


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[ ]:





# ### 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[ ]:


download_IERS_A()


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

# In[ ]:


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


# ### Setting up observatory location

# In[ ]:


# Let us manually add the location for the Liverpool Telescope
latitude = 28.76234 * u.deg
longitude = -17.87925 * u.deg
elevation = 2363 * u.m
location = EarthLocation.from_geodetic(longitude, latitude, elevation)

# You can optionally specify a time zone as well. 
# A list of timezones can be found at https://stackoverflow.com/questions/13866926/is-there-a-list-of-pytz-timezones
liverpool = Observer(location = location, timezone = 'Europe/Madrid',
                             name = "LT", description = "Liverpool Telescope")
liverpool


# In[ ]:


# If you know the name, you can set a location accordingly
location = EarthLocation.of_site("Roque de los Muchachos")


# 
# ### Sunset, Sunrise, Midnight

# In[ ]:


# Calculating the sunset, midnight and sunrise times for our observatory 
# What is astronomical twilight?
sunset_liverpool = liverpool.sun_set_time(now, which='nearest')
eve_twil_liverpool = liverpool.twilight_evening_astronomical(now, which='nearest')
midnight_liverpool = liverpool.midnight(now, which='next')
morn_twil_liverpool = liverpool.twilight_morning_astronomical(now, which='next')
sunrise_liverpool = liverpool.sun_rise_time(now, which='next')

print("Sunset at Liverpool will be at {0.iso} UTC".format(sunset_liverpool))
print("Astronomical evening twilight at Liverpool will be at {0.iso} UTC".format(eve_twil_liverpool))
print("Midnight at Liverpool will be at {0.iso} UTC".format(midnight_liverpool))
print("Astronomical morning twilight at Liverpool will be at {0.iso} UTC".format(morn_twil_liverpool))
print("Sunrise at Liverpool will be at {0.iso} UTC".format(sunrise_liverpool))


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

# In[ ]:





# ### Local Sidereal Time (LST)

# In[ ]:


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


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

# In[ ]:


# 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[ ]:


# 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[ ]:


# check if the target is up
print(liverpool.target_is_up(eve_twil_liverpool, m57, horizon=30*u.deg))
print(liverpool.target_is_up(midnight_liverpool, m57, horizon=30*u.deg))
print(liverpool.target_is_up(morn_twil_liverpool, m57, horizon=30*u.deg))


# In[ ]:


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


# Determine the time at which the target rises

# In[ ]:


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


# ### Dealing with moving targets

# In[ ]:


get_body('jupiter', now)


# In[ ]:


# get moon position at midnight 
get_moon(midnight_liverpool)


# In[ ]:


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


# In[ ]:


# We can turn solar system objects into 'pseudo-fixed' targets to plan observations
saturn_midnight = FixedTarget(name = 'Saturn', coord = get_body('saturn', midnight_liverpool))
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 Liverpool Telescope

# In[ ]:


#Is the target up at Liverpool Telescope at midnight?
liverpool.target_is_up(midnight_liverpool, target)


# In[ ]:


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


# That is about midway up from the horizon, not too bad.

# In[ ]:


#Find the airmass
target_altaz.secz


# ### Plots to help planning

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

# In[ ]:


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


# In[ ]:


#position of target at midnight
plot_sky(target, liverpool, midnight_liverpool);


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

# In[ ]:


t_start = eve_twil_liverpool
t_end = morn_twil_liverpool
t_observe = t_start + (t_end - t_start) * np.linspace(0.0, 1.0, 20)
plot_sky(target, liverpool, t_observe);


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

# In[ ]:


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


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

# ### Finder Charts

# In[ ]:


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


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

# In[ ]:


# field of view corresponding to the IO:O instrument on the Liverpool Telescope
fov = 10*u.arcmin

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


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

# In[ ]:


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[ ]:


liverpool.target_is_up(now, targets)


# Which of these targets are up at local midnight?

# In[ ]:


liverpool.target_is_up(midnight_liverpool, targets)


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

# In[ ]:





# How high is Vega above the horizion now?

# In[ ]:


liverpool.altaz(now, targets[0])


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

# In[ ]:


times = (t_start - 0.5 * u.h) + (t_end - t_start + 1 * u.h) * np.linspace(0.0, 1.0, 40)
elevations = liverpool.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[ ]:





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

# In[ ]:





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

# In[ ]:





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

# In[ ]:


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[ ]:


from astroplan import is_observable, is_always_observable, months_observable

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

# Are targets always observable in the time range?
always_observable = is_always_observable(constraints, liverpool, 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[ ]:


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[ ]:


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


# In[ ]:


# During what months are the targets ever observable?
months_observable(constraints, liverpool, 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[ ]: