Документ взят из кэша поисковой машины. Адрес оригинального документа : http://star.arm.ac.uk/~csj/idl/FAST/xcorrelate.pro Дата изменения: Tue Oct 20 18:43:14 2015 Дата индексирования: Sun Apr 10 05:47:00 2016 Кодировка: Поисковые слова: reflection nebula

PRO xcorrelate, star, temp_file, temp_v, temp_b, starname, wavwind

;; cross-correlates individual spectra with means

;; star ; select type of template (self, model, self-model)
;; temp_file ; file containing model spectrum
;; temp_v ; shift model by this velocity
;; temp_b ; broaden model by this width

; UNUSED....
; Apply filter function
; fft_filter, template, dw, g_cen, g_sig, tempfilt


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

;; Announce programme and options

print, 'XCORRELATE, star, gplot'
print
print, 'Cross-correlate a series of spectra with a template'
print
IF (star eq 0) THEN $
print, 'star =', star, '; template = mean spectrum'
IF (star eq 1) THEN $
print, 'star =', star, '; template = model'
IF (star eq 2) THEN $
print, 'star =', star, '; template = model - mean'


;; Define holders for CCFs

nwave=n_elements(wave)
nccf = 131072L
IF nwave LT 65537L THEN nccf = 65536L
IF nwave LT 32769L THEN nccf = 32768L
IF nwave LT 16385L THEN nccf = 16384L
IF nwave LT 8193L THEN nccf = 8192L
IF nwave LT 4097L THEN nccf = 4096L
IF nwave LT 2049L THEN nccf = 2048L
IF nwave LT 1025L THEN nccf = 1024L
ccf2n = fltarr(nccf, nspec)
ccf2y = fltarr(nccf, nspec)
ccfy = fltarr(nccf)

; x-axis for ccf
; mode; all;
ix1 = 0;
ix2 = n_elements(wave)-1;


;; 'WAVWIND' defines wavelength region to use in ccf
;; only works in this version for real wavelengths (not log)
IF ( n_params() GT 5 ) THEN BEGIN
IF (n_elements(wavwind) EQ 2 ) THEN BEGIN
ax1 = where(wave GE wavwind(0))
ax2 = where(wave LE wavwind(1))
ix1 = ax1(0)
ix2 = ax2(n_elements(ax2)-1)
ENDIF
ENDIF


;; wavelength
IF ( wave(0) gt 10. ) THEN BEGIN
dw = (wave(ix2) - wave(ix1)) / (ix2-ix1)
mw = wave(nwave/2)
ccfx = (findgen(nccf) - nccf/2) * dw / mw * 2.99792458e5

; log wavelength
ENDIF ELSE BEGIN
dw = (wave(ix2) - wave(ix1)) / (ix2-ix1)
ccfx = (findgen(nccf) - nccf/2) * dw
ccfx = (10.^(-ccfx) - 1.) * 2.99792458e5

ENDELSE



;; TEMPLATE definition

; Extract the template from the mean spectrum (specm)
template = fltarr(nccf)
template(0:n_elements(specm)-1) = specm - 1.0

colors = GetColor(/Load)
IF (strupcase(!d.name) EQ 'X') THEN BEGIN
window,10,title="xcorrelate - template"
plot, wave,template; ,yrange=[-0.6,0.2]
ENDIF

; do we want to redefine the template ?
IF (star gt 0) THEN BEGIN

; read model as template
bluestar = fltarr(nccf)
sp2rd,temp_file,tw,tf,sp_header
vcorr,tw,tf,temp_v/2.9979E5,1
tf = qqsm ( tw,tf,temp_b/2.534 )
twlog = alog( tw )
tflog = tf
dw = (wave(nwave-1)-wave(0))/(nwave-1)
resample,twlog,tflog,wave(0),wave(nwave-1),dw
tf = qsmooth(wave,tflog,0.00002)
bluestar(0:nwave-1) = tflog - 1.0

; define red star spectrum as difference
redstar = template - bluestar


; redefine template
IF (star eq 1) THEN template = bluestar
IF (star eq 2) THEN template = redstar
IF (star eq 3) THEN template = helium

ENDIF

IF (strupcase(!d.name) EQ 'X') THEN BEGIN
IF (star gt 0) THEN BEGIN
oplot,wave,bluestar-0.1,color=colors.blue
oplot,wave,redstar+0.1,color=colors.red
ENDIF
oplot,wave,template,color=colors.green

ENDIF



;; AUTOCORRELATION FUNCTION

; Construct acf for subtracting from ccf
fftemp = fft(template)
ffxff = fftemp * conj(fftemp) * (nccf/8)
acfy = float (fft ( ffxff, +1 ))

; Swap order around
work = acfy (nccf/2:nccf-1)
acfy(nccf/2:nccf-1) = acfy(0:nccf/2-1)
acfy(0:nccf/2-1)=work



testspec = fltarr(nccf)
IF (ix1 gt 1) THEN testspec (0:ix1-1) = 0
IF (ix2 lt nccf-2) THEN testspec (ix2+1:nccf-1) = 0




; Loop over spectra
IF (strupcase(!d.name) EQ 'X') THEN window,7,title='ccf spectra'
FOR it = 0,nspec-1 DO BEGIN

; construct test spectrum
testspec(ix1:ix2) = spec(ix1:ix2, it) - 1.0

;IF (strupcase(!d.name) EQ 'X') THEN BEGIN
; !p.multi=[0,0,2,0,0]
; window,11,title="xcorrelate - individual spectra"
;plot, wave,template,yrange=[-0.6,0.2]
;oplot, wave,testspec,color=cols.red
;ENDIF

; check data is valid
ccfy = 0
IF ( stddev(testspec(ix1:ix2)) gt 1.e-6 ) THEN BEGIN

; construct fast fourier transform
fftest = fft(testspec)

; convolve transforms
; ffxff = fftemp * conj(fftest) * nccf/8
ffxff = fftest * conj(fftemp) * nccf/8

; ccf is real part of inverse transform
ccfy = float (fft ( ffxff, +1 ))

; swap order around
work = ccfy (nccf/2:nccf-1)
ccfy(nccf/2:nccf-1) = ccfy(0:nccf/2-1)
ccfy(0:nccf/2-1)=work

; save residual ccf in 2d array
ccf2n(*,it) = ccfy - acfy

ENDIF

; save full ccf in 2d array
ccf2y(*,it) = ccfy

IF (strupcase(!d.name) EQ 'X') THEN BEGIN
plot,ccfy
ENDIF

ENDFOR


!p.multi=0

END