function sol = fndrc(tau1,tau2,zeta,lc)
%
% Purpose : computation of the zeros of the function F.
% See "Mathematical modeling and numerical simulation of magnetic
%      susceptibility artifacts in Magnetic Resonance Imaging" 
%      S. Balac and G. Caloz
%      Computer Methods in Biomechanics and Biomedical Engineering (3:335-349 (2000)) 
% 
% Arguments in :
%  see article
% 
% Arguments out :
% the zeros of the function F.
%
% ARAMIS - (c) S. Balac - ENSSAT - Université de Rennes 1 - 2007 (version 2.0)
%

global g
global r0
global B0
global Xm

rac = [];P = [];
P   = polnom(tau1,tau2,zeta,lc); % Computation of the polynome coefficients
rac = roots(P);                  % Computation of its 12 roots

so1=[];
% Test for real roots
for i=1:length(rac)
  if isreal(rac(i)), so1=[so1 , rac(i)]; end
end
% Test for roots satisfying: F(xs) = 0 
so2=[];
for i=1:length(so1)
  xs = so1(i);	
  xr = tau1 + (xs-zeta);
  xp = tau2;
  r  = sqrt(xr^2 + xp^2 + (xs+lc)^2);
  res = F(xr,xp,xs,zeta,lc);
  if ( abs(res) < 1e-10 & r > 1 ), so2=[so2 , so1(i)]; end
end
% Test for multiple roots
if length(so2) < 2; sol=so2; return; end;
so2 = sort(so2);  % Sort so2 in ascending order
sol = [so2(1)];
for i=2:length(so2)
  if abs(so2(i)-so2(i-1)) > 1e-3; sol = [sol , so2(i)]; end;
end
% Test for roots coming from the inside magnetic field expression
xs = zeta - 2*B0*Xm/(3*g*r0);
xr = tau1 + (xs-zeta);
xp = tau2;
r  = sqrt(xr^2 + xp^2 + (xs+lc)^2);
if ( r<1 ) sol = [sol xs]; end
