;-----------------------------------------------------------------------
function lwscoadd, filename, header, SILENT=silent
;-----------------------------------------------------------------------
;+
; NAME:
;       LWSCOADD
;
; PURPOSE:
;       Read an FITS file generated by the W. M. Keck Observatory's
;       Long Wavelength Spectrometer (LWS) into IDL data and header
;       variables. 
;
; CALLING SEQUENCE:
;       Result = LWSCOADD( Filename,[ Header, /SILENT] )
;
; 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 data array constructed from designated record.
;                If the specified file was not found, then Result = -1
;
; OPTIONAL OUTPUT:
;       Header = String array containing the header from the FITS file.
;
; OPTIONAL INPUT KEYWORDS:
;
;       /SILENT - Normally, LWSCOADD will display the size the array at the
;                 terminal.  The SILENT keyword will suppress this
;
; EXAMPLE:
;       Read an LWS FITS file TEST.FITS into an IDL image array, IM and FITS 
;       header array, H.
;
;              IDL> im = LWSCOADD( 'TEST.FITS', h)
;
;       To read in a file that has been compressed:
;
;              IDL> im = LWSCOADD( 'test.fits.gz',h )
;
; ERROR HANDLING:
;       If an error is encountered reading the FITS file, then 
;               (1) the system variable !ERROR is set (via the MESSAGE facility)
;               (2) the error message is displayed (unless /SILENT is set),
;                   and the message is also stored in !ERR_STRING
;               (3) LWSCOADD returns with a value of -1
;
; 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
;-
;-----------------------------------------------------------------------

silent = keyword_set( SILENT)
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

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

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

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

; verify dimensions...
n_dim = size( image, /n_dimensions)
if n_dim lt 4 then begin
    message, noprint=silent, /info,  $
      'image has too few dimensions (' + filename + ')'
    return, -1
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, noprint=silent, /info, 'image has wrong size (' + filename + ')'
    return, -1
endif
if n_dim ge chop_beam then begin
    if im_len[chop_beam] gt 2 then begin
        message, noprint=silent, /info,  $
          'image has more than 2 chopbeams (' + filename + ')'
        return, -1
    endif
endif
if n_dim ge nod_beam then begin
    if im_len[nod_beam] gt 2 then begin
        message, noprint=silent, /info,  $
          'image has more than 2 nodbeams (' + filename + ')'
        return, -1
    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, noprint=silent, /info, 'Wrong number of dimensions'
                    return, -1
                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

return, data
end