#################################################################################
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)+'/'
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('-----------------------------------------------------------')
#################################################################################

#################################################################################
# Plot temp vs. time:
templist = []
timelist = []
amblist=[]
sinklist=[]
explist = []
ratio=[]
wndlist=[]
for ii in range(0, N_images):
    im = fits.getdata(image_list[ii])
    hd = fits.getheader(image_list[ii])
    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)

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='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.png',dpi=100, bbox_inches = 'tight', transparent=False, facecolor='whitesmoke')
if obscode != '*':
    plt.savefig(str(daynum)+'.'+str(obscode)+'.plots.png',dpi=100, bbox_inches = 'tight', transparent=False, facecolor='whitesmoke')
#################################################################################

#################################################################################
# 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.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
                  if hd['FILTER'] != 6:
                      minval = np.mean(im[20:4075,20:4075])-np.std(im[20:4075,20:4075])*0.125
                      maxval = np.mean(im[20:4075,20:4075])+np.std(im[20:4075,20:4075])*0.75
                      meanval=np.mean(im[20:4075,20:4075])
                      stdval=np.std(im[20:4075,20:4075])
                  if hd['FILTER'] == 6:
                      minval = np.mean(im[20:4075,20:4075])-np.std(im[20:4075,20:4075])*0.125
                      maxval = 3000.
                      meanval=np.mean(im[20:4075,20:4075])
                      stdval=np.std(im[20:4075,20:4075])                      
                  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 = None)
                  #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['OBJECT'])+', '+str(hd['FILTER'])+ ' filter)', xy=(100,200), fontsize=10, color='white')
                  ax1.annotate('UTC  '+str.split(hd['DATE-OBS'],'T')[0]+'  '+str.split(hd['DATE-OBS'],'T')[1],xy=(100,300), fontsize=10, color='white')
                  ax1.annotate('Mean = '+str("{:.3f}".format(meanval)), xy=(100,400), fontsize=10, color='white')
                  ax1.annotate(r'$\sigma$ = '+str("{:.3f}".format(stdval)), xy=(100,500), fontsize=10, color='white')
                  ax1.annotate(r'Exposure Time ='+str(round(hd['EXPTIME'],3))+' s', xy=(100,600), fontsize=10, color='white')
                  ax1.annotate(hd['READOUTM'],xy=(100,700), fontsize=10, color='white')
                  #ax1.annotate(r'min = mean - 1 $\sigma$ = '+str(round(minval,3)),xy=(100,800), fontsize=10, color='white')
                  #ax1.annotate(r'max = mean + 2.5 $\sigma$ = '+str(round(maxval,3)),xy=(100,900), fontsize=10, color='white')
                  ax1.annotate(r'min = mean - 0.125 $\sigma$ = '+str(round(minval,3)),xy=(100,800), fontsize=10, color='white')
                  ax1.annotate(r'max = mean + 0.75 $\sigma$ = '+str(round(maxval,3)),xy=(100,900), fontsize=10, color='white')
                  ax1.annotate(r'Moon angle ='+str(round(hd['MOONANGL'],2))+'$^{\circ}$',xy=(100,1000), fontsize=10, color='white')
                  ax1.annotate(r'Moon phase ='+str(round(hd['MOONPHAS'],2))+'%',xy=(100,1100), fontsize=10, color='white')
                  ax1.annotate(r'Azimuth ='+str(round(azimdeg,2))+'$^{\circ}$',xy=(100,1200), fontsize=10, color='white')
                  ax1.annotate(r'Elevation ='+str(round(elevdeg,2))+'$^{\circ}$',xy=(100,1300), fontsize=10, color='white')
                  ax1.annotate(r'Hour Angle ='+str(ha)+'$^{\circ}$',xy=(100,1400), 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=(100,1500), fontsize=10, color='white')
                  ax1.annotate(r'Wind Direction ='+str(round(hd['WXWNDDIR'],2))+'$^{\circ}$ E of N',xy=(100,1600), fontsize=10, color='white')
                  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.pdf successfully created')
if obscode != '*':
    shutil.move(out_dir+str(daynum)+'.images.pdf',out_dir+str(daynum)+'.'+str(obscode)+'.images.pdf')
    print('File '+out_dir+str(daynum)+'.'+str(obscode)+'.images.pdf successfully created')
#plt.close()