#################################################################################
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-test2/'
if obscode != '*':
    image_dir = '/home/cannon/mount-test2/'
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('-----------------------------------------------------------')
#################################################################################

#################################################################################
# Create two-page PDF with plots vs. time:

templist = []
timelist = []
amblist=[]
sinklist=[]
explist = []
ratio=[]
wndlist=[]
azimlist=[]
elevlist=[]
halist=[]

for ii in range(0, N_images):
    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.
    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
    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
    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)

with PdfPages(out_dir+str(daynum)+'.plots.wcs.pdf') as pdf:
    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, templist, s=12, color='navy', label=r'CCD Temp ($^{\rm o}$C)')
    plt.axhline(y = -20.0, color = 'r', linestyle = '-', linewidth = 1)
    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, sinklist, s=12, color='saddlebrown', label=r'Heatsink Temp ($^{\rm o}$C)')
    plt.legend(facecolor='gainsboro',fontsize = 12)
    plt.ylabel(r"Temp ($^{\rm o}$C)",fontsize = 16)
        
    ax3 = plt.subplot(413, sharex=ax1)
    plt.scatter(timelist, amblist, s=12, color='darkcyan', label=r'Ambient Temp ($^{\rm o}$C)')
    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"Temp ($^{\rm o}$C)",fontsize = 16)

    ax4 = plt.subplot(414, sharex=ax1)
    plt.scatter(timelist, wndlist, s=12, color='tomato', 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)

    pdf.savefig(dpi=100)
    plt.close()

    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"Hour Angle (hours)",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='tomato', 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)
    
    pdf.savefig(dpi=100)
    plt.close()
#################################################################################


#################################################################################
# Sort the list by time: 
astlist = [x for _,x in sorted(zip(timelist,image_list))]
#################################################################################

#%matplotlib inline
#daynum = str.split(str.split(astlist[i],'/')[-1],'.')[0][3:6]
print('Processing '+str(N_images)+' images from day number '+str(daynum))
with PdfPages(out_dir+str(daynum)+'.images.wcs.pdf') as pdf:
      #for i in range(0,4):
      for i in range(len(image_list)):
            name = str.split(str.split(astlist[i],'/')[-1],'.')[0]
            with fits.open(str(astlist[i])) as hdul:
                  im = hdul[0].data
                  hd = hdul[0].header
                  wcs=WCS(hd)
                  if hd['FILTER'] != 6:
                      minval = np.mean(im)-np.std(im)*0.125
                      maxval = np.mean(im)+np.std(im)*0.75
                      meanval=np.mean(im)
                      stdval=np.std(im)
                  if hd['FILTER'] == 6:
                      minval = np.mean(im)-np.std(im)*0.125
                      maxval = 3000.
                      meanval=np.mean(im)
                      stdval=np.std(im)                      
                  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.
                  fig = plt.figure(figsize=(4,4))
                  ax1 = fig.add_subplot(projection=wcs)
                  #im0 = ax1.imshow(im[2500:4095,2500:4095],vmin=minval,vmax=maxval,cmap='cividis')
                  im0 = ax1.imshow(im,vmin=minval,vmax=maxval,cmap='cividis')
                  #ax1.annotate(str(name)+'  ('+str(hd['FILTER'])+ ' filter)', xy=(35,300), fontsize=7, color='white')               
                  #ax1.annotate(str(name)+'  ('+str(hd['OBJECT'])+', '+str(hd['FILTER'])+ ' filter)', xy=(100,200), fontsize=7, color='white')
                  ax1.annotate('UTC  '+str.split(hd['DATE-OBS'],'T')[0]+'  '+str.split(hd['DATE-OBS'],'T')[1],xy=(35,300), fontsize=7, color='white')
                  #ax1.annotate(r'Mean='+str("{:.2f}".format(meanval))+',  $\sigma$='+str("{:.2f}".format(stdval)), xy=(35,260), fontsize=7, color='white')
                  ax1.annotate(r'Mean='+str("{:.2f}".format(meanval))+',  $\sigma$='+str("{:.2f}".format(stdval)), xy=(35,260), fontsize=7, color='white')

                  #ax1.annotate(r'$\sigma$='+str("{:.2f}".format(stdval)), xy=(35,300), fontsize=7, color='white')
                  #ax1.annotate(r'Exposure Time ='+str(round(hd['EXPTIME'],3))+' s', xy=(50,600), fontsize=7, color='white')
                  #ax1.annotate(hd['READOUTM'],xy=(50,700), fontsize=7, color='white')
                  #ax1.annotate(r'min = mean - 1 $\sigma$ = '+str(round(minval,3)),xy=(50,800), fontsize=7, color='white')
                  #ax1.annotate(r'max = mean + 2.5 $\sigma$ = '+str(round(maxval,3)),xy=(50,900), fontsize=7, color='white')
                  ax1.annotate(r'min=mean-0.125$\sigma$='+str(round(minval,2)),xy=(35,220), fontsize=7, color='white')
                  ax1.annotate(r'max=mean+0.75$\sigma$='+str(round(maxval,2)),xy=(35,180), fontsize=7, color='white')
                  #ax1.annotate(r'Moon angle ='+str(round(hd['MOONANGL'],2))+'$^{\circ}$',xy=(50,1000), fontsize=7, color='white')
                  #ax1.annotate(r'Moon phase ='+str(round(hd['MOONPHAS'],2))+'%',xy=(50,1100), fontsize=7, color='white')
                  ax1.annotate(r'Azim.='+str(round(azimdeg,2))+'$^{\circ}$, '+'Elev.='+str(round(elevdeg,2))+'$^{\circ}$',xy=(35,140), fontsize=7, color='white')

                  #ax1.annotate(r'Azimuth ='+str(round(azimdeg,2))+'$^{\circ}$',xy=(35,180), fontsize=7, color='white')
                  #ax1.annotate(r'Elevation ='+str(round(elevdeg,2))+'$^{\circ}$',xy=(35,140), fontsize=7, color='white')
                  ax1.annotate(r'Hour Angle ='+str(ha)+'$^{\circ}$',xy=(35,100), fontsize=7, color='white')
                  ax1.annotate(r'Wind speed ='+str(round(hd['WXWNDSPD'],2))+' km/hr = '+str(round((hd['WXWNDSPD']/1.609344),2))+' miles/hr',xy=(35,60), fontsize=7, color='white')
                  ax1.annotate(r'Wind Direction ='+str(round(hd['WXWNDDIR'],2))+'$^{\circ}$ E of N',xy=(35,20), fontsize=7, color='white')
                  ax1.coords[0].set_auto_axislabel(False)
                  ax1.coords[1].set_auto_axislabel(False)
                  plt.xlabel("Right Ascension (J2000)",fontsize = 10)
                  plt.ylabel("Declination (J2000)",fontsize = 10)
                  plt.title(str(name)+'  ('+str(hd['FILTER'])+ ' filter)')
                  plt.tight_layout()
                  pdf.savefig(dpi=100)
                  plt.close()
                  print('Image #'+str(i+1)+' processed, '+str(round(100.*((i+1)/len(astlist)),1))+'% complete')

if obscode == '*':
    print('File '+out_dir+str(daynum)+'.images.wcs.pdf successfully created')
if obscode != '*':
    shutil.move(out_dir+str(daynum)+'.images.wcs.pdf',out_dir+str(daynum)+'.'+str(obscode)+'.images.wcs.pdf')
    print('File '+out_dir+str(daynum)+'.'+str(obscode)+'.images.wcs.pdf successfully created')
#plt.close()