Viewing contents of file '../idllib/iuedac/iuelib/pro/crscor.pro'
;**********************************************************************
;+
;*NAME:
;
; CRSCOR JUNE,1984
;
;*CLASS:
;
; Cross-correlation
;
;*CATEGORY:
;
;*PURPOSE:
;
; To derive the normalized cross correlation function for two
; spectra (using either linear wavelengths or velocity units), and
; to calculate its maximum by Gaussian fitting or finite differences.
; The location of the maximum represents the apparent shift of spectrum
; 1 with respect to spectrum 2.
;
;*CALLING SEQUENCE:
;
; CRSCOR,W1,F1,W2,F2,BEGL,ENDL,delta,vinc=vdel,ccf=ccf
;
;*PARAMETERS:
;
; W1,F1 (REQ) (I) (0) (S)
; Required wavelength and flux vectors to be analyzed for
; first spectrum.
;
; W2,F2 (REQ) (I) (0) (S)
; Required wavelength and flux vectors to be analyzed for
; reference spectrum.
; Calculated shifts are with respect to these wavelengths.
;
; BEGL (REQ) (I) (0) (F)
; Beginning wavelength for the cross-correlation, in the
; same units as the files to be correlated (default Angstroms)
;
; ENDL (REQ) (I) (0) (F)
; Ending wavelength for the cross-correlation, in the same
; units as BEGL.
;
; DELTA (OPT) (O) (0) (F)
; Output parameter giving the apparent shift between
; the two spectra in the same units as the wavelength scale.
; Note, if more than one Gaussian component is specified, only
; the first center value is output to DELTA although all
; values are listed on the display.
;
; VINC (KEY) (I) (0) (F)
; If specified, cross correlation is performed in velocity
; space using the specified value for the incremental velocity
; spacing. (3 - 10 km/sec is recommended for IUE high dispersion).
; If not specified, linear wavelengths are used.
;
; CCF (KEY) (O) (1) (F)
; optional keyword for outputting the normalized cross correlation
; function.
;
;*EXAMPLES:
;
; CRSCOR,W1,F1,W2,F2,2680,2760,DEL,vinc=10
; cross correlates order 85 of spectrum 1 with spectrum 2
; using the wavelength range from 2680 to 2760 angstroms, and
; producing the cross correlation function every 10 km/sec.
;
; CRSCOR,W1,F1,W2,F2,2400,2600,del
; correlation function is derived for wavelength region from
; 2400 to 2600 angstroms using linear wavelengths.
; (typical use for IUE low dispersion spectra).
;
; N1 = INDGEN(N_ELEMENTS(F1))
; N2 = INDGEN(N_ELEMENTS(F2))
; CRSCOR,N1,F1,N2,F2,0,N_ELEMENTS(F1)-1,CCF=CROSS
; cross correlates two arrays using indices instead of
; wavelengths, and outputs cross correlation function as cross.
;
;
;*INTERACTIVE INPUT:
;
; The programs first presents the user with a graph of
; the two input files. It then computes the
; normalized cross correlation function and displays 31 pts. over
; the interval +/- 15*VDEL (or +/- 15*delw). The user is prompted
; with the following options:
; 1) if a well-defined maximum is present the
; program can locate it automatically and determine its position
; by means of interpolation to zero in the first differences
; (finding the point of zero slope); 2) if it appears that a
; maximum exists outside the searched range, the range can be
; doubled and the function replotted; 3,4,..) if Gaussian fitting
; is desired, (e.g. if more than one real maximum appears present)
; then the cross-correlation function may be fit interactively
; with 1,2,.. components (as determined by GAUSSFIT). GAUSSFIT
; removes a constant baseline before fitting the Gaussians using
; the region outside the region specified by the user as the
; feature to be fitted. The user indicates with the cursor the
; approximate location and height of each component to be fit.
;
;*SYSTEM VARIABLES USED:
;
; !d.x_ch_size
; !d.y_ch_size
; !d.y_size
;
;*FILES USED:
;
; none
;
;*SUBROUTINES CALLED:
;
; PARCHECK
; IUEFETCH
; XYREADS
; IUETERP - resamples spectrum 1 to agree with spectrum 2 (if required)
; CRSPROD - computes a normalized cross-correlation function
; CRSMAX - determines maximum of cross-correlation function
;
;*PROCEDURE:
;
; While for most IUE applications the range (ENDL - BEGL) is
; much less than the wavelengths involved, it is more precise
; in high dispersion to find the shift in log lambda space instead
; of lambda directly (since dL/L = cnst). If velocity spacing is
; requested, the two arrays of wavelengths and fluxes are
; interpolated using QUADTERP to a common set of equally spaced
; log lambda points, with spacing corresponding to a velocity VDEL
; at the midpoint wavelength (BEGL+ENDL)/2. In either case, the mean
; flux for each spectrum is subtracted from that flux array, and the
; resulting arrays of positive and negative values are shifted with
; respect to each other to find the sum of the products at each shift.
; The array of product sums vs. shift is normalized to become the
; 'normalized cross-correlation function'. If the user wants to
; fit several Gaussians to the cross correlation function,
; GAUSSFITS fits the heights, widths, and central velocities for the
; specified number of components. It also computes the standard
; deviations of these parameters using the value computed from the
; baseline portion of the correlation function as the standard
; deviation of a single point. For correlation functions with
; a single well-defined maximum, the velocity shift can be
; calculated automatically by a weighted least squares fit to the
; 1st differences and interpolating to the point of zero slope.
; The weighting gives higher weight to points close to the midpoint
; of the correlation function.
; The output includes a plot of the cross correlation function,
; and the velocity shift determined from the maximum of the cross
; corrrelation function. If Gaussian fitting is used, the standard
; deviation of a single point is printed as are the values for
; the Gaussian parameters and their errors. A plot shows the data,
; the individual Gaussians, the sum of all components, and the
; deviations.
;
;*NOTES:
;
; - When a Gaussian fit is made to the cross correlation
; function, GAUSSFITS uses the following symbols for output:
; - histogram: data points
; - solid line: the total fit, the sum of all Gaussian
; components and the baseline
; - dots: individual Gaussian components
; - plus signs: O-C deviations from the fit
;
; - Before execution one should check IUE high dispersion
; spectra to be sure that the heliocentric velocity correction
; has been applied.
; For images processed before 10 Nov. 1981, the user should
; determine the orbital velocity corrections (see IUEVEL).
; Corrections may also be needed for velocity-like shifts due to
; IUE camera effects (see ASSESS). In any case, users may want to
; run IUEVEL for the most accurate velocity corrections.
;
; - Smoothing of the data might improve the fits.
;
; - Execution time for CRSCOR depends on the size of
; the arrays. It is recommended that no more than half of the
; wavelength interval covered by 3 IUE high dispersion orders be
; correlated at one time.
;
; - If a reasonable amount of correlation exists between the 2
; spaectra, a well-defined maximum will be present, often
; Gaussian in character but often with asymmetrical wings and
; secondary inflections which may or may not have significance.
; Note: if a spectrum has a strong slope, this can distort the
; function relative to what would be obtained from spectral
; lines with a flat continuum; a procedure such as NORM should
; be used in this case before running crscor.
; The precise maximum of the function can be found either by
; fitting a Gaussian function in the vicinity of the peak or by
; locating the shift value for which the slope of the function
; is zero. It was found that the procedure using the
; first differences of the function, locates this value with
; a precision better than 0.1*VDEL or about 1 km/sec generally.
;
; Program assumes velocity units are km/sec and wavelength units
; are in angstroms.
;
; tested with IDL Version 2.1.0 (sunos sparc) 25 Jun 91
; tested with IDL Version 2.1.0 (ultrix mispel) N/A
; tested with IDL Version 3.0 (vms vax) 09 Jul 93
;
;*MODIFICATION HISTORY:
;
; Sept. 1981 S. Parsons CSC initial program, adapted from
; a FORTRAN program by M. Slovak
; of the University of Texas.
; Apr. 1982 F.H. Schiffer GSFC program brought into agreement
; with RDAF standards
; Jun. 1984 N.R. Evans GSFC to test, document, alter output,
; permit fitting of the correlation
; function by Gaussian(s), and to
; determine the errors for fitting
; parameters.
; 18 Jul. 1984 RWT GSFC replaced subroutine PLPARM with
; call to PLTPARM, use LINFIT in
; first difference calculation, and
; update documentation.
; 8 Aug. 1984 RWT GSFC modified to use new version of
; LINFIT.
; 8 Nov. 1985 RWT GSFC DIDL modifications, NELEMENTS,
; FINDGEN, and indirect compilation
; 14 Mar. 1987 RWT GSFC VAX mods: add PARCHECK, use
; SET_XY, remove some EXTRACs
; 23 Sept. 1987 RWT GSFC compile RDAF version of GAUSSFIT,
; not IDL user's library procedure
; of the same name, and add listing
; of procedure call.
; 8 Oct. 1987 CAG GSFC Alter GAUSSFIT call to GAUSSFITS,
; to reflect change of procedure
; name. Convert BEGL and ENDL to
; be floating point, and correct one
; format statement.
; 14 Oct. 1987 RWT GSFC remove subroutine NETVEL, and
; compilation for PARCHECK.
; 16 Nov. 1987 RWT GSFC improve plot annotation.
; 19 Jan. 1988 RWT GSFC correct error for image numbers
; larger than 32767.
; 5 May 1988 CAG GSFC add VAX RDAF-style prolog.
; 23 Aug 1989 RWT GSFC Unix mods.
; 25 Jun. 1991 PJL GSFC cleaned up; tested on SUN and VAX;
; updated prolog
; 23 Feb. 1993 LLT GSFC Removed style and range keywords
; in OPLOT (not allowed in V. 3)
; 7/15/93 RWT GSFC remove ctstrim subroutine, allow linear wavelengths,
; make vdel an optional parameters, and use vectors
; as input rather than SAV files.
; 7/21/93 RWT GSFC add IUETERP to resample data in linear wavelength mode,
; and add ccf keyword
; 10/25/93 RWT GSFC write user prompts to text window instead of plot
; window (for x-window sessions), allow < 15 points
; in correlated wavelength region, and test for number
; of points available when search area is increased.
;-
;**********************************************************************
pro crscor,iw1,if1,iw2,if2,begl,endl,del,vinc=vdel,ccf=cross
;
npar = n_params()
if npar eq 0 then begin
print,' CRSCOR,W1,F1,W2,F2,BEGL,ENDL,delta,vinc=vdel,ccf=ccf'
retall
endif ; npar
parcheck,npar,[6,7],'CRSCOR'
if (strlowcase(!d.name) eq 'tek') then xw=0 else xw=1
;
; read data from .sav files
;
file1='spectrum 1'
file2='spectrum 2'
w1=iw1
f1=if1
w2=iw2
f2=if2
deltaw = w2(1) - w2(0) ; wavelength increment
print,' '
;
; check whether linear or logarithmic wavelengths are requested
; if logarithmic, set km/s resolution
;
if (keyword_set(vdel)) then begin
print,' Generating plot of normalized interpolated flux vs. log lambda'
ws = 1
mi = fix(599586. * (endl - begl) / (begl + endl) / vdel +.5) + 1
rwlog = findgen(mi) * (alog10(endl) - alog10(begl)) / (mi-1)
endif else begin
print,' Generating plot of normalized interpolated flux vs. lambda'
ws = 0
vdel = w2(1) - w2(0)
endelse
;
; check wavelengths limits
;
jump=0
begl=float(begl)
endl=float(endl)
n1 = n_elements(w1)
n2 = n_elements(w2)
if (w1(0) gt begl) or (w1(n1-1) lt endl) or $
(w2(0) gt begl) or (w2(n2-1) lt endl) then begin
print,' specified wavelength limits outside input wavelength array(s)'
jump = 1
endif
;
; extract and if requested, interpolate data to log-lambda scale
;
tabinv,w1,begl,j
tabinv,w1,endl,k
j=fix(j)-1>0
m1=fix(k+1)<(n1-1)
f1 = f1(j:m1)
w1 = w1(j:m1)
if (ws) then quadterp,alog10(w1)-alog10(begl),f1,rwlog,ffir $
else begin
ffir = f1
rwlog = w1
mi = m1 - j + 1
endelse
tabinv,w2,begl,j
tabinv,w2,endl,k
j=fix(j)-1>0
m2=fix(k+1)<(n2-1)
f2 = f2(j:m2)
w2 = w2(j:m2)
if (ws) then quadterp,alog10(w2)-alog10(begl),f2,rwlog,fsec $
else begin
fsec = f2
rwlog = w2
mi = mi < (m2-j+1)
if (total(w1(0:1)-w2(0:1)) ne 0.0) then begin
iueterp,w2,w1,f1
print,'warning: 1st spectrum resampled to match 2nd spectrum'
endif
endelse
nelfir = n_elements(ffir)
nelsec = n_elements(fsec)
print,' '
;
if jump eq 0 then begin ;data satisfactory
;
; plot the results of the data extraction
;
plot,rwlog,ffir/max(ffir)+0.6,psym=0,xstyle=1,ystyle=1, $
xrange=[rwlog(0),rwlog(mi-1)],yrange=[0.0,2.0]
oplot,rwlog,fsec/max(fsec)+0.1,psym=0
hpos = !d.x_ch_size
vpos = !d.y_size - !d.y_ch_size
vspace = 1.1 * !d.y_ch_size
st = ' Cross Correlation of spectrum 1 (top) vs. spectrum 2 (bottom)'
xyouts,hpos,vpos,font=0,/device,st
if (ws) then begin
stf = '(t10,"(abscissa is log lambda -",f7.4,")")'
st = string(format=stf,alog10(begl))
xyouts,hpos,vpos-vspace,font=0,/device,st
stf = '(1x,i4," points extracted from ",i4," (file2), interpolated to")'
st = string(format=stf,m2,n2)
xyouts,hpos,vpos-2*vspace,font=0,/device,st
stf = '(1x,i4," points in log lambda with spacing",e11.4)'
st = string(format=stf,mi,rwlog(1))
vpos = vpos - 3*vspace
xyouts,hpos,vpos,font=0,/device,st
endif else begin
stf = '(1x,i4," points from spectrum ",i1," extracted from ",i4)'
sn = 1b
st = string(format=stf,nelfir,sn,n1)
xyouts,hpos,vpos-2*vspace,font=0,/device,st
vpos = vpos - 3*vspace
sn = 2b
st = string(format=stf,nelsec,sn,n2)
xyouts,hpos,vpos,font=0,/device,st
endelse
;
; cross-correlate the data arrays
;
vpos = vpos - vspace
fst = '-- type 1 (one) to proceed, 0 (zero) to end --'
if (xw) then begin
print,fst
read,inp
endif else xyreads,hpos,vpos,inp,/device,fst
if inp eq 1 then begin ;correlation
nspr = 15 < (nelfir-1) ;initial size of correlation
repeat begin ;satisfactory correlation
if inp eq 2 then begin ;double range
nspr=nspr+nspr
if (nspr gt (nelfir-1)) then begin
if (xw) then begin
print,' '
print,' Requested search area is outside wavelength range'
print,' Returning from crscor'
endif else begin
xyouts,hpos,vpos-4*vspace,font=0,/device $
,'Requested search area is outside wavelength range'
xyouts,hpos,vpos-5*vspace,font=0,/device $
,'Returning from crscor'
endelse
retall
endif else begin
if (xw) then print,' doubling range' else $
xyouts,hpos,vpos-vspace,font=0,/device,' doubling range'
endelse
endif ; double range
crsprod,ffir,fsec,nspr,cross,crmin,crmax
if (ws) then vspac=(indgen(n_elements(cross))-nspr) * 599586. * $
(endl-begl)/(endl+begl)/(mi-1) else $
vspac = deltaw * (indgen(n_elements(cross)) - nspr)
;
; plot results of correlation
;
plot,vspac,cross,psym=6,xrange=[min(vspac),max(vspac)], $
yrange=[-0.10,1.15]
vpos = !d.y_size - !d.y_ch_size
xyouts,hpos,vpos,font=0,/device, $
' Normalized Cross Correlation Function '+file1+' vs.'+file2
vpos = vpos - vspace
stf = '(" min",e10.3,", max",e10.3)'
stf2 = '(" wavelengths:",f7.1," to",f7.1," A")'
st = string(format=stf,crmin,crmax) + string(format=stf2,begl,endl)
xyouts,hpos,vpos,font=0,/device,st
;
; get user's fitting preference
;
st1 = ' -- type 1 (one) to proceed if maximum well defined -'
st2 = ' -- type 2 to double the search area -'
st3 = ' -- type (2+n) to fit n gaussian components -'
if (xw) then begin
print,' '
print,st1
print,st2
print,st3
read,'?',inp
endif else begin
vpos = vpos - vspace
xyouts,hpos,vpos,font=0,/device,st1
xyouts,hpos,vpos-vspace,font=0,/device,st2
xyouts,hpos,vpos-2*vspace,font=0,/device,st3
xyouts,hpos,vpos-3*vspace,font=0,/device,'?'
xyreads,hpos,vpos-3*vspace,inp,/device
endelse
end until inp ne 2 ;satisfactory correlaton
;
; find maximum
;
if inp ge 1 then begin
vpos = vpos - 4*vspace
if (ws) then crsmax,cross,vspac,inp-2,del,sigy,hpos,vpos $
else crsmax,cross,vspac,inp-2,del,sigy,hpos,vpos,/lin
if inp lt 4 then begin
if sigy ne 0 then begin
if (ws) then st = string(format='(1x,a,f7.2," km/s")', $
'standard deviation of a single point is',sigy) else $
st = string(format='(1x,a,f7.2," angstroms")', $
'standard deviation of a single point is',sigy)
vpos = vpos + vspace
xyouts,hpos,vpos,font=0,/device,st
endif ; sigy
st = string(format='(1x,a,f8.3," angstroms")', $
'corresponding wavelength shift is',del*(begl+endl)/599586.)
if (ws) then xyouts,hpos,vpos-vspace,font=0,/device,st
endif ; inp lt 4
endif ; find maximum
endif ; correlation
endif ; data satisfactory
return
end ; crscor