#################################################################################
from pathlib import Path
import os
import sys
from astropy.stats import mad_std
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename
import glob, os
from astropy import stats
import matplotlib.pyplot as plt
import numpy as np
from datetime import date
today = date.today()
from astropy.convolution import convolve, Box1DKernel
from astropy.convolution import Box2DKernel
date = today.strftime("%b-%d-%Y")
from astropy.wcs import WCS
import warnings
warnings.filterwarnings('ignore')
from matplotlib.backends.backend_pdf import PdfPages
import shutil
#################################################################################

# Read in day number from command line:  

n = len(sys.argv)
if n == 0 or n == 1 or n > 3:
    print('Usage:  python process_movie.py ### $$$' )
    print('   where ### is a 3-digit identifier for the day number')
    print('   and $$$ is an (optional) 3-digit identifier for the observer code')
    # end program:
    sys.exit(0)

if n == 2:  
    daynum = sys.argv[1]
    obscode = '*'
if n == 3:  
    daynum = sys.argv[1]
    obscode = sys.argv[2]

#################################################################################
# Use glob to create list:
current_dir = os.getcwd()
if obscode == '*':
    image_dir = '/home/cannon/mount-test/'
if obscode != '*':
    image_dir = '/home/cannon/mount-test/'
out_dir = current_dir+'/'
os.makedirs(out_dir, exist_ok = True)
image_list = sorted(glob.glob(image_dir+'*.fits',recursive=True))
N_images = len(image_list)
print('-----------------------------------------------------------')
if obscode == '*':
    print('    --> '+str(N_images)+' images found for day '+str(daynum)+' from all observers')
if obscode != '*':
    print('    --> '+str(N_images)+' images found for day '+str(daynum)+' from observer '+str(obscode))
print('-----------------------------------------------------------')
#################################################################################

for ii in range(100, 101):
    hd = fits.getheader(image_list[ii])
    ha = hd['HA']
    print('--------------------')
    print(ha)
#    if str.split(ha,':')[0][0] == '-':
#        hahours = float(str.split(ha,':')[0]) + float(str.split(ha,':')[1])/60. + float(str.split(ha,':')[2])/3600.
#        hahours=-hahours
#        print('Hour angle = '+str(hahours))
#    if str.split(ha,':')[0][0] != '-':
#        hahours = float(str.split(ha,':')[0]) + float(str.split(ha,':')[1])/60. + float(str.split(ha,':')[2])/3600.
#        hahours = hahours
#        print('Hour angle = '+str(hahours))

    if str.split(ha,':')[0][0] == '-':
        hahours = -float(str.split(ha,':')[0]) + float(str.split(ha,':')[1])/60. + float(str.split(ha,':')[2])/3600.                                         
        print(ha, (str.split(ha,':')[0]), str.split(ha,':')[1], str.split(ha,':')[2])
        hahoursval=-hahours
        print('Hour angle = '+str(hahoursval))                                                  
    if str.split(ha,':')[0][0] != '-':
        hahours = float(str.split(ha,':')[0]) + float(str.split(ha,':')[1])/60. + float(str.split(ha,':')[2])/3600.
        hahoursval = hahours
#        print('Hour angle = '+str(hahoursval))


print()



