;==============================================================================
;+
; rivm_wrfits - Write reduced data to a fits file at the end of rivm from
;	    data, mask, and noise arrays, or can read from existing rivm
;	    output binary array files and write fits.  Also writes oblk
;	    structure (see ot_header for fields), as a byte array extension.
;
; CALLING SEQUENCE:
;	rivm_wrfits, oblk, rivmdat, rivmfile, mask, noise, lcarr, verbose
;
; DESCRIPTION:
;     Produce FITS output given rivm output arrays or files using these steps:
;	1. Convert oblk to fits header
;	2. Add header values for data size/shape/history
;	3. Write data array with header to fits file
;	4. Add extension containing mask array
;	5. Add extension containing noise array
;	6. Add extension containing oblk array
;	7. Add extension containing line-center array
;
; INPUT:
;	oblk = Input oblk structure (can be null if oblkfile is set)
;	rivmdat = Input data cube (can be null if rivmfile is set)
;	rivmfile = Rivm-output filename (as propagated through rivm)
;	mask = Mask array (can be null if maskfile is set)
;	noise = Noise array (can be null if noisefile is set)
;	lcarr = Line center array
;	verbose = 0=none, 1=minimal
;
; OPTIONAL KEYWORDS:
;	fitsfile  = Name of FITS file to write to [I/O: will be set if not input]
;	oblkfile  = Name of oblk file to read from
;	maskfile  = Name of mask file to read from
;	noisefile = Name of noise file to read from
;
; OUPUT:  [fits file is written to file specified by 'fitsfile']
;
; EXAMPLE:
;    wdir = '/data/SOLMAG/DATA/IVM/ARCHIVE_REDUCTION/2001/20010524/test/'
;    rivmfile = wdir+'IMG0_20010524.1630_df_rf_pl_sc_ds_sub_av_st_ca_ct_unb'
;    oblkfile = wdir+'oblk_20010524.1630'
;    maskfile = wdir+'mask_20010524.1630'
;    noisefile = wdir+'noise_20010524.1630'
;    lcarr = fltarr(512,512)
;    rivm_wrfits, oblk, rivmdat, rivmfile, mask, noise, lcarr, 2, oblkf=oblkfile, maskf=maskfile, noisef=noisefile
;
;    To Read the data array and extensions back in (or use rivm_rdfits.pro):
; 	frivmdat = readfits(fitsfile, fhead)
; 	fmask = readfits(fitsfile, fmaskhead, exten=1)
; 	fnoise = readfits(fitsfile, fnoisehead, exten=2)
; 	oblk_bytes = readfits(fitsfile, oblkhead, exten=3)
; 	flc = readfits(fitsfile, flchead, exten=4)
;
; HISTORY:
;    Written by Eric L. Wagner 2012.04.16, based on earlier test script
;    2012.04.27 ELW - Write oblk as third extension
;    2012.05.03 ELW - IMGS_ output filename, improve oblk output
;    2012.08.20 ELW - Added LC extension
;-
;==============================================================================

