; SIRS_READAT - Read Spitzer IRS data from a pipeline reduced file.
;               result = sirs_readat(filnam, /nocrop)
;
; SIRS_COMPRE - Compress the overlapping orders of an IRS spectrum.
;               result = sirs_compre(data)


;+
; NAME:
;       SIRS_READAT
;
; PURPOSE:
;       Read Spitzer IRS data from a pipeline reduced file.
;
; CATEGORY:
;       Spitzer IRS data.
;
; CALLING SEQUENCE:
;       result = sirs_readat(filnam [, /nocrop])
;
; REQUIRED INPUTS:
;       filnam      name of the file to read the data from
;
; OPTIONAL INPUTS:
;       None.
;
; KEYWORDS:
;       nocrop      select all spectral elements, not only the good ones
;
; OUTPUTS:
;       result      IDL structure containing the spectrum
;
; DESCRIPTION AND EXAMPLE:
;       This function reads a file with reduced Spitzer IRS spectra into an 
;       IDL structure. The input file must be the final product of the Spitzer 
;       IRS data reduction pipeline, e.g. a *bksub.tbl file. The resulting 
;       structure contains all fields present in the table file. A fiew fields 
;       are renamed to be consistent with the corresponding fields in other 
;       structures. If the keyword NOCROP is set, all spectral elements will 
;       be read in and not only those which are flagged as good frames (column 
;       'bit_flag' in the file). The function is used for example by:
;            files = file_search('NGC1234', 'SPITZER*bksub.tbl')
;            data  = [sirs_readat(files[0]), sirs_readat(files[1])]
;
; CALLED BY:
;       None.
;
; CALLING:
;       None.
;
; MODIFICATION HISTORY:
;       2008-10-04 Written by Konrad R. Tristram
;       2009-07-02 Konrad R. Tristram: Added description of to the function.
;

FUNCTION SIRS_READAT, FILNAM, NOCROP=NOCROP

; CHECK IF FILE EXISTS, OTHERWISE RETURN TO CALLER
;-------------------------------------------------------------------------------
if not file_test(filnam) then begin
	print, 'The file provided does not exist.'
	return, -1
endif

; INITIALISE A FEW PARAMETERS
;-------------------------------------------------------------------------------
lskipp = 0
header = 0
numdat = 0
datlin = '\'

; OPEN THE IRS TABLE FILE TO BE READ AND SKIP THE HEADER
;-------------------------------------------------------------------------------
openr, filuni, filnam, /get_lun
while not eof(filuni) and strmid(datlin, 0, 1) eq '\' do readf, filuni, datlin

; NOW THE COLUMN NAMES SHOULD BE NEXT IN THE FILE; IF NOT THEN QUIT
;-------------------------------------------------------------------------------
if not (strmid(datlin, 0, 1) eq '|') then begin
	print, 'The file provided does not seem to contain an IRS table.'
	return, -1
endif

; GET THE NAMES AND THE NUMBER OF COLUMNS AS IN THE FILE
;-------------------------------------------------------------------------------
headli = strtrim(strsplit(datlin, '|', /extract))
colnum = n_elements(headli)

; REPLACE A FEW COLUMN NAMES BY STANDARDISED ONES AND MAKE THEM IDL COMPATIBLE
;-------------------------------------------------------------------------------
colpos = where(headli eq 'wavelength')
if colpos ge 0 then headli[colpos] = 'wavele'
colpos = where(headli eq 'flux')
if colpos ge 0 then headli[colpos] = 'flxval'
colpos = where(headli eq 'error')
if colpos ge 0 then headli[colpos] = 'flxerr'
headli = idl_validname(headli, /convert_all)

; READ THE NEXT LINE AND DETERMINE THE DATA TYPES AND POSITIONS OF THE COLUMNS
;-------------------------------------------------------------------------------
readf, filuni, datlin
typstr = strtrim(strsplit(datlin, '|', /extract))

; CREATE A STRUCT WITH THE CORRECT DATA TYPES FOR THE COLUMNS
;-------------------------------------------------------------------------------
exestr = 'struct = create_struct(headli'
for i=0,colnum-1 do begin
	case typstr[i] of
		'int':	exestr = exestr + ', 0S'
		'real':	exestr = exestr + ', 0.0'
		else:	exestr = exestr + ', ""'
	endcase
endfor
exestr = exestr + ')'
noerrs = execute(exestr)
result = struct

