help-octave
[Top][All Lists]

## jybess.m (new)

 From: Eyal Doron Subject: jybess.m (new) Date: Tue, 17 Oct 1995 15:21:16 +0100 (MET)

Hi guys,

Here is a new version of jybess.m . It fixes some bugs, and adds the
possibility to calculate J_n and Y_n for complex orders as well as
complex arguments. Unfortunately, complex orders don't work well in
some parts of the complex plane :( .

This mail contains the following functions:

jybess.m - Bessel and Newmann functions for real/complex arguments and
orders.
cgamma.m - Gamma function for complex arguments. Required in jybess.m
if complex orders are needed.

and some other goodies:

cgammaln.m - log-gamma function for complex arguments.
cdigamma.m - Digamma (Psi) function for complex arguments.
sphbess.m  - Spherical j and y Bessel functions for complex arguments.

These routines should be regarded as beta, so let me know what you
think.

Eyal Doron

--------------------------jybess.m----------------------------
function [Jn,Yn] = jybess(n,x,ToDo);
% FUNCTION [Jn,Yn] = jybess(n,x,ToDo):
% Returns the Bessel and Newmann functions for complex arguments and order.
%
% n    - Maximal order, real or complex scalar.
% x    - Real or complex argument vector.
% ToDo - Optional string of command characters (case insensitive):
%        'S' - return both negative and positive orders (see below).
%        'Y' - if nargout<2, return Y_n(x) instead of J_n(x).
%        'L' - return only the highest computed order.
% -----------------------------------------
% Jn,Yn - Bessel functions of the first and (optionally) second kind.
%
% Jn, Yn are length(x) X (# orders) matrices. The orders returned are:
%   -- i*ni + [frac(nr):sign(nr):nr], or
%   -- i*ni + [-|nr|:-frac(|nr|), frac(|nr|):|nr|] if 'S'.
% where nr=real(n) and ni=imag(n).
%
% Note: Support for complex n is only partial; The algorithm doesn't work very
% well for large imag(n) and |x/n|>~1 .

% Algorithm:
% ----------
% 1) Calculate J_n by backwards recursion from a calculated starting point.
%    Minimal order is \nu=|n|-floor(|n|).
% 2) Normalization: Normalize J_\nu(x) by:
%    a) for integer n and |\Im(x)|<log(10), use
%                         1 + 2 sum_{k=1}^\infty J_{2k}(x) = 1
%    b) for integer n and |\Im(x)|>log(10), use
%                         1 + 2 sum_{k=1}^\infty (-)^k J_{2k}(x) = cos(x)
%    c) for half-integer n, use explicit formula for the spherical Bessels.
%    d) for complex n and |\Im(x)|<5, use
%       sum_{k=0}^\infty (\nu+2k)\Gamma(\nu+k)/k! J_{\nu+2k}(x) = (x/2)^\nu
%    e) for complex n and |\Im(x)|>5, use
%       e^{i\phi\nu} \sum_{k=0}^\infty [r/2*(1-e^{2i\phi})]^k/k! J_{\nu+k}(r)
%                   = J_\nu(r e^i\phi)
%       This requires a recursive call to obtain the J_{\nu+k}(r).
% 3) If Y_n, or J_n for negative non-integer n, are required, calculate
%    the Y_n(x) as follows:
% 4) Evaluate Y_\nu:
%    a) for nu=0, use
%      2/\pi(\log(x/2)+\gamma)J_0(x) - 4/\pi\sum_{k=1}^\infty (-)^k J_{2k}(x)/k
%                  = Y_0(x)
%    b) for nu=0.5, use explicit formula for the spherical Bessels.
%    c) otherwise, do the following:
%        i) Calculate J_{1-\nu}(x) and J_{2-\nu}(x) using a recursive call.
%       ii) Calculate J_{-\nu}(x) using one backwards recursion step.
%      iii) Use [J_\nu(x)\cos(\nu\pi)-J_{-\nu}(x)]/\sin(\nu\pi) = Y_\nu(x)
% 5) For real arguments, obtain the rest of the Y_{\nu+k}(x) using forward
%    recurrence, Y_{\nu+1}(x)={2\nu\over x}Y_\nu(x) - Y_{\nu-1}(x) .
% 6) Obtain the rest of the Y_{\nu+k}(x) using the Wronskian relation
%    J_{\nu+1}(x)Y_\nu(x) - J_\nu(x)Y_{\nu+1}(x) = 2/(\pi x)
%    This turns out to be much more stable for complex arguments than forward
%    recurrence technique, but fails at zeros of the J_\nu(x) .
% 7) If J for negative orders is required:
%    a) for nu=0, use J_n(x) = (-)^n J_{-n}(x)
%    b) otherwise, use 4.c.iii
% 8) If Y for negative orders is required:
%    a) for nu=0, use Y_n(x) = (-)^n Y_{-n}(x)
%    b) otherwise, use 4.c.iii
%
% For complex orders, jybess also requires the complex Gamma function cgamma.

