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