PRO TOMOGRAPHY,niters,delta,obsspec,recspec,pshift,ratio,simspec ; IDL version of tomography algorithm based on Wiemker's thesis ; Input: ; niters = # of iterations ; delta = convergence acceleration (0 to 1) ; obsspec = observed composite spectra ; recspec = reconstructed spectra (current estimate) ; pshift = pixel shifts for each spectrum, component ; ratios = flux ratio for each spectrum, component ; (sum to unity) ; Output: ; recspec = reconstructed spectra (revised estimate) ; simspec = composite spectra for same shifts as observed np=n_elements(recspec(*,0)) ; # Pixels 8 nc=n_elements(recspec(0,*)) ; # Components 2 j ns=n_elements(obsspec(0,*)) ; # composite Spectra 3 i ; Weighting factors is1=intarr(nc,ns) ; integer shifts is2=intarr(nc,ns) w1=fltarr(nc,ns) ; fractional shifts w2=fltarr(nc,ns) wt1=fltarr(nc,ns) ; weights wt2=fltarr(nc,ns) for i=0,ns-1 do begin ; for each spectrum for j=0,nc-1 do begin ; for each component is1(j,i)=fix(pshift(j,i)) ;integer shifts if (pshift(j,i) lt 0.) then is2(j,i)=is1(j,i)-1 else \$ is2(j,i)=is1(j,i)+1 w2(j,i)=abs(pshift(j,i)-is1(j,i));fraction pix shift coeff., rvfrac1,rvfrac2 w1(j,i)=1.-w2(j,i) endfor wtsum=total(w1(*,i)^2+w2(*,i)^2) wt1(0,i)=delta/ratio(*,i)*w1(*,i)/wtsum wt2(0,i)=delta/ratio(*,i)*w2(*,i)/wtsum w1(0,i)=ratio(*,i)*w1(*,i) w2(0,i)=ratio(*,i)*w2(*,i) endfor ; Main iteration ... npp=np+100;108 ;pad 50 pixel on each side, 50np50 rs=fltarr(npp,nc)+1. ; pad both sides with unity cs=fltarr(npp,nc) resid=fltarr(npp,ns) ; pad both sides with zeroes cr=fltarr(npp,ns) top=np+49 for i=0,niters-1 do begin for k=0,nc-1 do rs(50,k)=recspec(*,k) ; form matrix of simulated spectra simspec=0.*obsspec for j=0,ns-1 do begin ; compute each spectrum by------------------ for k=0,nc-1 do \$ ; adding contributions from each component cs(0,k)=w1(k,j)*shift(rs(*,k),is1(k,j)) + \$ w2(k,j)*shift(rs(*,k),is2(k,j));component flux pixel css=rebin(cs,npp,1)*nc ;compress cs to npp*1 array=css;essential =add two components spectra to a single spectra=composite spectra simspec(0,j)=css(50:top);css is 50np50=108pixel; top=np+49; 50:np+49 is just the middle np pixels ; calculate residuals resid(50,j)=obsspec(*,j)-simspec(*,j) ;50resid, obsspec(*,j) ;IDL> a=[1,2,3,4] ;IDL> a[1]=[5,6,7] ;IDL> print,a ; 1 5 6 7 endfor ; END/compute each spectrum by----------------------- ; calculate corrections matrix for j=0,nc-1 do begin ; for each component---------------- ; get corrections from each spectrum for k=0,ns-1 do \$ cr(0,k)=wt1(j,k)*shift(resid(*,k),-is1(j,k)) + \$ wt2(j,k)*shift(resid(*,k),-is2(j,k)) ;why resid here need to shift back?? recspec(0,j)=recspec(*,j)+rebin(cr(50:top,*),np,1) endfor ; enforce fluxes > 0 recspec=recspec > 0. endfor; END/ for each component----------------------- ; final matrix of simulated spectra simspec=0.*obsspec for k=0,nc-1 do rs(50,k)=recspec(*,k) for j=0,ns-1 do begin ; compute each spectrum by for k=0,nc-1 do \$ ; adding contributions from each component cs(0,k)=w1(k,j)*shift(rs(*,k),is1(k,j)) + \$ w2(k,j)*shift(rs(*,k),is2(k,j)) css=rebin(cs,npp,1)*nc simspec(0,j)=css(50:top) endfor ; write to disk ;writefits,'recspec.fits',recspec return end