if nargin==0
help jybess
return
end
if nargin<2
error('Not enough input parameters!');
end
if max(max(size(n)))>1
error('Max order must be a scalar!')
end
sym=0; Return_Y=(nargout==2); Return_Last=0;
if nargin>2
sym=any(ToDo==1) | any(ToDo=='s') | any(ToDo=='S');
Return_Y=Return_Y | any(ToDo=='y' | ToDo=='Y');
Return_Last=any(ToDo=='l' | ToDo=='L');
end
if sym & real(n)<0, n=-real(n)+i*imag(n); end
orig_n=n;
n_is_negative=(real(n)<0);
if n_is_negative, n=-n; end
nu=n-floor(real(n)); n=floor(real(n));
if nargin<3, sym=0; end
sym=(sym~=0);
calc_Y=(Return_Y | sym | n_is_negative);
Acc_Jnorm=(nu~=0.5);
Acc_Ynorm=calc_Y & nu==0;

xlen=prod(size(x)); x=reshape(x,1,xlen);
% x=x(:).'; xlen=length(x);
z = x==0; x = x + z;            % Temporarily replace x=0 with x=1

tiny = 16^(-250);

%                             Starting index for backwards recurrence.
c = [ 0.9507,    1.4208,   14.1850;
0.9507,    1.4208,   14.1850;
0.9507,    1.4208,   14.1850;
0.7629,    1.4222,   13.9554;
0.7369,    1.4289,   13.1756;
0.7674,    1.4311,   12.4523;
0.8216,    1.4354,   11.2121;
0.8624,    1.4397,    9.9718;
0.8798,    1.4410,    8.9217;
0.9129,    1.4360,    7.8214;
0.9438,    1.5387,    6.5014;
0.9609,    1.5216,    5.7256;
0.9693,    1.5377,    5.3565;
0.9823,    1.5220,    4.5659;
0.9934,    1.5049,    3.7902;
0.9985,    1.4831,    3.2100;
1.0006,    1.4474,    3.0239;
0.9989,    1.4137,    2.8604;
0.9959,    1.3777,    2.7760;
1.0005,    1.3500,    2.3099]';
j = 1+min(abs(n),19);
m = c(1,j).*max(3,j) + c(2,j).*(max(1,abs(x))-1) + ...
c(3,j)./(1-log(min(1,abs(x))));
m=max([m; abs(n)+10+0*m]);
m = 4*ceil(m/4);                      % Make sure rem(m,4)=0.

%                                       Prevent underflow
logtiny=log(tiny); logx=log(abs(x)/2)+1;
II=1:xlen;
while ~isempty(II)
est=-log(2*pi*m(II))/2+m(II).*(logx(II)-log(m(II)))<logtiny;
II=II(find(est)); m(II)=m(II)-4;
end;
mm = max(m(:));

% Normalization summation coefficients

if nu==0
Nrm=[1; 2*ones(ceil(mm/2),1)];
UseCos=find(abs(imag(x))>log(10));
if ~isempty(UseCos)                  % Use norm. to cos(z) instead of 1.
Nrm=Nrm(:,ones(xlen,1));
[r,c]=size(Nrm);
Nrm(2:2:r,UseCos)=-Nrm(2:2:r,UseCos);
end
elseif nu~=0.5               % nu=0.5 doesn't use normalizations
k=(1:ceil(mm/2)).';
if imag(nu)==0
Nrm=cumprod([nu*gamma(nu); (nu+2*k).*(nu+k-1)./k./(nu+2*k-2)]);
else                                    %     Needs complex gamma function
Nrm=cumprod([nu*cgamma(nu); (nu+2*k).*(nu+k-1)./k./(nu+2*k-2)]);
end
Nrm=Nrm(:,ones(xlen,1));
end

