#!/usr/bin/octave -qf
#
# variables which starts with uppercase char denotes matrix (or vector),
# while the ones with lowercase denotes scalar.
#
rSTEP = 0.1;			# measurement step in centimeter
mSIZE = 256;			# must be even, preferably 2**n.
global rMAX = 5;                # Radial Intensity must be zero for r > rMAX

wSTEP = 1 / (rSTEP * mSIZE);	# given in inverse centimeter

##############
# I(r) definition
##############
function y = RadialIntensity(r)
	y = 0.0;
	if (r < 2.0)
		y = 1;
	endif
endfunction

##############
# Make Pt()
##############
function y = Integrant(x)
	global tproj
	r = sqrt(x^2 + tproj^2);
	y = RadialIntensity(r);
endfunction

function InitInteg(t)
	global tproj
	tproj = t;
endfunction

function y = LineProjection(t)
	global rMAX;
	InitInteg(t);
	[y, ier, nfun, err] = quad("Integrant", -rMAX, rMAX);
endfunction

i=0;
for tm = 0 : rSTEP : rMAX
     i = i + 1;
     Pt(i) = LineProjection(tm);
endfor
pSIZE=rows(Pt);

##############
# Replication
##############
for j = 2 : pSIZE
	Pt(mSIZE - j + 2) = Pt(j);
endfor

##############
# Zero padding
##############
for j = pSIZE + 1 : mSIZE - (pSIZE - 1)
	Pt(j) = 0;
endfor

##############
#construct h(nt)
##############
h = zeros(256,1);
h(1) = 0.25;
for i = 2 : 2: 128
	h(i) = - 1 / (pi * (i-1))^2;
	h(i+1) = 0;
endfor
for i = 256:-1:129
	h(i) = h(256 - i + 2);
endfor
h = h ./ (rSTEP^2);
Hf = rSTEP .* abs( fft(h) );

##############
# Int in Rev. space
##############
St = rSTEP * fft(Pt);				# FFT
Wa = [0.25, 1:128, 127:-1:1 ]' .* wSTEP;	# |w|
W  = [0:128, -127:-1] .* wSTEP;			# w

Hamming = 0.54 + 0.46 * cos (2 * pi * (0:mSIZE-1)' / mSIZE);

for i = 1:pSIZE
	r = (i-1) * rSTEP;
	B = besselj(0, 2 * pi * r * W);
	RWa = abs (B * (Wa .* St .* Hamming) * pi * wSTEP);
	RHf = abs (B * (Hf .* St .* Hamming) * pi * wSTEP);
        Radial(i) = r;
        Iorig = RadialIntensity(r);
        dWa(i) = RWa - Iorig;
        dHf(i) = RHf - Iorig;
#	fprintf(stdout, "%7.4f\t%20.10f\n", r, Rp);
endfor

data = [Radial dWa dHf]
# gset term png color;
# gset output "test.png";
gset term tek40xx
#gplot data using 1:2 with lines, data using 1:3 with lines

