[dnu1,nudiff,mu1,mu2,theta,p,q,r,ip,iq,ir,c1,ic1,c2,ic2,ierror]=xcsrt(dnu1,nudiff,mu1,mu2,theta,p,q,r,ip,iq,ir,c1,ic1,c2,ic2,ierror);

function [dnu1,nudiff,mu1,mu2,theta,p,q,r,ip,iq,ir,c1,ic1,c2,ic2,ierror]=xcsrt(dnu1,nudiff,mu1,mu2,theta,p,q,r,ip,iq,ir,c1,ic1,c2,ic2,ierror);
persistent dmu dmu1 i ix1 ix2 j k l lm1 mu nu sx x1 x2 ; 

if isempty(i), i=0; end;
if isempty(ix1), ix1=0; end;
if isempty(ix2), ix2=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(lm1), lm1=0; end;
if isempty(mu), mu=0; end;
if isempty(dmu), dmu=0; end;
if isempty(dmu1), dmu1=0; end;
if isempty(nu), nu=0; end;
if isempty(sx), sx=0; end;
if isempty(x1), x1=0; end;
if isempty(x2), x2=0; end;
p_shape=size(p);p=reshape(p,1,[]);
ip_shape=size(ip);ip=reshape(ip,1,[]);
q_shape=size(q);q=reshape(q,1,[]);
iq_shape=size(iq);iq=reshape(iq,1,[]);
r_shape=size(r);r=reshape(r,1,[]);
ir_shape=size(ir);ir=reshape(ir,1,[]);
c1_shape=size(c1);c1=reshape(c1,1,[]);
ic1_shape=size(ic1);ic1=reshape(ic1,1,[]);
c2_shape=size(c2);c2=reshape(c2,1,[]);
ic2_shape=size(ic2);ic2=reshape(ic2,1,[]);
ierror = 0;
l = fix(nudiff +(mu2-mu1) + 1);
lm1 = fix(l - 1);
for i = 1 : l;
[p(i),ip(i),ierror]=xadj(p(i),ip(i),ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
[q(i),iq(i),ierror]=xadj(q(i),iq(i),ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
[r(i),ir(i),ierror]=xadj(r(i),ir(i),ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
end; i = fix(l+1);
if( mu2>mu1 )
sx = sin(theta);
for i = 1 : lm1;
c1(i) = 0.;
c2(i) = 0.;
mu = fix(mu1 + i - 1);
dmu = mu;
x1 = p(i+1).*q(i);
ix1 = fix(ip(i+1) + iq(i));
[x1,ix1,ierror]=xadj(x1,ix1,ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
x2 = p(i).*q(i+1);
ix2 = fix(ip(i) + iq(i+1));
[x2,ix2,ierror]=xadj(x2,ix2,ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
x1 =(dmu+dnu1+1.).*(dmu-dnu1).*x1;
[x1,ix1,dumvar3,ix2,c1(i),ic1(i),ierror]=xadd(x1,ix1,-x2,ix2,c1(i),ic1(i),ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
c1(i) = sx.*c1(i).*(-1).^mu;
[c1(i),ic1(i),ierror]=xadj(c1(i),ic1(i),ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
if( dmu<dnu1+1. || rem(dnu1,1.)~=0. )
x1 = r(i+1).*q(i);
ix1 = fix(ir(i+1) + iq(i));
[x1,ix1,ierror]=xadj(x1,ix1,ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
x2 = r(i).*q(i+1);
ix2 = fix(ir(i) + iq(i+1));
[x2,ix2,ierror]=xadj(x2,ix2,ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
[x1,ix1,dumvar3,ix2,c2(i),ic2(i),ierror]=xadd(x1,ix1,-x2,ix2,c2(i),ic2(i),ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
c2(i) = c2(i).*sx;
k = fix(2.*mu);
if( k>0 )
for j = 1 : k;
c2(i) = c2(i)./(dnu1+dmu+1.-j);
[c2(i),ic2(i),ierror]=xadj(c2(i),ic2(i),ierror);
end; j = fix(k+1);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
end;
end;
end; i = fix(lm1+1);
else;
dmu1 = mu1;
for i = 1 : lm1;
c1(i) = 0.;
c2(i) = 0.;
nu = dnu1 + i - 1.;
x1 = p(i+1).*q(i);
ix1 = fix(ip(i+1) + iq(i));
[x1,ix1,ierror]=xadj(x1,ix1,ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
x2 = p(i).*q(i+1);
ix2 = fix(ip(i) + iq(i+1));
[x2,ix2,ierror]=xadj(x2,ix2,ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
x1 =(dmu1+nu+1.).*x1;
x2 =(dmu1-nu-1.).*x2;
[x1,ix1,x2,ix2,c1(i),ic1(i),ierror]=xadd(x1,ix1,x2,ix2,c1(i),ic1(i),ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
[c1(i),ic1(i),ierror]=xadj(c1(i),ic1(i),ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
c1(i) = c1(i).*(-1).^mu1;
if( dmu1>=nu+1. && rem(nu,1.)==0. )
c2(i) = 0.;
ic2(i) = 0;
else;
x1 = r(i+1).*q(i);
ix1 = fix(ir(i+1) + iq(i));
[x1,ix1,ierror]=xadj(x1,ix1,ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
x2 = r(i).*q(i+1);
ix2 = fix(ir(i) + iq(i+1));
[x2,ix2,ierror]=xadj(x2,ix2,ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
[x1,ix1,dumvar3,ix2,c2(i),ic2(i),ierror]=xadd(x1,ix1,-x2,ix2,c2(i),ic2(i),ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
k = fix(2.*mu1 - 1);
for j = 1 : k;
if( k>0 )
c2(i) = c2(i)./(nu+dmu1+1.-j);
end;
[c2(i),ic2(i),ierror]=xadj(c2(i),ic2(i),ierror);
end; j = fix(k+1);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
if( k<=0 )
c2(i) =(nu+1.).*c2(i);
end;
end;
end; i = fix(lm1+1);
end;
for i = 1 : lm1;
[c1(i),ic1(i),ierror]=xred(c1(i),ic1(i),ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
[c2(i),ic2(i),ierror]=xred(c2(i),ic2(i),ierror);
if( ierror~=0 )
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
return;
end;
end; i = fix(lm1+1);
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
iq_shape=zeros(iq_shape);iq_shape(:)=iq(1:numel(iq_shape));iq=iq_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
c1_shape=zeros(c1_shape);c1_shape(:)=c1(1:numel(c1_shape));c1=c1_shape;
ic1_shape=zeros(ic1_shape);ic1_shape(:)=ic1(1:numel(ic1_shape));ic1=ic1_shape;
c2_shape=zeros(c2_shape);c2_shape(:)=c2(1:numel(c2_shape));c2=c2_shape;
ic2_shape=zeros(ic2_shape);ic2_shape(:)=ic2(1:numel(ic2_shape));ic2=ic2_shape;
end %subroutine xcsrt