pro tomoft,tol,obsspec,recspec,pshift,ratio,simspec ;a similar pro to tomograph.pro, but in FT domain,3times faster ;NO regularization ;Inputs:-------------- ;tol: SVD singular value tolerance, smaller than tol will be treated as zero; ; usually 1d-4,doesn't matter really ;obsspec: npix*nrv ;recspec: npix*2; set inital values are prior spectra ;pshift: 2*nrv ;ratio: 2*nrv ;(nstar=2 for now) ;Outputs:--------------- ;recspec: npix*2 ;simspec: npix*nrv yobs=obsspec;dblarr(npix,nrv) ydim=size(yobs,/dimension) npix=ydim[0] nrv=ydim[1] pdim=size(pshift,/dimension) nstar=pdim[0] ym=dcomplexarr(npix,nrv);FT of the obsspec for i=0,nrv-1 do begin yobs(*,i)=obsspec(*,i)-mean(obsspec(*,i)) ym(*,i)=fft(yobs(*,i),/double) endfor ;------------------------- ;x1m0=fft(x1,/double) ;x2m0=fft(x2,/double); npix*1 array,each element is a complex number ;bet=dblarr(5,2) bet=dblarr(nrv,nstar) bet=transpose(pshift) ;rv1arr=[ -4.32491d,4.42757d,4.29319d,4.04728d,3.69769d];,\$ ;rv2arr=[ 4.05250d,-4.36323d,-4.25879d,-4.05420d,-3.75160d];,\$ ;bet(*,0)=rv1arr ;bet(*,1)=rv2arr ;ratio=dblarr(2,nrv) ;ratio(0,*)=f1 ;ratio(1,*)=f2 ;L=dblarr(2,nrv) L=ratio ;L=[f1,f2] Farr=dcomplexarr(nstar,nrv,npix) ;nrv obs spec, 2 stars,npix pixels for i=1,nrv do begin for j=1,nstar do begin Lj=L(j-1,i-1) col=dblarr(npix) ;npix=n_elements(col) temp=abs(bet(i-1,j-1)) low=fix(temp) high=low+1.0 fhigh=temp-low flow=1.0-fhigh if(bet(i-1,j-1) lt 0.)then begin if(low eq 0.)then begin;speciall cases when RV within[-1,0],2 nonzero elements are col(npix-1)and col(0) col(0)=flow*Lj col(npix-high)=fhigh*Lj endif else begin col(npix-low)=flow*Lj col(npix-high)=fhigh*Lj endelse endif else begin col(low)=flow*Lj col(high)=fhigh*Lj endelse ;print,i,j,'col=',col Farr(j-1,i-1,*)=fft(col,/double) endfor endfor ;print,Farr(0,0,*)*npix;Fm11 compare with Fm results in tomo9.pro ;print,'-------' ;print,Farr(1,0,*)*npix;Fm12 Farr=Farr*npix x_inv=dcomplexarr(npix,nstar);ft of component spectra ysim=dcomplexarr(npix,nrv);ft of simspec for mm=0,npix-1 do begin ;y=[y1m[mm],y2m[mm],y3m[mm],y4m[mm],y5m[mm]] y=ym(mm,0) for iy=1,nrv-1 do y=[y,ym(mm,iy)];form the FT(yobs)at freq=m ;x=[x1m0[mm],x2m0[mm]] Fm=Farr(*,*,mm);the Fm matrix for freq=mm ;----svd solution for each mm la_svd,Fm,w,u,v wfilter=w for i=0,n_elements(w)-1 do begin if(w[i] le tol)then begin wfilter[i]=0.0d endif else begin wfilter[i]=1.0/w[i] endelse endfor ;for each freq=mm,x_inv is a 2 elements array x_inv(mm,*)=v##diag_matrix(wfilter)##transpose(conj(u))##y ysim(mm,*)=Fm##x_inv(mm,*);nrv*nstar nstar*1 =nrv*1,all complex arrays ;print,mm,ysim(mm,*); endfor ;x1_inv=x_inv[0,*] ;x2_inv=x_inv[1,*] ;x1orig=fft(x1_inv,/inverse,/double) ;x2orig=fft(x2_inv,/inverse,/double) ;recspec=dblarr(npix,2) for i=0,nstar-1 do begin recspec(*,i)=real_part(fft(x_inv[*,i],/inverse,/double)) endfor ;ysim(npix,nrv) simspec=dblarr(npix,nrv) for i=0,nrv-1 do begin simspec(*,i)=real_part(fft(ysim(*,i),/inverse,/double)) simspec(*,i)=simspec(*,i)+mean(obsspec(*,i)) ;dblarr(npix,nrv) endfor end