[nperod,n,m,a,bb,c,idimq,q,b,b2,b3,w,w2,w3,d,tcos,p]=postg2(nperod,n,m,a,bb,c,idimq,q,b,b2,b3,w,w2,w3,d,tcos,p);

function [nperod,n,m,a,bb,c,idimq,q,b,b2,b3,w,w2,w3,d,tcos,p]=postg2(nperod,n,m,a,bb,c,idimq,q,b,b2,b3,w,w2,w3,d,tcos,p);
persistent fi fnum fnum2 i i2r i2rby2 ii ijump ip ipstor j jm1 jm2 jm3 jp1 jp2 jp3 jr jstart jstep jstop k k1 k2 k3 k4 kr lr mr nlast nlastp np nr nrod nrodpr t ; 

if isempty(fi), fi=0; end;
if isempty(fnum), fnum=0; end;
if isempty(fnum2), fnum2=0; end;
if isempty(t), t=0; end;
if isempty(i), i=0; end;
if isempty(i2r), i2r=0; end;
if isempty(i2rby2), i2rby2=0; end;
if isempty(ii), ii=0; end;
if isempty(ijump), ijump=0; end;
if isempty(ip), ip=0; end;
if isempty(ipstor), ipstor=0; end;
if isempty(j), j=0; end;
if isempty(jm1), jm1=0; end;
if isempty(jm2), jm2=0; end;
if isempty(jm3), jm3=0; end;
if isempty(jp1), jp1=0; end;
if isempty(jp2), jp2=0; end;
if isempty(jp3), jp3=0; end;
if isempty(jr), jr=0; end;
if isempty(jstart), jstart=0; end;
if isempty(jstep), jstep=0; end;
if isempty(jstop), jstop=0; end;
if isempty(k), k=zeros(1,4); end;
if isempty(k1), k1=0; end;
if isempty(k2), k2=0; end;
if isempty(k3), k3=0; end;
if isempty(k4), k4=0; end;
if isempty(kr), kr=0; end;
if isempty(lr), lr=0; end;
if isempty(mr), mr=0; end;
if isempty(nlast), nlast=0; end;
if isempty(nlastp), nlastp=0; end;
if isempty(np), np=0; end;
if isempty(nr), nr=0; end;
if isempty(nrod), nrod=0; end;
if isempty(nrodpr), nrodpr=0; end;
%***BEGIN PROLOGUE  POSTG2
%***SUBSIDIARY
%***PURPOSE  Subsidiary to POISTG
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (POSTG2-S)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     subroutine to solve Poisson's equation on a staggered grid.
%
%***SEE ALSO  POISTG
%***ROUTINES CALLED  COSGEN, S1MERG, TRI3, TRIX
%***REVISION HISTORY  (YYMMDD)
%   801001  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%   920130  Modified to use merge routine S1MERG rather than deleted
%           routine MERGE.  (WRB)
%***end PROLOGUE  POSTG2
%
a_shape=size(a);a=reshape(a,1,[]);
bb_shape=size(bb);bb=reshape(bb,1,[]);
c_shape=size(c);c=reshape(c,1,[]);
q_shape=size(q);q=reshape([q(:).',zeros(1,ceil(numel(q)./prod([idimq])).*prod([idimq])-numel(q))],idimq,[]);
b_shape=size(b);b=reshape(b,1,[]);
b2_shape=size(b2);b2=reshape(b2,1,[]);
b3_shape=size(b3);b3=reshape(b3,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
w2_shape=size(w2);w2=reshape(w2,1,[]);
w3_shape=size(w3);w3=reshape(w3,1,[]);
d_shape=size(d);d=reshape(d,1,[]);
tcos_shape=size(tcos);tcos=reshape(tcos,1,[]);
p_shape=size(p);p=reshape(p,1,[]);
% equivalence(k(1),k1) ::;
% equivalence(k(2),k2) ::;
% equivalence(k(3),k3) ::;
% equivalence(k(4),k4) ::;
%***FIRST EXECUTABLE STATEMENT  POSTG2
np = fix(nperod);
fnum = (0.5.*(fix(np./3)));
fnum2 = (0.5.*(fix(np./2)));
mr = fix(m);
ip = fix(-mr);
ipstor = 0;
i2r = 1;
jr = 2;
nr = fix(n);
nlast = fix(n);
kr = 1;
lr = 0;
if( nr>3 )
while( true );
jr = fix(2.*i2r);
nrod = 1;
if(((fix(nr./2)).*2)==nr )
nrod = 0;
end;
jstart = 1;
jstop = fix(nlast - jr);
if( nrod==0 )
jstop = fix(jstop - i2r);
end;
i2rby2 = fix(fix(i2r./2));
if( jstop>=jstart )
%
%     REGULAR REDUCTION.
%
ijump = 1;
for j = jstart : jr: jstop ;
jp1 = fix(j + i2rby2);
jp2 = fix(j + i2r);
jp3 = fix(jp2 + i2rby2);
jm1 = fix(j - i2rby2);
jm2 = fix(j - i2r);
jm3 = fix(jm2 - i2rby2);
if( j==1 )
[i2r,dumvar2,fnum,dumvar4,tcos]=cosgen(i2r,1,fnum,0.5,tcos);
if( i2r==1 )
for i = 1 : mr;
b(i) = q(i,1);
q(i,1) = q(i,2);
end; i = fix(mr+1);
else;
for i = 1 : mr;
b(i) = q(i,1) + 0.5.*(q(i,jp2)-q(i,jp1)-q(i,jp3));
q(i,1) = q(i,jp2) + q(i,1) - q(i,jp1);
end; i = fix(mr+1);
end;
else;
if( ijump~=2 )
ijump = 2;
[i2r,dumvar2,dumvar3,dumvar4,tcos]=cosgen(i2r,1,0.5,0.0,tcos);
end;
if( i2r==1 )
for i = 1 : mr;
b(i) = 2..*q(i,j);
q(i,j) = q(i,jm2) + q(i,jp2);
end; i = fix(mr+1);
else;
for i = 1 : mr;
fi = q(i,j);
q(i,j) = q(i,j) - q(i,jm1) - q(i,jp1) + q(i,jm2) + q(i,jp2);
b(i) = fi + q(i,j) - q(i,jm3) - q(i,jp3);
end; i = fix(mr+1);
end;
end;
[i2r,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(i2r,0,mr,a,bb,c,b,tcos,d,w);
for i = 1 : mr;
%
%     end OF REDUCTION FOR REGULAR UNKNOWNS.
%
q(i,j) = q(i,j) + b(i);
end; i = fix(mr+1);
end; j = fix(jstop +1);
%
%     BEGIN SPECIAL REDUCTION FOR LAST UNKNOWN.
%
j = fix(jstop + jr);
else;
j = fix(jr);
end;
nlast = fix(j);
jm1 = fix(j - i2rby2);
jm2 = fix(j - i2r);
jm3 = fix(jm2 - i2rby2);
if( nrod==0 )
%
%     EVEN NUMBER OF UNKNOWNS
%
jp1 = fix(j + i2rby2);
jp2 = fix(j + i2r);
if( i2r==1 )
for i = 1 : mr;
b(i) = q(i,j);
end; i = fix(mr+1);
tcos(1) = 0.;
[dumvar1,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(1,0,mr,a,bb,c,b,tcos,d,w);
ip = 0;
ipstor = fix(mr);
for i = 1 : mr;
p(i) = b(i);
b(i) = b(i) + q(i,n);
end; i = fix(mr+1);
tcos(1) = -1. + (2.*(fix(np./2)));
tcos(2) = 0.;
[dumvar1,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(1,1,mr,a,bb,c,b,tcos,d,w);
for i = 1 : mr;
q(i,j) = q(i,jm2) + p(i) + b(i);
end; i = fix(mr+1);
else;
for i = 1 : mr;
b(i) = q(i,j) + .5.*(q(i,jm2)-q(i,jm1)-q(i,jm3));
end; i = fix(mr+1);
if( nrodpr==0 )
for i = 1 : mr;
ii = fix(ip + i);
b(i) = b(i) + p(ii);
end; i = fix(mr+1);
else;
for i = 1 : mr;
b(i) = b(i) + q(i,jp2) - q(i,jp1);
end; i = fix(mr+1);
end;
[i2r,dumvar2,dumvar3,dumvar4,tcos]=cosgen(i2r,1,0.5,0.0,tcos);
[i2r,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(i2r,0,mr,a,bb,c,b,tcos,d,w);
ip = fix(ip + mr);
ipstor = fix(max(ipstor,ip+mr));
for i = 1 : mr;
ii = fix(ip + i);
p(ii) = b(i) + .5.*(q(i,j)-q(i,jm1)-q(i,jp1));
b(i) = p(ii) + q(i,jp2);
end; i = fix(mr+1);
if( lr==0 )
for i = 1 : i2r;
ii = fix(kr + i);
tcos(ii) = tcos(i);
end; i = fix(i2r+1);
else;
[lr,dumvar2,fnum2,dumvar4,tcos(sub2ind(size(tcos),max(i2r+1,1)):end)]=cosgen(lr,1,fnum2,0.5,tcos(sub2ind(size(tcos),max(i2r+1,1)):end));
i2r_orig=i2r;    [tcos,dumvar2,i2r,dumvar4,lr,kr]=s1merg(tcos,0,i2r,i2r,lr,kr);    i2r(dumvar4~=i2r_orig)=dumvar4(dumvar4~=i2r_orig);
end;
[kr,dumvar2,fnum2,dumvar4,tcos]=cosgen(kr,1,fnum2,0.5,tcos);
kr_orig=kr;    [kr,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(kr,kr,mr,a,bb,c,b,tcos,d,w);    kr(dumvar2~=kr_orig)=dumvar2(dumvar2~=kr_orig);
for i = 1 : mr;
ii = fix(ip + i);
q(i,j) = q(i,jm2) + p(ii) + b(i);
end; i = fix(mr+1);
end;
lr = fix(kr);
kr = fix(kr + jr);
else;
%
%     ODD NUMBER OF UNKNOWNS
%
if( i2r==1 )
for i = 1 : mr;
b(i) = q(i,j);
q(i,j) = q(i,jm2);
end; i = fix(mr+1);
else;
for i = 1 : mr;
b(i) = q(i,j) + .5.*(q(i,jm2)-q(i,jm1)-q(i,jm3));
end; i = fix(mr+1);
if( nrodpr==0 )
for i = 1 : mr;
ii = fix(ip + i);
q(i,j) = q(i,jm2) + p(ii);
end; i = fix(mr+1);
ip = fix(ip - mr);
else;
for i = 1 : mr;
q(i,j) = q(i,j) - q(i,jm1) + q(i,jm2);
end; i = fix(mr+1);
end;
if( lr~=0 )
[lr,dumvar2,fnum2,dumvar4,tcos(sub2ind(size(tcos),max(kr+1,1)):end)]=cosgen(lr,1,fnum2,0.5,tcos(sub2ind(size(tcos),max(kr+1,1)):end));
end;
end;
[kr,dumvar2,fnum2,dumvar4,tcos]=cosgen(kr,1,fnum2,0.5,tcos);
[kr,lr,mr,a,bb,c,b,tcos,d,w]=trix(kr,lr,mr,a,bb,c,b,tcos,d,w);
for i = 1 : mr;
q(i,j) = q(i,j) + b(i);
end; i = fix(mr+1);
kr = fix(kr + i2r);
end;
nr =fix(fix((nlast-1)./jr) + 1);
if( nr<=3 )
break;
end;
i2r = fix(jr);
nrodpr = fix(nrod);
end;
end;
%
%      BEGIN SOLUTION
%
j = fix(1 + jr);
jm1 = fix(j - i2r);
jp1 = fix(j + i2r);
jm2 = fix(nlast - i2r);
if( nr==2 )
%
%     CASE OF GENERAL N AND NR = 2 .
%
for i = 1 : mr;
ii = fix(ip + i);
b3(i) = 0.;
b(i) = q(i,1) + p(ii);
q(i,1) = q(i,1) - q(i,jm1);
b2(i) = q(i,1) + q(i,nlast);
end; i = fix(mr+1);
k1 = fix(kr + jr);
k2 = fix(k1 + jr);
[dumvar1,dumvar2,dumvar3,dumvar4,tcos(sub2ind(size(tcos),max(k1+1,1)):end)]=cosgen(jr-1,1,0.0,1.0,tcos(sub2ind(size(tcos),max(k1+1,1)):end));
if( np==2 )
[dumvar1,dumvar2,dumvar3,dumvar4,tcos(sub2ind(size(tcos),max(k2,1)):end)]=cosgen(kr+1,1,0.5,0.0,tcos(sub2ind(size(tcos),max(k2,1)):end));
else;
tcos(k2) = 2.*np - 4;
[kr,dumvar2,dumvar3,dumvar4,tcos(sub2ind(size(tcos),max(k2+1,1)):end)]=cosgen(kr,1,0.0,1.0,tcos(sub2ind(size(tcos),max(k2+1,1)):end));
end;
k4 = fix(1 - fix(np./3));
[tcos,k1]=s1merg(tcos,k1,jr-k4,k2-k4,kr+k4,0);
if( np==3 )
k1 = fix(k1 - 1);
end;
k2 = fix(kr);
[kr,dumvar2,fnum2,dumvar4,tcos(sub2ind(size(tcos),max(k1+1,1)):end)]=cosgen(kr,1,fnum2,0.5,tcos(sub2ind(size(tcos),max(k1+1,1)):end));
k4 = fix(k1 + kr);
[lr,dumvar2,fnum2,dumvar4,tcos(sub2ind(size(tcos),max(k4+1,1)):end)]=cosgen(lr,1,fnum2,0.5,tcos(sub2ind(size(tcos),max(k4+1,1)):end));
k3 = fix(lr);
k4 = 0;
[mr,a,bb,c,k,b,b2,b3,tcos,d,w,w2,w3]=tri3(mr,a,bb,c,k,b,b2,b3,tcos,d,w,w2,w3);
for i = 1 : mr;
b(i) = b(i) + b2(i);
end; i = fix(mr+1);
if( np==3 )
tcos(1) = 2.;
[dumvar1,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(1,0,mr,a,bb,c,b,tcos,d,w);
end;
for i = 1 : mr;
q(i,1) = q(i,1) + b(i);
end; i = fix(mr+1);
elseif( lr~=0 ) ;
%
%     CASE OF GENERAL N WITH NR = 3 .
%
for i = 1 : mr;
b(i) = q(i,1) - q(i,jm1) + q(i,j);
end; i = fix(mr+1);
if( nrod==0 )
for i = 1 : mr;
ii = fix(ip + i);
b(i) = b(i) + p(ii);
end; i = fix(mr+1);
else;
for i = 1 : mr;
b(i) = b(i) + q(i,nlast) - q(i,jm2);
end; i = fix(mr+1);
end;
for i = 1 : mr;
t = .5.*(q(i,j)-q(i,jm1)-q(i,jp1));
q(i,j) = t;
b2(i) = q(i,nlast) + t;
b3(i) = q(i,1) + t;
end; i = fix(mr+1);
k1 = fix(kr + 2.*jr);
[dumvar1,dumvar2,dumvar3,dumvar4,tcos(sub2ind(size(tcos),max(k1+1,1)):end)]=cosgen(jr-1,1,0.0,1.0,tcos(sub2ind(size(tcos),max(k1+1,1)):end));
k2 = fix(k1 + jr);
tcos(k2) = 2.*np - 4;
k4 =fix((np-1).*(3-np));
k3 = fix(k2 + 1 - k4);
[dumvar1,dumvar2,dumvar3,dumvar4,tcos(sub2ind(size(tcos),max(k3,1)):end)]=cosgen(kr+jr+k4,1,k4./2.,1.-k4,tcos(sub2ind(size(tcos),max(k3,1)):end));
k4 = fix(1 - fix(np./3));
[tcos,k1]=s1merg(tcos,k1,jr-k4,k2-k4,kr+jr+k4,0);
if( np==3 )
k1 = fix(k1 - 1);
end;
k2 = fix(kr + jr);
k4 = fix(k1 + k2);
[kr,dumvar2,fnum2,dumvar4,tcos(sub2ind(size(tcos),max(k4+1,1)):end)]=cosgen(kr,1,fnum2,0.5,tcos(sub2ind(size(tcos),max(k4+1,1)):end));
k3 = fix(k4 + kr);
[jr,dumvar2,fnum,dumvar4,tcos(sub2ind(size(tcos),max(k3+1,1)):end)]=cosgen(jr,1,fnum,0.5,tcos(sub2ind(size(tcos),max(k3+1,1)):end));
[tcos,k4,kr,k3,jr,k1]=s1merg(tcos,k4,kr,k3,jr,k1);
k4 = fix(k3 + jr);
[lr,dumvar2,fnum2,dumvar4,tcos(sub2ind(size(tcos),max(k4+1,1)):end)]=cosgen(lr,1,fnum2,0.5,tcos(sub2ind(size(tcos),max(k4+1,1)):end));
[tcos,k3,jr,k4,lr]=s1merg(tcos,k3,jr,k4,lr,k1+k2);
[kr,dumvar2,fnum2,dumvar4,tcos(sub2ind(size(tcos),max(k3+1,1)):end)]=cosgen(kr,1,fnum2,0.5,tcos(sub2ind(size(tcos),max(k3+1,1)):end));
k3 = fix(kr);
k4 = fix(kr);
[mr,a,bb,c,k,b,b2,b3,tcos,d,w,w2,w3]=tri3(mr,a,bb,c,k,b,b2,b3,tcos,d,w,w2,w3);
for i = 1 : mr;
b(i) = b(i) + b2(i) + b3(i);
end; i = fix(mr+1);
if( np==3 )
tcos(1) = 2.;
[dumvar1,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(1,0,mr,a,bb,c,b,tcos,d,w);
end;
for i = 1 : mr;
q(i,j) = q(i,j) + b(i);
b(i) = q(i,1) + q(i,j);
end; i = fix(mr+1);
[jr,dumvar2,fnum,dumvar4,tcos]=cosgen(jr,1,fnum,0.5,tcos);
[jr,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(jr,0,mr,a,bb,c,b,tcos,d,w);
if( jr==1 )
for i = 1 : mr;
q(i,1) = b(i);
end; i = fix(mr+1);
else;
for i = 1 : mr;
q(i,1) = q(i,1) - q(i,jm1) + b(i);
end; i = fix(mr+1);
end;
elseif( n==3 ) ;
%
%     CASE N = 3.
%
if( np==2 )
for i = 1 : mr;
b(i) = q(i,2);
b2(i) = q(i,3);
b3(i) = q(i,1);
end; i = fix(mr+1);
[dumvar1,dumvar2,dumvar3,dumvar4,tcos]=cosgen(3,1,0.5,0.0,tcos);
tcos(4) = -1.;
tcos(5) = 1.;
tcos(6) = -1.;
tcos(7) = 1.;
k1 = 3;
k2 = 2;
k3 = 1;
k4 = 1;
else;
for i = 1 : mr;
b(i) = q(i,2);
b2(i) = q(i,1) + q(i,3);
b3(i) = 0.;
end; i = fix(mr+1);
if( np==1 || np==2 )
tcos(1) = -2.;
tcos(2) = 1.;
tcos(3) = -1.;
k1 = 2;
else;
tcos(1) = -1.;
tcos(2) = 1.;
k1 = 1;
end;
k2 = 1;
k3 = 0;
k4 = 0;
end;
[mr,a,bb,c,k,b,b2,b3,tcos,d,w,w2,w3]=tri3(mr,a,bb,c,k,b,b2,b3,tcos,d,w,w2,w3);
for i = 1 : mr;
b(i) = b(i) + b2(i) + b3(i);
end; i = fix(mr+1);
if( np~=1 && np~=2 )
tcos(1) = 2.;
[dumvar1,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(1,0,mr,a,bb,c,b,tcos,d,w);
end;
for i = 1 : mr;
q(i,2) = b(i);
b(i) = q(i,1) + b(i);
end; i = fix(mr+1);
tcos(1) = -1. + 4..*fnum;
[dumvar1,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(1,0,mr,a,bb,c,b,tcos,d,w);
for i = 1 : mr;
q(i,1) = b(i);
end; i = fix(mr+1);
jr = 1;
i2r = 0;
else;
%
%     CASE N = 2**P+1
%
for i = 1 : mr;
b(i) = q(i,j) + q(i,1) - q(i,jm1) + q(i,nlast) - q(i,jm2);
end; i = fix(mr+1);
if( np==2 )
for i = 1 : mr;
fi =(q(i,j)-q(i,jm1)-q(i,jp1))./2.;
b2(i) = q(i,1) + fi;
b3(i) = q(i,nlast) + fi;
end; i = fix(mr+1);
k1 = fix(nlast + jr - 1);
k2 = fix(k1 + jr - 1);
[dumvar1,dumvar2,dumvar3,dumvar4,tcos(sub2ind(size(tcos),max(k1+1,1)):end)]=cosgen(jr-1,1,0.0,1.0,tcos(sub2ind(size(tcos),max(k1+1,1)):end));
[nlast,dumvar2,dumvar3,dumvar4,tcos(sub2ind(size(tcos),max(k2+1,1)):end)]=cosgen(nlast,1,0.5,0.0,tcos(sub2ind(size(tcos),max(k2+1,1)):end));
[tcos,k1,dumvar3,k2,nlast]=s1merg(tcos,k1,jr-1,k2,nlast,0);
k3 = fix(k1 + nlast - 1);
k4 = fix(k3 + jr);
[jr,dumvar2,dumvar3,dumvar4,tcos(sub2ind(size(tcos),max(k3+1,1)):end)]=cosgen(jr,1,0.5,0.5,tcos(sub2ind(size(tcos),max(k3+1,1)):end));
[jr,dumvar2,dumvar3,dumvar4,tcos(sub2ind(size(tcos),max(k4+1,1)):end)]=cosgen(jr,1,0.0,0.5,tcos(sub2ind(size(tcos),max(k4+1,1)):end));
jr_orig=jr;    [tcos,k3,jr,k4,dumvar5,k1]=s1merg(tcos,k3,jr,k4,jr,k1);    jr(dumvar5~=jr_orig)=dumvar5(dumvar5~=jr_orig);
k2 = fix(nlast - 1);
k3 = fix(jr);
k4 = fix(jr);
else;
for i = 1 : mr;
b2(i) = q(i,1) + q(i,nlast) + q(i,j) - q(i,jm1) - q(i,jp1);
b3(i) = 0.;
end; i = fix(mr+1);
k1 = fix(nlast - 1);
k2 = fix(nlast + jr - 1);
[dumvar1,dumvar2,dumvar3,dumvar4,tcos(sub2ind(size(tcos),max(nlast,1)):end)]=cosgen(jr-1,1,0.0,1.0,tcos(sub2ind(size(tcos),max(nlast,1)):end));
tcos(k2) = 2.*np - 4;
[jr,dumvar2,dumvar3,dumvar4,tcos(sub2ind(size(tcos),max(k2+1,1)):end)]=cosgen(jr,1,0.5-fnum,0.5,tcos(sub2ind(size(tcos),max(k2+1,1)):end));
k3 =fix(fix((3-np)./2));
[tcos,k1]=s1merg(tcos,k1,jr-k3,k2-k3,jr+k3,0);
k1 = fix(k1 - 1 + k3);
[jr,dumvar2,fnum,dumvar4,tcos(sub2ind(size(tcos),max(k1+1,1)):end)]=cosgen(jr,1,fnum,0.5,tcos(sub2ind(size(tcos),max(k1+1,1)):end));
k2 = fix(jr);
k3 = 0;
k4 = 0;
end;
[mr,a,bb,c,k,b,b2,b3,tcos,d,w,w2,w3]=tri3(mr,a,bb,c,k,b,b2,b3,tcos,d,w,w2,w3);
for i = 1 : mr;
b(i) = b(i) + b2(i) + b3(i);
end; i = fix(mr+1);
if( np==3 )
tcos(1) = 2.;
[dumvar1,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(1,0,mr,a,bb,c,b,tcos,d,w);
end;
for i = 1 : mr;
q(i,j) = b(i) + .5.*(q(i,j)-q(i,jm1)-q(i,jp1));
b(i) = q(i,j) + q(i,1);
end; i = fix(mr+1);
[jr,dumvar2,fnum,dumvar4,tcos]=cosgen(jr,1,fnum,0.5,tcos);
[jr,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(jr,0,mr,a,bb,c,b,tcos,d,w);
for i = 1 : mr;
q(i,1) = q(i,1) - q(i,jm1) + b(i);
end; i = fix(mr+1);
end;
%
%     START BACK SUBSTITUTION.
%
while( true );
j = fix(nlast - jr);
for i = 1 : mr;
b(i) = q(i,nlast) + q(i,j);
end; i = fix(mr+1);
jm2 = fix(nlast - i2r);
if( jr==1 )
for i = 1 : mr;
q(i,nlast) = 0.;
end; i = fix(mr+1);
elseif( nrod==0 ) ;
for i = 1 : mr;
ii = fix(ip + i);
q(i,nlast) = p(ii);
end; i = fix(mr+1);
ip = fix(ip - mr);
else;
for i = 1 : mr;
q(i,nlast) = q(i,nlast) - q(i,jm2);
end; i = fix(mr+1);
end;
[kr,dumvar2,fnum2,dumvar4,tcos]=cosgen(kr,1,fnum2,0.5,tcos);
[lr,dumvar2,fnum2,dumvar4,tcos(sub2ind(size(tcos),max(kr+1,1)):end)]=cosgen(lr,1,fnum2,0.5,tcos(sub2ind(size(tcos),max(kr+1,1)):end));
[kr,lr,mr,a,bb,c,b,tcos,d,w]=trix(kr,lr,mr,a,bb,c,b,tcos,d,w);
for i = 1 : mr;
q(i,nlast) = q(i,nlast) + b(i);
end; i = fix(mr+1);
nlastp = fix(nlast);
while( true );
jstep = fix(jr);
jr = fix(i2r);
i2r = fix(fix(i2r./2));
if( jr==0 )
w(1) = ipstor;




a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
bb_shape=zeros(bb_shape);bb_shape(:)=bb(1:numel(bb_shape));bb=bb_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
b2_shape=zeros(b2_shape);b2_shape(:)=b2(1:numel(b2_shape));b2=b2_shape;
b3_shape=zeros(b3_shape);b3_shape(:)=b3(1:numel(b3_shape));b3=b3_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
w2_shape=zeros(w2_shape);w2_shape(:)=w2(1:numel(w2_shape));w2=w2_shape;
w3_shape=zeros(w3_shape);w3_shape(:)=w3(1:numel(w3_shape));w3=w3_shape;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
tcos_shape=zeros(tcos_shape);tcos_shape(:)=tcos(1:numel(tcos_shape));tcos=tcos_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
return;
end;
jstart = fix(1 + jr);
kr = fix(kr - jr);
if( nlast+jr>n )
jstop = fix(nlast - jr);
else;
kr = fix(kr - jr);
nlast = fix(nlast + jr);
jstop = fix(nlast - jstep);
end;
lr = fix(kr - jr);
[jr,dumvar2,dumvar3,dumvar4,tcos]=cosgen(jr,1,0.5,0.0,tcos);
for j = jstart : jstep: jstop ;
jm2 = fix(j - jr);
jp2 = fix(j + jr);
if( j==jr )
for i = 1 : mr;
b(i) = q(i,j) + q(i,jp2);
end; i = fix(mr+1);
else;
for i = 1 : mr;
b(i) = q(i,j) + q(i,jm2) + q(i,jp2);
end; i = fix(mr+1);
end;
if( jr==1 )
for i = 1 : mr;
q(i,j) = 0.;
end; i = fix(mr+1);
else;
jm1 = fix(j - i2r);
jp1 = fix(j + i2r);
for i = 1 : mr;
q(i,j) = .5.*(q(i,j)-q(i,jm1)-q(i,jp1));
end; i = fix(mr+1);
end;
[jr,dumvar2,mr,a,bb,c,b,tcos,d,w]=trix(jr,0,mr,a,bb,c,b,tcos,d,w);
for i = 1 : mr;
q(i,j) = q(i,j) + b(i);
end; i = fix(mr+1);
end; j = fix(jstop +1);
nrod = 1;
if( nlast+i2r<=n )
nrod = 0;
end;
if( nlastp~=nlast )
break;
end;
end;
end;
%
%     RETURN STORAGE REQUIREMENTS FOR P VECTORS.
%
w(1) = ipstor;




a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
bb_shape=zeros(bb_shape);bb_shape(:)=bb(1:numel(bb_shape));bb=bb_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
b2_shape=zeros(b2_shape);b2_shape(:)=b2(1:numel(b2_shape));b2=b2_shape;
b3_shape=zeros(b3_shape);b3_shape(:)=b3(1:numel(b3_shape));b3=b3_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
w2_shape=zeros(w2_shape);w2_shape(:)=w2(1:numel(w2_shape));w2=w2_shape;
w3_shape=zeros(w3_shape);w3_shape(:)=w3(1:numel(w3_shape));w3=w3_shape;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
tcos_shape=zeros(tcos_shape);tcos_shape(:)=tcos(1:numel(tcos_shape));tcos=tcos_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
end %subroutine postg2
%DECK PPADD