% MEHLER The Mehler conical Legendre functions P_{1/2+tau*i}^m(x) for x >=1, m=0,1,2,...,M
%
function [P,nu0,nu]=mehler(x,tau,M,eps0)
numax=1000; lambda=zeros(numax,1);
nu0=0; nu=0;
P=zeros(M+1,1); Pa=ones(M+1,1); 
if M>0, R=zeros(M,1); end
if x==1
  P(1)=1; return
end
t=tau^2; x1=sqrt(x^2-1); x2=x+x1;
sum0=cos(tau*log(x2))/sqrt(x2); x1=2*x/x1;
lambda(1)=1/(.25+t); lambda(2)=(3-4*t)/((.25+t)*(2.25+t));
for m=2:numax-1
  lambda(m+1)=(1+1/m)*(2*lambda(m)-lambda(m-1)) ...
    /((1+.5/m)^2+(tau/m)^2);
end
nu0=30+floor((1+(.14+.0246*tau)*(x-1))*M);
nu=nu0;
while any(abs(P-Pa)>eps0*abs(P))
  Pa=P;
  nu=nu+40; 
  if nu>1000 
    fprintf('P may be inaccurate\n'), 
    break
  end
  r=0; s=0;
  for m=nu:-1:1
    r=-((1-.5/m)^2+(tau/m)^2)/(x1+(1+1/m)*r);
    s=r*(lambda(m)+s);
    if m<=M, R(m)=r; end
  end
  P(1)=sum0/(1+s);
  if M>0
    for m=1:M
      P(m+1)=R(m)*P(m);
    end
  end
end
fac=1;
if M>0
  for m=1:M
    fac=m*fac;
    P(m+1)=fac*P(m+1);
  end
end