Документ взят из кэша поисковой машины. Адрес оригинального документа : http://star.arm.ac.uk/~csj/idl/FAST/xcorr_fits.pro Дата изменения: Tue Oct 20 19:50:06 2015 Дата индексирования: Sun Apr 10 05:46:58 2016 Кодировка: Поисковые слова: п п п п п п п п п п п

PRO xcorr_fits, dfl, mode, gplot, xsave


;; calculates centroids of cross-correlation functions
;; from 2d array

;; dfl is size of window to fit peak (in pixels / 2)
;; mode = 0, fit Gaussian,
;; = 1, fit parabola
;; = 2, compute moments
;; gplot = 0, no plots,
;; = 1, plot ccfs and fits
;; xsave = 0, no save
;; = 1, save last ccf and fit



COMMON fitsblock
COMMON xcorblock, ccfx, ccfy, ccf2y, ccf2n, nccf, afits, efits, template
COMMON colorblock


;; Announce programme and options

print, 'XCORR_FITS, dfl, model, gplot, xsave'
print
print, 'Analyse a series of CCFs'
print
print, 'dfl =',dfl, '; half-width (pixels) of region used for fit '
IF (mode eq 0) THEN $
print, 'mode =', mode, '; fit Gaussian'
IF (mode eq 1) THEN $
print, 'mode =', mode, '; fit parabola'
IF (mode eq 2) THEN $
print, 'mode =', mode, '; compute moments 0 - 3: 1 = centroid'
IF (gplot eq 0) THEN $
print, 'gplot =', gplot, '; no graphics'
IF (gplot eq 1) THEN $
print, 'gplot =', gplot, '; graphics on'
IF (xsave eq 0) THEN $
print, 'xsave =', xsave, '; no saves'
IF (xsave eq 1) THEN $
print, 'xsave =', xsave, '; save last xcorr and fit'
print


; Define window for fitting ccf
plot_factor=20
nfl=(nccf/2)-dfl
nfu=(nccf/2)+dfl
pfl=plot_factor*dfl
wfl=(nccf/2)-pfl
wfu=(nccf/2)+pfl

efits=afits
; Loop over spectra
FOR it = 0,nspec-1 DO BEGIN
afit = [0.0,0.0,1.0,0.0,0.0,1.0E-8]
afits(*,it) = afit

; copy 1d ccf to work array and check valid
ccfy = ccf2y(*,it)
IF ( stddev(ccfy) gt 1.e-6 ) THEN BEGIN

ccfpeak=ccfy(wfl:wfu)
xxfmax=max(ccfpeak,maxccf_x)
maxccf=maxccf_x+wfl
kfl=maxccf-dfl
kfu=maxccf+dfl

; select central fraction of ccf
v = ccfx(kfl:kfu)
p = ccfy(kfl:kfu)

; fit gaussian to peak
IF (mode eq 0) THEN BEGIN
gsb = p(0) ; get estimates for Gaussian
gst = p(dfl)
a0 = [gst-gsb,0.0,dfl*30,gsb,0.0,0.0]
ccfit = gaussfit(v,p,afit,estimates=a0,nterms=6)

IF (afit(1) GT dfl*30) THEN afit = [0.0,0.0,1.0,0.0,0.0,1.0E-6]
IF (afit(0) LT 0.01) THEN afit = [0.0,0.0,1.0,0.0,0.0,1.0E-6]
IF TOTAL(FINITE(afit)) NE N_ELEMENTS(afit) THEN $
afit = [0.0,0.0,1.0,0.0,0.0,1.0E-8]
IF (ABS(afit(5)) LT 1.E-8) THEN afit(5) = 1.0E-8
afits(*,it) = afit
ENDIF

; fit parablola to peak
IF (mode eq 1) THEN BEGIN
prbola, v,p,a,as,t,ts,is,pbfit
afits(0:1,it) = t(1:0:-1)
efits(0:1,it) = ts(1:0:-1)
print,t(0),'+/-',ts(0), 'min/max=', is
ENDIF

; compute moments of line profile
IF (mode eq 2) THEN BEGIN
mom = fltarr(4)
lp_moment, v,p,0.,mom
afits(0:3,it) = mom
ENDIF

; plotting the individual ccfs - for debugging only
IF (gplot gt 0) THEN BEGIN
window,20,title="ccf fits"
ccf_min=ccfx(kfl-pfl)
ccf_max=ccfx(kfu+pfl)
plot, ccfx,ccfy,xrange=[ccf_min,ccf_max]
oplot, v,p,psym=1
IF (mode eq 0) THEN oplot, v,ccfit, color=cols.red
IF (mode eq 1) THEN oplot, v,pbfit, color=cols.green
ENDIF

wait,0.2

ENDIF

ENDFOR

; IF (gplot gt 0) THEN window,0
; IF (gplot gt 0) THEN shade_surf, ccf2y(nfl:nfu,*)

IF (xsave eq 1) THEN BEGIN
openw,lxfit,'xfit.data', /get_lun
FOR ix=0,n_elements(v)-1 DO BEGIN
IF (mode eq 0) THEN printf, lxfit, v(ix),p(ix),ccfit(ix), format='(3f10.3)'
IF (mode eq 1) THEN printf, lxfit, v(ix),p(ix),pbfit(ix), format='(3f10.3)'
ENDFOR
free_lun, lxfit
ENDIF


END