;+
; NAME:
;  spe2fits
; PURPOSE:   (one line only)
;  Convert WinView SPE files to FITS format files
; DESCRIPTION:
;  This routine is used to process one directory of files at a time and
;    convert the data from SPE format to FITS format.  The SPE data do not
;    include all the information needed for the header.  Therefore, additional
;    information must be provided to support the converstion.
;  All files that also have supporting information will be written out if
;    they don't already exist.
; CATEGORY:
;  CCD data processing
; CALLING SEQUENCE:
;  spe2fits
; INPUTS:
; OPTIONAL INPUT PARAMETERS:
; KEYWORD INPUT PARAMETERS:
;  SPEDIR - String, name of directory to process.  Default is current directory.
; OUTPUTS:
; KEYWORD OUTPUT PARAMETERS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
;  The SPE files need to be constructed in the following way:
;    YYYYMMDD_N.SPE where YYYY is the year, MM is the numeric month, DD
;    is the day, and N is the file number (variable number of digits).
;
;  The output file names will be reconstructed like this:
;    YYYYMMDD.NNN.fits where NNN ends up being zero padded to a 3-digit number.
;
;  The key element of the incoming file name is the underscore.  This serves
;    to delineate the root of the file name from the file number.  If the
;    underscore is not found the program will quit.
;
; PROCEDURE:
;  A list of SPE files is generated based on the files found in the input
;    directory.  This list is built from all files ending in .SPE.  Note that
;    the case of the file name is ignored.  These are the candidate files
;    for processing.
;
;  A file named, header.toc, must exist and by default is in the directory with the data.  This file
;    contains supplemental information that will be added to the image headers.
;    Only those files with entries in this file will be processed.
;    The file contains the following columns (space-delimited):
;      1 - File number (integer) of the file to be read (file name suffix)
;      2 - Object name (string).  Embedded blanks are not allowed.
;             The name can be a NAIF/geteph object code which is useful
;             if you need to calculate the position to populate the header.
;             There is a special case supported.  The name can be two fields
;             separated by a colon (:).  The first field is a NAIF/geteph
;             object code (ie., P301 for the Moon) and the second field is
;             some location on the body (ie., Tycho for the name of the crater).
;             The object code is used (if desired) to generate the target
;             position.  The NAIF code is converted to a common name before
;             writing to the header.
;      3 - dt (float)  This is the inter-pulse period as programmed on the
;             GPS-based triggering device.  This is in units of seconds.
;             If you are not using GPS triggering this value should be -1.0
;      4 - UT date of observation.  YYYY/MM/DD format.  YYYY-MM-DD is also
;             allowed.  This must always be provided.
;      5 - UT time.  HH:MM:SS.sss format.  Fractional seconds are optional.
;             If the time string is "spe" (no quotes) then the UT time from
;             the header is used.  In this case, the start time is effectively
;             derived from the system clock on the computer taking the data.
;      6 - Right ascension of image/telescope, expected to be somewhere near
;             the center of the image.
;      7 - Declination of image/telescope.
;             If ra or dec have a value of "calculate" (no quotes), then
;             an ephemeris of the object
;      8 - Equinox of coordinates in columns 6 and 7.
;      9 - Name of filter (string, no embedded blanks)
;
;
;  Next, a file named, site.dat, must be found and then loaded.  This file
;    contains the static information for the night.  This file can have
;    embedded blanks, where appropriate.  The following lines must be in
;    the file, in this order:
;      Line 1: Name of the observatory
;      Line 2: Latitute of observatory (N or S)DD:MM:SS.sss, WGS84 datum
;      Line 3: Longitude of observatory (E or W)DDD:MM:SS.sss
;      Line 4: Observatory elevation in meters
;      Line 5: Name(s) of observer(s)
;      Line 6: Name of camera
;      Line 7: Name of telescope (should include aperture)
;
;  If a file is found having the same root name with a suffix of .hdr, this
;    file is loaded and added to the FITS header before writing.  It will
;    appear after the SPE header data.  Note that these header lines are
;    expected to be pre-formatted FITS header lines and will be transferred
;    with no processing at all.  There must be no END statement in this file
;    and it should not duplicate keywords already present in the data as
;    created and manipulated by this program.
;
;  Each file is then processed to combine the embedded SPE header information
;    plus the extra information and then written out as a FITS file.
;
; MODIFICATION HISTORY:
;  Written by Marc W. Buie, Southwest Research Institute, 2009/10/13
;  2009/10/19, MWB, added optional header support
;-
pro spe2fits,SPEDIR=spedir,OUTDIR=outdir,AUXDIR=auxdir

   self='SPE2FITS: '
   version = 4
   if badpar(spedir,[0,7],0,caller=self+'(SPEDIR) ',default='') then return
   if spedir ne '' then begin
      if not exists(spedir) then begin
         print,self,'SPEdir, [',spedir,'] does not exist.  Aborting.'
         return
      endif
      spedir=addslash(spedir)
   endif

   if badpar(spedir,[0,7],0,caller=self+'(AUXDIR) ',default='') then return
   if auxdir ne '' then begin
      if not exists(auxdir) then begin
         print,self,'AUXdir, [',auxdir,'] does not exist.  Aborting.'
         return
      endif
      auxdir=addslash(auxdir)
   endif else begin
      auxdir = spedir
   endelse

   if badpar(outdir,[0,7],0,caller=self+'(OUTDIR) ',default='') then return
   if outdir ne '' then begin
      if not exists(outdir) then begin
         ans=''
         print,outdir
         read,prompt='Output directory does not exists.  Create and continue?',ans
         if strlowcase(strmid(ans,0,1)) eq 'y' then begin
            file_mkdir,outdir
         endif else begin
            return
         endelse
      endif
      outdir=addslash(outdir)
   endif

   pushd,spedir
   fnlist = file_search('*.spe',count=nfiles,/fold_case)
   popd

   print,'Found ',strn(nfiles),' files to process in directory ',spedir

   ; Load the toc file
   if not exists(auxdir+'header.toc') then begin
      print,self,' Header information file could not be found.'
      return
   endif
   readcol,auxdir+'header.toc',num,object,dt,date,time,ra,dec,epoch,filter, $
      format='i,a,f,a,a,a,a,f,a'

   ; Load the site information file
   if not exists(auxdir+'site.dat') then begin
      print,self,' Site information file could not be found.'
      return
   endif
   openr,lun,auxdir+'site.dat',/get_lun
   obsname=''
   lats=''
   lons=''
   alt=''
   observers=''
   camera=''
   telescope=''
   readf,lun,obsname,format='(a)'
   readf,lun,lats,format='(a)'
   readf,lun,lons,format='(a)'

   ; This is a trick to ensure the lat/lon is properly formatted (despite
   ;   documentation header information).
   cvtsixty,lats,-0.5d0*!dpi,0.5d0*!dpi,0,['N','S'],lat
   cvtsixty,lat,-0.5d0*!dpi,0.5d0*!dpi,0,['N','S'],lats

   cvtsixty,lons,-1.0d0*!dpi,1.0d0*!dpi,1,['W','E'],lon
   cvtsixty,lon,-1.0d0*!dpi,1.0d0*!dpi,1,['W','E'],lons

   readf,lun,alt,format='(a)'
   alt = float(alt)
   readf,lun,observers,format='(a)'
   readf,lun,camera,format='(a)'
   readf,lun,telescope,format='(a)'
   free_lun,lun

   ; eye candy
   print,obsname,' ',lats,' ',lons,' ',alt

   for i=0,nfiles-1 do begin
      words=strsplit(fnlist[i],'.',/extract)
      words=strsplit(words[0],'_',/extract)
      if n_elements(words) ne 2 then begin
         print,self,fnlist[i],' file name is not in standard form.  Skipping.'
         continue
      endif
      root = words[0]
      fno  = fix(words[1])

      outname = root+'.'+string(fno,format='(i3.3)')+'.fits'
      if exists(outdir+outname) then continue
      print,fnlist[i],'-->',outname

      ; need to provide the following keywords for the new FITS header
      ; DATE-OBS  (YYYY-MM-DD)    - from toc or spe header
      ; UT                        - from toc or spe header
      ; OBJECT                    - from toc
      ; RA                        - from toc
      ; DEC                       - from toc
      ; EQUINOX                   - from toc
      ; FILTER                    - from toc
      ; EXPTIME                   - from spe header
      ; EXPDELTA                  - from spe header?
      ; FILENAME                  - internal
      ; AIRMASS                   - Calculated

      z=where(fno eq num,count)
      if count ne 1 then begin
         print,self,'No header information found for ',fnlist[i]
         continue
      endif
      z=z[0]

      ; Read the data file
      read_princeton_gz,spedir+fnlist[i],data,header=header,compress=-1
      mkhdr,hdr,data
      hdr=spe_head2fits(header,infitsh=hdr)

      exptime = header.exp_sec
      expdelta = exptime + header.delaytime

      ; SPE date is contained within DATE as DDmonYYYY, mon is name
      ;  however, this is the LOCAL time and is not trustworthy.  Do not use.

      ; SPE time is contained within EXTIMUTC as hhmmss string (no colons)
      if time[z] eq 'spe' then begin
         timestr = sxpar(hdr,'EXTIMUTC')
         time[z] = strmid(timestr,0,2)+':'+ $
                   strmid(timestr,2,2)+':'+ $
                   strmid(timestr,4,2)
      endif

      jd = jdparse(date[z]+' '+time[z])

      jdstr,jd,0,jds,timesep='-'
      words=strsplit(jds,' ',/extract)
      date[z] = words[0]
      time[z] = words[1]

      namfld = strsplit(object[z],':',/extract)

      if ra[z] eq 'calculate' or dec[z] eq 'calculate' then begin