; LOOP OVER THE REMAINING LINES IN THE FILE UNTIL END OF FILE IS REACHED
;-------------------------------------------------------------------------------
while not eof(filuni) do begin
	; READ THE CURRENT LINE
	;-----------------------------------------------------------------------
	datlin = ''
	readf, filuni, datlin
	; CHECK IF THE CURRENT LINE CONTAINS ANY FURTHER COMMENTS OR HEADERS
	;-----------------------------------------------------------------------
	if total(strmid(datlin, 0, 1) eq ['|', '\']) eq 0 then begin
		; ADD AN ENTRY TO THE RESULT STRUCT
		;---------------------------------------------------------------
		result = [result, struct]
		; PLIT THE STRING AND THEN COPY IT INTO THE RESULT STRUCT
		;---------------------------------------------------------------
		datlin = strtrim(strsplit(datlin, /extract))
		for i=0,colnum-1 do result[numdat].(i) = datlin[i]
		; INCREASE THE COUNTER FOR THE NUMBER OF DATA SETS
		;---------------------------------------------------------------
		numdat += 1
	endif
endwhile

; CLOSE THE FILE AND FREE THE LUN
;-------------------------------------------------------------------------------
close, filuni
free_lun, filuni

; CHECK IF ANY LINES WERE READ
;-------------------------------------------------------------------------------
if numdat eq 0 then begin
	print, 'There were no lines with data in the table.'
	return, -1
endif

; CROP THE LAST ENTRY
;-------------------------------------------------------------------------------
result = result[0:numdat-1]

; SELECT ONLY LINES WHICH ARE NOT FLAGGED
;-------------------------------------------------------------------------------
if (~ keyword_set(nocrop)) and (total(headli eq 'bit_flag') eq 1) then begin
	select = where(result.bit_flag eq 0)
	result = result[select]
endif

; RETURN THE RESULT
;-------------------------------------------------------------------------------
return, result

END


;+
; NAME:
;       SIRS_COMPRE
;
; PURPOSE:
;       Compress the overlapping orders of an IRS spectrum into one spectrum.
;
; CATEGORY:
;       Spitzer IRS data.
;
; CALLING SEQUENCE:
;       result = sirs_compre(data)
;
; REQUIRED INPUTS:
;       data        structure containing the raw spectrum with overlaps
;
; OPTIONAL INPUTS:
;       None.
;
; KEYWORDS:
;       None.
;
; OUTPUTS:
;       result      the spectrum without overlaps in wavelength as a structure
;
; DESCRIPTION AND EXAMPLE:
;       This function removes all overlaps in wavelengths for a spectrum 
;       composed of several overlapping segments. For each overlapping 
;       portion, the middle elements of the two segments are determined 
;       and a new spectrum is constructed by joining the segments at this 
;       position. For the programme to work correctly, the segments must be 
;       sorted with increasing wavelengths and in every segment the elements 
;       must be ordered with increasing wavelengths. The function is thought 
;       to be used in conjunction with SIRS_READAT:
;            files = file_search('NGC1234', 'SPITZER*bksub.tbl')
;            data  = [sirs_readat(files[0]), sirs_readat(files[1])]
;            data  = sirs_compre(data)
;
; CALLED BY:
;       None.
;
; CALLING:
;       None.
;
; MODIFICATION HISTORY:
;       2008-10-04 Written by Konrad R. Tristram
;       2009-07-02 Konrad R. Tristram: Added description to the function.
;

FUNCTION SIRS_COMPRE, DATA

; GET THE NUMBER OF ELEMENTS IN THE SPECTRUM AND INITIALISE A FEW PARAMETERS
;-------------------------------------------------------------------------------
numele = n_elements(data)
k = 1L
j = 1L

; LOOP OVER ALL ELEMENTS OF THE SPECTRUM
;-------------------------------------------------------------------------------
for i=0,numele-2 do begin
	; CHECK IF THERE IS A JUMP BACKWARDS IN WAVELENGTH OR THE END OF DATA
	;-----------------------------------------------------------------------
	if (data[i].wavele gt data[i+1].wavele) or (i eq (numele-2)) then begin
		; FIND THE RANGE OF VALUES TO RETAIN
		;---------------------------------------------------------------
		differ = min(data[j:i].wavele-data[i+1].wavele, idxtop, /absol)
		differ = min(data[j:i].wavele-data[j-1].wavele, idxlow, /absol)
		idxtop = (idxtop + i + j) / 2
		idxlow = idxlow / 2 + j
		; COPY THE DATA RANGE INTO THE RESULT STRUCT
		;---------------------------------------------------------------
		if j eq 1 then result = data[idxlow:idxtop] $
		          else result = [result, data[idxlow:idxtop]]
		j = i + 1
	endif
endfor

; RETURN THE RESULT
;-------------------------------------------------------------------------------
return, result

END