pro tomoftregu2,tol,bound,nalpha,prior,obsspec,recspec2,pshift,ratio,simspec,resi_arr,x_xprior,alpha_arr,alphaBest
;tomoft with regularization,use tomoftergu and find best alpha
;this is the same as in gou81tomo3.pro but for more genreate use.

;tol: SVD singular value tolerance, smaller than tol will be treated as zero
;bound: two elements array,lower and upper bounds of alpha,e.g[0.01,10]
;nalpha: # of regular parameter calculated 
;prior: npix*2
;obsspec:	npix*nrv
;recspec2: npix*2,inital value is going to be set as prior
;pshift: 2*nrv
;ratio:  2*nrv
;----------Outputs:
;recspec2:		npix*2; final recspec
;simspec:		npix*nrv
;resi_arr: nalpha elements array, store the |Gx-d|^2
;x_xprior:	nalpha elements array, store the |x-xprior|^2
;alpha_arr: nalpha elements array,store the regular parameters alpha
;bestalpha: the alpha adopted for final reconstruction
;
;yobs=obsspec;dblarr(npix,nrv)
ydim=size(obsspec,/dimension)
npix=ydim[0]
nrv=ydim[1]

;nalpha=20.
low=alog10(bound[0])
up=alog10(bound[1])

alpha0=dindgen(nalpha)/nalpha*(up-low)+low;low->up
alpha_arr=10.0^alpha0;

resi_arr=dblarr(nalpha)
x_xprior=dblarr(nalpha)

for al=0,nalpha-1 do begin
;alpha=0.1
alpha=alpha_arr[al]

recspec2=prior;need to reset prior, recspec2 initial values=prior
tomoftregu,tol,alpha,obsspec,recspec2,pshift,ratio,simspec
rec1=recspec2(*,0)
rec2=recspec2(*,1)

x_xprior[al]=total((rec1-prior(*,0))^2)+total((rec2-prior(*,1))^2)

sum=0.0d
for j=0,nrv-1 do begin
sum=sum+total((obsspec(*,j)-simspec(*,j))^2)
endfor
resi_arr[al]=sum

print,x_xprior[al],resi_arr[al],alpha
;alpha=-2->2, resi,xxprior,alpha saved in xyz.txt,and used by test_findcorner
endfor

window,18,xsize=1200,ysize=600
x=alog10(resi_arr)
y=alog10(x_xprior)
plot,x,y,psym=4,yrange=[min(y),max(y)],xrange=[min(x),max(x)],xstyle=1,ystyle=1,$
title='log(resi) vs log(x_xpior)',xtitle='resi',ytitle='x_xprior'

;use	findcorner
findcorner,resi_arr,x_xprior,alpha_arr,alphaBest,g,kap
print,'best alpha=',alphaBest,' index=',g
oplot,[x(g)],[y(g)],psym=1,color=cgcolor('red')
for i=0,nalpha-1 do begin
	str=strtrim(string(alpha_arr[i],format='(f6.3)'),2)
	xyouts,[x(i)],[y(i)],str
endfor

;calculate	curvature	numerially
dy1=deriv(x,y)
dy2=deriv(x,dy1)
curv=abs(dy2)/(1.0+dy1^2)^1.5; Curvature
gg=where(curv eq max(curv))
;print,gg,alpha(gg)
oplot,x[gg],y[gg],psym=4,color=cgcolor('red')
window,19
plot,alpha_arr,curv,psym=1,title='alpha vs curvature'
oplot,alpha_arr,kap,psym=4,color=cgcolor('red');curvature from findcorner.pro approx.method
;--------------Lcurve end-----------------------------also see test_findcorner.pro
;use the best alpha to do final tomo again
recspec2=prior
tomoftregu,tol,alphaBest,obsspec,recspec2,pshift,ratio,simspec

end