#! /bin/csh
# steep.sh  Steepest-Descent Iteration on Radon Inversion
# All programs here in /quake/data/rg
# Set Number of Iterations
set niter=3
# Data already computed from True Model test.ds

# Set Initial Residual to Data test.dt , R = Y
cp test.dt resid.dt
# Set Initial Model to Zero
fac 22500 0 <test.ds >temp.ds

# Iterate, Two Transforms per Iteration
set i=0
while ($i < $niter )

	# Find Gradient g = (A transpose)R
	radon par=islpar i=resid.dt o=grad.ds >&/dev/null

	# Condition Gradient by Smoothing
	smooth i=grad.ds o=grad.ds nx=150 nt=150 xw=3 tw=3

	# Find Conjugate Gradient G = (A)g
	radon par=slpar i=grad.ds o=grad.dt >&/dev/null

	# Find Residual Size R.R
	set RR=`dot 22500 resid.dt resid.dt `

	# Compute alpha = (R.G)/(G.G)
	set RG=`dot 22500 resid.dt grad.dt `
	set GG=`dot 22500 grad.dt grad.dt `
	set alpha=`divide $RG $GG `
	echo $i $RR >>steep.res
	echo Iteration $i : residual $RR ,  alpha $alpha

	# Estimate Next Model x <- x + alpha.g
	fac 22500 $alpha <grad.ds >agrad.ds
	#rm grad.ds
	add 22500 agrad.ds temp.ds >new.ds
	rm agrad.ds
	mv new.ds temp.ds

	# Estimate Next Residual R <- R - alpha.G
	set malpha=`product -1 $alpha `
	fac 22500 $malpha <grad.dt >magrad.dt
	#rm grad.dt
	add 22500 resid.dt magrad.dt >new.dt
	rm magrad.dt
	mv new.dt resid.dt

	#cat new.ds >>all.ds
	@ i = ( $i + 1)
end
mv temp.ds steep.3.ds