Viewing contents of file '../idllib/iuedac/iuelib/pro/coadd.pro'
;******************************************************************************
;+
;*NAME:
;
; COADD (RDAF General Production Library) July, 1986
;
;*CLASS:
;
; correction
;
;*CATEGORY:
;
; IUESIPS NEWSIPS
;
;*PURPOSE:
;
; To coadd IUESIPS RDAF .SAV formatted files OR NEWSIPS fits formatted
; files. For IUESIPS data, a .SAV file is generated as output containing
; the coadded results. For NEWSIPS data, a fits file is generated as
; output containing the coadded results. See FILES USED: category below.
;
;*CALLING SEQUENCE:
;
; COADD,SPEC,ELIM,WOPT,CNAME,/newsips
;
;*PARAMETERS:
;
; SPEC (REQ) (I) (0 1) (I S)
; Number of spectra to be coadded, 0 = user is prompted for
; input. Can also be an array of IUESIPS .SAV filenames OR
; NEWSIPS fits filenames, or an argument for FINDFILE.
;
; ELIM (REQ) (I) (0) (I)
; If IUESIPS data:
; Cutoff epsilon value for weighting purposes. If ELIM < 0
; only points with epsilons > ELIM will be considered good
; points.
; < 0 = largest EPS value to be treated as bad data
; 0 = user is prompted for input
; > 0 = all points are considered good (none excluded)
;
; If NEWSIPS data:
; Nu flag value(s) for weighting purposes. Multiple nu flags
; may be set. If ELIM < 0 points with the absolute value nu
; flags bit settings that match the absolute value of ELIM's
; bit settings are considered bad points.
; < 0 = absolute value bit settings used to determine bad
; points
; 0 = user is prompted for input
; > 0 = all points are considered good (none excluded)
;
; WOPT (REQ) (I) (0) (I)
; Exposure weighting option
; 0 = prompt user for input
; 1 = use SQRT(exposure time) for weighting
; 2 = use exposure time directly for weighting
; 3 = ignore exposure times
; 4 = user will specify a constant for weighting
; 5 = use sigmas - NEWSIPS only (can not be used when
; interpolation is necessary)
;
; CNAME (REQ) (I) (1) (S)
; Name of output file; may include path information.
; (user is prompted for name if null string or blanks are
; specified).
;
; NEWSIPS (KEY) (I) (0) (I)
; If the number of spectra to be coadded is given, the files
; involved will be assumed to be IUESIPS RDAF .SAV files unless
; the NEWSIPS flag is set to indicate that the files are NEWSIPS
; fits files.
;
;*EXAMPLES:
;
; coadd,0,0,0,' ' ; For IUESIPS data. Prompt user for all
; ; necessary input.
;
; coadd,0,0,0,' ',/newsips ; For NEWSIPS data. Prompt user for all
; ; necessary input.
;
; coadd,2,-350,2,'favg' ; For IUESIPS data. Coadd 2 spectra, delete
; ; points with EPS values < or equal to -350.
; ; Weight by exposure times and save results
; ; in file favg.sav
;
; coadd,2,-4160,2,'favg',/newsips
; ; For NEWSIPS data. Coadd 2 spectra, delete
; ; points with nu flags values of -4096 and
; ; -64 (thus -4160). Weight by exposure
; ; times and save results in file favg.fit
;
; coadd,5,1,1,'favg5' ; For IUESIPS data. Coadd 5 spectra, do not
; ; exclude any points, weight by the square
; ; root of the exposure time and store
; ; results as favg5.sav
;
; coadd,5,1,1,'favg5',/newsips
; ; For NEWSIPS data. Coadd 5 spectra, do not
; ; exclude any points, weight by the square
; ; root of the exposure time and store
; ; results as favg5.fit
;
; file = findfile('swp*.sav')
; COADD,file,1,3,'swpcoadd'
; ; For IUESIPS data. Create an array of .sav
; ; file names. Coadd all SWP .sav files in
; ; the user's current directory with equal
; ; weighting and including all points.
; ; Results will be saved in swpcoadd.sav
;
; file=findfile('swp*.mxlo')
; COADD,file,1,3,'swpcoadd',/newsips
; ; For NEWSIPS data. Create an array of
; ; .mxlo file names. Coadd all SWP .mxlo
; ; files in the user's current directory with
; ; equal weighting and including all points.
; ; Results will be saved in swpcoadd.fit
;
; COADD,'swp*.sav',1,3,'add' ; For IUESIPS data. Find all files of the
; ; form swp*.sav and coadd them with equal
; ; weighting and including all points.
; ; Results will be saved in add.sav
;
; COADD,'swp*.mxlo',1,3,'add',/newsips
; ; For NEWSIPS data. Find all files of the
; ; form swp*.mxlo and coadd them with equal
; ; weighting and including all points.
; ; Results will be saved in add.fit
;
;*SYSTEM VARIABLES USED:
;
; !iuer.dat
; !iuer.inf
;
;*INTERACTIVE INPUT:
;
; Besides above input prompting, user is prompted for names of IUESIPS
; RDAF .sav OR NEWSIPS fits files to be coadded, unless the user has
; entered an array of IUESIPS RDAF .sav OR NEWSIPS fits filenames to use.
; The filenames may also include path specifications. For IUESIPS data,
; wavelengths from the first file entered will be used as a reference for
; resampling all subsequent spectra.
;
;*SUBROUTINES CALLED:
;
; PARCHECK
; PCHECK
; CHKFITS
; IUETERP
; READMX
; IUESAVE
; IUEFETCH
; INTIME
; NUFLAGS
; SHOWTXT
;
;*FILES USED:
;
; !iuer.inf EPSILON.DOC (I)
; File describing IUESIPS data quality flag assignments (displayed if
; ELIM = 0).
;
; !iuer.inf NUFLAGS.DOC (I)
; File describing IUE NEWSIPS data quality flag assignments
; (displayed if ELIM = 0).
;
; CNAME.SAV or CNAME.FIT (O)
; RDAF SAV or fits file of coadded spectrum generated from the
; resulting header, wavelength, average weighted flux, and weight
; arrays described above. The E vector describes the fraction of
; points used at each wavelength in the coadded spectrum. If no good
; points were found at a particular wavelength, the E vector entry is
; set to -1000.
;
;*SIDE EFFECTS:
;
;*RESTRICTIONS:
;
; Weighting option 5 (NEWSIPS only) can not be used if the vectors must
; be interpolated.
;
;*NOTES:
;
; COADD assumes that the user has corrected for wavelengths shifts
; between spectra (e.g. see CRSCOR.PRO).
;
; If ELIM = 0, the files are IUESIPS RDAF .sav files, or the NEWSIPS
; flag is not set, the disk file IUER_INF EPSILON.DOC or
; IUER_INF NUFLAGS.DOC is displayed.
;
;*PROCEDURE:
;
; COADD begins by prompting the user to input each of the first 3
; parameters which were specified as 0 in the procedure call. Once these
; options are set, the user specifies the name of each of the IUESIPS
; RDAF sav or fits files to be coadded (this step is skipped if the user
; entered a string or an array of filenames).
;
; The wavelengths in the first file named will be used as a reference to
; which all subsequent spectra will be compared. If the wavelength
; sampling for any succeeding file differs from that of the first
; spectrum, IUETERP will be called to resample the wave, flux, and
; IUESIPS eps or NEWSIPS nu flags vectors.
;
; Depending on the weighting option specified, an array WGT is generated
; in which 0s designate bad points and 1s designate good points. The sum
; of the WGT arrays (called NGOOD),is used to describe the fraction of
; good data points used at each wavelength of the coadded spectrum. If
; no points are used at a particular wavelength, the entry in NGOOD is
; set to -1000, and the flux value is set to zero. The WGT array is then
; scaled according to the users choice of WOPT and the calculated
; exposure times (if so requested). The average weighted flux is
; calculated by dividing the sum of the weighted fluxes (i.e. FLUX*WGT),
; by the sum of the weights. A plot is displayed of average weighted
; flux versus wavelength with X's drawn at points where no good points
; were found.
;
; The final vectors H,W,FAVG, & NGOOD are stored in a RDAF .sav file or
; a fits file.
;
;*I_HELP nn
;
;*MODIFICATION HISTORY:
;
; Written by J. K. Feggans for Dr. Bergstahl 7/85.
; Updated for Dr. J. Neff 10/85.
; Updated for use by Dr. Nousek 11/1/85
; Documentation C. Grady 12/85
; Further revisions J.K. Feggans and C. Grady 1/86
; 2-28-86 RWT remove EOPT parameter, update prolog, and improved
; parameter definitions
; Optimized JKF 3/86.
; 3-10-86 RWT use NGOOD to generate final array of weights
; 4-15-86 RWT prompt user for ELIM & WOPT before SAV file names
; 5-13-86 RWT add PCHECK for CNAME and comments about quotes
; 3-14-87 RWT VAX mods: add PARCHECK, & use N_ELEMENTS
; 10- 7-87 RWT add procedure call listing, replace IUESTOR with new
; IUESAVE, and remove EPSWGT & EXPWGT subroutines
; 2-19-88 RWT make GETEXP a separate routine called CO_GETEXP
; 4-18-88 RWT rename COADD_GETEXP
; 5-11-88 HAA add RDAF Prolog
; 1-03-90 RWT UNIX mods: convert to lower case, change logical name
; references, remove routine COADD_GETEXP, change format
; statements
; 4-10-91 KBC spawn appropriate system call based on operating system
; type for compatibility on SUN/DEC/VAX
; 27 Jun 91 LLT add call to GRIDTERP, changed logical; clean up;
; tested on VAX; updated prolog
; 12 Jul 91 PJL corrected error with wopt = 5; tested on SUN and VAX;
; updated prolog
; 21 Nov 91 GRA removed IUER_USERDATA and 16 character name limit
; from IUESAVE, updated prolog; tested.
; 30 Jan 92 LLT GRIDTERP when wavelength vectors are of unequal size,
; even if the existing elements all match; add ability to
; input an array of filenames.
; 4 Feb 92 LLT add call to FINDFILE when a scalar string is given
; (SPEC).
; 22 Jun 93 EER (JPL) to put total exposure time into header.
; 30 Jun 93 PJL changed !version.os if to case; replaced GRIDTERP with
; IUETERP
; 26 Jul 93 PJL compatibility with NEWSIPS - including READMX and
; NUFLAGS procedures; NEWSIPS keyword
; 5 Aug 93 PJL use IUETERP on NEWSIPS data when necessary - except for
; weighting option 5
; 14 Mar 94 LLT if newsips keyword not set and no extension in filename
; assume .SAV
; 15 Jun 94 RWT allow respecifying file name, add print message so
; "files not found" is distinguished from "incorrect
; file type", correct telescop keyword value, store
; NEWSIPS wavelengths as single-precision rather than
; double-precision, allow nu flag description to be
; read before screen is erased, and allow cname to include
; extension.
; 23 Jun 94 PJL replaced !version with PLATFORM
; 18 Oct 94 PJL added code for cat command equals 'NA'
; 21 Oct 94 PJL replace environment variables with system structure
; 24 Jan 95 RWT replace spawn with call to showtxt routine
;-
;******************************************************************************
pro coadd,spec,elim,wopt,cname,newsips=newsips
;
npar = n_params(0)
if (npar eq 0) then begin
print,'COADD,SPEC,ELIM,WOPT,CNAME,/newsips'
retall
endif ; npar eq 0
parcheck,npar,4,'COADD'
;pcheck,nspec,1,110,1010 ; integer scalar or string vector
pcheck,elim ,2,100,0010 ; scalar
pcheck,wopt ,3,100,0010 ; scalar
pcheck,cname,4,100,1000 ; string
wsum = 0. ; total weight summary
fsum = 0. ; total weighted fluxes
ngood = 0. ; # of good points
dashes = string(bytarr(79) + 45b) ; dashed line
ftype = ['RDAF .SAV','fits']
;
; determine input files
;
snspec = size(spec) ; Is spec a scalar number or an array of filenames?
type = snspec(n_elements(snspec)-2) ; Type code for spec - 7 if string
if (type eq 7) then begin
nspec = n_elements(spec) ; Number of files to be coadded.
;
; If only one file, or a scalar string call FINDFILE.
;
if (nspec eq 1) then begin
spec = findfile(spec)
;
; Number of files found by FINDFILE to be coadded.
;
nspec = n_elements(spec)
endif ; nspec eq 1
if (nspec le 1) then begin ; If less than two were found, retall.
if (nspec eq 1) then print,'Only one file was found: ',spec else $
print,'No files were found.'
print,'ACTION: retall'
retall
endif ; nspec le 1
;
; determine if IUESIPS (newsips = 0) or NEWSIPS (newsips = 1)
;
if (keyword_set(newsips)) then newsips = 1 else newsips = 0
sipschk = intarr(nspec)
for i = 0,nspec-1 do begin
chkfits,spec(i),res,/silent
sipschk(i) = res
endfor ; i
;
; If the NEWSIPS keyword has been set, then all the files must be fits files.
; If the NEWSIPS keyword has not been set, then all the files must be fits
; or all the files must be RDAF .sav.
;
if (newsips) then begin
temp = where(sipschk ne newsips,ct)
if (ct gt 0) then begin
print,' '
print,'NEWSIPS flag set.'
print,strtrim(ct,2) + ' of the ' + strtrim(nspec,2) + $
' files are not fits files.'
print,'ACTION: retall'
retall
endif ; ct gt 0
newsips = sipschk(0)
endif else begin
temp = where(sipschk ne sipschk(0),ct)
if (ct gt 0) then begin
print,' '
print,'Of the ' + strtrim(nspec,2) + ' files, ' + strtrim(ct,2) + $
' are ' + ftype(sipschk(0)) + ' files and'
print,strtrim(nspec-ct,2) + ' are ' + ftype(sipschk(temp(0))) + $
' files.'
print,'ACTION: retall'
retall
endif ; ct gt 0
newsips = sipschk(0)
endelse ; newsips
;
print,'The following ',strtrim(string(nspec),2),' files will be coadded:'
print,spec
if (elim eq 0) then begin
print,' '
temp = ''
read,'Please press RETURN to continue.',temp
endif ; elim eq 0
endif else begin
nspec = spec
;
; number of spectra to be coadded option.
;
if (nspec le 0) then read,'Enter number of spectra to be coadded: ',nspec
if (nspec le 0) then begin ; abort for invalid number of spectra
print,'Invalid coadd parameter...number of spectra = ',nspec
print,'ACTION: Returning
retall
endif ; end of abort loop - nspec le 0
;
; determine if IUESIPS (newsips = 0) or NEWSIPS (newsips = 1)
;
if (keyword_set(newsips)) then newsips = 1 else begin
newsips = 0
print,' '
print,'Assuming IUESIPS RDAF .SAV files as input.'
print,'For NEWSIPS files, please use the calling sequence:'
print,' coadd,spec,elim,wopt,cname,/newsips'
print,' '
if (elim eq 0) then begin
print,' '
temp = ''
read,'Please press RETURN to continue.',temp
endif ; elim eq 0
endelse ; keyword_set(newsips)
endelse ; type eq 7
;
; select epsilon or nu flag weighting and exposure time weighting
;
if (elim eq 0) then begin ; let user specify elim
print,' '
print,dashes
;
if (newsips) then begin
print,' IUE data quality flags for selecting nu flag weighting:'
da = !iuer.inf + 'nuflags.doc'
endif else begin
print,' IUE data quality flags for selecting epsilon weighting:'
da = !iuer.inf + 'epsilon.doc'
endelse ; newsips
;
; display text file
;
showtxt,da
;
print,dashes
;
if (newsips) then begin
temp = ''
read,'Please press RETURN to continue.',temp
nuflags,elim,'Please select nu flags to be considered bad points'
endif else begin
print,'Please enter the largest epsilon value to be considered ' + $
'a bad point'
read,' (Enter a value > 0 to include all data points) ',elim
endelse ; newsips
endif ; elim eq 0
;
if (wopt eq 0) then begin
print,' '
print,dashes
print,'Exposure time weighting options:'
print,' enter 1 to use square root of exposure times.'
print,' enter 2 to use exposure times directly.'
print,' enter 3 to ignore exposure times...set to unity.'
print,' enter 4 to specify constant factor.'
if (newsips) then print,' enter 5 to use sigmas.'
print,' '
print,' enter 0 to abort'
print,' '
print,dashes
read,' ',wopt
endif ; wopt eq 0
;
if (newsips) then begin
if ((wopt le 0) or (wopt gt 5)) then begin
print,' Exposure time option specified...option = ',wopt
print,' ACTION: Returning
retall
endif ; ((wopt le 0) or (wopt gt 5))
endif else begin
if ((wopt le 0) or (wopt ge 5)) then begin
print,' Exposure time option specified...option = ',wopt
if (wopt eq 5) then print,'Sigma option only for use with NEWSIPS data.'
print,' ACTION: Returning
retall
endif ; ((wopt le 0) or (wopt ge 5))
endelse ; newsips
;
print,' '
;
; fetch each save file, weight according to epsilons and exposure times.
;
if (not(newsips)) then begin
print,'WARNING:'
print,'Unless the wavelength sampling is already the same for each ' + $
'spectrum,'
print,'it will be interpolated to match that of the FIRST file given.'
print,' '
endif ; not(newsips)
print,' '
;
net_exps = 0. ; net exposure time for coadd (EER)
; gives zero for wopt > 2
for i = 0,nspec-1 do begin
;
; user must input each file name
;
if (type ne 7) then begin
print,' '
print,format='(1x,a,i3)','Enter ' + ftype(newsips) + $
' file name for entry: ',i+1
savfil = ' ' ; .sav string file name
read,'(without quotes) ',savfil ; read name
if not newsips then begin
decompose,savfil,sdisk,spath,sname,sext,svers
if sext eq '' then savfil=sdisk+spath+sname+'.sav'+svers
endif ;not newsips
dum = findfile(savfil,count=cnt)
if (cnt eq 0) then begin
print,'specified file: ',savfil,' not found'
read,'please respecify full file name: ',savfil ; read one more time
endif
chkfits,savfil,res,/silent
if (res eq 0) and (newsips eq 1) then begin
print,'NEWSIPS flag is set but '+savfil+' is an RDAF .SAV file.'
print,'ACTION: retall'
retall
endif
if (res eq 1) and (newsips eq 0) then begin
print,'IUESIPS file assumed, but '+savfil+' is a fits file.'
print,'To specify NEWSIPS, please set the NEWSIPS keyword'
print,' (i.e., coadd,spec,elim,wopt,cname,/newsips)'
print,'ACTION: retall'
retall
endif
if (res eq 2) then begin
print,savfil,' still not found'
print,'ACTION: retall'
retall
endif
endif else savfil = spec(i)
;
; read the file
;
print,' '
print,'Reading file: ' + savfil
if (newsips) then readmx,savfil,h,wave,flux,eps,sigma else $
iuefetch,savfil,h,wave,flux,eps
;
; check the wavelength vector
;
if (i eq 0) then w0 = wave else begin
s = size(where(w0 ne wave) ne -1)
if (s(0) ne 0) or (n_elements(w0) ne n_elements(wave)) then begin
;;; gridterp,w0,wave,flux,eps,a,b
if (newsips) then begin
if (wopt eq 5) then begin
print,' '
print,'Weighting option 5, using the sigmas, can _not_ ' + $
'be used when interpolation is necessary.'
print,'ACTION: retall'
print,' '
retall
endif else iueterp,w0,wave,flux,nu=eps,header=h
endif else iueterp,w0,wave,flux,eps=eps
endif ; (s(0) ne 0) or (n_elements(w0) ne n_elements(wave))
endelse ; i eq 0
;
; update record zero or create/update fits header
;
if (newsips) then begin
if (i eq 0) then begin
bline = string(bytarr(80) + 32b) ; a blank line
fil = !iuer.dat + 'genhd.txt'
a = bline
openr,lu1,fil,/get_lun
hsum = strarr(36)
for j = 0,2 do readf,lu1,a ; skip 1st 3 lines
j = 0
while not eof(lu1) do begin
temp = bline
readf,lu1,a
strput,temp,a
hsum(j) = temp
j = j + 1
endwhile ; not eof(lu1)
free_lun,lu1
delpar,hsum,'COMMENT'
iuedaf,hsum,temp,/override
addpar,hsum,'TELESCOP','IUE ', $
' International Ultraviolet Explorer','IUEDAC'
addpar,hsum,'NSPEC',nspec,' number of spectra COADDed','IUEDAC'
if (elim eq 0) then addpar,hsum,'HISTORY', $
'All points considered good',' ','IUEDAC' else $
addpar,hsum,'HISTORY','Nu flags in ' + strtrim(elim,2) + $
' considered bad points in COADD.','','IUEDAC'
endif ; i eq 0
;
; information about spectra being coadded
;
camname = ''
imagenum = ''
aper = ''
dispersn = ''
exps = ''
;
; camera name
;
stpar,h,'camera',camname,err1
if (err1 ne 0) then begin
print,'No camera name found.'
camname = ''
endif ; err1 ne 0
;
; image number
;
stpar,h,'image',imagenum,err2
if (err2 ne 0) then begin
print,'No image number found.'
imagenum = ''
endif ; err2 ne 0
;
; aperture
;
stpar,h,'extaper',aper,err3
if (err3 ne 0) then begin
stpar,h,'aperture',aper,err4
if (err4 ne 0) then begin
print,' '
print,'Unable to determine the aperture.'
print,'Please select one of the following options:'
print,' '
print,' 1 - large aperture'
print,' 2 - small aperture'
print,' 3 - exit'
print,' '
aperch = ''
read,'Option: ',aperch
;
case aperch of
'1': aper = 'LARGE'
'2': aper = 'SMALL'
'3': retall
else: begin
print,'Invalid option.'
print,'ACTION: Returning.'
retall
end ; else
endcase ; aperch
endif ; err4 ne 0
err3 = 0
endif ; err3 ne 0
;
if (err3 ne 0) then begin
print,'No aperture found.'
aper = ''
endif else begin
if (aper eq 'LARGE') then aper = 'LG'
if (aper eq 'SMALL') then aper = 'SM'
endelse ; err3 ne 0
;
if ( (err1 eq 0) and (err2 eq 0) and (err3 eq 0) ) then $
addpar,hsum,'SPEC' + strtrim(i+1,2) + 'N',strtrim(camname,2) + $
strtrim(imagenum,2) + aper, ' name of spectra ' + strtrim(i+1,2) + $
' COADDed','IUEDAC'
;
; dispersion
;
stpar,h,'dispersn',dispersn,err
if (err ne 0) then begin
print,'No dispersion found.'
disp = ''
endif else addpar,hsum,'SPEC' + strtrim(i+1,2) + 'D',dispersn, $
' dispersion of spectra ' + strtrim(i+1,2) + ' COADDed','IUEDAC'
;
; expsoure time
;
stpar,h,strmid(strtrim(aper,2),0,1)+'EXPTIME',exps,errexp
if (errexp ne 0) then exps = '' else addpar,hsum,'SPEC' + $
strtrim(i+1,2) + 'E',exps,' exposure time (seconds) of spectra ' + $
strtrim(i+1,2) + ' COADDed','IUEDAC'
;
endif else begin
if (i eq 0) then begin
hsum = h ; header vector
hsum([4,70]) = [0,nspec] ; remove isn & store spec. #
endif ; i eq 0
hsum(71+i) = h(4) ; store isn's of spectra used in coadd
endelse ; newsips
;
; set wgt = 0 for bad pts, wgt = 1 for good pts
;
if (newsips) then begin
if (elim le 0) then begin
wgt = (abs(eps) and abs(elim))
tempnz = where(wgt ne 0,ctnz) ; "bad" points
tempz = where(wgt eq 0,ctz) ; "good" points
if (ctnz gt 0) then wgt(tempnz) = 0
if (ctz gt 0) then wgt(tempz) = 1
endif else wgt = eps * 0. + 1.
endif else begin
if (elim le 0) then wgt = (eps gt elim) ; remove bad eps or nu values
if (elim gt 0) then wgt = eps * 0. + 1. ; no epsilon or nu weighting
endelse ; newsips
ngood = ngood + wgt ; total good points
;
; for wopt < 3 or eq 5 (for NEWSIPS), get exposure time out of IUESIPS header
; record and convert to seconds, or use the exposure time from the NEWSIPS
; fits keyword xEXPTIME (where x is L or S). If missing, user may either
; specify time or abort.
;
if ( (wopt lt 3) or (wopt eq 5) ) then begin
if (newsips) then begin
;
; unknown NEWSIPS exposure time
;
if (errexp ne 0) then begin
print,' '
print,'Unable to determine exposure time.'
exps = ''
print,'Please enter exposure time in seconds (0 to exit): ',exps
numornot,exps,ans,1
if (ans) then exps = float(exps) else begin
print,' '
print,'Invalid exposure time.'
print,'ACTION: Retall'
retall
endelse ; ans
if (exps eq 0) then retall
endif ; errexp ne 0
;
if (exps le 0) then begin
print,'Invalid exposure time - less than or equal to 0 seconds.'
print,'ACTION: retall'
retall
endif ; exps le 0
endif else begin
;
; IUESIPS exposure time
;
if (total(h(39:41)) le 0) then begin
print,dashes
print,' Exposure time missing from header'
print,' enter 1 to input exposure time'
print,' enter 2 to abort'
print,dashes
read,' ',option
print,' '
if (option ne 1) then begin
print,' ACTION: Returning'
retall
endif else intime,h ; option ne 1
endif ; exps le 0
exps = (60. * h(39)) + h(40) + (h(41) / 1000.)
endelse ; newsips
net_exps = net_exps + exps ; store total exposure time (EER)
if (wopt eq 1) then exps = sqrt(exps) ; compute sqrt of exp. time
endif ; (wopt lt 3) or (wopt eq 5)
;
; compute weightings
;
if (wopt eq 5) then begin
;
; user selection 5 - NEWSIPS sigmas
;
wgt = 1.0/(sigma*sigma)
wsum = wsum + wgt * exps
fsum = fsum + (wgt * exps * flux)
endif else begin
;
; compute weighting according to user selections 3 or 4
; note weighting directly by exposure time is assumed for selection 2
;
if (wopt eq 3) then exps = 1
if (wopt eq 4) then begin
print,' '
read,' Enter exposure time weighting factor for ' + savfil + $
': ',fact
exps = fact
net_exps = net_exps + exps
endif ; wopt eq 4
;
wgt = wgt * exps ; exposure weighted
wsum = wsum + wgt ; weight sum
fsum = fsum + (wgt * flux) ; sum weighted fluxes
endelse ; wopt eq 5
endfor ; end file loop
;
; compute average weighted flux & array describing number of good data pts
; used. for wavelengths with no valid points set ngood = -1000.
;
favg = fsum / (wsum > 1e-36) ; weighted flux
ngood = (ngood / nspec) + ( -1000 * (ngood eq 0) ) ; weight array
;
; allow user to specify resulting file name...max. 16 char.'s
;
cname = strtrim(cname,2)
if (cname eq '') then $
read,'Enter name for coadded spectrum (no quotes) ',cname
decompose,cname,cdi,cpa,cnam,cext,cver
if (newsips) and (cext eq '') then cname = cdi+cpa+cnam + '.fit'+cver
;
; display weighted flux vs wavelength, overplot x's on bad points.
;
wave = float(wave)
stf = '(1x," Weighted Average for",i3," Spectra")'
st = string(format=stf,fix(nspec))
plot,title=st,ytitle='Weighted Flux',xtitle='Wavelength',wave,favg
bd = where(ngood lt 0,count) ; flag bad points
if (count gt 1) then oplot,psym=7,wave(bd),favg(bd) else $
if (count eq 1) then oplot,psym=7,[wave(bd),wave(bd)],[favg(bd),favg(bd)]
;
if (newsips) then begin
;
; weighting option used
;
case wopt of
1: addpar,hsum,'HISTORY','COADD used SQRT(exposure time) for ' + $
'weighting','','IUEDAC'
2: addpar,hsum,'HISTORY','COADD used exposure time directly for ' + $
'weighting','','IUEDAC'
3: addpar,hsum,'HISTORY','COADD ignored exposure times','','IUEDAC'
4: addpar,hsum,'HISTORY','COADD used user specified constant for ' + $
'weighting','','IUEDAC'
5: addpar,hsum,'HISTORY','COADD used sigmas','','IUEDAC'
else: print,'This should not happen.'
endcase ; wopt
addpar,hsum,'HISTORY','COADD ' + !stime,'','IUEDAC'
;
; update header to reflect added spectrum rather than first spectrum (EER)
; gives zero if wopt=3, sum of weighting factors if wopt=4
;
addpar,hsum,'NEXPTIME',net_exps,' net exposure time in seconds','IUEDAC'
ifitswrt,hsum,wave,favg,ngood,ofn=cname,p1t='WAVELENGTH', $
p1u='ANGSTROM',p2t='FLUX',p2u='ERG/CM2/S/A',p3t='WEIGHT'
print,' '
print,'Output file may be read with IFITSRD'
endif else begin
;
; update header to reflect added spectrum rather than first spectrum (EER)
; gives zero if wopt=3, sum of weighting factors if wopt=4
;
hsum(39) = net_exps/60. ; minutes
hsum(40) = fix(net_exps) - hsum(39)*60 ; seconds
hsum(41) = (net_exps - fix(net_exps))*1000. ; milliseconds
;
; save header, wavelength, average weighted flux, and weight vectors
;
iuesave,cname,hsum,wave,favg,ngood,1
endelse ; newsips
;
return
end ; coadd