function [A,B] = abmat(NW,depw,rhow,grhow,cwnsw,cgwnsw,HW,...
                         NB,depb,rhob,grhob,cwnsb,cgwnsb,HB,...
                         NL,enl,fnl,gnl,hnl,M,FMAXWB)
% The following function computes the complex symetric matrix A
% (not Hermitian) and positive definite symetric matrix B needed 
% in the Galerkin method to form the generalized eigenvalue problem 
% that is solved to find eigenvalued and eigenfunctions, in terms of 
% Legerdre polynomials. NL is the over all number of Legendre polynomials
% used to form a double basis of dimension M=2*NL-2. The coefficients
% enl, fnl, gnl and hnl are chosen so that the functions in the
% combined basis satisfy the boundary and interface conditions.
% The number of Gauss-Legendre quadrature points is set to NL+2,
% in each layer.
%
NWM1=NW-1; % # of layers in the water
NBM1=NB-1;% # of layers in the bottom
%
NLM1=NL-1;
NLP1=NL+1;
NLP2=NL+2;
%
J=1:NLP2; % Vectorize the index in Gauss-Legendre quadrature
 %
hwo2=HW/2.0;
hbmhwo2=(HB-HW)/2.0;
hbphwo2=(HB+HW)/2.0;
%
% Find zeros and weights for Gauss-Legendre integration in each layer
% in the water and bottom. Define the reciprocal of the density and the
% complex wave number squared at the zeros. There are NL+2 zeros.
%
zw=zeros(NLP2,NWM1);
ww=zeros(NLP2,NWM1);
xw=zeros(NLP2,NWM1);
rzw=zeros(NLP2,NWM1);
qzw=zeros(NLP2,NWM1)+1i*zeros(NLP2,NWM1);
%
for L=1:NWM1
z1=depw(L);
z2=depw(L+1);
D=[z1 z2];
[z,w]=legpts(NLP2,D);
%
zw(J,L)=z(J);
ww(J,L)=w(J);
%
   xw(J,L)=(zw(J,L)-hwo2)/hwo2;
   rzw(J,L)=exp(-grhow(L)*(zw(J,L)-z1))/rhow(L);
   qzw(J,L)=FMAXWB-(cwnsw(L)+cgwnsw(L)*zw(J,L));
%
end
%
zb=zeros(NLP2,NBM1);
wb=zeros(NLP2,NBM1);
xb=zeros(NLP2,NBM1);
rzb=zeros(NLP2,NBM1);
qzb=zeros(NLP2,NWM1)+1i*zeros(NLP2,NWM1);
%
for L=1:NBM1
z1=depb(L);
z2=depb(L+1);
D=[z1 z2];
[z,w]=legpts(NLP2,D);
zb(J,L)=z(J);
wb(J,L)=w(J);
%
   xb(J,L)=(zb(J,L)-hbphwo2)/hbmhwo2;
   rzb(J,L)=exp(-grhob(L)*(zb(J,L)-z1))/rhob(L);
   qzb(J,L)=FMAXWB-(cwnsb(L)+cgwnsb(L)*zb(J,L));
%
end
%
% Compute the matrices A and B
%
A=zeros(M,M)+1i*zeros(M,M);
B=zeros(M,M);
pw1=zeros(NLP2,1);
pw2=zeros(NLP2,1);
pb1=zeros(NLP2,1);
pb2=zeros(NLP2,1);
dw1=zeros(NLP2,1);
dw2=zeros(NLP2,1);
db1=zeros(NLP2,1);
db2=zeros(NLP2,1);
%
% Build A and B entry by entry. This is a vestige Fortran.
%
for I=1:M % DO 500 
    for II=I:M % DO 400 with II >= I
