[ts,tf,m,mbdcnd,bdts,bdtf,ps,pf,n,nbdcnd,bdps,bdpf,elmbda,f,idimf,pertrb,am,bm,cm,sn,ss,sint,d]=hwsss1(ts,tf,m,mbdcnd,bdts,bdtf,ps,pf,n,nbdcnd,bdps,bdpf,elmbda,f,idimf,pertrb,am,bm,cm,sn,ss,sint,d);

function [ts,tf,m,mbdcnd,bdts,bdtf,ps,pf,n,nbdcnd,bdps,bdpf,elmbda,f,idimf,pertrb,am,bm,cm,sn,ss,sint,d]=hwsss1(ts,tf,m,mbdcnd,bdts,bdtf,ps,pf,n,nbdcnd,bdps,bdpf,elmbda,f,idimf,pertrb,am,bm,cm,sn,ss,sint,d);
persistent at cf cnp cp csp ct den dfn dfs dnn dns dphi dphi2 dsn dss dth dth2 fim1 fjj fm fn hdth hld hne i ierror ii iid inp ising isp itf itfm its itsp j jpf jpfm jps jpsp mbr mp1 munk nbr np1 nunk rtn rts sum1 sum2 summlv t1 tdp tdt theta wp wpf wps wtf wts yhld ; 

