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

PRO xvclean , tvec, dvec, evec, gplot, flat, tclean, dclean



;; Announce programme and options

IF (gplot ne 2) THEN BEGIN

print, 'XVCLEAN, tvec, dvec, gplot, flat, tclean, dclean'
print
print, 'clean variable star data by removing outliers and background trends'
print, 'Input:'
print, 'tvec, dvec : arrays of time and data values '
print, 'gplot = 0 (no graphics), 1 (X or PS), 2 (X window), >2 (more) '
print, 'flat = 0 (no background), 1 (sine wave), 2 (polynomial) '
print
print, 'Output:'
print, 'tclean, dclean : cleaned data '
print

IF (gplot eq 0) THEN $
print, 'gplot =', gplot, '; no graphics'
IF (gplot eq 1) THEN $
print, 'gplot =', gplot, '; graphics on / X or PS'
IF (strupcase(!d.name) EQ 'X') THEN window,3,title='xvclean'
IF (gplot eq 2) THEN $
print, 'gplot =', gplot, '; graphics on / X window'
IF (strupcase(!d.name) EQ 'X') THEN window,3,title='xvclean'
IF (gplot gt 2) THEN $
print, 'gplot =', gplot, '; diagnostic graphics'
IF (flat eq 1) THEN $
print, 'flat =', flat, '; sine wave background'
IF (flat eq 2) THEN $
print, 'flat =', flat, '; polynomial background'
print

ENDIF


;; clean the data

; Initialize
tclean = tvec
dclean = dvec

; Compute sigma, omitting far out data
nspec = N_ELEMENTS(dvec)
svec = STDDEV(dvec)
mvec = TOTAL(dvec)/nspec
FOR it=0,nspec-1 DO BEGIN
IF (ABS(dvec(it)-mvec) gt 2*svec) THEN BEGIN
dclean(it) = mvec
ENDIF
ENDFOR
svec = stddev(dclean)

; Select data within +/- 3 sigma
ic = 0
FOR it=0,nspec-1 DO BEGIN
IF (abs(dvec(it)-mvec) lt 3*svec) THEN BEGIN
tclean(ic) = tvec(it)
dclean(ic) = dvec(it)
ic = ic + 1
ENDIF
ENDFOR

IF (ic eq 0) THEN ic = nspec
tclean = tclean(0:ic-1)
dclean = dclean(0:ic-1)
mclean = TOTAL(dclean)/ic
dclean = dclean - mclean


; plot the raw data
IF (gplot gt 2) THEN BEGIN
print, mvec, mclean, svec
plot,tvec,dvec, $
psym=1, $
/xstyle, xtitle='Time', $
ytitle='cleaned data', $
yrange=[mclean-3*svec,mclean+3*svec], $
title=ptitle, $
charsize=1.4
ENDIF



;; Subtract low-frequency background


; initialize background fit
bg_fit = fltarr( ic )
bg_fit = 0.


; Polynomial background
IF ( flat eq 2 ) THEN BEGIN
bg_coeffs = poly_fit (tclean, dclean, 2, bg_fit)
dclean = dclean - bg_fit
print, "Background coefficients"
print, "bg = a0 + a1 * x + a2 * x^2"
print, bg_coeffs, format='(3f10.2)'
ENDIF


; Sine wave background
IF ( flat eq 1 ) THEN BEGIN
bg_coeffs = [0.,5.,0.55,1.]
wvec = fltarr( ic )
wvec = [replicate(1.0, ic)]
bg_fit = curvefit (tclean, dclean, wvec, $
bg_coeffs, function_name='sinefit', $
bg_sigma, noderivative=1 )
dclean = dclean - bg_fit
print, "Background coefficients"
print, "bg = a0 + a1 * sin ( 2.pi.(x-a2) / P )
print, bg_coeffs
print, "+/-", bg_sigma
ENDIF


;; Plot the data and the background fit

IF (gplot gt 0) THEN BEGIN
plot,tclean,dclean+bg_fit, $
psym=+1, $
/xstyle, xtitle='Time', $
ytitle='Velocity (km/s)', $
$ ; yrange=[-50,+50], $
title=ptitle, $
charsize=1.4
oploterr,tclean,dclean+bg_fit,evec
IF (flat gt 0 ) THEN oplot,tclean,bg_fit,color=45
ENDIF


END