pro rivm_wrfits, oblk, rivmdat, rivmfile, mask, noise, lcarr, verbose, $
		fitsfile=fitsfile, oblkfile=oblkfile, $
		maskfile=maskfile, noisefile=noisefile

  compile_opt idl2
  pn = 'rivm_wrfits: '

  ;============================================================================
  ; Get data parameters from oblk array or file
  ;============================================================================
  if (keyword_set(oblkfile)) then begin
      oblk=ivm_rd_oblk(oblkfile, err=err)
      if (err) then begin
          print, pn, 'Warning!  Unable to read oblk.  Stopping.'
	  stop
      endif
  endif

  nx = long(oblk.spix)
  ny = long(oblk.ppix)
  np = oblk.nframes
  nw = oblk.n_image_sets

  if (n_elements(rivmdat) lt 120) then begin
      if (verbose gt 1) then print, pn, 'Reading rivm output from file: ',rivmfile
      rivmdat = fltarr(nx, ny, np, nw)
      openr, lun, rivmfile, /get
      readu, lun, rivmdat
      close, lun & free_lun, lun
  endif

  ;============================================================================
  ; Make filename for fits output if 'fitsfile' not set
  ;============================================================================
  if (~keyword_set(fitsfile)) then begin
      iloc = strpos(rivmfile, 'IMG0_', /reverse_search)
    ; fitsfile = strmid(rivmfile, 0, iloc) + 'IMGS_' + strmid(rivmfile, iloc+5, 13) + '.fits'
      ; ELW20121220:  Adding _bn for binned data for comparison:
      fitsfile = strmid(rivmfile, 0, iloc) + 'IMGS_' + strmid(rivmfile, iloc+5, 13) + ((nx eq 256) ? '_bn.fits' : '.fits')
      if (verbose gt 1) then print, pn, 'FITS filename set to: ',fitsfile
  endif
  rivm_basefile = file_basename(rivmfile)

  ;============================================================================
  ; Convert oblk, add fields, write data array
  ;============================================================================
  head = ivm_oblk2fits(oblk)
  sxaddpar, head, 'SIMPLE',  'T'
  sxaddpar, head, 'NAXIS',  4, 'NUMBER OF DIMENSIONS', AFTER='SIMPLE'
  sxaddpar, head, 'NAXIS1', nx, 'X PIXELS IN SERIAL DIRECTION (AXIS 1)', AFTER='NAXIS'
  sxaddpar, head, 'NAXIS2', ny, 'Y PIXELS IN PARALLEL DIRECTION (AXIS 2)', AFTER='NAXIS1'
  sxaddpar, head, 'NAXIS3', np, 'STOKES POLARIZATION STATES (AXIS 3)', AFTER='NAXIS2'
  sxaddpar, head, 'NAXIS4', nw, 'WAVELENGTHS (AXIS 4)', AFTER='NAXIS3'
  sxaddpar, head, 'EXTEND', 'T', 'FITS FILE CONTAINS EXTENSIONS', AFTER='NAXIS4'
  sxaddpar, head, 'INSTRUME', 'IVM', 'Instrument name'
  sxaddpar, head, 'COMMENT', 'Final wavelength point is at continuum.'
  sxaddpar, head, 'HISTORY', oblk.comments
  sxaddpar, head, 'HISTORY', 'Reduced with Rivm V5.0 on ' + systime()
  sxaddpar, head, 'HISTORY', 'Rivm filename: ' + rivm_basefile
  if (verbose gt 0) then print, pn, 'Writing FITS file: ', fitsfile
  writefits, fitsfile, rivmdat, head

  ;============================================================================
  ; Extension #1:  Mask - read from file if maskfile set
  ;============================================================================
  if (keyword_set(maskfile)) then begin
      maskSet = ivm_rd_mask(nx, maskFile, mask)
      if (~maskSet) then begin
          print, pn, 'Warning!  Mask not set.  Stopping.'
	  stop
      endif
  endif
  sxaddpar, maskhead, 'XTENSION',  'IMAGE   ', 'MASK IMAGE EXTENSION'
  sxaddpar, maskhead, 'NAXIS',  2, 'NUMBER OF DIMENSIONS', AFTER='XTENSION'
  sxaddpar, maskhead, 'NAXIS1', nx, 'X PIXELS IN SERIAL DIRECTION (AXIS 1)', AFTER='NAXIS'
  sxaddpar, maskhead, 'NAXIS2', ny, 'Y PIXELS IN PARALLEL DIRECTION (AXIS 2)', AFTER='NAXIS1'
  sxaddpar, maskhead, 'PCOUNT', 0, '', AFTER='NAXIS2'
  sxaddpar, maskhead, 'GCOUNT', 1, '', AFTER='PCOUNT'
  writefits, fitsfile, mask, maskhead, /append

  ;============================================================================
  ; Extension #2:  Noise - read from file if noisefile set
  ;============================================================================
  if (keyword_set(noisefile)) then begin
      noise=fltarr(nx,ny,np)
      openr, lun, noisefile, error=err, /get
      if (err) then begin
          print, pn, 'Warning!  Unable to open noise file.  Stopping.'
	  stop
      endif
      readu, lun, noise
      close, lun & free_lun, lun
  endif
  sxaddpar, noisehead, 'XTENSION',  'IMAGE   ', 'COMPUTED NOISE IMAGE EXTENSION'
  sxaddpar, noisehead, 'NAXIS',  3, 'NUMBER OF DIMENSIONS', AFTER='XTENSION'
  sxaddpar, noisehead, 'NAXIS1', nx, 'X PIXELS IN SERIAL DIRECTION (AXIS 1)', AFTER='NAXIS'
  sxaddpar, noisehead, 'NAXIS2', ny, 'Y PIXELS IN PARALLEL DIRECTION (AXIS 2)', AFTER='NAXIS1'
  sxaddpar, noisehead, 'NAXIS3', np, 'STOKES POLARIZATION STATES (AXIS 3)', AFTER='NAXIS2'
  sxaddpar, noisehead, 'PCOUNT', 0, '', AFTER='NAXIS3'
  sxaddpar, noisehead, 'GCOUNT', 1, '', AFTER='PCOUNT'
  writefits, fitsfile, noise, noisehead, /append

  ;============================================================================
  ; Extension #3:  Set up and write an extension for oblk
  ;============================================================================
  obytes = ivm_oblk2raw(oblk)
  sxaddpar, oblkhead, 'XTENSION',   'IMAGE   ', '2048 BYTE OBLK STRUCTURE'
  sxaddpar, oblkhead, 'NAXIS',  1,  'NUMBER OF DIMENSIONS', AFTER='XTENSION'
  sxaddpar, oblkhead, 'NAXIS1', 2048, 'BYTES IN IVM OBLK STRUCTURE', AFTER='NAXIS'
  sxaddpar, oblkhead, 'PCOUNT', 0, '', AFTER='NAXIS1'
  sxaddpar, oblkhead, 'GCOUNT', 1, '', AFTER='PCOUNT'
  writefits, fitsfile, obytes, oblkhead, /append

  ;============================================================================
  ; Extension #4:  Set up and write an extension for line center
  ;============================================================================
  if (n_elements(lcarr) eq nx*ny) then begin
      sxaddpar, lchead, 'XTENSION',  'IMAGE   ', 'COMPUTED LINE CENTER IMAGE EXTENSION'
      sxaddpar, lchead, 'NAXIS',  2, 'NUMBER OF DIMENSIONS', AFTER='XTENSION'
      sxaddpar, lchead, 'NAXIS1', nx, 'X PIXELS IN SERIAL DIRECTION (AXIS 1)', AFTER='NAXIS'
      sxaddpar, lchead, 'NAXIS2', ny, 'Y PIXELS IN PARALLEL DIRECTION (AXIS 2)', AFTER='NAXIS1'
      sxaddpar, lchead, 'PCOUNT', 0, '', AFTER='NAXIS2'
      sxaddpar, lchead, 'GCOUNT', 1, '', AFTER='PCOUNT'
      writefits, fitsfile, lcarr, lchead, /append
  endif

end