stop
         ephem,jd,500,2,namfld[0],eph
         rastr,eph[0],1,rah
         decstr,eph[1],0,dech
      endif else begin
         rah = ra[z]
         dech = dec[z]
      endelse

      objname = naifname(namfld[0])
      if n_elements(namfld) gt 1 then begin
         objname = objname+':'+namfld[1]
      endif

      ; compute the airmass
      era  = raparse(rah)
      edec = decparse(dech)
      am=airmass(jd,era,edec,lat,lon)
      am=am[0]

      sxaddpar,hdr,'SPE2FITS', version, before='DATE', 'SPE2FITS version'
      sxaddpar,hdr,'DATE-OBS',date[z]+'T'+time[z],before='DATE', $
         ' UT date and time of start'
      sxaddpar,hdr,'UT',time[z],before='DATE', $
         ' UT time of start'
      sxaddpar,hdr,'OBJECT',objname,before='DATE'
      sxaddpar,hdr,'RA',rah,before='DATE'
      sxaddpar,hdr,'DEC',dech,before='DATE'
      sxaddpar,hdr,'EQUINOX',epoch[z],before='DATE'
      sxaddpar,hdr,'FILTER',filter[z],before='DATE'
      sxaddpar,hdr,'EXPTIME',exptime,before='DATE', $
         ' commanded exposure time in seconds'
      if size(data,/n_dimensions) eq 3 then begin
         if dt[z] gt 0 then begin
            sxaddpar,hdr,'EXPDELTA',dt[z],before='DATE', $
               ' interval between successive images (sec)'
         endif else begin
            readt = sxpar(hdr,'READOUTT')
            delt  = sxpar(hdr,'EXP_SEC')
            if readt gt delt then delt=readt
            sxaddpar,hdr,'EXPDELTA',delt,before='DATE', $
               ' interval between successive images (sec)'
         endelse
      endif
      sxaddpar,hdr,'FILENAME',outname,before='DATE'
      sxaddpar,hdr,'AIRMASS',am,before='DATE'
      sxaddpar,hdr,'OBSLAT',lats,before='DATE',' Location, WGS84 datum'
      sxaddpar,hdr,'OBSLON',lons,before='DATE'
      sxaddpar,hdr,'OBSALT',alt,before='DATE',' in meters'
      sxaddpar,hdr,'OBSERVER',observers,before='DATE'
      sxaddpar,hdr,'CAMERA',camera,before='DATE'
      sxaddpar,hdr,'TELESCOP',telescope,before='DATE'
      sxaddpar,hdr,'BSCALE',1.0,before='DATE'
      sxaddpar,hdr,'BZERO',32768,before='DATE'

      ; Look for extra header file, if found load it and add it to header
      hdrnam=strmid(fnlist[i],0,strlen(fnlist[i])-4)+'.hdr'
      if exists(spedir+hdrnam) then begin
         nextra=file_lines(spedir+hdrnam)
         line=''
         xhdr=strarr(nextra)
         openr,lun,spedir+hdrnam,/get_lun
         for j=0,nextra-1 do begin
            readf,lun,line,format='(a)'
            xhdr[j] = line
         endfor
         free_lun,lun
         ; scan header and look for END
         for j=n_elements(hdr)-1,0,-1 do begin
            if strmid(hdr[j],0,20) eq 'END                 ' then break
         endfor
         if j gt 0 then hdr = [hdr[0:j-1],xhdr,hdr[j:*]]
      endif

      writefits,outdir+outname,data,hdr

   endfor

end