#################################################################################
# Plot temp vs. time:
templist = []
timelist = []
amblist=[]
sinklist=[]
explist = []
ratio=[]
wndlist=[]
azimlist=[]
elevlist=[]
halist=[]
#for ii in range(100,110):
for ii in range(0, N_images):
    #im = fits.getdata(image_list[ii])                                                                                                                      
    hd = fits.getheader(image_list[ii])
    elev = hd['ELEVATIO']
    azim = hd['AZIMUTH']
    ha = hd['HA']
    azimdeg = float(str.split(azim,':')[0]) + float(str.split(azim,':')[1])/60. + float(str.split(azim,':')[2])/3600.
    elevdeg = float(str.split(elev,':')[0]) + float(str.split(elev,':')[1])/60. + float(str.split(elev,':')[2])/3600.
    #hahours = float(str.split(ha,':')[0]) + float(str.split(ha,':')[1])/60. + float(str.split(ha,':')[2])/3600.
    if str.split(ha,':')[0][0] == '-' or str.split(ha,':')[0][1] == '-':
        hahours = -float(str.split(ha,':')[0]) + float(str.split(ha,':')[1])/60. + float(str.split(ha,':')[2])/3600.
        hahoursval=-hahours
        #print('Hour angle = '+str(hahours))
    if str.split(ha,':')[0][0] != '-':
        hahours = float(str.split(ha,':')[0]) + float(str.split(ha,':')[1])/60. + float(str.split(ha,':')[2])/3600.
        hahoursval = hahours
        #print('Hour angle = '+str(hahours))
    #if str.split(ha,':')[0][0] == '-':
    #    hahours = float(str.split(ha,':')[0]) + float(str.split(ha,':')[1])/60. + float(str.split(ha,':')[2])/3600.
    #    #print('Hour angle = '+str(hahours))
    #if str.split(ha,':')[0][0] != '-':
    #    hahours = -(float(str.split(ha,':')[0]) + float(str.split(ha,':')[1])/60. + float(str.split(ha,':')[2])/3600.)
        #print('Hour angle = '+str(hahours))
    #if np.abs(float(str.split(ha,':')[0])) > 0.:
    #    hahours = float(str.split(ha,':')[0]) + float(str.split(ha,':')[1])/60. + float(str.split(ha,':')[2])/3600.
    #if np.abs(float(str.split(ha,':')[0])) < 0.:
    #    hahours = -float(str.split(ha,':')[0]) - float(str.split(ha,':')[1])/60. - float(str.split(ha,':')[2])/3600.
    exptime = hd['EXPTIME']
    logexptime= np.log(exptime)/np.log(2)
    DATE=str.split(hd['DATE-OBS'],'T')
    day = DATE[0]
    hour    = float(str.split(DATE[1],':')[0])
    minute  = float(str.split(DATE[1],':')[1])
    second  = float(str.split(DATE[1],':')[2])    
    time=hour + (minute/60.) + (second/3600.)
    temp=hd['CCD-TEMP']
    hsink=hd['HSINKT']
    ambient=hd['WXTEMP'] 
    windspd=hd['WXWNDSPD']
    timelist.append(time)
    templist.append(temp)
    explist.append(logexptime)
    amblist.append(ambient)
    sinklist.append(hsink)
    ratio.append(ambient/hsink)
    wndlist.append(windspd)
    elevlist.append(elevdeg)
    azimlist.append(azimdeg)
    halist.append(hahoursval)
    #print(hahoursval)

fig = plt.figure(figsize=(9,12))
fig.subplots_adjust(hspace=0)
ax1 = plt.subplot(411)
plt.title("RLMT "+str(DATE[0]),fontsize = 16)
plt.scatter(timelist, halist, s=12, color='navy', label=r'Hour Angle (hours)')
plt.grid(color = 'grey', linestyle = '--', linewidth = 0.5)
plt.tick_params('x', labelbottom=False)
plt.legend(facecolor='gainsboro',fontsize = 12)
plt.ylabel(r"Temp ($^{\rm o}$C)",fontsize = 16)

ax2 = plt.subplot(412, sharex=ax1)
plt.grid(color = 'grey', linestyle = '--', linewidth = 0.5)
plt.scatter(timelist, azimlist, s=12, color='saddlebrown', label=r'Telescope Azimuth ($^{\rm o}$)')
plt.legend(facecolor='gainsboro',fontsize = 12)
plt.ylabel(r"Azimuth ($^{\rm o}$)",fontsize = 16)

ax3 = plt.subplot(413, sharex=ax1)
plt.scatter(timelist, elevlist, s=12, color='darkcyan', label=r'Source Elevation ($^{\rm o}$)')
plt.tick_params('x')
plt.tick_params('y')
plt.grid(color = 'grey', linestyle = '--', linewidth = 0.5)
plt.legend(facecolor='gainsboro',fontsize = 12)
plt.xlabel("Time (UT hours)",fontsize = 16)
plt.ylabel(r"Elevation ($^{\rm o}$)",fontsize = 16)

ax4 = plt.subplot(414, sharex=ax1)
plt.scatter(timelist, wndlist, s=12, color='darkcyan', label=r'Wind Speed (kph)')
plt.tick_params('x')
plt.tick_params('y')
plt.grid(color = 'grey', linestyle = '--', linewidth = 0.5)
plt.legend(facecolor='gainsboro',fontsize = 12)
plt.xlabel("Time (UT hours)",fontsize = 16)
plt.ylabel(r"Wind Speed (kph)",fontsize = 16)
if obscode == '*':
    plt.savefig(str(daynum)+'.plots.v2.png',dpi=100, bbox_inches = 'tight', transparent=False, facecolor='whitesmoke')
if obscode != '*':
    plt.savefig(str(daynum)+'.'+str(obscode)+'.plots.v2.png',dpi=100, bbox_inches = 'tight', transparent=False, facecolor='whitesmoke')
#################################################################################