| [a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w]=hstplr(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]=hstplr(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w);
persistent a1 a2 deltar deltht dlrsq dlthsq i ierr1 isw iwb iwc iwr j k lp mb np ;
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(dlrsq), dlrsq=0; end;
if isempty(dlthsq), dlthsq=0; end;
%***BEGIN PROLOGUE HSTPLR
%***PURPOSE Solve the standard five-point finite difference
% approximation on a staggered grid to the Helmholtz equation
% in polar coordinates.
%***LIBRARY SLATEC (FISHPACK)
%***CATEGORY I2B1A1A
%***TYPE SINGLE PRECISION (HSTPLR-S)
%***KEYWORDS ELLIPTIC, FISHPACK, HELMHOLTZ, PDE, POLAR
%***AUTHOR Adams, J., (NCAR)
% Swarztrauber, P. N., (NCAR)
% Sweet, R., (NCAR)
%***DESCRIPTION
%
% HSTPLR solves the standard five-point finite difference
% approximation on a staggered grid to the Helmholtz equation in
% polar coordinates
%
% (1/R,d/DR,R(dU/DR)) + (1/R**2,d/dTHETA,dU/dTHETA)
%
% + LAMBDA*U = F(R,THETA)
%
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%
% * * * * * * * * 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 grid points in the interval (A,B). The grid points
% in the R-direction are given by R(I) = A + (I-0.5)DR for
% I=1,2,...,M where DR =(B-A)/M. M must be greater than 2.
%
% 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.
% (see note 1 below)
%
% = 3 If the derivative of the solution with respect to R is
% specified at R = A (see note 2 below) and R = B.
%
% = 4 If the derivative of the solution with respect to R is
% specified at R = A (see note 2 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 1: If A = 0, MBDCND = 2, and NBDCND = 0 or 3, the system of
% equations to be solved is singular. The unique solution
% is determined by extrapolation to the specification of
% U(0,THETA(1)). But in this case the right side of the
% system will be perturbed by the constant PERTRB.
%
% NOTE 2: 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 that specifies the boundary
% values (if any) of the solution at R = A. When MBDCND = 1 or 2,
%
% BDA(J) = U(A,THETA(J)) , J=1,2,...,N.
%
% When MBDCND = 3 or 4,
%
% BDA(J) = (d/dR)U(A,THETA(J)) , J=1,2,...,N.
%
% When MBDCND = 5 or 6, BDA is a dummy variable.
%
% BDB
% A one-dimensional array of length N that specifies the boundary
% values of the solution at R = B. When MBDCND = 1,4, or 5,
%
% BDB(J) = U(B,THETA(J)) , J=1,2,...,N.
%
% When MBDCND = 2,3, or 6,
%
% BDB(J) = (d/dR)U(B,THETA(J)) , J=1,2,...,N.
%
% C,D
% The range of THETA, i.e. C .LE. THETA .LE. D. C must be less
% than D.
%
% N
% The number of unknowns in the interval (C,D). The unknowns in
% the THETA-direction are given by THETA(J) = C + (J-0.5)DT,
% J=1,2,...,N, where DT = (D-C)/N. N must be greater than 2.
%
% NBDCND
% Indicates the type of boundary conditions at THETA = C
% and THETA = D.
%
% = 0 If the solution is periodic in THETA, i.e.
% U(I,J) = U(I,N+J).
%
% = 1 If the solution is specified at THETA = C and THETA = D
% (see note below).
%
% = 2 If the solution is specified at THETA = C and the derivative
% of the solution with respect to THETA is specified at
% THETA = D (see note below).
%
% = 3 If the derivative of the solution with respect to THETA is
% specified at THETA = C and THETA = D.
%
% = 4 If the derivative of the solution with respect to THETA is
% specified at THETA = C and the solution is specified at
% THETA = d (see note below).
%
% NOTE: When NBDCND = 1, 2, or 4, do not use MBDCND = 5 or 6 (the
% former indicates that the solution is specified at R = 0; the
% latter indicates the solution is unspecified at R = 0). Use
% instead MBDCND = 1 or 2.
%
% BDC
% A one dimensional array of length M that specifies the boundary
% values of the solution at THETA = C. When NBDCND = 1 or 2,
%
% BDC(I) = U(R(I),C) , I=1,2,...,M.
%
% When NBDCND = 3 or 4,
%
% BDC(I) = (d/dTHETA)U(R(I),C), I=1,2,...,M.
%
% When NBDCND = 0, BDC is a dummy variable.
%
% BDD
% A one-dimensional array of length M that specifies the boundary
% values of the solution at THETA = D. When NBDCND = 1 or 4,
%
% BDD(I) = U(R(I),D) , I=1,2,...,M.
%
% When NBDCND = 2 or 3,
%
% BDD(I) = (d/dTHETA)U(R(I),D) , I=1,2,...,M.
%
% When NBDCND = 0, BDD is a dummy variable.
%
% ELMBDA
% The constant LAMBDA in the Helmholtz equation. If LAMBDA is
% greater than 0, a solution may not exist. However, HSTPLR will
% attempt to find a solution.
%
% F
% A two-dimensional array that specifies the values of the right
% side of the Helmholtz equation. For I=1,2,...,M and J=1,2,...,N
%
% F(I,J) = F(R(I),THETA(J)) .
%
% F must be dimensioned at least M X N.
%
% IDIMF
% The row (or first) dimension of the array F as it appears in the
% program calling HSTPLR. This parameter is used to specify the
% variable dimension of F. IDIMF must be at least M.
%
% W
% A one-dimensional array that must be provided by the user for
% work space. W may require up to 13M + 4N + M*INT(log2(N))
% locations. The actual number of locations used is computed by
% HSTPLR and is returned in the location W(1).
%
%
% * * * * * * On Output * * * * * *
%
% F
% contains the solution U(I,J) of the finite difference
% approximation for the grid point (R(I),THETA(J)) for
% I=1,2,...,M, J=1,2,...,N.
%
% PERTRB
% If a combination of periodic, derivative, or unspecified
% boundary conditions is specified for a Poisson equation
% (LAMBDA = 0), a solution may not exist. PERTRB is a con-
% stant, calculated and subtracted from F, which ensures
% that a solution exists. HSTPLR 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 that indicates invalid input parameters.
% Except for numbers 0 and 11, a solution is not attempted.
%
% = 0 No error
if isempty(i), i=0; end;
if isempty(ierr1), ierr1=0; end;
if isempty(isw), isw=0; end;
if isempty(iwb), iwb=0; end;
if isempty(iwc), iwc=0; end;
if isempty(iwr), iwr=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(lp), lp=0; end;
if isempty(mb), mb=0; end;
if isempty(np), np=0; end;
%
% = 1 A .LT. 0
%
% = 2 A .GE. B
%
% = 3 MBDCND .LT. 1 or MBDCND .GT. 6
%
% = 4 C .GE. D
%
% = 5 N .LE. 2
%
% = 6 NBDCND .LT. 0 or NBDCND .GT. 4
%
% = 7 A = 0 and MBDCND = 3 or 4
%
% = 8 A .GT. 0 and MBDCND .GE. 5
%
% = 9 MBDCND .GE. 5 and NBDCND .NE. 0 or 3
%
% = 10 IDIMF .LT. M
%
% = 11 LAMBDA .GT. 0
%
% = 12 M .LE. 2
%
% Since this is the only means of indicating a possibly
% incorrect call to HSTPLR, 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),BDB(N),BDC(M),BDD(M),F(IDIMF,N),
% Arguments W(see ARGUMENT LIST)
%
% Latest June 1, 1977
% Revision
%
% Subprograms HSTPLR,POISTG,POSTG2,GENBUN,POISD2,POISN2,POISP2,
% Required COSGEN,MERGE,TRIX,TRI3,PIMACH
%
% Special NONE
% Conditions
%
% Common NONE
% Blocks
%
% I/O NONE
%
% Precision Single
%
% Specialist Roland Sweet
%
% Language FORTRAN
%
% History Written by Roland Sweet at NCAR in February, 1977
%
% Algorithm This subroutine defines the finite-difference
% equations, incorporates boundary data, adjusts the
% right side when the system is singular and calls
% either POISTG or GENBUN which solves the linear
% system of equations.
%
% Space 8265(decimal) = 20111(octal) LOCATIONS ON THE
% Required NCAR Control Data 7600
%
% Timing and The execution time T on the NCAR Control Data
% Accuracy 7600 for subroutine HSTPLR is roughly proportional
% to M*N*log2(N). Some typical values are listed in
% the table below.
% The solution process employed results in a loss
% of no more than four significant digits for N and M
% as large as 64. More detailed information about
% accuracy can be found in the documentation for
% subroutine POISTG which is the routine that
% actually solves the finite difference equations.
%
%
% M(=N) MBDCND NBDCND T(MSECS)
% ----- ------ ------ --------
%
% 32 1-6 1-4 56
% 64 1-6 1-4 230
%
% Portability American National Standards Institute Fortran.
% The machine dependent constant PI is defined in
% function PIMACH.
%
% Required COS
% Resident
% Routines
%
% Reference Schumann, U. and R. Sweet,'A Direct Method For
% The Solution Of Poisson's Equation With Neumann
% Boundary Conditions On A Staggered Grid of
% Arbitrary Size,' J. Comp. Phys. 20(1976),
% pp. 171-182.
%
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%
%***REFERENCES U. Schumann and R. Sweet, A direct method for the
% solution of Poisson's equation with Neumann boundary
% conditions on a staggered grid of arbitrary size,
% Journal of Computational Physics 20, (1976),
% pp. 171-182.
%***ROUTINES CALLED GENBUN, POISTG
%***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 HSTPLR
%
%
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 HSTPLR
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<=2 )
ierror = 5;
end;
if( nbdcnd<0 || nbdcnd>=5 )
ierror = 6;
end;
if( a==0. &&(mbdcnd==3 || mbdcnd==4) )
ierror = 7;
end;
if( a>0. && mbdcnd>=5 )
ierror = 8;
end;
if( mbdcnd>=5 && nbdcnd~=0 && nbdcnd~=3 )
ierror = 9;
end;
if( idimf<m )
ierror = 10;
end;
if( m<=2 )
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;
deltar =(b-a)./m;
dlrsq = deltar.^2;
deltht =(d-c)./n;
dlthsq = deltht.^2;
np = fix(nbdcnd + 1);
isw = 1;
mb = fix(mbdcnd);
if( a==0. && mbdcnd==2 )
mb = 6;
end;
%
% DEFINE A,B,C COEFFICIENTS IN W-ARRAY.
%
iwb = fix(m);
iwc = fix(iwb + m);
iwr = fix(iwc + m);
for i = 1 : m;
j = fix(iwr + i);
w(j) = a +(i-0.5).*deltar;
w(i) =(a+(i-1).*deltar)./dlrsq;
k = fix(iwc + i);
w(k) =(a+i.*deltar)./dlrsq;
k = fix(iwb + i);
w(k) =(elmbda-2../dlrsq).*w(j);
end; i = fix(m+1);
for i = 1 : m;
j = fix(iwr + i);
a1 = w(j);
for j = 1 : n;
f(i,j) = a1.*f(i,j);
end; j = fix(n+1);
end; i = fix(m+1);
%
% ENTER BOUNDARY DATA FOR R-BOUNDARIES.
%
if( mb==3 || mb==4 )
a1 = deltar.*w(1);
w(iwb+1) = w(iwb+1) + w(1);
for j = 1 : n;
f(1,j) = f(1,j) + a1.*bda(j);
end; j = fix(n+1);
elseif( mb~=5 && mb~=6 ) ;
a1 = 2..*w(1);
w(iwb+1) = w(iwb+1) - w(1);
for j = 1 : n;
f(1,j) = f(1,j) - a1.*bda(j);
end; j = fix(n+1);
end;
if( mb==2 || mb==3 || mb==6 )
a1 = deltar.*w(iwr);
w(iwc) = w(iwc) + w(iwr);
for j = 1 : n;
f(m,j) = f(m,j) - a1.*bdb(j);
end; j = fix(n+1);
else;
a1 = 2..*w(iwr);
w(iwc) = w(iwc) - w(iwr);
for j = 1 : n;
f(m,j) = f(m,j) - a1.*bdb(j);
end; j = fix(n+1);
end;
%
% ENTER BOUNDARY DATA FOR THETA-BOUNDARIES.
%
a1 = 2../dlthsq;
if( np~=1 )
if( np==4 || np==5 )
a1 = 1../deltht;
for i = 1 : m;
j = fix(iwr + i);
f(i,1) = f(i,1) + a1.*bdc(i)./w(j);
end; i = fix(m+1);
else;
for i = 1 : m;
j = fix(iwr + i);
f(i,1) = f(i,1) - a1.*bdc(i)./w(j);
end; i = fix(m+1);
end;
a1 = 2../dlthsq;
if( np~=1 )
if( np==3 || np==4 )
a1 = 1../deltht;
for i = 1 : m;
j = fix(iwr + i);
f(i,n) = f(i,n) - a1.*bdd(i)./w(j);
end; i = fix(m+1);
else;
for i = 1 : m;
j = fix(iwr + i);
f(i,n) = f(i,n) - a1.*bdd(i)./w(j);
end; i = fix(m+1);
end;
end;
end;
%
% ADJUST RIGHT SIDE OF SINGULAR PROBLEMS TO INSURE EXISTENCE OF A
% SOLUTION.
%
pertrb = 0.;
if( elmbda>=0 )
if( elmbda~=0 )
ierror = 11;
elseif( mb~=1 && mb~=2 && mb~=4 && mb~=5 ) ;
if( np~=2 && np~=3 && np~=5 )
isw = 2;
for j = 1 : n;
for i = 1 : m;
pertrb = pertrb + f(i,j);
end; i = fix(m+1);
end; j = fix(n+1);
pertrb = pertrb./(m.*n.*0.5.*(a+b));
for i = 1 : m;
j = fix(iwr + i);
a1 = pertrb.*w(j);
for j = 1 : n;
f(i,j) = f(i,j) - a1;
end; j = fix(n+1);
end; i = fix(m+1);
a2 = 0.;
for j = 1 : n;
a2 = a2 + f(1,j);
end; j = fix(n+1);
a2 = a2./w(iwr+1);
end;
end;
end;
%
% MULTIPLY I-TH EQUATION THROUGH BY R(I)*DELTHT**2
%
for i = 1 : m;
j = fix(iwr + i);
a1 = dlthsq.*w(j);
w(i) = a1.*w(i);
j = fix(iwc + i);
w(j) = a1.*w(j);
j = fix(iwb + i);
w(j) = a1.*w(j);
for j = 1 : n;
f(i,j) = a1.*f(i,j);
end; j = fix(n+1);
end; i = fix(m+1);
lp = fix(nbdcnd);
w(1) = 0.;
w(iwr) = 0.;
%
% CALL POISTG OR GENBUN TO SOLVE THE SYSTEM OF EQUATIONS.
%
if( lp==0 )
[lp,n,dumvar3,m,w,dumvar6,dumvar7,idimf,f,ierr1,dumvar11]=genbun(lp,n,1,m,w,w(sub2ind(size(w),max(iwb+1,1)):end),w(sub2ind(size(w),max(iwc+1,1)):end),idimf,f,ierr1,w(sub2ind(size(w),max(iwr+1,1)):end)); dumvar6i=find((w(sub2ind(size(w),max(iwb+1,1)):end))~=(dumvar6));dumvar7i=find((w(sub2ind(size(w),max(iwc+1,1)):end))~=(dumvar7));dumvar11i=find((w(sub2ind(size(w),max(iwr+1,1)):end))~=(dumvar11)); w(iwb+1-1+dumvar6i)=dumvar6(dumvar6i); w(iwc+1-1+dumvar7i)=dumvar7(dumvar7i); w(iwr+1-1+dumvar11i)=dumvar11(dumvar11i);
else;
[lp,n,dumvar3,m,w,dumvar6,dumvar7,idimf,f,ierr1,dumvar11]=poistg(lp,n,1,m,w,w(sub2ind(size(w),max(iwb+1,1)):end),w(sub2ind(size(w),max(iwc+1,1)):end),idimf,f,ierr1,w(sub2ind(size(w),max(iwr+1,1)):end)); dumvar6i=find((w(sub2ind(size(w),max(iwb+1,1)):end))~=(dumvar6));dumvar7i=find((w(sub2ind(size(w),max(iwc+1,1)):end))~=(dumvar7));dumvar11i=find((w(sub2ind(size(w),max(iwr+1,1)):end))~=(dumvar11)); w(iwb+1-1+dumvar6i)=dumvar6(dumvar6i); w(iwc+1-1+dumvar7i)=dumvar7(dumvar7i); w(iwr+1-1+dumvar11i)=dumvar11(dumvar11i);
end;
w(1) = w(iwr+1) + 3.*m;
if( a==0. && mbdcnd==2 && isw==2 )
a1 = 0.;
for j = 1 : n;
a1 = a1 + f(1,j);
end; j = fix(n+1);
a1 =(a1-dlrsq.*a2./16.)./n;
if( nbdcnd==3 )
a1 = a1 +(bdd(1)-bdc(1))./(d-c);
end;
a1 = bda(1) - a1;
for i = 1 : m;
for j = 1 : n;
f(i,j) = f(i,j) + a1;
end; j = fix(n+1);
end; i = fix(m+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 hstplr
%DECK HSTSSP
|
|
