| [a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w]=hwscyl(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w); |
function [a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w]=hwscyl(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w);
persistent a1 a2 deltar deltht dlrby2 dlrsq dlthsq i id2 id3 id4 id5 id6 ierr1 ij istart j k l mp1 mstart mstop munk np np1 nsp1 nstart nstm1 nstop nunk r s s1 s2 ;
if isempty(a1), a1=0; end;
if isempty(a2), a2=0; end;
if isempty(deltar), deltar=0; end;
if isempty(deltht), deltht=0; end;
if isempty(dlrby2), dlrby2=0; end;
if isempty(dlrsq), dlrsq=0; end;
if isempty(dlthsq), dlthsq=0; end;
if isempty(r), r=0; end;
if isempty(s), s=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(i), i=0; end;
if isempty(id2), id2=0; end;
if isempty(id3), id3=0; end;
if isempty(id4), id4=0; end;
if isempty(id5), id5=0; end;
if isempty(id6), id6=0; end;
if isempty(ierr1), ierr1=0; end;
if isempty(ij), ij=0; end;
if isempty(istart), istart=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(mp1), mp1=0; end;
if isempty(mstart), mstart=0; end;
if isempty(mstop), mstop=0; end;
if isempty(munk), munk=0; end;
%***BEGIN PROLOGUE HWSCYL
%***PURPOSE Solve a standard finite difference approximation
% to the Helmholtz equation in cylindrical coordinates.
%***LIBRARY SLATEC (FISHPACK)
%***CATEGORY I2B1A1A
%***TYPE SINGLE PRECISION (HWSCYL-S)
%***KEYWORDS CYLINDRICAL, ELLIPTIC, FISHPACK, HELMHOLTZ, PDE
%***AUTHOR Adams, J., (NCAR)
% Swarztrauber, P. N., (NCAR)
% Sweet, R., (NCAR)
%***DESCRIPTION
%
% subroutine HWSCYL solves a finite difference approximation to the
% Helmholtz equation in cylindrical coordinates:
%
% (1/R,d/dR,R(dU/dR)) + (d/dZ,dU/dZ)
%
% + (LAMBDA/R**2)U = F(R,Z)
%
% This modified Helmholtz equation results from the Fourier
% transform of the three-dimensional Poisson equation.
%
% * * * * * * * * Parameter Description * * * * * * * * * *
%
% * * * * * * On Input * * * * * *
%
% A,B
% The range of R, i.e., A .LE. R .LE. B. A must be less than B
% and A must be non-negative.
%
% M
% The number of panels into which the interval (A,B) is
% subdivided. Hence, there will be M+1 grid points in the
% R-direction given by R(I) = A+(I-1)DR, for I = 1,2,...,M+1,
% where DR = (B-A)/M is the panel width. M must be greater than 3.
%
% MBDCND
% Indicates the type of boundary conditions at R = A and R = B.
%
% = 1 If the solution is specified at R = A and R = B.
% = 2 If the solution is specified at R = A and the derivative of
% the solution with respect to R is specified at R = B.
% = 3 If the derivative of the solution with respect to R is
% specified at R = A (see note below) and R = B.
% = 4 If the derivative of the solution with respect to R is
% specified at R = A (see note below) and the solution is
% specified at R = B.
% = 5 If the solution is unspecified at R = A = 0 and the
% solution is specified at R = B.
% = 6 If the solution is unspecified at R = A = 0 and the
% derivative of the solution with respect to R is specified
% at R = B.
%
% NOTE: If A = 0, do not use MBDCND = 3 or 4, but instead use
% MBDCND = 1,2,5, or 6 .
%
% BDA
% A one-dimensional array of length N+1 that specifies the values
% of the derivative of the solution with respect to R at R = A.
% When MBDCND = 3 or 4,
%
% BDA(J) = (d/dR)U(A,Z(J)), J = 1,2,...,N+1 .
%
% When MBDCND has any other value, BDA is a dummy variable.
%
% BDB
% A one-dimensional array of length N+1 that specifies the values
% of the derivative of the solution with respect to R at R = B.
% When MBDCND = 2,3, or 6,
%
% BDB(J) = (d/dR)U(B,Z(J)), J = 1,2,...,N+1 .
%
% When MBDCND has any other value, BDB is a dummy variable.
%
% C,D
% The range of Z, i.e., C .LE. Z .LE. D. C must be less than D.
%
% N
% The number of panels into which the interval (C,D) is
% subdivided. Hence, there will be N+1 grid points in the
% Z-direction given by Z(J) = C+(J-1)DZ, for J = 1,2,...,N+1,
% where DZ = (D-C)/N is the panel width. N must be greater than 3.
%
% NBDCND
% Indicates the type of boundary conditions at Z = C and Z = D.
%
% = 0 If the solution is periodic in Z, i.e., U(I,1) = U(I,N+1).
% = 1 If the solution is specified at Z = C and Z = D.
% = 2 If the solution is specified at Z = C and the derivative of
% the solution with respect to Z is specified at Z = D.
% = 3 If the derivative of the solution with respect to Z is
% specified at Z = C and Z = D.
% = 4 If the derivative of the solution with respect to Z is
% specified at Z = C and the solution is specified at Z = D.
%
% BDC
% A one-dimensional array of length M+1 that specifies the values
% of the derivative of the solution with respect to Z at Z = C.
% When NBDCND = 3 or 4,
%
% BDC(I) = (d/dZ)U(R(I),C), I = 1,2,...,M+1 .
%
% When NBDCND has any other value, BDC is a dummy variable.
%
% BDD
% A one-dimensional array of length M+1 that specifies the values
% of the derivative of the solution with respect to Z at Z = D.
% When NBDCND = 2 or 3,
%
% BDD(I) = (d/dZ)U(R(I),D), I = 1,2,...,M+1 .
%
% When NBDCND has any other value, BDD is a dummy variable.
%
% ELMBDA
% The constant LAMBDA in the Helmholtz equation. If
% LAMBDA .GT. 0, a solution may not exist. However, HWSCYL will
% attempt to find a solution. LAMBDA must be zero when
% MBDCND = 5 or 6 .
%
% F
% A two-dimensional array that specifies the values of the right
% side of the Helmholtz equation and boundary data (if any). For
% I = 2,3,...,M and J = 2,3,...,N
%
% F(I,J) = F(R(I),Z(J)).
%
% On the boundaries F is defined by
%
% MBDCND F(1,J) F(M+1,J)
% ------ --------- ---------
%
% 1 U(A,Z(J)) U(B,Z(J))
% 2 U(A,Z(J)) F(B,Z(J))
% 3 F(A,Z(J)) F(B,Z(J)) J = 1,2,...,N+1
% 4 F(A,Z(J)) U(B,Z(J))
% 5 F(0,Z(J)) U(B,Z(J))
% 6 F(0,Z(J)) F(B,Z(J))
%
% NBDCND F(I,1) F(I,N+1)
% ------ --------- ---------
%
% 0 F(R(I),C) F(R(I),C)
% 1 U(R(I),C) U(R(I),D)
% 2 U(R(I),C) F(R(I),D) I = 1,2,...,M+1
% 3 F(R(I),C) F(R(I),D)
% 4 F(R(I),C) U(R(I),D)
%
% F must be dimensioned at least (M+1)*(N+1).
%
% NOTE
%
% If the table calls for both the solution U and the right side F
% at a corner then the solution must be specified.
%
% IDIMF
% The row (or first) dimension of the array F as it appears in the
% program calling HWSCYL. This parameter is used to specify the
% variable dimension of F. IDIMF must be at least M+1 .
%
% W
% A one-dimensional array that must be provided by the user for
% work space. W may require up to 4*(N+1) +
% (13 + INT(log2(N+1)))*(M+1) locations. The actual number of
% locations used is computed by HWSCYL and is returned in location
% W(1).
%
%
% * * * * * * On Output * * * * * *
%
% F
% contains the solution U(I,J) of the finite difference
% approximation for the grid point (R(I),Z(J)), I = 1,2,...,M+1,
% J = 1,2,...,N+1 .
%
% PERTRB
% If one specifies a combination of periodic, derivative, and
% unspecified boundary conditions for a Poisson equation
% (LAMBDA = 0), a solution may not exist. PERTRB is a constant,
% calculated and subtracted from F, which ensures that a solution
% exists. HWSCYL then computes this solution, which is a least
% squares solution to the original approximation. This solution
% plus any constant is also a solution. Hence, the solution is
% not unique. The value of PERTRB should be small compared to the
% right side F. Otherwise, a solution is obtained to an
% essentially different problem. This comparison should always
% be made to insure that a meaningful solution has been obtained.
%
% IERROR
% An error flag which indicates invalid input parameters. Except
% for numbers 0 and 11, a solution is not attempted.
%
% = 0 No error.
% = 1 A .LT. 0 .
% = 2 A .GE. B.
% = 3 MBDCND .LT. 1 or MBDCND .GT. 6 .
% = 4 C .GE. D.
% = 5 N .LE. 3
% = 6 NBDCND .LT. 0 or NBDCND .GT. 4 .
% = 7 A = 0, MBDCND = 3 or 4 .
% = 8 A .GT. 0, MBDCND .GE. 5 .
if isempty(np), np=0; end;
if isempty(np1), np1=0; end;
if isempty(nsp1), nsp1=0; end;
if isempty(nstart), nstart=0; end;
if isempty(nstm1), nstm1=0; end;
if isempty(nstop), nstop=0; end;
if isempty(nunk), nunk=0; end;
% = 9 A = 0, LAMBDA .NE. 0, MBDCND .GE. 5 .
% = 10 IDIMF .LT. M+1 .
% = 11 LAMBDA .GT. 0 .
% = 12 M .LE. 3
%
% Since this is the only means of indicating a possibly incorrect
% call to HWSCYL, the user should test IERROR after the call.
%
% W
% W(1) contains the required length of W.
%
% *Long Description:
%
% * * * * * * * program Specifications * * * * * * * * * * * *
%
% Dimension of BDA(N+1),BDB(N+1),BDC(M+1),BDD(M+1),F(IDIMF,N+1),
% Arguments W(see argument list)
%
% Latest June 1, 1976
% Revision
%
% Subprograms HWSCYL,GENBUN,POISD2,POISN2,POISP2,COSGEN,MERGE,
% Required TRIX,TRI3,PIMACH
%
% Special NONE
% Conditions
%
% Common NONE
% Blocks
%
% I/O NONE
%
% Precision Single
%
% Specialist Roland Sweet
%
% Language FORTRAN
%
% History Standardized September 1, 1973
% Revised April 1, 1976
%
% Algorithm The routine defines the finite difference
% equations, incorporates boundary data, and adjusts
% the right side of singular systems and then calls
% GENBUN to solve the system.
%
% Space 5818(decimal) = 13272(octal) locations on the NCAR
% Required Control Data 7600
%
% Timing and The execution time T on the NCAR Control Data
% Accuracy 7600 for subroutine HWSCYL is roughly proportional
% to M*N*log2(N), but also depends on the input
% parameters NBDCND and MBDCND. Some typical values
% are listed in the table below.
% The solution process employed results in a loss
% of no more than three significant digits for N and
% M as large as 64. More detailed information about
% accuracy can be found in the documentation for
% subroutine GENBUN which is the routine that
% solves the finite difference equations.
%
%
% M(=N) MBDCND NBDCND T(MSECS)
% ----- ------ ------ --------
%
% 32 1 0 31
% 32 1 1 23
% 32 3 3 36
% 64 1 0 128
% 64 1 1 96
% 64 3 3 142
%
% Portability American National Standards Institute FORTRAN.
% The machine dependent constant PI is defined in
% function PIMACH.
%
% Required COS
% Resident
% Routines
%
% Reference Swarztrauber, P. and R. Sweet, 'Efficient FORTRAN
% Subprograms for the Solution of Elliptic Equations'
% NCAR TN/IA-109, July, 1975, 138 pp.
%
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%
%***REFERENCES P. N. Swarztrauber and R. Sweet, Efficient Fortran
% subprograms for the solution of elliptic equations,
% NCAR TN/IA-109, July 1975, 138 pp.
%***ROUTINES CALLED GENBUN
%***REVISION HISTORY (YYMMDD)
% 801001 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890531 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE HWSCYL
%
%
f_shape=size(f);f=reshape([f(:).',zeros(1,ceil(numel(f)./prod([idimf])).*prod([idimf])-numel(f))],idimf,[]);
bda_shape=size(bda);bda=reshape(bda,1,[]);
bdb_shape=size(bdb);bdb=reshape(bdb,1,[]);
bdc_shape=size(bdc);bdc=reshape(bdc,1,[]);
bdd_shape=size(bdd);bdd=reshape(bdd,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
%***FIRST EXECUTABLE STATEMENT HWSCYL
ierror = 0;
if( a<0. )
ierror = 1;
end;
if( a>=b )
ierror = 2;
end;
if( mbdcnd<=0 || mbdcnd>=7 )
ierror = 3;
end;
if( c>=d )
ierror = 4;
end;
if( n<=3 )
ierror = 5;
end;
if( nbdcnd<=-1 || nbdcnd>=5 )
ierror = 6;
end;
if( a==0. &&(mbdcnd==3 || mbdcnd==4) )
ierror = 7;
end;
if( a>0. && mbdcnd>=5 )
ierror = 8;
end;
if( a==0. && elmbda~=0. && mbdcnd>=5 )
ierror = 9;
end;
if( idimf<m+1 )
ierror = 10;
end;
if( m<=3 )
ierror = 12;
end;
if( ierror~=0 )
f_shape=zeros(f_shape);f_shape(:)=f(1:numel(f_shape));f=f_shape;
bda_shape=zeros(bda_shape);bda_shape(:)=bda(1:numel(bda_shape));bda=bda_shape;
bdb_shape=zeros(bdb_shape);bdb_shape(:)=bdb(1:numel(bdb_shape));bdb=bdb_shape;
bdc_shape=zeros(bdc_shape);bdc_shape(:)=bdc(1:numel(bdc_shape));bdc=bdc_shape;
bdd_shape=zeros(bdd_shape);bdd_shape(:)=bdd(1:numel(bdd_shape));bdd=bdd_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
return;
end;
mp1 = fix(m + 1);
deltar =(b-a)./m;
dlrby2 = deltar./2.;
dlrsq = deltar.^2;
np1 = fix(n + 1);
deltht =(d-c)./n;
dlthsq = deltht.^2;
np = fix(nbdcnd + 1);
%
% DEFINE RANGE OF INDICES I AND J FOR UNKNOWNS U(I,J).
%
mstart = 2;
mstop = fix(m);
if( mbdcnd~=1 )
if( mbdcnd==2 )
mstop = fix(mp1);
elseif( mbdcnd==3 || mbdcnd==6 ) ;
mstart = 1;
mstop = fix(mp1);
else;
mstart = 1;
end;
end;
munk = fix(mstop - mstart + 1);
nstart = 1;
nstop = fix(n);
if( np~=1 && np~=5 )
if( np==3 )
nstart = 2;
nstop = fix(np1);
elseif( np==4 ) ;
nstop = fix(np1);
else;
nstart = 2;
end;
end;
nunk = fix(nstop - nstart + 1);
%
% DEFINE A,B,C COEFFICIENTS IN W-ARRAY.
%
id2 = fix(munk);
id3 = fix(id2 + munk);
id4 = fix(id3 + munk);
id5 = fix(id4 + munk);
id6 = fix(id5 + munk);
istart = 1;
a1 = 2../dlrsq;
ij = 0;
if( mbdcnd==3 || mbdcnd==4 )
ij = 1;
end;
if( mbdcnd>4 )
w(1) = 0.;
w(id2+1) = -2..*a1;
w(id3+1) = 2..*a1;
istart = 2;
ij = 1;
end;
for i = istart : munk;
r = a +(i-ij).*deltar;
j = fix(id5 + i);
w(j) = r;
j = fix(id6 + i);
w(j) = 1../r.^2;
w(i) =(r-dlrby2)./(r.*dlrsq);
j = fix(id3 + i);
w(j) =(r+dlrby2)./(r.*dlrsq);
k = fix(id6 + i);
j = fix(id2 + i);
w(j) = -a1 + elmbda.*w(k);
end; i = fix(munk+1);
if( mbdcnd~=1 && mbdcnd~=5 )
if( mbdcnd==3 || mbdcnd==6 )
w(id2) = a1;
w(id3+1) = a1.*istart;
elseif( mbdcnd==4 ) ;
w(id3+1) = a1.*istart;
else;
w(id2) = a1;
end;
end;
%
% ENTER BOUNDARY DATA FOR R-BOUNDARIES.
%
if( mbdcnd==3 || mbdcnd==4 )
a1 = 2..*deltar.*w(1);
for j = nstart : nstop;
f(1,j) = f(1,j) + a1.*bda(j);
end; j = fix(nstop+1);
elseif( mbdcnd~=5 && mbdcnd~=6 ) ;
a1 = w(1);
for j = nstart : nstop;
f(2,j) = f(2,j) - a1.*f(1,j);
end; j = fix(nstop+1);
end;
if( mbdcnd==2 || mbdcnd==3 || mbdcnd==6 )
a1 = 2..*deltar.*w(id4);
for j = nstart : nstop;
f(mp1,j) = f(mp1,j) - a1.*bdb(j);
end; j = fix(nstop+1);
else;
a1 = w(id4);
for j = nstart : nstop;
f(m,j) = f(m,j) - a1.*f(mp1,j);
end; j = fix(nstop+1);
end;
%
% ENTER BOUNDARY DATA FOR Z-BOUNDARIES.
%
a1 = 1../dlthsq;
l = fix(id5 - mstart + 1);
if( np~=1 )
if( np==4 || np==5 )
a1 = 2../deltht;
for i = mstart : mstop;
f(i,1) = f(i,1) + a1.*bdc(i);
end; i = fix(mstop+1);
else;
for i = mstart : mstop;
f(i,2) = f(i,2) - a1.*f(i,1);
end; i = fix(mstop+1);
end;
a1 = 1../dlthsq;
if( np~=1 )
if( np==3 || np==4 )
a1 = 2../deltht;
for i = mstart : mstop;
f(i,np1) = f(i,np1) - a1.*bdd(i);
end; i = fix(mstop+1);
else;
for i = mstart : mstop;
f(i,n) = f(i,n) - a1.*f(i,np1);
end; i = fix(mstop+1);
end;
end;
end;
%
% ADJUST RIGHT SIDE OF SINGULAR PROBLEMS TO INSURE EXISTENCE OF A
% SOLUTION.
%
pertrb = 0.;
while (1);
if( elmbda>=0 )
if( elmbda==0 )
w(id5+1) = .5.*(w(id5+2)-dlrby2);
if( mbdcnd~=1 && mbdcnd~=2 && mbdcnd~=4 && mbdcnd~=5 )
if( mbdcnd~=3 )
w(id5+1) = .5.*w(id5+1);
end;
if( np==1 )
a2 = 1.;
elseif( np==2 || np==3 || np==5 ) ;
break;
else;
a2 = 2.;
end;
k = fix(id5 + munk);
w(k) = .5.*(w(k-1)+dlrby2);
s = 0.;
for i = mstart : mstop;
s1 = 0.;
nsp1 = fix(nstart + 1);
nstm1 = fix(nstop - 1);
for j = nsp1 : nstm1;
s1 = s1 + f(i,j);
end; j = fix(nstm1+1);
k = fix(i + l);
s = s +(a2.*s1+f(i,nstart)+f(i,nstop)).*w(k);
end; i = fix(mstop+1);
s2 = m.*a +(.75+(m-1).*(m+1)).*dlrby2;
if( mbdcnd==3 )
s2 = s2 + .25.*dlrby2;
end;
s1 =(2.+a2.*(nunk-2)).*s2;
pertrb = s./s1;
for i = mstart : mstop;
for j = nstart : nstop;
f(i,j) = f(i,j) - pertrb;
end; j = fix(nstop+1);
end; i = fix(mstop+1);
end;
else;
ierror = 11;
end;
end;
break;
end;
%
% MULTIPLY I-TH EQUATION THROUGH BY DELTHT**2 TO PUT EQUATION INTO
% CORRECT FORM FOR SUBROUTINE GENBUN.
%
for i = mstart : mstop;
k = fix(i - mstart + 1);
w(k) = w(k).*dlthsq;
j = fix(id2 + k);
w(j) = w(j).*dlthsq;
j = fix(id3 + k);
w(j) = w(j).*dlthsq;
for j = nstart : nstop;
f(i,j) = f(i,j).*dlthsq;
end; j = fix(nstop+1);
end; i = fix(mstop+1);
w(1) = 0.;
w(id4) = 0.;
%
% CALL GENBUN TO SOLVE THE SYSTEM OF EQUATIONS.
%
[nbdcnd,nunk,dumvar3,munk,dumvar5,dumvar6,dumvar7,idimf,f(sub2ind(size(f),mstart,nstart):end),ierr1,dumvar11]=genbun(nbdcnd,nunk,1,munk,w(sub2ind(size(w),max(1,1)):end),w(sub2ind(size(w),max(id2+1,1)):end),w(sub2ind(size(w),max(id3+1,1)):end),idimf,f(sub2ind(size(f),mstart,nstart):end),ierr1,w(sub2ind(size(w),max(id4+1,1)):end)); dumvar5i=find((w(sub2ind(size(w),max(1,1)):end))~=(dumvar5));dumvar6i=find((w(sub2ind(size(w),max(id2+1,1)):end))~=(dumvar6));dumvar7i=find((w(sub2ind(size(w),max(id3+1,1)):end))~=(dumvar7));dumvar11i=find((w(sub2ind(size(w),max(id4+1,1)):end))~=(dumvar11)); w(1-1+dumvar5i)=dumvar5(dumvar5i); w(id2+1-1+dumvar6i)=dumvar6(dumvar6i); w(id3+1-1+dumvar7i)=dumvar7(dumvar7i); w(id4+1-1+dumvar11i)=dumvar11(dumvar11i);
w(1) = w(id4+1) + 3.*munk;
if( nbdcnd==0 )
for i = mstart : mstop;
f(i,np1) = f(i,1);
end; i = fix(mstop+1);
end;
f_shape=zeros(f_shape);f_shape(:)=f(1:numel(f_shape));f=f_shape;
bda_shape=zeros(bda_shape);bda_shape(:)=bda(1:numel(bda_shape));bda=bda_shape;
bdb_shape=zeros(bdb_shape);bdb_shape(:)=bdb(1:numel(bdb_shape));bdb=bdb_shape;
bdc_shape=zeros(bdc_shape);bdc_shape(:)=bdc(1:numel(bdc_shape));bdc=bdc_shape;
bdd_shape=zeros(bdd_shape);bdd_shape(:)=bdd(1:numel(bdd_shape));bdd=bdd_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
end %subroutine hwscyl
%DECK HWSPLR
|
|
