import matplotlib.pyplot as plt
import numpy as np

from astropy.io import fits

#
# read HLSP
#
hlsp_file = 'TWA_21_2010-02-18_NaCo_H_diva_hlsp.fits'

hdu = fits.open(hlsp_file)
hdr = hdu[0].header

#
# list extensions
#
nextend = hdr['NEXTEND']

print()
print('File has {0} extensions:'.format(nextend))
extensions = {}
for i in range(1, nextend+1):
    extname = hdr['EXT{0}NAME'.format(i)]
    exttype = hdr['EXT{0}TYPE'.format(i)]

    print(' * {0}; {1}'.format(extname, exttype))

    # save for later
    extensions[extname] = i

# retrieve basic information
surv   = hdr['SURVNAME']
survpi = hdr['SURVPI']
instru = hdr['INSTRUME']
filt   = hdr['FILTER']
target = hdr['TARGNAME']

#
# data information
#
print()
print('Data information:')

idx = extensions.get('DATA_INFORMATION')
data_info = hdu[idx].columns

# list columns
print(data_info)

# retrieve obs information
date = data_info['UT_Midpoint_Date_of_Observation'].array[0]
Texp = data_info['Exposure_Time'].array[0]
FoV_rot = data_info['Combined_rotation_angle'].array[0]

#
# reduced image
#
idx = extensions.get('REDUCED_DATA')
if (idx is not None):
    img = hdu[idx].data

    # display
    fig = plt.figure(0, figsize=(10, 10))
    plt.clf()
    plt.imshow(img, vmin=-0.005*np.nanmax(img), vmax=0.005*np.nanmax(img), aspect='equal')
    plt.figtext(0.07, 0.9, 'Texp = {0:.1f} min \nFoV rot = {1:.1f} deg'.format(Texp/60, FoV_rot), fontsize=15)
    plt.title('{0} / {1} {2} / {3} / {4} ({5})'.format(target, instru, filt, date, surv, survpi))
    plt.tight_layout()

#
# detection limit
#
idx = extensions.get('DETECTION_LIMIT')
if (idx is not None):
    dat = hdu[idx].data.columns

    # threshold
    nsigma = hdu[idx].header['NSIGMA']

    # values
    sep  = dat['Radius'].array
    dlim = dat['Detection_limit'].array

    # units
    sep_u  = dat['Radius'].unit
    dlim_u = dat['Detection_limit'].unit

    # plot
    fig = plt.figure(1, figsize=(10, 7))
    plt.clf()
    plt.plot(sep, dlim, linestyle='-', color='b')
    plt.gca().invert_yaxis()
    plt.xlabel(r'Angular separation [{0}]'.format(sep_u))
    plt.ylabel(r'Sensitivity limit @ {0:.0f}$\sigma$ [{1}]'.format(nsigma, dlim_u))
    plt.title(r'{0} / {1} {2} / {3} / {4} ({5})'.format(target, instru, filt, date, surv, survpi))
    plt.tight_layout()

#
# candidates
#
print()
print('Detected candidates:')

candnum = hdr['CANDNUM']
if (candnum == -1):
    print(' * this survey has no information on detected candidates')
elif (candnum > 0):
    idx = extensions.get('SOURCE_DETECTION')
    dat = hdu[idx].data.columns

    # list properties
    print('idx  dRA        err_dRA   dDEC      err_dDEC   dmag    err_dmag')
    print('---------------------------------------------------------------')
    for c in range(candnum):
        print('{:2d}{:10.4f}{:10.4f}{:10.4f}{:10.4f}{:10.2f}{:7.2f}'.
              format(c, dat['dRA'].array[c], dat['err_dRA'].array[c],
                        dat['dDEC'].array[c], dat['err_dDEC'].array[c],
                        dat['Contrast'].array[c], dat['err_Contrast'].array[c]))
      

;;
;; read HLSP
;;
hlsp_file = 'TWA_21_2010-02-18_NaCo_H_diva_hlsp.fits'

hdu = mrdfits(hlsp_file, 0, hdr, /silent)

;;
;; list extensions
;;
nextend = sxpar(hdr, 'NEXTEND')

print
print,'File has '+string(nextend, format='(I0)')+' extensions:'
extensions = ['']
for i=1,nextend do begin
   extname = sxpar(hdr, string(i, format='("EXT",I0,"NAME")'))
   exttype = sxpar(hdr, string(i, format='("EXT",I0,"TYPE")'))

   print,' * '+extname+'; '+exttype

   ;; save for later
   extensions = [extensions, extname]
endfor

;; retrieve basic information
surv   = strtrim(sxpar(hdr, 'SURVNAME'))
survpi = strtrim(sxpar(hdr, 'SURVPI'))
instru = strtrim(sxpar(hdr, 'INSTRUME'))
filt   = strtrim(sxpar(hdr, 'FILTER'))
target = strtrim(sxpar(hdr, 'TARGNAME'))

;;
;; data information
;;
print
print,'Data information:'

idx = where(extensions eq 'DATA_INFORMATION')
data_info = mrdfits(hlsp_file, idx[0], /silent)

;; list colmuns
help,/str,data_info

;; retrieve obs information
date = data_info.UT_Midpoint_Date_of_Observation[0]
Texp = data_info.Exposure_Time[0]
FoV_rot = data_info.Combined_rotation_angle[0]

;;
;; reduced image
;;
idx = where(extensions eq 'REDUCED_DATA')
if (idx ne -1) then begin
   img = mrdfits(hlsp_file, idx[0], /silent)

   img[where(finite(img) eq 0)] = 0
   img <=  0.005*max(img)
   img >= -0.005*max(img)

   window, 0, xs=600, ys=600
   tvscl, congrid(img, 600, 600)
   xyouts, 10, 10, target+' / '+instru+' '+filt+' / '+date+' / '+surv+' ('+survpi+')', $
           color=255, /device, size=1.8
endif

;;
;; detection limit
;;
idx = where(extensions eq 'DETECTION_LIMIT')
if (idx ne -1) then begin
   dat = mrdfits(hlsp_file, idx[0], chdr, /silent)

   ;; threshold
   nsigma = sxpar(chdr, 'NSIGMA')

   ;; values
   sep  = dat.Radius
   dlim = dat.Detection_limit

   ;; plot
   window, 1, xs=640, ys=480
   plot,sep,dlim,yr=[max(dlim, /nan), min(dlim, /nan)], xtitle='Angular separation', $
        ytitle='Sensitivity limite @ '+string(nsigma, format='(I0)')+'-sigma', $
        title=target+' / '+instru+' '+filt+' / '+date+' / '+surv+' ('+survpi+')'
endif

;;
;; candidates
;;
print
print,'Detected candidates:'

candnum = sxpar(hdr, 'CANDNUM')
if (candnum eq -1) then begin
   print,' * this survey has no information on detected candidates'
endif else if (candnum gt 0) then begin
   idx = where(extensions eq 'SOURCE_DETECTION')
   dat = mrdfits(hlsp_file, idx[0], /silent)

   ;; list properties
   print,'idx  dRA        err_dRA   dDEC      err_dDEC   dmag    err_dmag'
   print,'---------------------------------------------------------------'
   for c=0,candnum-1 do begin
      print, c, dat[c].dRA, dat[c].err_dRA, dat[c].dDEC, dat[c].err_dDEC, $
             dat[c].Contrast, dat[c].err_Contrast,format='(I2,F10.4,F10.4,F10.4,F10.4,F10.2,F7.2)'
   endfor

endif

end