Документ взят из кэша поисковой машины. Адрес оригинального документа : http://star.arm.ac.uk/~csj/idl/UCLES/ucles_merge.pro Дата изменения: Fri Oct 9 20:35:52 2015 Дата индексирования: Sun Apr 10 05:50:22 2016 Кодировка:

PRO ucles_merge, nw, wb,wi, fn,fs, wm,fm, thresh, notweak=notwk


; Procedure to merge orders from a 2d echelle frame
; into a single spectrum.

; Input
; nw n wavelength / order
; no n orders
; wb nw x no wavelengths
; wi no number of valid wavelengths in ech order
; fs nw x no resampled rectification function
; fn nw x no resampled normalized spectrum

; thresh cutoff low fluxes at order ends

; Output
; wm merged wavelength grid
; fm merged fluxes


PRINT,'UCLES_MERGE: Merging echelle orders...'
colors = GetColor(/Load)

no = n_elements(wb)/nw

;; do not use data where local mean is small
;; fs is now the weighting function to be used in
;; the order merging process.

FOR j = 0,no-1 DO BEGIN
sf = fs(0:wi(j),j)
bar = MAX ( sf )
msk = WHERE ( sf LT thresh * bar, cnt )
IF ( cnt GT 0 ) THEN fs(msk,j) = 0.
ENDFOR


; Find overlaps and tweak rectification

wm = reform(wb(0:wi(no-1),no-1))
fm = reform(fn(0:wi(no-1),no-1))

iovr_lo = fltarr(no)
iovr_hi = fltarr(no)
iovr_lo(no-1) = 0
iovr_hi(0) = wi(0)


; work from low wavelength to high wavelength
FOR j = no-2,0,-1 DO BEGIN

print,'Order:',j,wb(0,j),wb(wi(j),j)

im = n_elements(wm) - 1

k = j + 1
ik = 0
WHILE ( wb(ik,k) LT wb(0,j)+0.001 AND ik LT wi(k) ) DO ik = ik + 1
iovr_lo(k) = ik


ij = 0
WHILE ( wb(ij,j) LT wb(wi(k),k)+0.001 AND ij LT wi(j) ) DO ij = ij + 1
iovr_hi(j) = ij

window,2,title="order overlaps"
plot,wb(0:wi(k),k),fn(0:wi(k),k),xrange=[wb(0,k),wb(wi(j),j)]
oplot,wb(0:wi(j),j),fn(0:wi(j),j),color=colors.yellow
wait, 0.2

print,k,ik,wb(ik,k),wb(wi(k),k),j,ij,wb(0,j),wb(ij,j)

ENDFOR


; Mean values in overlap regions

wbar_lo = fltarr(no)
wbar_hi = fltarr(no)
fbar_lo = fltarr(no)
fbar_hi = fltarr(no)
wbar_hi(0) = wb(wi(0),0 )
fbar_hi(0) = 1.0
wbar_lo(no-1) = wb(0,no-1)
fbar_lo(no-1) = 1.0

print,0,wbar_hi(0),fbar_hi(0),no-1,wbar_lo(no-1),fbar_lo(no-1)


;; investigated a slow convergence approach -- but do-it-in-one was as good.
maxit = 0
FOR niter = 0,maxit DO BEGIN


FOR j = no-2,0,-1 DO BEGIN

k = j+1
ilk = iovr_lo(k)
iuk = wi(k)
IF ( ilk NE 0 ) THEN BEGIN
wbar = moment ( wb(ilk:iuk,k) )
fbar = moment ( fn(ilk:iuk,k) )
wbar_hi(k) = wbar(0)
fbar_hi(k) = fbar(0)
ENDIF

ilj = 0
iuj = iovr_hi(j)
IF ( iuj NE 0 ) THEN BEGIN
wbar = moment ( wb(ilj:iuj,j) )
fbar = moment ( fn(ilj:iuj,j) )
wbar_lo(j) = wbar(0)
fbar_lo(j) = fbar(0)
ENDIF

print,k,wbar_hi(k),fbar_hi(k),j,wbar_lo(j),fbar_lo(j)

ENDFOR


; now regrade the orders so that overlaps merge

FOR j = no-1,0,-1 DO BEGIN

