function u=om_spher_pot_aniso(x,q,p,sigmas,xis,rk) ;
% Calculates the electrical potential at point x on the surface of a layered
% conductive sphere induced by a current dipole of strength q at point
% p. The radii of the spheres are given by a vector rk in ascending order,
% the radial (resp. tangential) conductivities are in sigmas (resp. xis) SI units.
% The norm of x must be equal to the last rk.
%
% Usage: u=om_spher_pot_aniso(x,q,p,sigmas,xis,rk) ;
%
% If called without arguments, makes a demo
%
% All vectors (x,q,p) should be in 3D.
%
% Based on formulas (1I),(2I) in Zhi Zhang: "A fast method to compute
% surface potential generated by dipoles within multilayer anisotropic
% spheres", in Phys. Med. Biol. 40 (1995), pages 335-349.
%
% PROBLEM: Getting imaginary numbers solved ad hoc by real()

% Copyright (C) 2010-2017, OpenMEEG developers

[n,m] = size(x);
if m~=3,
    error('x should be a set of row 3D vectors.')
end

q=q(:) ; p=p(:) ; % make vectors
q = q';p = p';
matp = ones(n,1)*p;
matq = ones(n,1)*q;

N = length(sigmas) ; % number of levels

% calculate angles
alpha = om_vangle(matp,matq);
cosgamma = cos(om_vangle(matp,x));
beta = om_vangle(cross(matp,matq),cross(matp,x));
relerr = 1e-20; % keep adding terms as longer as the term to be added is greater

n = 1;
sum = 0;

while 1,
    term=(2*n+1)/n*(norm(p)/rk(N))^(n-1) ;
    v1 = (sqrt(1+4*n*(n+1)*xis(1)./sigmas(1))-1)./2;
    v = v1;
    % initialize matrix M
    rsigma=sigmas(1)/sigmas(2) ;
    rerk = rk(N)/rk(1);
    M = eye(2);
    rkprod = (norm(p)./rk(1)).^v1;
    for k=2:N,
        % calculate v_k
        vold = v;
        % calculate v_(k+1)
        v = (sqrt(1+4*n*(n+1)*xis(k)./sigmas(k))-1)./2;
        rsigma=sigmas(k-1)/sigmas(k) ;
        rerk=rk(N)/rk(k-1) ;
        % calculate matrix M
        M=M*[ v+(vold+1)*rsigma ((vold+1)*rsigma-(v+1))*rerk.^(2*v+1) ; ...
              (vold*rsigma-v)/(rerk.^(2*vold+1)) ...
              ((v+1)+vold*rsigma)*rerk.^(2*(v-vold)) ] ./ (2*vold+1);
        rkprod = rkprod*(rk(k-1)/rk(k)).^v;
    end
    g = n*(2*v+1)/(2*v1+1)/(v*M(2,2)+(1+v)*M(2,1)) *(norm(p)/rk(N))^(2*v1-n)/ rkprod;
    f = v1*g./n;
    P = legendre(n,cosgamma);
    % Attention! Matlab's definition of Legendre includes (-1)^m, unlike
    % the definition supposed for the following formula
    term=term*(f*n*cos(alpha).*P(1,:)'-g*cos(beta).*sin(alpha).*P(2,:)');
    %   disp(['Adding term n=' num2str(n) ' value=' num2str(term) ]);
    sum = sum + term;
    % finish if contribution is small or if something goes wrong
    if (prod(abs(term)<abs(sum)*relerr) & n>20) | n>200,
        break;
    end
    n = n+1;
end
u = sum*norm(q)/(4*pi*sigmas(N)*rk(N)^2);
u = real(u); % to cancel very small imaginary parts
