#################################################################################
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 = '/mnt/imagesbucket/*/*/2023/'+str(daynum)+'/'
if obscode != '*':
    image_dir = '/mnt/imagesbucket/*/'+str(obscode)+'/2023/'+str(daynum)+'/'
print(image_dir)
out_dir = current_dir+'/'
#os.makedirs(out_dir, exist_ok = True)
image_list = sorted(glob.glob(image_dir+'*.fts',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[256:3840,256:3840])-np.std(im[256:3840,256:3840])*0.5
                      maxval = np.mean(im[256:3840,256:3840])+np.std(im[256:3840,256:3840])*2
                      meanval=np.mean(im[256:3840,256:3840])
                      stdval=np.std(im[256:3840,256:3840])
                  if hd['FILTER'] == 6:
                      minval = np.mean(im[256:3840,256:3840])-np.std(im[256:3840,256:3840])*0.125
                      maxval = 3000.
                      meanval=np.mean(im[256:3840,256:3840])
                      stdval=np.std(im[256:3840,256:3840])                      
                  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=(8,8))
                  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(r'($\alpha$,$\delta$) = ('+str(hd['OBJRA'])+','+str(hd['OBJDEC'])+')',xy=(35,1600), fontsize=12, color='white')
                  ax1.annotate('UTC  '+str.split(hd['DATE-OBS'],'T')[0]+'  '+str.split(hd['DATE-OBS'],'T')[1],xy=(35,1450), fontsize=10, color='white')
                  ax1.annotate(r'Mean='+str("{:.2f}".format(meanval))+',  $\sigma$='+str("{:.2f}".format(stdval)), xy=(35,1300), fontsize=10, color='white')
                  ax1.annotate(r'min=mean-0.5$\sigma$='+str(round(minval,2)),xy=(35,1150), fontsize=10, color='white')
                  ax1.annotate(r'max=mean+2$\sigma$='+str(round(maxval,2)),xy=(35,1000), fontsize=10, color='white')
                  ax1.annotate(r'Moon angle ='+str(round(hd['MOONANGL'],2))+'$^{\circ}$',xy=(35,850), fontsize=10, color='white')
                  ax1.annotate(r'Moon phase ='+str(round(hd['MOONPHAS'],2))+'%',xy=(35,700), fontsize=10, color='white')
                  ax1.annotate(r'Azim.='+str(round(azimdeg,2))+'$^{\circ}$, '+'Elev.='+str(round(elevdeg,2))+'$^{\circ}$',xy=(35,550), fontsize=10, color='white')
                  ax1.annotate(r'Hour Angle ='+str(ha)+'$^{\circ}$',xy=(35,400), fontsize=10, 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,250), fontsize=10, color='white')
                  ax1.annotate(r'Wind Direction ='+str(round(hd['WXWNDDIR'],2))+'$^{\circ}$ E of N',xy=(35,100), fontsize=10, 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()