function sol = fndrc(tau1,tau2,zeta,r1,r2,adim)
%
% Computes the roots of F.
% All quantities in this function are dimensionless.
%
% (c) S. Balac - MAPLY - INSA Lyon - mars 2001
%

global g   % parameters of the experiment
global mu
global mu0
global CA
global CB
global CC
global CD
global nr
global np

dmax = 20; % maximal allowed distorsion = (dmax * adim) meters
dmin = 0.1; % minimal distance between 2 "distincts" solution = (dmin* adim) meters
           % This is used to sort the "physical solutions"
	
% Domain 3
% ========
alpha = 0;
beta = mu(3)*mu0*CB;
rac = [];
so1 = [];
P = [];                          % Computation of the polynom coefficients
P = polnom(tau1,tau2,zeta,g(1),g(2),g(3),nr(1),np(1),alpha,beta,adim);
rac = roots(P);                  % Computation of its 5 roots
% Test for real roots
for i=1:length(rac)
  if isreal(rac(i)), so1=[so1 , rac(i)]; end
end
% Test if the point belongs to domain 3
rac3=[];
for i=1:length(so1)
  xs = so1(i);
  xr = tau1 + (xs-zeta);
  xp = tau2;
  daxe2 = xr^2+xp^2; 
  if ((daxe2 > r2^2)&(abs(xs)<dmax)) rac3=[rac3 so1(i)];end;
end
% Domain 2
% ========
alpha = -mu(2)*mu0*CC;
beta = mu(2)*mu0*CD;
rac = [];
so2 = [];
P = [];                          % Computation of the polynom coefficients
P = polnom(tau1,tau2,zeta,g(1),g(2),g(3),nr(1),np(1),alpha,beta,adim);
rac = roots(P);                  % Computation of its 5 roots
% Test for real roots
for i=1:length(rac)
  if isreal(rac(i)), so2=[so2 , rac(i)]; end
end
% Test if the point belongs to domain 2
rac2=[];
for i=1:length(so2)
  xs = so2(i);
  xr = tau1 + (xs-zeta);
  xp = tau2;
  daxe2 = xr^2+xp^2; 
  if ((daxe2 >= r1^2) & (daxe2 <= r2^2) & (abs(xs)<dmax)) ...
     rac2=[rac2 so2(i)];end;
end
% Domain 1
% ========
alpha = -mu(1)*mu0*CA;
beta = 0;
rac = [];
so3 = [];
P = [];                          % Computation of the polynom coefficients
P = polnom(tau1,tau2,zeta,g(1),g(2),g(3),nr(1),np(1),alpha,beta,adim);
rac = roots(P);                  % Computation of its 5 roots
% Test for real roots
for i=1:length(rac)
  if isreal(rac(i)), so3=[so3 , rac(i)]; end
end
% Test if the point belongs to domain 1
rac1=[];
for i=1:length(so3)
  xs = so3(i);
  xr = tau1 + (xs-zeta);
  xp = tau2;
  daxe2 = xr^2+xp^2; 
  if ((daxe2 < r1^2) & (abs(xs)<dmax)) rac1=[rac1 so3(i)];end;
end


soo = [rac1,rac2,rac3];
% Test for multiple roots
if length(soo) < 2; sol=soo; return; end;
soo = sort(soo);  % Sort so1 in ascending order
sol = [soo(1)];
for i=2:length(soo)
  if abs(soo(i)-soo(i-1)) > dmin; sol = [sol , soo(i)]; end;
end

