I have problem with the following code for fan beam image reconstruction matlab gives the following error ??? Error using ==> interp1 at 165 The interpolation points XI should be real . Please help of you have any idea why its not working.thanks

                clear all; close all;
%figure(1);
%x=linspace(-1,1);
%y=x;
%[xx,yy]=meshgrid(x,y);
%rr=xx.^2+yy.^2;
%B=(rr<=9/16)-0.75*(rr<=1/4)+0.25*(rr<=1/16);
%imshow(B,[0 1]);
%pause;
%a is half the image diagonal dimension
%a=sqrt(size(B,1)^2 + size(B,2)^2)/2;
%a;
%a=71.4178;
%D is the distance from the fan-beam vertex to the center of rotation
%D should be a few pixels larger than a
D=80;
alpha=1;
gamam=20;
 t=-1:0.05:1;
gama=-gamam:alpha:gamam;
rb1=radon_bull(gama);
rb= radon_bull(t);
figure(3); hold on;
plot(t,rb,'k-','Linewidth',3);
%plot(x,rb,'r.','Markersize',10);
title('Radon transform of Bulls Eye');
pause;
r2_vals=r2(t);
figure(4); 
plot(t,r2_vals,'k-','Linewidth',3);
pause;
r3_vals=r2(gama);
 n=-gamam/alpha:gamam/alpha;
g2=g(n,alpha);
g2 = [g2(gamam+1:2*gamam+1) g2(1:gamam)];
h=convp(r3_vals,g2);%this is a convolution for one beta angle
figure(5)
plot(h,'k-','Linewidth',3);
pause;
s = -gamam:1/50:gamam;
Ih = interp1(n,h,s,'nearest');
figure(6); hold on;
plot(s,Ih,'k-','Linewidth',3);
title('Interpolation of h');
pause;
dx=1/100;
x=0:dx:2;
y=0:dx:2;
[xx yy] = meshgrid(x,y);
l=-gamam:gamam;
B=back3(l,h,xx,yy,D);
figure(7); hold on;
imshow(B);
axis image; 
axis xy;
title('weighted back projection')
%and radon_bull.m
function rb = radon_bull(t)
%
t = abs(t);
rb = 0.5*sqrt(9 - 16*t.^2).*(t <= 3/4) - 0.75*sqrt(1 - 4*t.^2).*(t <= 0.5) + 1/8*sqrt(1 - 16*t.^2).*(t <= 0.25);
%r2.m is
function [R]=r(x)
R=0.5*sqrt(9-16*abs(x).^2).*(abs(x)<=3/4)-0.75*sqrt(1-4*abs(x).^2).*(abs(x)<=0.5)+1/8*sqrt(1-16*abs(x).^2).*(abs(x)<=0.25);
end
%g.m
function [G]=g(n,alpha)
%G=(1/(4*pi))*(n/sin(n))^2*rlf(n);
%
%g(gama)=1/(4*pi)*((gama)/sin(gama))^2*rlf(gama)
%rlf is the ram-lak filter for L=10;
for i=1:length(n)
 G(i)=1/(4*pi)*(i*alpha/sin(i*alpha))^2*rlf(i);
 end
end
 %rlf.m
function [RLF]=rlf(n)
for i=1:length(n)
 if n(i)==0;
 RLF(i)=100./(2*pi);
  elseif mod(n(i),2)==0;
   RLF(i)=0;
  else
    RLF(i)=-200./(pi^3*(n(i).^2));
 end
end
 %l.m
function [L] = l(r,fi,beta,D)
L=sqrt((D+r*sin(beta-fi))^2+(r*cos(beta-fi))^2);
end
 back3.m
function B = back3(t,h,x,y,D)
%
%r=sqrt(x^2+y^2);
%fi=acos(x/r);
%B=zeros(M,M);
%=-1:1/100:1;
%y=-1:1/100:1;
%[fi,r] = cart2pol(x,y);
%gamam=20;
%[M M] = size(x);
%B = zeros(M,M);
[M M] = size(x);
B = zeros(M,M);
r=sqrt(x^2+y^2);
fi=acos(x/(sqrt(x^2+y^2)));
%for beta=0:36/sqrt(2):2*180
for beta=0:10:360
v = atan(r*cos(beta-fi)/(D+r*sin(beta-fi)));
  B(:) = B(:) + interp1(t,h,v(:),'nearest')*1/(l(r,fi,beta,D)^2);
%  B(:) = B(:) + interp1(t,h,v(:),'linear');
%x=r*cos(fi);
%y=r*sin(fi);
end;
B=B*10;
convp.m
function h = convp(f,g)
%CONVP Convolution of discrete periodic functions.
%   h = CONVP(f,g) convolves the discrete N-periodic functions f and g.
%   The resulting function is a discrete N-periodic function. 
N  =length(f);
for m=1:N
  h(m) = f(1:m)*g(m:-1:1)' + f(m+1:N)*g(N:-1:m+1)';
end;

[Geoff]

Could you please edit this question and format the code using the 'code' button in the toolbar? It's unreadable right now.

[Vera]

I did it