if isempty(at), at=0; end;
if isempty(cf), cf=0; end;
if isempty(cnp), cnp=0; end;
if isempty(cp), cp=0; end;
if isempty(csp), csp=0; end;
if isempty(ct), ct=0; end;
if isempty(den), den=0; end;
if isempty(dfn), dfn=0; end;
if isempty(dfs), dfs=0; end;
if isempty(dnn), dnn=0; end;
if isempty(dns), dns=0; end;
if isempty(dphi), dphi=0; end;
if isempty(dphi2), dphi2=0; end;
if isempty(dsn), dsn=0; end;
if isempty(dss), dss=0; end;
if isempty(dth), dth=0; end;
if isempty(dth2), dth2=0; end;
if isempty(fim1), fim1=0; end;
if isempty(fjj), fjj=0; end;
if isempty(fm), fm=0; end;
if isempty(fn), fn=0; end;
if isempty(hdth), hdth=0; end;
if isempty(hld), hld=0; end;
if isempty(hne), hne=0; end;
if isempty(rtn), rtn=0; end;
if isempty(rts), rts=0; end;
if isempty(summlv), summlv=0; end;
if isempty(sum1), sum1=0; end;
if isempty(sum2), sum2=0; end;
if isempty(t1), t1=0; end;
if isempty(tdp), tdp=0; end;
if isempty(tdt), tdt=0; end;
if isempty(theta), theta=0; end;
if isempty(wp), wp=0; end;
if isempty(wpf), wpf=0; end;
if isempty(wps), wps=0; end;
if isempty(wtf), wtf=0; end;
if isempty(wts), wts=0; end;
if isempty(yhld), yhld=0; end;
if isempty(i), i=0; end;
if isempty(ierror), ierror=0; end;
if isempty(ii), ii=0; end;
if isempty(iid), iid=0; end;
if isempty(inp), inp=0; end;
if isempty(ising), ising=0; end;
if isempty(isp), isp=0; end;
if isempty(itf), itf=0; end;
if isempty(itfm), itfm=0; end;
if isempty(its), its=0; end;
if isempty(itsp), itsp=0; end;
if isempty(j), j=0; end;
if isempty(jpf), jpf=0; end;
if isempty(jpfm), jpfm=0; end;
if isempty(jps), jps=0; end;
if isempty(jpsp), jpsp=0; end;
if isempty(mbr), mbr=0; end;
if isempty(mp1), mp1=0; end;
if isempty(munk), munk=0; end;
if isempty(nbr), nbr=0; end;
if isempty(np1), np1=0; end;
if isempty(nunk), nunk=0; end;
%***BEGIN PROLOGUE  HWSSS1
%***SUBSIDIARY
%***PURPOSE  Subsidiary to HWSSSP
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (HWSSS1-S)
%***AUTHOR  (UNKNOWN)
%***SEE ALSO  HWSSSP
%***ROUTINES CALLED  GENBUN
%***REVISION HISTORY  (YYMMDD)
%   801001  DATE WRITTEN
%   891009  Removed unreferenced variables.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  HWSSS1
f_shape=size(f);f=reshape([f(:).',zeros(1,ceil(numel(f)./prod([idimf])).*prod([idimf])-numel(f))],idimf,[]);
bdts_shape=size(bdts);bdts=reshape(bdts,1,[]);
bdtf_shape=size(bdtf);bdtf=reshape(bdtf,1,[]);
bdps_shape=size(bdps);bdps=reshape(bdps,1,[]);
bdpf_shape=size(bdpf);bdpf=reshape(bdpf,1,[]);
am_shape=size(am);am=reshape(am,1,[]);
bm_shape=size(bm);bm=reshape(bm,1,[]);
cm_shape=size(cm);cm=reshape(cm,1,[]);
ss_shape=size(ss);ss=reshape(ss,1,[]);
sn_shape=size(sn);sn=reshape(sn,1,[]);
d_shape=size(d);d=reshape(d,1,[]);
sint_shape=size(sint);sint=reshape(sint,1,[]);
%
%***FIRST EXECUTABLE STATEMENT  HWSSS1
mp1 = fix(m + 1);
np1 = fix(n + 1);
fn = n;
fm = m;
dth =(tf-ts)./fm;
hdth = dth./2.;
tdt = dth + dth;
dphi =(pf-ps)./fn;
tdp = dphi + dphi;
dphi2 = dphi.*dphi;
dth2 = dth.*dth;
cp = 4../(fn.*dth2);
wp = fn.*sin(hdth)./4.;
for i = 1 : mp1;
fim1 = i - 1;
theta = fim1.*dth + ts;
sint(i) = sin(theta);
if( sint(i)~=0 )
t1 = 1../(dth2.*sint(i));
am(i) = t1.*sin(theta-hdth);
cm(i) = t1.*sin(theta+hdth);
bm(i) = -am(i) - cm(i) + elmbda;
end;
end; i = fix(mp1+1);
inp = 0;
isp = 0;
%
% BOUNDARY CONDITION AT THETA=TS
%
mbr = fix(mbdcnd + 1);
if( mbr==2 || mbr==3 || mbr==8 )
at = am(2);
its = 2;
elseif( mbr==4 || mbr==5 || mbr==9 ) ;
at = am(1);
its = 1;
cm(1) = am(1) + cm(1);
elseif( mbr==6 || mbr==7 || mbr==10 ) ;
at = am(2);
inp = 1;
its = 2;
else;
its = 1;
end;
%
% BOUNDARY CONDITION THETA=TF
%
if( mbr==2 || mbr==5 || mbr==6 )
ct = cm(m);
itf = fix(m);
elseif( mbr==3 || mbr==4 || mbr==7 ) ;
ct = cm(m+1);
am(m+1) = am(m+1) + cm(m+1);
itf = fix(m + 1);
elseif( mbr==8 || mbr==9 || mbr==10 ) ;
itf = fix(m);
isp = 1;
ct = cm(m);
else;
itf = fix(m);
end;
%
% COMPUTE HOMOGENEOUS SOLUTION WITH SOLUTION AT POLE EQUAL TO ONE
%
itsp = fix(its + 1);
itfm = fix(itf - 1);
wts = sint(its+1).*am(its+1)./cm(its);
wtf = sint(itf-1).*cm(itf-1)./am(itf);
munk = fix(itf - its + 1);
if( isp>0 )
d(its) = cm(its)./bm(its);
for i = itsp : m;
d(i) = cm(i)./(bm(i)-am(i).*d(i-1));
end; i = fix(m+1);
ss(m) = -d(m);
iid = fix(m - its);
for ii = 1 : iid;
i = fix(m - ii);
ss(i) = -d(i).*ss(i+1);
end; ii = fix(iid+1);
ss(m+1) = 1.;
end;
if( inp>0 )
sn(1) = 1.;
d(itf) = am(itf)./bm(itf);
iid = fix(itf - 2);
for ii = 1 : iid;
i = fix(itf - ii);
d(i) = am(i)./(bm(i)-cm(i).*d(i+1));
end; ii = fix(iid+1);
sn(2) = -d(2);
for i = 3 : itf;
sn(i) = -d(i).*sn(i-1);
end; i = fix(itf+1);
end;
%
% BOUNDARY CONDITIONS AT PHI=PS
%
nbr = fix(nbdcnd + 1);
wps = 1.;
wpf = 1.;
if( nbr==2 || nbr==3 )
jps = 2;
elseif( nbr==4 || nbr==5 ) ;
jps = 1;
wps = .5;
else;
jps = 1;
end;
%
% BOUNDARY CONDITION AT PHI=PF
%
if( nbr==2 || nbr==5 )
jpf = fix(n);
elseif( nbr==3 || nbr==4 ) ;
wpf = .5;
jpf = fix(n + 1);
else;
jpf = fix(n);
end;
jpsp = fix(jps + 1);
jpfm = fix(jpf - 1);
nunk = fix(jpf - jps + 1);
fjj = jpfm - jpsp + 1;
%
% SCALE COEFFICIENTS FOR SUBROUTINE GENBUN
%
for i = its : itf;
cf = dphi2.*sint(i).*sint(i);
am(i) = cf.*am(i);
bm(i) = cf.*bm(i);
cm(i) = cf.*cm(i);
end; i = fix(itf+1);
am(its) = 0.;
cm(itf) = 0.;
ising = 0;
if( mbr~=2 && mbr~=3 && mbr~=5 && mbr~=6 && mbr~=8 )
if( nbr~=2 && nbr~=3 && nbr~=5 )
if( elmbda>=0 )
ising = 1;
summlv = wts.*wps + wts.*wpf + wtf.*wps + wtf.*wpf;
if( inp>0 )
summlv = summlv + wp;
end;
if( isp>0 )
summlv = summlv + wp;
end;
sum1 = 0.;
for i = itsp : itfm;
sum1 = sum1 + sint(i);
end; i = fix(itfm+1);
summlv = summlv + fjj.*(sum1+wts+wtf);
summlv = summlv +(wps+wpf).*sum1;
hne = summlv;
end;
end;
end;
if( mbr==1 )
elseif( mbr==2 || mbr==3 || mbr==8 ) ;
for j = jps : jpf;
f(2,j) = f(2,j) - at.*f(1,j);
end; j = fix(jpf+1);
elseif( mbr==4 || mbr==5 || mbr==9 ) ;
for j = jps : jpf;
f(1,j) = f(1,j) + tdt.*bdts(j).*at;
end; j = fix(jpf+1);
elseif( nbdcnd==3 ) ;
yhld = f(1,jps) - 4../(fn.*dphi.*dth2).*(bdpf(2)-bdps(2));
for j = 1 : np1;
f(1,j) = yhld;
end; j = fix(np1+1);
end;
if( mbr==1 )
elseif( mbr==2 || mbr==5 || mbr==6 ) ;
for j = jps : jpf;
f(m,j) = f(m,j) - ct.*f(m+1,j);
end; j = fix(jpf+1);
elseif( mbr==3 || mbr==4 || mbr==7 ) ;
for j = jps : jpf;
f(m+1,j) = f(m+1,j) - tdt.*bdtf(j).*ct;
end; j = fix(jpf+1);
elseif( nbdcnd==3 ) ;
yhld = f(m+1,jps) - 4../(fn.*dphi.*dth2).*(bdpf(m)-bdps(m));
for j = 1 : np1;
f(m+1,j) = yhld;
end; j = fix(np1+1);
end;
if( nbr==1 )
elseif( nbr==4 || nbr==5 ) ;
for i = its : itf;
f(i,1) = f(i,1) + tdp.*bdps(i)./(dphi2.*sint(i).*sint(i));
end; i = fix(itf+1);
else;
for i = its : itf;
f(i,2) = f(i,2) - f(i,1)./(dphi2.*sint(i).*sint(i));
end; i = fix(itf+1);
end;
if( nbr==1 )
elseif( nbr==3 || nbr==4 ) ;
for i = its : itf;
f(i,n+1) = f(i,n+1) - tdp.*bdpf(i)./(dphi2.*sint(i).*sint(i));
end; i = fix(itf+1);
else;
for i = its : itf;
f(i,n) = f(i,n) - f(i,n+1)./(dphi2.*sint(i).*sint(i));
end; i = fix(itf+1);
end;
pertrb = 0.;
if( ising~=0 )
summlv = wts.*wps.*f(its,jps) + wts.*wpf.*f(its,jpf)+ wtf.*wps.*f(itf,jps) + wtf.*wpf.*f(itf,jpf);
if( inp>0 )
summlv = summlv + wp.*f(1,jps);
end;
if( isp>0 )
summlv = summlv + wp.*f(m+1,jps);
end;
for i = itsp : itfm;
sum1 = 0.;
for j = jpsp : jpfm;
sum1 = sum1 + f(i,j);
end; j = fix(jpfm+1);
summlv = summlv + sint(i).*sum1;
end; i = fix(itfm+1);
sum1 = 0.;
sum2 = 0.;
for j = jpsp : jpfm;
sum1 = sum1 + f(its,j);
sum2 = sum2 + f(itf,j);
end; j = fix(jpfm+1);
summlv = summlv + wts.*sum1 + wtf.*sum2;
sum1 = 0.;
sum2 = 0.;
for i = itsp : itfm;
sum1 = sum1 + sint(i).*f(i,jps);
sum2 = sum2 + sint(i).*f(i,jpf);
end; i = fix(itfm+1);
summlv = summlv + wps.*sum1 + wpf.*sum2;
pertrb = summlv./hne;
for j = 1 : np1;
for i = 1 : mp1;
f(i,j) = f(i,j) - pertrb;
end; i = fix(mp1+1);
end; j = fix(np1+1);
end;
%
% SCALE RIGHT SIDE FOR SUBROUTINE GENBUN
%
for i = its : itf;
cf = dphi2.*sint(i).*sint(i);
for j = jps : jpf;
f(i,j) = cf.*f(i,j);
end; j = fix(jpf+1);
end; i = fix(itf+1);
[nbdcnd,nunk,dumvar3,munk,am(sub2ind(size(am),max(its,1)):end),bm(sub2ind(size(bm),max(its,1)):end),cm(sub2ind(size(cm),max(its,1)):end),idimf,f(sub2ind(size(f),its,jps):end),ierror,d]=genbun(nbdcnd,nunk,1,munk,am(sub2ind(size(am),max(its,1)):end),bm(sub2ind(size(bm),max(its,1)):end),cm(sub2ind(size(cm),max(its,1)):end),idimf,f(sub2ind(size(f),its,jps):end),ierror,d);
while (1);
if( ising>0 )
if( inp<=0 )
if( isp>0 )
for j = 1 : np1;
f(m+1,j) = 0.;
end; j = fix(np1+1);
break;
end;
elseif( isp<=0 ) ;
for j = 1 : np1;
f(1,j) = 0.;
end; j = fix(np1+1);
break;
end;
end;
if( inp>0 )
summlv = wps.*f(its,jps) + wpf.*f(its,jpf);
for j = jpsp : jpfm;
summlv = summlv + f(its,j);
end; j = fix(jpfm+1);
dfn = cp.*summlv;
dnn = cp.*((wps+wpf+fjj).*(sn(2)-1.)) + elmbda;
dsn = cp.*(wps+wpf+fjj).*sn(m);
if( isp<=0 )
cnp =(f(1,1)-dfn)./dnn;
for i = its : itf;
hld = cnp.*sn(i);
for j = jps : jpf;
f(i,j) = f(i,j) + hld;
end; j = fix(jpf+1);
end; i = fix(itf+1);
for j = 1 : np1;
f(1,j) = cnp;
end; j = fix(np1+1);
break;
end;
elseif( isp<=0 ) ;
break;
end;
summlv = wps.*f(itf,jps) + wpf.*f(itf,jpf);
for j = jpsp : jpfm;
summlv = summlv + f(itf,j);
end; j = fix(jpfm+1);
dfs = cp.*summlv;
dss = cp.*((wps+wpf+fjj).*(ss(m)-1.)) + elmbda;
dns = cp.*(wps+wpf+fjj).*ss(2);
if( inp<=0 )
csp =(f(m+1,1)-dfs)./dss;
for i = its : itf;
hld = csp.*ss(i);
for j = jps : jpf;
f(i,j) = f(i,j) + hld;
end; j = fix(jpf+1);
end; i = fix(itf+1);
for j = 1 : np1;
f(m+1,j) = csp;
end; j = fix(np1+1);
else;
rtn = f(1,1) - dfn;
rts = f(m+1,1) - dfs;
if( ising>0 )
csp = 0.;
cnp = rtn./dnn;
elseif( abs(dnn)<=abs(dsn) ) ;
den = dns - dss.*dnn./dsn;
rtn = rtn - rts.*dnn./dsn;
csp = rtn./den;
cnp =(rts-dss.*csp)./dsn;
else;
den = dss - dns.*dsn./dnn;
rts = rts - rtn.*dsn./dnn;
csp = rts./den;
cnp =(rtn-csp.*dns)./dnn;
end;
for i = its : itf;
hld = cnp.*sn(i) + csp.*ss(i);
for j = jps : jpf;
f(i,j) = f(i,j) + hld;
end; j = fix(jpf+1);
end; i = fix(itf+1);
for j = 1 : np1;
f(1,j) = cnp;
f(m+1,j) = csp;
end; j = fix(np1+1);
end;
break;
end;
if( nbdcnd==0 )
for i = 1 : mp1;
f(i,jpf+1) = f(i,jps);
end; i = fix(mp1+1);
end;
f_shape=zeros(f_shape);f_shape(:)=f(1:numel(f_shape));f=f_shape;
bdts_shape=zeros(bdts_shape);bdts_shape(:)=bdts(1:numel(bdts_shape));bdts=bdts_shape;
bdtf_shape=zeros(bdtf_shape);bdtf_shape(:)=bdtf(1:numel(bdtf_shape));bdtf=bdtf_shape;
bdps_shape=zeros(bdps_shape);bdps_shape(:)=bdps(1:numel(bdps_shape));bdps=bdps_shape;
bdpf_shape=zeros(bdpf_shape);bdpf_shape(:)=bdpf(1:numel(bdpf_shape));bdpf=bdpf_shape;
am_shape=zeros(am_shape);am_shape(:)=am(1:numel(am_shape));am=am_shape;
bm_shape=zeros(bm_shape);bm_shape(:)=bm(1:numel(bm_shape));bm=bm_shape;
cm_shape=zeros(cm_shape);cm_shape(:)=cm(1:numel(cm_shape));cm=cm_shape;
ss_shape=zeros(ss_shape);ss_shape(:)=ss(1:numel(ss_shape));ss=ss_shape;
sn_shape=zeros(sn_shape);sn_shape(:)=sn(1:numel(sn_shape));sn=sn_shape;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
sint_shape=zeros(sint_shape);sint_shape(:)=sint(1:numel(sint_shape));sint=sint_shape;
end %subroutine hwsss1
%DECK HWSSSP