%                           Backwards recursion for the Jn
Jn=zeros(n+1,xlen);
k = mm;
bkp1 = 0*x;
bk = tiny*(m==k);
if Acc_Jnorm, t = Nrm(k/2+1,:).*bk; end
if Acc_Ynorm, y = 2*bk/k; end
sig = 1;
for k = mm-1:-1:0
bkp2 = bkp1;
bkp1 = bk;
bk = 2*(k+1+nu)*bkp1./x - bkp2 + tiny*(m==k);
if k<=n
Jn(k+1,:)=bk;
end
if (floor(k/2)*2-k==0)
if Acc_Jnorm, t = t + Nrm(k/2+1,:).*bk; end
if Acc_Ynorm & k>0
sig = -sig;
y = y + sig*2*bk/k;
end
end
end
clear Nrm bkp2 sig

if nu==0.5               % Use explicit formulae for the spherical Bessels
si=sin(x);
II=(abs(si)>0.1);
F1=find(II); F2=find(1-II);
nrm=ones(1,xlen);
SpFact=sqrt(2*x/pi);
if ~isempty(F1)
nrm(F1)=SpFact(F1).*(si(F1)./x(F1))./Jn(1,F1);
end
if ~isempty(F2)
nrm(F2)=SpFact(F2).*(si(F2)./x(F2).^2-cos(x(F2))./x(F2))./bkp1(F2);
end
elseif nu==0
nrm=ones(size(x));
if ~isempty(UseCos)
nrm(UseCos)=cos(x(UseCos));
end
abst=abs(t); nrm=nrm./abst./(t./abst);
else
II=(abs(imag(x))<5);
F1=find(II); F2=find(1-II);
nrm=ones(1,xlen);
if ~isempty(F1)
nrm(F1)=(x(F1)/2).^nu ./t(F1);
end
if ~isempty(F2)               % Use mult. Theorem to calc J_0
r=abs(x(F2)); ph=angle(x(F2)); rn=max(ceil(r));
MaxOrder=ceil(18+1.25*rn);
Jtmp=jybess(nu+MaxOrder,r);
k=(0:MaxOrder).';
Kg=cumsum([0;log(1:MaxOrder).']);          % log(k!)
v=i*ph-i*pi/2+log(imag(x(F2)));

nrm(F2)=exp(i*nu*ph).*sum(exp(k*v-Kg(:,ones(length(F2),1))).*Jtmp.')./Jn(1,F2);
end
end
Jn=Jn.*nrm(ones(n+1,1),:);           % Normalizing condition.

%                Restore results for x = 0; j0(0) = 1, jn(0) = 0 for n > 0.
if any(z)
Ii=find(z); LI=max(size(Ii));
Jn(:,Ii)=[ones(1,LI);zeros(n,LI)];
end
Jn=Jn.';

if calc_Y                                    % Also Yn
%                         First, get Y_0(x)
if n==0
J1 = bkp1.*nrm;
else
J1=Jn(:,2).';
end
if nu==0                          % Use dependence on sums of J
gam = 0.57721566490153286;     % Euler number
lx=log(x/2);
II=find(imag(x)==0 & real(x)<0);
if ~isempty(II)
lx(II)=log(abs(x(II))/2);
end
y = 2/pi*(lx + gam).*bk - 4/pi*y;
Y0 = y.*nrm;
elseif nu==0.5                    % explicit formula
Y0=-cos(x)./x.*SpFact;
else                              % Use linear relations between J and Y
if real(nu)==0
Jminus=jybess(-nu,x);
else                            % One step of the recursion for the J
Jtmp=jybess(2-nu,x);
Jminus=2*(1-nu)./x.'.*Jtmp(:,1)-Jtmp(:,2);
end
Y0=(Jn(:,1).'*cos(pi*nu)-Jminus.')/sin(pi*nu);
end

Yn=[Y0.', zeros(xlen,n)];
if n>0
Yn(:,2)=(Jn(:,2).*Yn(:,1)-2/pi./x(:).')./Jn(:,1);
end
if n>1
II=(imag(x)==0) * (imag(nu)==0);
if any(II)                        % Iterate using forward recurrence
IN=find(II);
for l=2:n
Yn(IN,l+1)=2*(nu+l-1)./x(IN).'.*Yn(IN,l) - Yn(IN,l-1);
end
end
if ~all(II)                       % Iterate using Wronskian relation
IN=find(~II);
for l=2:n
Yn(IN,l+1)=(Jn(IN,l+1).*Yn(IN,l)-2/pi./x(IN).')./Jn(IN,l);
end
end
end

if any(z), Yn(find(z),:)=-inf*ones(size(Yn(find(z),:))); end
end

if nu==0
if n_is_negative                                 % negative order
Jn(:,2:2:n+1)=-Jn(:,2:2:n+1); Jn=fliplr(Jn);
if Return_Y
Yn(:,2:2:n+1)=-Yn(:,2:2:n+1); Yn=fliplr(Yn);
end
end
if sym & n>0                                     % symmetric output
Rev=n:-1:1; Rev=2*(floor(Rev/2)*2==Rev)-1;
Jn=[Jn(:,n+1:-1:2).*Rev(ones(xlen,1),:) Jn];
if Return_Y
Yn=[Yn(:,n+1:-1:2).*Rev(ones(xlen,1),:) Yn];
end
end
else
if n_is_negative | sym
ord=nu+(0:n); co=cos(pi*ord); so=sin(pi*ord);
Jminus=Jn.*co(ones(xlen,1),:)-Yn.*so(ones(xlen,1),:);
if Return_Y
ord=nu+(0:n);co=cos(pi*ord); so=sin(-pi*ord);
Yminus=(Jminus.*co(ones(xlen,1),:)-Jn)./so(ones(xlen,1),:);
end
end
if sym
if real(nu)==0
Jn=[fliplr(Jminus(:,2:n+1)), Jn];
else
Jn=[fliplr(Jminus), Jn];
end
if Return_Y
if real(nu)==0
Yn=[fliplr(Yminus(:,2:n+1)), Yn];
else
Yn=[fliplr(Yminus), Yn];
end
end
elseif n_is_negative
Jn=fliplr(Jminus);
if Return_Y, Yn=fliplr(Yminus); end
end
end
if Return_Last
[r,c]=size(Jn);
if n_is_negative, c=1; end
Jn=Jn(:,c);
if calc_Y, Yn=Yn(:,c); end
end
if nargout<2 & Return_Y, Jn=Yn; end

-----------------------------cgamma.m------------------------
function Y=cgamma(Z);
% function Y=cgamma(Z) returns the value of the Gamma function for complex
% argument matrix Z.

% Algorithm:
% ----------
% 1) for real Z, use the builtin gamma function.
% 2) If real(Z)<0, use the reflection formula for Z -> -Z.
% 3) For |Z|>6 use an order-15 Stirling's approximation.
% 4) For the rest, use gamma(z+1)=z*gamma(z) and evaluate gamma(z+N) such that
%    |z+N|>6 .

Asymp_coef=[
+1.0;
+0.83333333333333333333e-1;
+0.34722222222222222222e-2;
-0.26813271604938271605e-2;
-0.22947209362139917695e-3;
+0.78403922172006662747e-3;
+0.69728137583658577743e-4;
-0.59216643735369388286e-3;
-0.51717909082605921934e-4;
+0.83949872067208727999e-3;
+0.72048954160200105591e-4;
-0.19144384985654775265e-2;
-0.16251626278391581690e-3;
+0.64033628338080697948e-2;
+0.54016476789260451518e-3;
];

[r,c]=size(Z); Z=Z(:).'; Y=zeros(size(Z));

n_real=imag(Z)==0;
n_complex=~n_real;
if any(n_real)
IN=find(n_real);
Y(IN)=gamma(Z(IN));
end
if any(n_complex)
IN=find(n_complex);
z=Z(IN);
FN=find(real(z)<0);
z(FN)=-z(FN);
N=36-abs(z).^2;
II=find(N>0); I2=find(N<=0);
N(II)=ceil(sqrt(N(II)+real(z(II)).^2)-real(z(II)));
N(I2)=0*N(I2);
z(II)=z(II)+N(II);

y=sqrt(2*pi./z).*z.^z.*exp(-z).*(exp(-(0:length(Asymp_coef)-1)'*log(z)).'*Asymp_coef);
for l=II
y(l)=y(l)/prod(z(l)-(1:N(l)));    % Take back the added N
end
z(II)=z(II)-N(II);
y(FN)=-pi./(z(FN).*sin(pi*z(FN)).*y(FN));
Y(IN)=y;
end

Y=reshape(Y,r,c);

----------------------------cgammaln.m--------------------------------
function Y=cgammaln(Z);
% function Y=cgammaln(Z) returns the value of the log-gamma function for
% complex argument matrix Z.

% Algorithm:
% ----------
% 1) for real Z, use the builtin log-gamma function.
% 2) If real(Z)<0, use the reflection formula for Z -> -Z.
% 3) For |Z|>6 use an order-15 Stirling's approximation.
% 4) Near the zeros |Z+/-1|<=0.1 and |Z+/-2|<=0.1 use a series expansion.
% 5) For the rest, use gamma(z+1)=z*gamma(z) and evaluate gamma(z+N) such that
%    |z+N|>6 .

