Pro disp_2d,m_wave,m_pix,m_ord,m_flag,maxordx,maxordy,npix,nord,obase,wvc \$ ,resid_pix=resid_pix,resid_m_s=resid_m_s common thid_common, thid_order, thid_pixel, thid_maxordx, thid_maxordy ;Prepare line list for 2D fit nm = n_elements(m_wave) m_w = double(m_wave) m_p = double(m_pix) m_o = double(m_ord) m_f = double(m_flag) igd = where(m_flag eq 0, ngd) m_w = m_w(igd) m_p = m_p(igd) m_o = m_o(igd) m_f = m_f(igd) m_ml = m_o * m_w ;List of sp. orders in the correct sequence olist = obase + indgen(nord) ;list of orders if(m_ord(nm-1) lt m_ord(0)) then olist = rotate(olist, 5) ;Determine what order polynomials can be constrained. ordxok = intarr(maxordx+1) for i=0, n_elements(ordxok)-1 do begin bsiz = npix / float(i+1) h = histogram(m_p, min=0, max=npix-1, bin=bsiz) if min(h) gt 0 then ordxok(i) = 1 endfor nordxok = total(ordxok) ordyok = intarr(maxordy+1) for i=0, n_elements(ordyok)-1 do begin bsiz = nord / float(i+1) h = histogram(m_o, min=obase, max=obase+nord-1, bin=bsiz) if min(h) gt 0 then ordyok(i) = 1 endfor nordyok = total(ordyok) ;Choose polynomial order in X. if nordxok eq 1 then begin iwhr = where(ordxok eq 1) ordx = iwhr(0) print, 'Polynomial order in X: ', strtrim(ordx, 2) endif else begin list = '' for i=0, nordxok-1 do begin if ordxok(i) eq 1 then begin list = list + ',' + strtrim(i, 2) endif endfor ordx=max(where(ordxok eq 1)) ; print, 'Polynomial order in X: ', strtrim(ordx, 2) endelse ;Choose polynomial order in Y. if nordyok eq 1 then begin iwhr = where(ordyok eq 1) ordy = iwhr(0) ; print, 'Polynomial order in Y: ', strtrim(ordy, 2) endif else begin list = '' for i=0, nordyok-1 do begin if ordyok(i) eq 1 then begin list = list + ',' + strtrim(i, 2) endif endfor ordy=max(where(ordyok eq 1)) ; print, 'Polynomial order in Y: ', strtrim(ordy, 2) endelse ;Decide on cross-terms (use all that are OK). maxcross=6 xok = intarr(maxcross) xlist = '' nxok = 0 if ordx ge 1 and ordy ge 1 and ngd ge ordx+ordy+nxok+2 then begin xok(0) = 1 xlist = xlist + 'XY ' nxok = nxok + 1 endif if ordx ge 2 and ordy ge 1 and ngd ge ordx+ordy+nxok+2 then begin xok(1) = 1 xlist = xlist + 'XXY ' nxok = nxok + 1 endif if ordx ge 1 and ordy ge 2 and ngd ge ordx+ordy+nxok+2 then begin xok(2) = 1 xlist = xlist + 'XYY ' nxok = nxok + 1 endif if ordx ge 2 and ordy ge 2 and ngd ge ordx+ordy+nxok+2 then begin xok(3) = 1 xlist = xlist + 'XXYY ' nxok = nxok + 1 endif if ordx ge 3 and ordy ge 1 and ngd ge ordx+ordy+nxok+2 then begin xok(4) = 1 xlist = xlist + 'XXXY ' nxok = nxok + 1 endif if ordx ge 1 and ordy ge 3 and ngd ge ordx+ordy+nxok+2 then begin xok(5) = 1 xlist = xlist + 'XYYY ' nxok = nxok + 1 endif nxok = total(xok) ; if nxok gt 0 then print, 'Cross terms: ' + xlist ;Load common block used by fitting routine. thid_order = m_o / 1d2 thid_pixel = m_p / 1d3 thid_maxordx = maxordx thid_maxordy = maxordy ;Use Marquardt to fit marked lines. dummy = fltarr(ngd) sig = replicate(1.0, ngd) par = dblarr(maxordx + 1 + maxordy + maxcross) par(0:1) = poly_fit(m_p/1d3, m_ml, 1, /double) dpar = [ 1d0, replicate(1.d-2, maxordx) \$ , replicate(1.d-2, maxordy) \$ , replicate(1.d-2, maxcross) ] if ordx lt maxordx then dpar(ordx+1:maxordx) = 0.0 if ordy lt maxordy then dpar(ordy+maxordx+1:maxordx+maxordy) = 0.0 iwhr = where(xok eq 0, nwhr) if nwhr gt 0 then dpar(maxordx+maxordy+1+iwhr) = 0.0 mlfit = marq('thid_func', dummy, m_ml, sig, par, dpar, trace=0) ;Construct wavelength coefficients. vers = 2.5 fill = 0.0 wvc = [ vers, npix, nord, obase, fill, fill, fill \$ , maxcross, maxordx, maxordy, par ] ;Calculate wavelengths for every pixel. if(arg_present(resid_pix) or arg_present(resid_m_s)) then begin wave=dblarr(npix,nord) dummy = fltarr(npix) thid_pixel = dindgen(npix) / 1d3 for i=0, nord-1 do begin thid_order = replicate(olist(i), npix) / 1d2 wave(*,i) = thid_func(dummy, par) / olist(i) endfor resid_pix = fltarr(nm) resid_m_s = fltarr(nm) xpix = dindgen(npix) for i=0,n_elements(olist)-1 do begin order=olist[i] iwhr = where(m_ord eq order, nwhr) if nwhr gt 0 then begin resid_pix(iwhr) = m_pix(iwhr) - interpol(xpix, wave(*,order-obase), m_wave(iwhr)) resid_m_s(iwhr) = m_wave(iwhr) - interpol(wave(*,order-obase), xpix, m_pix(iwhr)) endif endfor resid_m_s=resid_m_s/m_wave*2.9979246d8 endif return end ; f=findfile('\nikolai\reduction\elodie_order*.sav') ; f=f(sort(f)) ; m_wave=0.d0 ; m_pix=0.d0 ; m_flag=0 ; m_ord=0. ; npix=1024 ; ord=0 ; for i=0,n_elements(f)-1 do begin ; order=fix(strmid(f[i],strlen(f[i])-7,3)) ; restore,f[i] ; nlines=n_elements(line) ; ord=[ord,order] ; m_wave=[m_wave,line.wll] ; m_pix=[m_pix,line.posm] ; m_flag=[m_flag,line.flag] ; m_ord=[m_ord,replicate(order,nlines)] ; endfor ; ord=ord[1:*] ; m_wave=m_wave[1:*] ; m_pix=m_pix[1:*] ; m_flag=m_flag[1:*] ; m_ord=m_ord[1:*] ; nm=n_elements(m_wave) ; ; obase=min(m_ord) ; nord=max(m_ord)-obase+1 ; maxordx=5 ; maxordy=5 ; ; SDtoFWHM=2.d0*sqrt(2.d0*alog(2.d0)) ; !p.region=0 & !p.position=0 ; !p.multi=0 ; device,decomposed=0 & colors ; iter=0 ; xmin=-0.5>min(resid) ; xmax= 0.5min(resid) ; xmax= 6.*a(2)/SDtoFWHM