# hack to get the path right
import sys
sys.path.append('..')

import ztf_sim
from astropy.time import Time
import pandas as pd
import numpy as np
import astropy.units as u
import pylab as plt

"""
Explanation: ztf_sim_introduction
This notebook illustrates basic use of the ztf_sim modules.
End of explanation
"""


ztf_sim.fields.generate_test_field_grid()

"""
Explanation: First we'll generate a test field grid.  You only need to do this the first time you run the simulator.
End of explanation
"""


f = ztf_sim.fields.Fields()

"""
Explanation: Let's load the Fields object with the default field grid.  Fields is a thin wrapper around a pandas DataFrame containing the field information.
End of explanation
"""


f.fields.head()

"""
Explanation: The raw fieldid and coordinates are stored as a pandas Dataframe in the .fields attribute:
End of explanation
"""


f.alt_az(Time.now()).head()

"""
Explanation: Now let's calculate their altitude and azimuth at a specific time using the astropy.time.Time object:
End of explanation
"""


f.fields.loc[853]

"""
Explanation: Demonstrating low-level access to fields by the fieldid index (usually not required):
End of explanation
"""


f.fields['dec'] > -30.

"""
Explanation: We can select fields with conditionals:
End of explanation
"""


cuts = f.select_fields(dec_range=[0,10],gridid=0,ecliptic_lat_range=[-5,5])
cuts.head()

"""
Explanation: It's easier to use the select_fields convenience function, though.  It returns a boolean Series indexed by fieldid that we can use to do calculations on subsets of the field grid.
End of explanation
"""


f.alt_az(Time.now(),cuts=cuts)

"""
Explanation: Calculate the current altitude and azimuth of the selected fields:
End of explanation
"""


f.overhead_time(853,Time.now())

f = ztf_sim.fields.Fields()
Exposure_time = 60*u.second
Night_length=9*u.h


time0 = Time('2015-09-10 20:00:00') + 7*u.h
time = time0
f.fields = f.fields.join(pd.DataFrame(np.zeros(len(f.fields)),columns=['observed']))
f.fields = f.fields.join(pd.DataFrame(np.zeros(len(f.fields)),columns=['possibleToObserve']))

def observe(f, nightStart):
    time=nightStart
    goodAltitude = f.alt_az(time)['alt'] > 20
    shouldObserve = f.fields['observed'] == 0
    good = goodAltitude & shouldObserve #& f.alt_az(time+1*u.h)['alt'] < 20 # start with a field which won't be observable later
    
    if np.all(good) is False:
        good = goodAltitude & shouldObserve
        
    fid = f.fields[good].iloc[0].name
    f.fields['observed'][fid]+=1
    f.fields['possibleToObserve'][goodAltitude] = 1
    time += Exposure_time

    while time < nightStart + Night_length:
        goodAltitude = f.alt_az(time)['alt'] > 20
        shouldObserve = f.fields['observed'] == 0
        good = goodAltitude & shouldObserve
        f.fields['possibleToObserve'][goodAltitude] = 1
        
        if not np.any(good):
            time += 60*u.s
            continue
            
        slewTime = f.overhead_time(fid,time)[good]
        fid = int(slewTime.idxmin())
        # print slewTime['overhead_time'][fid]
        time += Exposure_time + slewTime['overhead_time'][fid]*u.second
        f.fields['observed'][fid]+=1
        # print time-7*u.h
        

# First night
observe(f,time)
fieldsPossible = np.sum(f.fields['possibleToObserve'])
print fieldsPossible
fieldsObserved = np.sum(f.fields['observed'])
print fieldsObserved
meanTime = (Night_length.to(u.s)-fieldsObserved*Exposure_time)/(fieldsObserved-1)
print meanTime

# Second night
time=time0+24*u.h
observe(f,time)
fieldsPossible = np.sum(f.fields['possibleToObserve'])
print fieldsPossible
fieldsObserved = np.sum(f.fields['observed'])
print fieldsObserved
meanTime = (2*Night_length.to(u.s)-fieldsObserved*Exposure_time)/(fieldsObserved-1)
print meanTime

for dec in np.append(np.linspace(-90,90,10),0):
    ra=np.linspace(0, 360,1000)
    x,y = raDec2xy(ra,dec)
    plt.plot(x,y,'k')
    
for ra in np.linspace(0,360,10):
    dec=np.linspace(-90, 90,1000)
    x,y = raDec2xy(ra,dec)
    plt.plot(x,y,'k')
    
x,y = raDec2xy(f.fields['ra'],f.fields['dec'])
plt.plot(x,y,'o',color=(.8,.8,.8))    
plt.show()

def raDec2xy(ra,dec):
    # Using Aitoff projections (from Wiki) returns x-y coordinates on a plane of RA and Dec
    theta = np.deg2rad(dec)
    phi = np.deg2rad(ra)-np.pi #the range is [-pi,pi]
    alpha=np.arccos(np.cos(theta)*np.cos(phi/2))
    x=2*np.cos(theta)*np.sin(phi/2)/np.sinc(alpha/np.pi) # The python's sinc is normalazid, hence the /pi
    y=np.sin(theta)/np.sinc(alpha/np.pi)
    return x,y

"""
Explanation: Calculating the overhead time (max of ha, dec, dome slews and readout time):
End of explanation
"""


