;#> irsfringe.dc1
; Identifier irsfringe
;
; Purpose remove fringes from IRS data
;
; Synopsis struct = irsfringe( struct [, order = order $
; , modules = modules $
; , bcd = bcd $
; , slitpos = slitpos $
; , fmask = fmask $
; , wrange = wrange $
; , lowpass = lowpass $
; , suppress = suppress $
; , supwid = supwid $
; , status = status $
; , plot = plot $
; , /noplot $
; , /clearplots $
; MODFLAT specifics (when FLAT is provided)
; , /modflat
; , nbin = nbin $
; , edge = edge $
; , method = method $
; , noshift = noshift $
; , nosmooth = nosmooth $
; , fitorder = fitorder $
; , /sinefit $
; SINE defringing specifics
; , cycle = cycle $
; , ncycle = ncycle $
; , cystep = cystep $
; , /auto $
; , nfringes = nfringes $
; , tol = tol $
; , /weight $
; , fringelist = fringelist $
; , /apply_fringelist
; , evidence = evidence $
; , /gui ] )
;
;
; Arguments Name I/O Type: Description:
; ------------------------------------------------------------
; struct I struct IRS 1-D spectra
;
; data selection
; ---------------
; order I intarr list of IRS orders to be processed
; -1 : combine all orders in a module
; -(list) : combine selected orders
; NB in V 1.x this does not yet
; give optimum results
; see Comments
; [ default: all separate ]
;
; modules I intarr data from IRS modules to process
; 0 ... 5 --> SL-ap1, SL-ap2,
; LL-ap1, LL-ap2, SH, LH
; NB. - it is possible to process
; SL and LL data. However we do
; not yet know if it makes sense
; to do so,
; - for LH the defaults have not been
; verified for lack of testdata
; Any feedback is useful
; [ default: all ]
;
; bcd I intarr bcd's
; -1 : combine all bcd's
; -(list) : combine selected bcd's
; [ default: separate ]
;
; slitpos I intarr slit positions ( 0 .. 30 )
; -1 : combine all slit positions
; -(list) : combine selected sl.pos
; [ default: separate ]
;
; fmask I long select only data with
; ( flag AND fmask ) EQ 0
; [ default: '40000000'XL (== 2L^30)
; == BAD data ]
;
; wrange I fltarr[2] limit operations to wavelengths
; within wrange.
; [ default: defined by module and order ]
;
; lowpass I keyword only apply a lowpass filter for the
; fringe isolation. By default a bandpass
; with the high frequency at the
; sampling frequency (or 1/4 of the
; lower cutoff) is applied.
; [ default: bandpass ]
;
; suppress I fltarr set of wavelengths to disregard during
; continuum and fringe fitting
; (eg. at emission lines where the
; automatic masking fails)
;
; supwid I fltarr width of the suppression
; for each "suppress" wavelength or
; one for all
; [ default : instrumental fwhm ]
; general
; -------
; plot I int -1 : no plotting (== /noplot)
; 0 : plot the data and result
; 1-3: various plotting options
; see irsfringe_modflat and
; irsfringe_sine for details
; [ default: 0 ]
;
; noplot I keyword no plotting (same as plot=-1)
;
; clearplots I keyword remove plot windows on exit
;
; status O int exit status
; 0: successfully processed
; 1: error, input struct is returned
;
; MODFLAT specifics (not yet available, requires improved extraction tools)
; -----------------
; /modflat I -- Use the MODFLAT option, in this
; case no sine defringing is done
; unless /SINEFIT is set.
; nbin I int Number of wavelength regions into which
; to split each order for derivation of
; shift and smooth/enhance parameters.
; Give INTARR(1) if nbin same for all
; orders, or INTARR(nord) for each
; order separately.
; [default: 1]
;
; edge I flt What fraction of wavelength coverage
; of each order (in %) should be masked
; on the edges in shift/smooth/enhance
; determination because of, for example,
; weak signal at order edges or some
; edge anomaly.
; [default: 15% of wavelength range
; per order at each edge]
;
; method I str Method with which to determine shift
; and smooth/enhance factors.
; Possibilities are:
; 'LIN' : minimum chi square for
; linear range of parameters
; 'AMO' : amoeba minimization method,
; which has an elegant way of
; searching parameter space
; and is probably also faster
; than 'LIN'.
; Give STRARR(1) if method same for all
; orders, or STRARR(nord) for each order
; separately.
; [default: 'AMO' for all orders]
;
; noshift I intarr Do not determine nor apply phase shift
; correction to fringes in flat field
; spectrum. Either give 1 if you don't
; want shifting for any order, or an
; array with 1's and 0's for each
; echelle order.
; [default: 0, determine shift]
;
; nosmooth I intarr Do not determine nor apply amplitude
; smoothing or enhancing correction to
; fringes in flat field spectrum.
; Either give 1 if you don't want
; smoothing for any order, or an array
; with 1's and 0's for each Eachelle order.
; [default: 0, determine smooth factors]
;
; fitorder I int Fit the smooth/enhance and/or shift
; factors with a smooth polynomial as a
; function of order (as plotted)?
; -1: No, just apply the smooth/enhance
; and shift values derived for each
; order without doing polynomial fit.
; 0: Interactively go through the
; various options.
; 1: Use polynomial fit of smooth/enhance
; factors as function of order.
; Do NOT rederive shifts.
; 2: Use polynomial fit of smooth/enhance
; factors as function of order.
; Do rederive shifts with this.
; 3: Use polynomial fit of shifts as
; function of order.
; Do NOT rederive smooth/enhance
; factors.
; 4: Use polynomial fit of shifts as
; function of order.
; Do rederive smooth/enhance factors.
; 5: Use polynomial fits of shifts AND
; smooth/enhance factors as function
; of order.
; ORDER = -1 (combining all orders)
; implies FITORDER=5 unless FITORDER
; is also specified, in which case
; the specified fitorder is used.
; [default: fitorder=-1]
;
; noshift I int Do not determine nor apply phase
; shift correction to fringes in
; flat field spectrum. Either
; give 1 if you don't want shifting
; for any order, or an array with 1's
; and 0's for each echelle order.
; [default: 0, determine shift
; for each order]
;
; nosmooth I int Do not determine nor apply amplitude
; smoothing or enhancing correction
; to fringes in flat field spectrum.
; Either give 1 if you don't want
; smoothing for any order, or an array
; with 1's and 0's for each echelle order.
; [default: 0, determine smooth/enhance
; factors for each order]
;
; /iter I -- If method='LIN' then it may be
; necessary in some cases to iterate
; over shift and smooth/enhance
; determination to get the best solution.
; [default: do not iterate]
;
; filename I string File into which the fit parameters
; will be saved. If the filename is
; without extension it will get the
; default extension .xdr
; [default: modflat_params.xdr]
;
; /sinefit I -- When /MODFLAT is epcified an extra
; SINE defringing run is done
;
;
; SINE defringing specifics
; -------------------------
; cycle I real cycle of the fringe component
; searched for
; [ default SL : TBI
; LL : 650
; SH : [3500,8500]
; LH : [3000,8000] ]
;
; ncycle I int search window in cycles
; [ default SL : TBI
; LL : 600
; SH : 2000
; LH : 2500 ]
;
; cystep I real step in the cycles
; ( only useful when computer
; resources are too limited )
; [ default = 1 ]
;
; /auto I keyword automatic defringing.
;
; nfringes I int maximum number of fringes removed
; Specifying NFRINGES=-1 will
; not give any sine defringing.
; [ default: determined automatically ]
;
; tol I real remove fringes until the reduction
; in chisq is less than tol
; ( 0.01 == 1 percent )
; [ default: 0 == disabled ]
;
; weight I keyword set all weights to 1
; [ default: less weights to outliers ]
;
; fringelist O struct cumulative!! list of extracted fringes
;
; { module int IRS module
; order int order
; bcd int bcd number
; slitpos int slit position
; frin_nr int fringe number
; cycle float fringe cycle
; cosamp float amplitude of cosine component
; sinamp float amplitude of sine component
; chisq float chisq
; chi_red } float accumulative chisq reduction
;
; apply_fringelist I -- apply the input fringelist to the data
;
; evidence O flt(arr) evidence == log( probability)
; (undefined if not auto)
;
; /gui I -- Start a graphical user interface.
; First one 'default' run through
; IRSFRINGE is done after which the
; gui starts.
;
;
; Returns defringed data
;
; Description IRSFRINGE employs two methods to reduce the effect
; of fringes in IRS data. The first method is MODFLAT
; which tries to minimize the fringe residuals by
; matching the science data with the flatfield by
; shifting and smoothing the flatfield.
; The second method fits robust sine functions (which
; is the first order approximation of the formal
; optical interference function).
;
; Both methods can be used independently or together.
;
; Modifying flatfield
; -------------------
;
; Defringes the data using the correlation method, similar
; to the procedure resp_inter for ISO/SWS. The fringe phase
; is a function of source position in the slit, while the
; spectral resolution depends on the filling factor of
; the source in the slit. Straight division of observation
; and flat is thus expected to give quite strong residual
; fringes.
;
; To minimize fringe residuals, MODFLAT uses the isolated
; fringes from the astronomical observation and flat field.
; It shifts the fringe phase (i.e. wavelength scale) and
; matches the flat field spectral resolution to that of
; the observation by smoothing or artificially enhancing
; the fringe amplitudes. The modified flat field fringe
; spectrum is then multiplied back into the flat field
; 'continuum'.
;
; The astronomical spectrum is subsequently divided
; through the final modified flat field, giving a spectrum
; hopefully free of fringes and full of reliable exciting
; astronomical features.
;
; Resolved spectral features like emission lines or glitches
; are automatically detected and masked and subsequently given
; a weight zero. When there are still points that should
; not be taken along in the calculations SUPPRESS and SUPWID
; can be used to exclude these.
;
;
; Robust SINE fitting
; -------------------
;
; The sine fitting used in IRSFRINGE is the first order
; approximation of the full optical interference model.
;
; T^2 I
; It = ------- * -----------------------------
; (1-R)^2 1+(4R/(1-R)^2) * sin^2(wD!PI)
;
; where T is the transmission, R the reflectance,
; w the wavelength in wavenumbers and D the optical
; thickness of the interfering layer. For small amplitudes
; of the reflectance R this approximates to:
;
; It = A * cos(2!PIwD+a)
;
; which can be rewritten in the computationally more handy
; form used in irsfringe:
;
; It = A cos(2!PIwD) + B sin(2!PIwD)
;
; Fringes are detected in wavenumber (1/micron) space, where
; they have a constant `wavelength', in cycles per wavenumber
; (C/Wn [==micron]).
;
; To find fringes we fit a cosine function to the data
; f(n) = alfa.cos( p.n + phi ) == A.cos(p.n) + B.sin(p.n)
; A, B and p are parameters to be estimated. p the 'wavelength'
; of the fringe.
; For SH the main fringe is located at app. 3500
; cycle/wavenumber and for LH at app. 3000. A second
; component at app. 8000 is sometimes seen in very good
; S/N data.
; For a given p f(n) is a linear function which can easily
; be fitted to the data. So we search all wavelengths for
; e.g. SH from CYCLE (=2500) to CYCLE+NCYCLE*CYSTEP (=4500),
; fit A and B and find the minimum in chi-sq.
;
; In automatic mode we calculate the Bayesian evidence including
; every new fringe. We keep the fringe if this evidence is
; larger than the evidence carried by the previous set of
; fringes. The process of identifying new fringes continues
; as long as the evidence increases with each new fringe. For
; all fringes the amplitudes are determined and the total
; fringe-set is divided out.
;
; In non-automatic mode fringes are extracted until the reduction
; in chi-sq is less than TOL or the number of identified
; fringes is larger than NFRINGES. If the chi-sq increases
; the proccess stops even if the number of found fringes
; is less than NFRINGES.
;
; Before the fitting is started the fringes are isolated
; from the data using a bandpass filter for the main component
; and a lowpass filter for the high frequency component.
; It is possible to use a lowpass filter in all cases by
; using the /LOWPASS option.
;
; Resolved spectral features like emission lines or glitches
; are automatically detected and masked and subsequently given
; a weight zero. Large unmasked outliers in the background
; subtracted data are given less weight unless /WEIGHT is set.
; When there are still points that should not be taken along
; in the calculations SUPPRESS and SUPWID can be used to
; exclude these.
;
; Depending on the value of PLOT several sets of plots pop up.
; -1. (or /NOPLOT) no plotting.
; 0. The result of the fringe removal is plotted. Original data
; in white, smoothed background in red. The defringed data
; are displayed in green, offset somewhat for clarity.
; In the lower panel the background-removed data are plotted
; in white (wgt=1), red (wgt=0) or yellow (0<wgt<1).
; Overplotted in blue are the fringe complex which is found.
; 1. Three panels with the background subtracted data in white/
; yellow/red depending on its weights. The total of the
; fringe complex is overplotted in blue. All these are
; plotted against wavenumber.
; The top panel repeats the next plot, with all extracted
; fringes in green
; 2. The parameters A (white) and B (yellow) as a function
; of p (cycles per wavenumber).
; Chi-sq (scaled and red) is also plotted as function of p. In
; green is the wavelength of the selected fringe.
; Also a preliminary version of 1 is popping up.
; 3. Several plots of the background filtering algorithm
; (for debugging only)
;
; See also the Examples.
;
; Comment See:
; Kester et al. 2001 "SWS Fringes and Models"
; in ISO Legacy Conference Proceedings.
; Fred Lahuis and Adwin Boogert. 2002
; "How to get rid of Fringes in SIRTF-IRS data"
; in "Chemistry as a Diagnostic of Star Formation"
; August 2002, University of Waterloo, Canada
;
; For 3 reasons the robust SINE fitting is superior to
; the use of FFTs.
; 1. The wavelength of the fringe can be estimated much better
; than in an FFT where it is limited to the Nyquist frequency,
; which in turn is determined by the number of datapoints.
; 2. It works on unevenly spaced data just as well as on a
; regularly gridded data set. Even when there are large gaps
; in the data eg. as the result of line masking.
; 3. It automatically determines how many fringe components are
; present, avoiding possible overfitting and underfitting.
;
; The fringe pattern is extremely sensitive to small changes in
; the wavelength, changes so small that they may well be within
; the calibration accuracy of IRS. This is one of the concerns
; and studies to be undertaken once more in-flight data becomes
; available.
;
; Future improvements
; -------------------
; - incorporate knowledge of the detector dimensions
; and optical properties into the fringe model
;
;
; Example To defringe an IRS spectrum:
; irs = irsfringe( irs, order=11 ) ; defringe order 11
; irs = irsfringe( irs ) ; defringe all orders
;
; Use the gui
; irs = irsfringe( irs, /gui )
;
; Category UTIL
;
; Filename irsfringe.pro
;
; Author Fred Lahuis
;
; Version 3.3
;
; History 1.0 01-May-2002 FL, Beta release
; 1.1 12-Jun-2002 FL, properly running selection loops
; added gui call
; 1.1.1 03-Jul-2002 FL, improved continuum fitting using
; frin_filter
; 2.0 02-Oct-2002 FL, incorporated MODFLAT and extracted
; the SINE wave fitting into a
; separate function
; 2.1 10-Oct-2002 FL, follow last modifications to modflat
; at the end return stdev for each order
; 2.2 19-Jun-2003 FL, midcycle parameter
; 2.3 06-Aug-2003 FL, bandpass used in frin_filter
; /lowpass to suppress
; 2.4 07-Aug-2003 FL, redefined cycle,ncycle definition
; 2.5 11-Aug-2003 FL, make sure internal IRS struct is complete
; 2.5.1 18-Aug-2003 FL, fringe_test id added
; 3.0 01-Sep-2003 FL, launch delivery version
; 3.1 29-Oct-2003 FL, proper wrange, updated defaults
; 3.2 11-Dec-2003 FL, second component, apply input fringelist
; 3.3 18-Feb-2004 FL, SSC release updates
; Mar-2004 FL, updates
;#<