%
% Sum over water layers
%
      for L=1:NWM1 % DO 200
      xwJ=xw(J,L);
      [p,d]=legpoly(NL,xwJ);
      %
      if I <= NLM1 % I=1:NL-1
      %
      ik=I-1; % ik is the order of the polynomial
      pw1(J)=p(J,ik+1)+(enl(ik+1)+(-1)^(ik+NL))*p(J,NL)+ ...
             enl(ik+1)*p(J,NLP1);
      dw1(J)=(d(J,ik+1)+(enl(ik+1)+(-1)^(ik+NL))*d(J,NL)+ ...
              enl(ik+1)*d(J,NLP1))/hwo2;
      %
      else % I=NL:2*NL-2
      %
      ik=I-NL; % ik is the order of the polynomial
      pw1(J)=fnl(ik+1)*(p(J,NL)+p(J,NLP1));
      dw1(J)=fnl(ik+1)*(d(J,NL)+d(J,NLP1))/hwo2;
      %
      end
%
      if II <= NLM1 % II=1:NL-1
      %
      iik=II-1; % iik is the order of the polynomial
      pw2(J)=p(J,iik+1)+(enl(iik+1)+(-1)^(iik+NL))*p(J,NL)+ ...
             enl(iik+1)*p(J,NLP1);
      dw2(J)=(d(J,iik+1)+(enl(iik+1)+(-1)^(iik+NL))*d(J,NL)+ ...
              enl(iik+1)*d(J,NLP1))/hwo2;      
      %
      else % II=NL:2*NL-2
      %
      iik=II-NL; % iik is the order of the polynomial
      pw2(J)=fnl(iik+1)*(p(J,NL)+p(J,NLP1));
      dw2(J)=fnl(iik+1)*(d(J,NL)+d(J,NLP1))/hwo2;
      end
      %
%
      A(II,I)=A(II,I)+sum(ww(J,L).*rzw(J,L).*pw1(J).*pw2(J).*qzw(J,L));
      A(II,I)=A(II,I)+sum(ww(J,L).*rzw(J,L).*dw1(J).*dw2(J));
      B(II,I)=B(II,I)+sum(ww(J,L).*rzw(J,L).*pw1(J).*pw2(J));
%
      end % 200 CONTINUE
%
%     Sum over bottom layers
%
      for L=1:NBM1 % DO 300
      xbJ=xb(J,L);
      [p,d]=legpoly(NL,xbJ);
      %
      if I <= NLM1
      %
      ik=I-1;
      pb1(J)=gnl(ik+1)*(-p(J,NL)+p(J,NLP1));
      db1(J)=gnl(ik+1)*(-d(J,NL)+d(J,NLP1))/hbmhwo2;
      %
      else
      %
      ik=I-NL;
      pb1(J)=p(J,ik+1)-(hnl(ik+1)+1.0d0)*p(J,NL)+ ...
             hnl(ik+1)*p(J,NLP1);
      db1(J)=(d(J,ik+1)-(hnl(ik+1)+1.0d0)*d(J,NL)+ ...
              hnl(ik+1)*d(J,NLP1))/hbmhwo2;
      end
      %
      if II <= NLM1
      %
      iik=II-1;
      pb2(J)=gnl(iik+1)*(-p(J,NL)+p(J,NLP1));
      db2(J)=gnl(iik+1)*(-d(J,NL)+d(J,NLP1))/hbmhwo2;
      %
      else
      %
      iik=II-NL;
      pb2(J)=p(J,iik+1)-(hnl(iik+1)+1.0d0)*p(J,NL)+ ...
             hnl(iik+1)*p(J,NLP1);
      db2(J)=(d(J,iik+1)-(hnl(iik+1)+1.0d0)*d(J,NL)+ ...
              hnl(iik+1)*d(J,NLP1))/hbmhwo2;
      %
      end
%      
      A(II,I)=A(II,I)+sum(wb(J,L).*rzb(J,L).*pb1(J).*pb2(J).*qzb(J,L));
      A(II,I)=A(II,I)+sum(wb(J,L).*rzb(J,L).*db1(J).*db2(J));
      B(II,I)=B(II,I)+sum(wb(J,L).*rzb(J,L).*pb1(J).*pb2(J));
%
      end % 300 CONTINUE
%
%     Define the sub-diagonal entries by symetry
%
      if II ~= I 
          A(I,II)=A(II,I);
          B(I,II)=B(II,I);
      end
%
    end % 400 CONTINUE
end % 500 CONTINUE