One_coef=[
-0.57721566490153286061;
+0.82246703342411321826;
-0.40068563438653142846;
+0.27058080842778454789;
-0.20738555102867398526;
+0.16955717699740818997;
-0.14404989676884611811;
+0.12550966952474304243;
-0.11133426586956469049;
+0.10009945751278180854;
-0.90954017145829042236e-1;
+0.83353840546109004037e-1;
-0.76932516411352191469e-1;
+0.71432946295361336071e-1;
%    -0.66668705882420468034e-1;
%    +0.62500955141213040754e-1;
%    -0.58823978658684582341e-1;
%    +0.55555767627403611114e-1;
%    -0.52631679379616660732e-1;
];
Asymp_coef=[
0.08333333333333333;         %  1/12
-0.002777777777777778;        % -1/360
0.0007936507936507937;       %  1/1260
-0.0005952380952380953;       % -1/1680
0.0008417508417508417;       %  1/1188
-0.001917526917526918;        % -691/360360
0.00641025641025641;         %  1/156
];

[r,c]=size(Z); Z=Z(:).'; Y=zeros(size(Z));

n_real=imag(Z)==0;
n_one=~n_real & (abs(Z-1)<=0.1 | abs(Z+1)<=0.1);
n_two=~n_real & (abs(Z-2)<=0.1 | abs(Z+2)<=0.1);
n_complex=~(n_real | n_one | n_two);
if any(n_real)
IN=find(n_real);
if exist('OCTAVE_VERSION')
Y(IN)=lgamma(Z(IN));
else
Y(IN)=gammaln(Z(IN));
end
end
if any(n_one)
IN=find(n_one);
z=Z(IN);
FN=find(abs(z+1)<=0.1);
z(FN)=-z(FN);
y=exp((1:length(One_coef))'*log(z-1)).'*One_coef;
y(FN)=i*pi+log(pi)-log(z)-log(sin(pi*z(FN)))-y(FN);
Y(IN)=y;
end
if any(n_two)
IN=find(n_two);
z=Z(IN);
FN=find(abs(z+2)<=0.1);
z(FN)=-z(FN);
y=exp((1:length(One_coef))'*log(z-2)).'*One_coef+log(z-1);
y(FN)=i*pi+log(pi)-log(z)-log(sin(pi*z(FN)))-y(FN);
Y(IN)=y;
end
if any(n_complex)
IN=find(n_complex);
z=Z(IN);
FN=find(real(z)<0);
z(FN)=-z(FN);
N=36-abs(z).^2;
II=find(N>0); I2=find(N<=0);
N(II)=ceil(sqrt(N(II)+real(z(II)).^2)-real(z(II)));
N(I2)=0*N(I2);
z(II)=z(II)+N(II);                     % Make it asymptotic
lz=log(z);

y=(log(2*pi)-lz)/2+(lz-1).*z+exp(-(1:2:2*length(Asymp_coef)-1)'*lz).'*Asymp_coef;
for l=II
y(l)=y(l)-sum(log(z(l)-(1:N(l))));    % Take back the added N
end
z(II)=z(II)-N(II);
y(FN)=log(-pi./(z(FN).*sin(pi*z(FN))))-y(FN);
Y(IN)=y;
end

Y=reshape(Y,r,c);

---------------------------cdigamma.m------------------------------
function Y=cdigamma(Z);
% function Y=cdigamma(Z) returns the value of the digamma function for
% complex argument matrix Z.

% Algorithm:
% ----------
% 1) If real(Z)<0, use the reflection formula for Z -> -Z.
% 2) For |Z|>6 use an order-15 Stirling's approximation.
% 3) For the rest, use Psi(z+1)=1/z+Psi(z) and evaluate Psi(z+N) such that
%    |z+N|>6 .