f_hi = 1 & IF j LT no-1 THEN f_hi = fbar_hi(j+1)
f_lo = 1 & IF j GT 0 THEN f_lo = fbar_lo(j-1)
a_j = fbar_lo(j)
b_j = fbar_hi(j)
alph_j = 0.8 * a_j + 0.2 * f_hi
beta_j = 0.8 * b_j + 0.2 * f_lo
IF niter EQ maxit THEN alph_j = 0.5 * a_j + 0.5 * f_hi
IF niter EQ maxit THEN beta_j = 0.5 * b_j + 0.5 * f_lo

del_j = wbar_hi(j)-wbar_lo(j)
w_j = wb(*,j) -wbar_lo(j)
grad_j = (alph_j * del_j + w_j * (beta_j-alph_j)) / ( a_j * del_j + w_j * (b_j-a_j) )
if keyword_set(notwk) THEN grad_j=grad_j*0.+1.

print,j,f_hi,alph_j,a_j,b_j,beta_j,f_lo

IF j EQ 0 THEN !x.range=[wb(0,2),wb(wi(0),0)]
IF j EQ no-1 THEN !x.range=[wb(0,no-1),wb(wi(no-3),no-3)]
IF j LT no-1 AND j GT 0 THEN !x.range=[wb(0,j+1),wb(wi(j-1),j-1)]
plot,wb(0:wi(j),j),fn(0:wi(j),j)*grad_j(0:wi(j)),color=colors.red
IF j LT no-1 THEN oplot,wb(0:wi(j+1),j+1),fn(0:wi(j+1),j+1),color=colors.blue
IF j GT 0 THEN oplot,wb(0:wi(j-1),j-1),fn(0:wi(j-1),j-1),color=colors.blue
oplot,wb(0:wi(j),j),fn(0:wi(j),j),color=colors.yellow

oplot,[wbar_lo(j),wbar_hi(j)],[fbar_lo(j),fbar_hi(j)],color=colors.green
oplot,[wbar_lo(j),wbar_hi(j)],[alph_j,beta_j],color=colors.green

wait,0.5

fn(*,j) = fn(*,j) * grad_j



ENDFOR


ENDFOR
!x.range=0


; Find overlaps and merge

wm = reform(wb(0:wi(no-1),no-1))
fm = reform(fn(0:wi(no-1),no-1))


!y.range=[0,1.5]
FOR j = no-2,0,-1 DO BEGIN

im = n_elements(wm) - 1
plot,wm,fm,xrange=[wb(0,j+1),wb(wi(j),j)]
; print,wm(0),wm(im),im
; print,wb(0,j),wb(wi(j),j),j


ilm = 0L
WHILE ( wm(ilm) LT wb(0,j)+0.001 AND ilm LT im ) DO ilm = ilm + 1
ium = im

ilj = 0
iuj = im - ilm
; WHILE ( wb(iuj,j) GT wm(im)-0.001 ) DO iuj = iuj - 1

k = j+1
ilk = wi(k) - (ium-ilm)
iuk = wi(k)

IF (ium-ilm NE iuj-ilj OR ium-ilm NE iuk-ilk) THEN BEGIN
print,'overlap error',ilj,iuj-ilj,ilm,ium-ilm,iuk,iuk-ilk
print,wm(0),wm(ilm),wm(ium)
print,0., wb(ilk,k),wb(iuk,k),wb(wi(k),k)
print,0., wb(ilj,j),wb(iuj,j),wb(wi(j),j)
STOP
ENDIF

; check that overlap actually exsists before plotting it
IF ( ilj NE 0 OR iuj NE 0 ) THEN BEGIN

wdge = ( 0.5 + 0.5 * cos ( !pi * findgen(ium-ilm+1) / (ium-ilm) ) ) ^ 2


oplot,wb(ilj:iuj,j),fn(ilj:iuj,j),color=colors.red
oplot,wb(ilj:iuj,j),(1.-wdge),color=colors.green
oplot,wb(iuj:wi(j),j),fn(iuj:wi(j),j),color=colors.yellow
oplot,wb(ilj:iuj,j),fs(ilj:iuj,j),color=colors.blue
wait,0.5

fm(ilm:ium) = ( fm(ilm:ium) * fs(ilk:iuk,k)^2 * wdge $
+ fn(ilj:iuj,j) * fs(ilj:iuj,j)^2 * (1.-wdge) ) $
/ ( fs(ilk:iuk,k)^2 * wdge + fs(ilj:iuj,j)^2 * (1.-wdge) )

ENDIF

wm = [wm , wb(iuj+1:wi(j),j)]
fm = [fm , fn(iuj+1:wi(j),j)]


ENDFOR

!y.range=0

END