#! /bin/csh

# Set Number of Iterations
set niter=200

# Data already computed from True Model test.ds
# Estimate Initial Model from Data test.dt by (A transpose)Y
radon par=islpar i=test.dt o=temp.ds >&/dev/null
#cat temp.ds >all.ds

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

	# Find Estimated Data (A)x
	radon par=slpar i=temp.ds o=new.dt >&/dev/null

	# Find Residual R = Y - (A)x
	neg 22500 <new.dt >nnew.dt
	rm new.dt
	add 22500 nnew.dt test.dt >resid.dt
	rm nnew.dt

	# Find Residual Size R.R
	set res=`dot 22500 resid.dt resid.dt `
	echo Iteration $i : residual $res
	echo $i $res >>gradst.res

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

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

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

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