Asymp_coef=[
-0.08333333333333333;    %    -1/12
0.008333333333333333;   %    +1/120
-0.003968253968253968;   %    -1/252
0.004166666666666667;   %    +1/240
-0.007575757575757576;   %    -1/132
0.0210927960927961;     %    +691/32760
-0.08333333333333333;    %    -1/12
];

[r,c]=size(Z); Z=Z(:).'; Y=zeros(size(Z));

FN=find(real(Z)<0);
Z(FN)=-Z(FN);
N=36-abs(Z).^2;
II=find(N>0); I2=find(N<=0);
N(II)=ceil(sqrt(N(II)+real(Z(II)).^2)-real(Z(II)));
N(I2)=0*N(I2);
Z(II)=Z(II)+N(II);                     % Make it asymptotic
lZ=log(Z);
Y=lZ-1/Z/2+exp(-(2:2:2*length(Asymp_coef))'*lZ).'*Asymp_coef;
for l=II
Y(l)=Y(l)-sum(1./(Z(l)-(1:N(l))));    % Take back the added N
end
Z(II)=Z(II)-N(II);
II=(imag(Z(FN))==0) & (abs(ceil(Z(FN)-0.5)-(Z(FN)-0.5))<=eps);
F1=FN(find(II)); F2=FN(find(~II));
Y(F1)=Y(F1)+1./Z(F1);
Y(F2)=Y(F2)+1./Z(F2)+pi./tan(pi*Z(F2));

Y=reshape(Y,r,c);

----------------------------------sphbess.m----------------------
function [jn,yn] = sphbess(n,X,ToDo);
% [jn,yn]=SPHBESS(n,X,ToDo):
% Returns the spherical Bessel and Newmann functions for complex arguments.
%
% n    - Maximal order, integer scalar.
% x    - Real or complex argument vector.
% ToDo - Optional string of command characters (case not important):
%        'S' - return orders -n:n rather than 0:n.
%        'Y' - if nargout<2, return y_n(x) instead of j_n(x).
%        'L' - return only the highest computed order.
% -----------------------------------------
% jn,yn - Spherical Bessel functions of the 1st and (optionally) 2nd kind,
%         length(X) x (n+1) or length(X) x (2n+1) matrices containing orders
%         0:n for n>=0, n:0 for n<0, or -|n|:|n| for 'S'.

% Algorithm:
% ----------
% sphbess calls jybess(n+1/2,X), and multiplies by sqrt(pi/(2X)). The point
% X=0 is taken care of separately.

if nargin<1
help sphbess
return
end
if nargin<2, error('Not enough input arguments!'), end
if max(size(n))>1 | any(ceil(n)~=n),
error('"n" must be an integer scalar!')
end
if nargin<3, ToDo=''; end
if ~isstr(ToDo), error('"ToDo" must be a string!'), end

NewToDo=''; do_J=1; do_Y=(nargout>1); sym=0;
if any(ToDo=='Y') | any(ToDo=='y')
do_J=0; do_Y=1; NewToDo=[NewToDo, 'Y'];
end
OnlyLast=any(ToDo=='L') | any(ToDo=='l');
if any(ToDo=='S') | any(ToDo=='s')
if n~=0
sym=1; NewToDo=[NewToDo, 'S']; n=abs(n);
end
Orders=-n:n;
else
Orders=0:n;
end
if n<0
n_negative=1; NewToDo=[NewToDo, 'S']; n=-n;
end

X=X(:);

%           Take care of X=0
z=(X==0); X=X+z; z=find(z);

%                  Calculate Bessel functions
if do_J & do_Y
[Jn,Yn]=jybess(n+1/2,X,NewToDo);
elseif do_J
Jn=jybess(n+1/2,X,NewToDo);
else
Yn=jybess(n+1/2,X,NewToDo);
end

fac=sqrt(pi./(2*X));
ZeroOrder=find(Orders==0);
IntOrder=find(Orders~=0);
if do_J
if n_negative, Jn=Jn(:,2:2+n); end
if sym, Jn=Jn(:,2:2*n+2); end
[r,c]=size(Jn);
jn=Jn.*fac(:,ones(1,c));
if ~isempty(z)
jn(z,ZeroOrder)=ones(size(jn(z,ZeroOrder)));
if ~isempty(IntOrder), jn(z,IntOrder)=zeros(size(jn(z,IntOrder))); end
end
if OnlyLast
if n_negative
jn=jn(:,1);
else
jn=jn(:,c);
end
end
end
if do_Y
if n_negative, Yn=Yn(:,2:2+n); end
if sym, Yn=Yn(:,2:2*n+2); end
[r,c]=size(Jn);
yn=Yn.*fac(:,ones(1,c));
if ~isempty(z), yn(z,:)=yn(z,:)-Inf; end
if OnlyLast
if n_negative
yn=yn(:,1);
else
yn=yn(:,c);
end
end
if nargout<2, jn=yn; end
end