;-----------------------------------------------------------------------
pro lwscoadd_koa, filename, outfile
;-----------------------------------------------------------------------
;+
; NAME:
;       LWSCOADD_KOA
;
; PURPOSE:
;       Read an FITS file generated by the W. M. Keck Observatory's
;       Long Wavelength Spectrometer (LWS) into IDL data and header
;       variables, coadds extensions, and writes result to output FITS file 
;;
; INPUTS:
;       FILENAME = Scalar string containing the name of the FITS file  
;                 (including extension) to be read.   If the filename has
;                  a *.gz extension, it will be treated as a gzip compressed
;                  file.   If it has a .Z extension, it will be treated as a
;                  Unix compressed file.
;
; OUTPUTS:
;       Result = Co-added (2-D) FITS file constructed from designated record.
;;
; EXAMPLE:
;       Read an LWS FITS file TEST.FITS into an IDL image array, coadd the
;       array and write to TEST_coadded.fits  FITS file in same directory
;
;              IDL> lwscoadd, TEST.FITS, TEST_coadded.fits
;
; ERROR HANDLING:
;       If an error is encountered reading the FITS file, then 
;       the error message is displayed 
;
; PROCEDURES USED:
;       Functions:   READFITS(), SXPAR(), WHERENAN()
;       Procedures:  IEEE_TO_HOST, SXADDPAR, SXDELPAR
;
; REQUIREMENTS:
; Needs routines listed above from the IDL Astronomy Library,
; http://idlastro.gsfc.nasa.gov/homepage.html
;
; NOTE:
;       Under Unix, LWSCOADD() can also read gzip or Unix compressed FITS files.
;       See http://idlastro.gsfc.nasa.gov/fitsio.html for other ways of
;       reading FITS files with IDL.
;
; AUTHOR:
; Gregory D. Wirth, W. M. Keck Observatory
;
; MODIFICATION HISTORY:
; 19-Jan-2000  GDW	Original version
; 22-Aug-2014  JH       Modifications for KOA: 
;                        -> Removed messages/returned-numbers
;                        -> added writefits
;
;-----------------------------------------------------------------------
message, /reset                 ; clear the error state

; define constants...
sz_lws    = 128                 ; define lws detector size
chop_beam = 3                   ; dimension representing chop beams
nod_beam  = 5                   ; dimension representing nod beams
chop_on   = 1                   ; index number of the chop "on" beam
chop_off  = 0                   ; index number of the chop "off" beam
nod_on    = 0                   ; index number of the nod "on" beam
nod_off   = 1                   ; index number of the nod "off" beam

header = headfits(filename, ext=0,/silent)

; read the FITS image...
image = readfits( filename, header)

; check for error...
if !error_state.code ne 0 then return

; verify datatype...
type = size( image, /type)
if type ne 3 then begin
    message, /info, 'image data type is not LONG (' + filename + ')'
    return
endif

; verify dimensions...
n_dim = size( image, /n_dimensions)
if n_dim lt 4 then begin
    message, /info,  $
      'image has too few dimensions (' + filename + ')'
    return
endif

; verify image size...
im_len = size( image)
if im_len[1] ne sz_lws or im_len[2] ne sz_lws then begin
    message, /info, 'image has wrong size (' + filename + ')'
    return
endif
if n_dim ge chop_beam then begin
    if im_len[chop_beam] gt 2 then begin
        message, /info,  $
          'image has more than 2 chopbeams (' + filename + ')'
        return
    endif
endif
if n_dim ge nod_beam then begin
    if im_len[nod_beam] gt 2 then begin
        message,  /info,  $
          'image has more than 2 nodbeams (' + filename + ')'
        return
    endif
endif

; create output array...
data = lonarr( im_len[1], im_len[2])

; mash nodsets...
if n_dim ge 3 then chpbeams = im_len[3] else chpbeams = 1
if n_dim ge 4 then savesets = im_len[4] else savesets = 1
if n_dim ge 5 then nodbeams = im_len[5] else nodbeams = 1
if n_dim ge 6 then nodsets  = im_len[6] else nodsets  = 1

; co-add the data...
for i=0,chpbeams-1 do begin      
    for j=0,savesets-1 do begin  
        for k=0,nodbeams-1 do begin 
            for l=0,nodsets-1 do begin
                ; extract the relevant data...
                if n_dim eq 3 then begin
                    buf = image[*,*,i]
                endif else if n_dim eq 4 then begin
                    buf = image[*,*,i,j]
                endif else if n_dim eq 5 then begin
                    buf = image[*,*,i,j,k]
                endif else if n_dim eq 6 then begin
                    buf = image[*,*,i,j,k,l]
                endif else begin
                    message, /info, 'Wrong number of dimensions'
                    return
                endelse

                ; presume line is to be added by default...
                pm = long(1)
                
                ; determine whether the chop beam is to be added or subtracted...
                if i eq CHOP_OFF then pm = -pm 
                
                ; determine whether the nod beam is to be added or subtracted...
                if k eq NOD_OFF then pm = -pm
                
                ; add this line into the final array...
                data = data + pm*buf
            endfor
        endfor
    endfor
endfor

writefits, outfile, data, header

return
end