[nu1,nu2,mu,theta,id,pqa,ipqa,ierror]=dxpqnu(nu1,nu2,mu,theta,id,pqa,ipqa,ierror);

function [nu1,nu2,mu,theta,id,pqa,ipqa,ierror]=dxpqnu(nu1,nu2,mu,theta,id,pqa,ipqa,ierror);
persistent a di dmu factmu flok i ia ifmlv ipq ipq1 ipq2 ipsik ipsix ix1 ixs j j0 k nu pq pq1 pq2 r w x x1 x2 xs y z ; 

if isempty(i), i=0; end;
if isempty(ia), ia=0; end;
if isempty(ifmlv), ifmlv=0; end;
if isempty(ipq), ipq=0; end;
if isempty(ipq1), ipq1=0; end;
if isempty(ipq2), ipq2=0; end;
if isempty(ipsik), ipsik=0; end;
if isempty(ipsix), ipsix=0; end;
if isempty(ix1), ix1=0; end;
if isempty(ixs), ixs=0; end;
if isempty(j), j=0; end;
if isempty(j0), j0=0; end;
if isempty(k), k=0; end;
global dxblk1_1; if isempty(dxblk1_1), dxblk1_1=0; end;
%***BEGIN PROLOGUE  DXPQNU
%***SUBSIDIARY
%***PURPOSE  To compute the values of Legendre functions for DXLEGF.
%            This subroutine calculates initial values of P or Q using
%            power series, then performs forward nu-wise recurrence to
%            obtain P(-MU,NU,X), Q(0,NU,X), or Q(1,NU,X). The nu-wise
%            recurrence is stable for P for all mu and for Q for mu=0,1.
%***LIBRARY   SLATEC
%***CATEGORY  C3A2, C9
%***TYPE      doubleprecision (XPQNU-S, DXPQNU-D)
%***KEYWORDS  LEGENDRE FUNCTIONS
%***AUTHOR  Smith, John M., (NBS and George Mason University)
%***ROUTINES CALLED  DXADD, DXADJ, DXPSI
%***COMMON BLOCKS    DXBLK1
%***REVISION HISTORY  (YYMMDD)
%   820728  DATE WRITTEN
%   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
%   901019  Revisions to prologue.  (DWL and WRB)
%   901106  Changed all specific intrinsics to generic.  (WRB)
%           Corrected order of sections in prologue and added TYPE
%           section.  (WRB)
%   920127  Revised PURPOSE section of prologue.  (DWL)
%***end PROLOGUE  DXPQNU
if isempty(a), a=0; end;
if isempty(nu), nu=0; end;
if isempty(pq), pq=0; end;
if isempty(r), r=0; end;
if isempty(w), w=0; end;
if isempty(x), x=0; end;
if isempty(x1), x1=0; end;
if isempty(x2), x2=0; end;
if isempty(xs), xs=0; end;
if isempty(y), y=0; end;
if isempty(z), z=0; end;
if isempty(di), di=0; end;
if isempty(dmu), dmu=0; end;
if isempty(pq1), pq1=0; end;
if isempty(pq2), pq2=0; end;
if isempty(factmu), factmu=0; end;
if isempty(flok), flok=0; end;
pqa_shape=size(pqa);pqa=reshape(pqa,1,[]);
ipqa_shape=size(ipqa);ipqa=reshape(ipqa,1,[]);
% common :: ;
%% common /dxblk1/ nbitsf;
%% common /dxblk1/ dxblk1_1;
% save :: ;
% save ::       ;
%
%        J0, IPSIK, AND IPSIX ARE INITIALIZED IN THIS SUBROUTINE.
%        J0 IS THE NUMBER OF TERMS USED IN SERIES EXPANSION
%        IN SUBROUTINE DXPQNU.
%        IPSIK, IPSIX ARE VALUES OF K AND X RESPECTIVELY
%        USED IN THE CALCULATION OF THE DXPSI FUNCTION.
%
%***FIRST EXECUTABLE STATEMENT  DXPQNU
ierror = 0;
j0 = fix(dxblk1_1);
ipsik = fix(1 +(fix(dxblk1_1./10)));
ipsix = fix(5.*ipsik);
ipq = 0;
%        FIND NU IN INTERVAL [-.5,.5) IF ID=2  ( CALCULATION OF Q )
nu = rem(nu1,1.0d0);
if( nu>=.5d0 )
nu = nu - 1.0d0;
end;
%        FIND NU IN INTERVAL (-1.5,-.5] IF ID=1,3, OR 4  ( CALC. OF P )
if( id~=2 && nu>-.5d0 )
nu = nu - 1.0d0;
end;
%        CALCULATE MU FACTORIAL
k = fix(mu);
dmu = mu;
if( mu>0 )
factmu = 1.0d0;
ifmlv = 0;
for i = 1 : k;
factmu = factmu.*i;
[factmu,ifmlv,ierror]=dxadj(factmu,ifmlv,ierror);
end; i = fix(k+1);
if( ierror~=0 )
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
return;
end;
end;
if( k==0 )
factmu = 1.0d0;
end;
if( k==0 )
ifmlv = 0;
end;
%
%        X=COS(THETA)
%        Y=SIN(THETA/2)**2=(1-X)/2=.5-.5*X
%        R=TAN(THETA/2)=SQRT((1-X)/(1+X)
%
x = cos(theta);
y = sin(theta./2.0d0).^2;
r = tan(theta./2.0d0);
%
%        use ASCENDING SERIES TO CALCULATE TWO VALUES OF P OR Q
%        FOR USE AS STARTING VALUES IN RECURRENCE RELATION.
%
pq2 = 0.0d0;
for j = 1 : 2;
ipq1 = 0;
if( id==2 )
%
%        Z=-LN(R)=.5*LN((1+X)/(1-X))
%
z = -log(r);
[w ,dumvar2,ipsik,ipsix]=dxpsi(nu+1.0d0,ipsik,ipsix);
xs = 1.0d0./sin(theta);
%
%        SERIES SUMMATION FOR Q ( ID = 2 )
%        Q(0,NU,X)=SUM(FROM 0 TO J0-1,(.5*LN((1+X)/(1-X))
%    +DXPSI(J+1,IPSIK,IPSIX)-DXPSI(NU+1,IPSIK,IPSIX)))*A(J)*(.5-.5*X)**J
%
%        Q(1,NU,X)=-SQRT(1./(1.-X**2))+SQRT((1-X)/(1+X))
%             *SUM(FROM 0 T0 J0-1,-NU*(NU+1)/2*LN((1+X)/(1-X))
%                 +(J-NU)*(J+NU+1)/(2*(J+1))+NU*(NU+1)*
%     (DXPSI(NU+1,IPSIK,IPSIX)-DXPSI(J+1,IPSIK,IPSIX))*A(J)*(.5-.5*X)**J
%
%        NOTE, IN THIS LOOP K=J+1
%
pq = 0.0d0;
ipq = 0;
ia = 0;
a = 1.0d0;
for k = 1 : j0;
flok = k;
if( k~=1 )
a = a.*y.*(flok-2.0d0-nu).*(flok-1.0d0+nu)./((flok-1.0d0+dmu).*(flok-1.0d0));
[a,ia,ierror]=dxadj(a,ia,ierror);
if( ierror~=0 )
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
return;
end;
end;
if( mu>=1 )
x1 =(nu.*(nu+1.0d0).*(z-w+dxpsi(flok,ipsik,ipsix))+(nu-flok+1.0d0).*(nu+flok)./(2.0d0.*flok)).*a;
ix1 = fix(ia);
pq_orig=pq; ipq_orig=ipq;    [pq,ipq,x1,ix1,dumvar5,dumvar6,ierror]=dxadd(pq,ipq,x1,ix1,pq,ipq,ierror);    ipq(dumvar6~=ipq_orig)=dumvar6(dumvar6~=ipq_orig); pq(dumvar5~=pq_orig)=dumvar5(dumvar5~=pq_orig);
if( ierror~=0 )
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
return;
end;
else;
x1 =(dxpsi(flok,ipsik,ipsix)-w+z).*a;
ix1 = fix(ia);
pq_orig=pq; ipq_orig=ipq;    [pq,ipq,x1,ix1,dumvar5,dumvar6,ierror]=dxadd(pq,ipq,x1,ix1,pq,ipq,ierror);    ipq(dumvar6~=ipq_orig)=dumvar6(dumvar6~=ipq_orig); pq(dumvar5~=pq_orig)=dumvar5(dumvar5~=pq_orig);
if( ierror~=0 )
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
return;
end;
end;
end; k = fix(j0+1);
if( mu>=1 )
pq = -r.*pq;
end;
ixs = 0;
if( mu>=1 )
pq_orig=pq; ipq_orig=ipq;    [pq,ipq,dumvar3,ixs,dumvar5,dumvar6,ierror]=dxadd(pq,ipq,-xs,ixs,pq,ipq,ierror);    ipq(dumvar6~=ipq_orig)=dumvar6(dumvar6~=ipq_orig); pq(dumvar5~=pq_orig)=dumvar5(dumvar5~=pq_orig);
end;
if( ierror~=0 )
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
return;
end;
if( j==2 )
mu = fix(-mu);
end;
if( j==2 )
dmu = -dmu;
end;
else;
%
%        SERIES FOR P ( ID = 1, 3, OR 4 )
%        P(-MU,NU,X)=1./FACTORIAL(MU)*SQRT(((1.-X)/(1.+X))**MU)
%                *SUM(FROM 0 TO J0-1)A(J)*(.5-.5*X)**J
%
ipq = 0;
pq = 1.0d0;
a = 1.0d0;
ia = 0;
for i = 2 : j0;
di = i;
a = a.*y.*(di-2.0d0-nu).*(di-1.0d0+nu)./((di-1.0d0+dmu).*(di-1.0d0));
[a,ia,ierror]=dxadj(a,ia,ierror);
if( ierror~=0 )
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
return;
end;
if( a==0.0d0 )
break;
end;
pq_orig=pq; ipq_orig=ipq;    [pq,ipq,a,ia,dumvar5,dumvar6,ierror]=dxadd(pq,ipq,a,ia,pq,ipq,ierror);    ipq(dumvar6~=ipq_orig)=dumvar6(dumvar6~=ipq_orig); pq(dumvar5~=pq_orig)=dumvar5(dumvar5~=pq_orig);
if( ierror~=0 )
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
return;
end;
end;
if( mu>0 )
x2 = r;
x1 = pq;
k = fix(mu);
for i = 1 : k;
x1 = x1.*x2;
[x1,ipq,ierror]=dxadj(x1,ipq,ierror);
end; i = fix(k+1);
if( ierror~=0 )
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
return;
end;
pq = x1./factmu;
ipq = fix(ipq - ifmlv);
[pq,ipq,ierror]=dxadj(pq,ipq,ierror);
if( ierror~=0 )
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
return;
end;
end;
end;
if( j==1 )
pq2 = pq;
end;
if( j==1 )
ipq2 = fix(ipq);
end;
nu = nu + 1.0d0;
end;
k = 0;
if( nu-1.5d0>=nu1 )
k = fix(k + 1);
pqa(k) = pq2;
ipqa(k) = fix(ipq2);
if( nu>nu2+.5d0 )
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
return;
end;
end;
while( true );
pq1 = pq;
ipq1 = fix(ipq);
if( nu>=nu1+.5d0 )
k = fix(k + 1);
pqa(k) = pq;
ipqa(k) = fix(ipq);
if( nu>nu2+.5d0 )
break;
end;
end;
%
%        FORWARD NU-WISE RECURRENCE FOR F(MU,NU,X) FOR FIXED MU
%        USING
%        (NU+MU+1)*F(MU,NU,X)=(2.*NU+1)*F(MU,NU,X)-(NU-MU)*F(MU,NU-1,X)
%        WHERE F(MU,NU,X) MAY BE P(-MU,NU,X) OR IF MU IS REPLACED
%        BY -MU THEN F(MU,NU,X) MAY BE Q(MU,NU,X).
%        NOTE, IN THIS LOOP, NU=NU+1
%
x1 =(2.0d0.*nu-1.0d0)./(nu+dmu).*x.*pq1;
x2 =(nu-1.0d0-dmu)./(nu+dmu).*pq2;
[x1,ipq1,dumvar3,ipq2,pq,ipq,ierror]=dxadd(x1,ipq1,-x2,ipq2,pq,ipq,ierror);
if( ierror~=0 )
break;
end;
[pq,ipq,ierror]=dxadj(pq,ipq,ierror);
if( ierror~=0 )
break;
end;
nu = nu + 1.0d0;
pq2 = pq1;
ipq2 = fix(ipq1);
end;
%
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
end %subroutine dxpqnu
%DECK DXPSI