function [zr,cdfr,rkm,kh,cpr,NRNG] = ...
    rimLG(FREQ,ZS,NL,MS,NDEP,ZMIN,ZINC,RMIN,RMAX,RINC, ...
    HW,HB,depw,svpw,dbpwlw,rhow,NW,depb,svpb,dbpwlb,rhob,NB)
%
% Run this script by typing "run_rimLG" at the Matlab prompt. The
% scipt run_rimLG.m automates the process of running this script.
%
% Compute the range independent modal (rim) acoustic field using the 
% Legendre-Galerkin method: Setup index of refraction table, compute
% the matrices A and B, solve the generalized eigenvlaue problem
% and compute the acoustic field.  
% Calls table_gen, abmat(legpts, legpoly), csrt, basic_depth_fct 
% and the Matlab function 'eig.'
%
% See Legendre_Galerkin.pdf for documentation on the computational
% method. (R. B. Evans, N. Stonington, CT, Jan. 5, 2016)
%
% Compute an piecewise linear representation of the complex wave
% number squared (intercepts and gradients). Also find the logarithmic
% gradients in the density
%
% Water column
[cwnsw,cgwnsw,grhow,ROHSW]=table_gen(FREQ,ZS,depw,svpw,dbpwlw,rhow,NW);
% Define the density at the bottom or the water column
ROHW=rhow(NW); % A spelling error from Couple is retained and captilized
FMAXW=max(real(cwnsw));
%
%Bottom sediment
[cwnsb,cgwnsb,grhob,ROHSB]=table_gen(FREQ,ZS,depb,svpb,dbpwlb,rhob,NB);
% Define the density at the top of the bottom sediment
ROHB=rhob(1);
FMAXB=max(real(cwnsb));
%
FMAXWB=max([FMAXW FMAXB]);
%
% Define the density at the source
if (depw(1) < ZS) && (ZS <= depw(NW))
    ROHS=ROHSW;
elseif (depb(1) < ZS) && (ZS <= depb(NB))
    ROHS=ROHSB;
else
    fprintf('%s\n','ZS not in water or bottom')
    fprintf('%s\n','Type "dbquit" to quit:')
    keyboard
end
%
% Compute coefficient in the dual basis of Legrendre polynomials, which
% enforce the boundary conditions and the continuity conditions at the 
% water sediment interface
%
RW=1.0D0/(ROHW*HW);
RB=1.0D0/(ROHB*(HB-HW));
TWORS=2.0D0*(RW+RB);
NLM1=NL-1;
NLSQ=NL^2;
enl=zeros(NLM1,1);
fnl=zeros(NLM1,1);
gnl=zeros(NLM1,1);
hnl=zeros(NLM1,1);
% The coefficients depends on whether NL is even or odd.
for II=1:NLM1
    I=II-1;
    J=(-1)^I;
    L=(-1)^NL; % sign change for NL odd
    enl(II)=(-L*RB*(J+L) - (RW/NLSQ)*( NLM1*NL*J*L + I*(I+1) ) )/TWORS;
    fnl(II)= (  RB*(J+L) - (RB/NLSQ)*( NLM1*NL*L + I*(I+1)*J ) )/TWORS;
    gnl(II)= (  RW*(J+L) - (RW/NLSQ)*( NLM1*NL*J + I*(I+1)*L ) )/TWORS;
    hnl(II)=(-L*RW*(J+L) - (RB/NLSQ)*( NLM1*NL + I*(I+1)*J*L ) )/TWORS;
end
%
% The number of functions in the dual basis is M = 2*NL-2 = 2*(NL-1).
% Compute the MxM matrices A and B in the generalize matrix eigenvalue
% problem. A is a complex symmetric matrix (A = A.'). A is not a complex 
% conjugate symmetric matrix (A ~= A'). B is a real symmetric positive
% definite matrix.
%
M=2*NL-2;
[A,B]=abmat(NW,depw,rhow,grhow,cwnsw,cgwnsw,HW, ...
            NB,depb,rhob,grhob,cwnsb,cgwnsb,HB, ...
            NL,enl,fnl,gnl,hnl,M,FMAXWB);
% Find generalized eigenvalues and eigenvectors
[V,D]=eig(A,B,'chol');
%[V,D]=eig(A,B);
egvl=diag(D);
clear D;
% Test plot of eigenvalues commented out
%figure
%plot(real(egvl),imag(egvl),'bo')
%xlim([0.0 1.0]);
%xlabel('real(egvl)');
%ylabel('imag(egvl)');
%title('Matrix Eigenvalues');
%
% Make eigenvector orthonormal so that V is orthogonal in the sense
% that V.'*B*V=I (the energy inner product). This gives the correct
% normalization to the modes
d=diag(V.'*B*V);
sqrtd=sqrt(d);
rsqrtd=1.0./sqrtd;
D=diag(rsqrtd);
V=V*D;
% Compute horozontal wave numbers. Use a square root with a positive
% imaginary part. The offset is the same as in abmat.m
kh=csrt(FMAXWB-egvl);
% Sort by descending real parts of the horizontal wave numbers. This is 
% the uaual ordering of normal modes in underwater acoustics.
[s_rekh,J]=sort(real(kh),1,'descend');
kh=kh(J);
% Reorganize the columns of V in the same order as kh
V=V(1:M,J);
% Construct the basic depth functions at the source and receivers
[bdfs]=basic_depth_fct(ZS,HW,HB,NL,enl,fnl,gnl,hnl,M);
bdfr=zeros(M,NDEP);
zr=zeros(NDEP,1);
for j=1:NDEP
    zr(j)=ZMIN+(j-1)*ZINC;
    [bdf]=basic_depth_fct(zr(j),HW,HB,NL,enl,fnl,gnl,hnl,M);
     bdfr(1:M,j)=bdf;
end
% Define the complex depth functions at the source and receivers
cdfs=bdfs'*V/ROHS; % Divide by the density at ZS
cdfr=bdfr'*V;
% Define ranges to calculate complex pressure
rkm=RMIN:RINC:RMAX;
NRNG=length(rkm);
%
% Compute the complex field and transmission loss by summing over MS
% modes at each range and depth.
%
cpr=zeros(NDEP,NRNG);
%
for n=1:NRNG
rng=1000.0*rkm(n);
crf=sqrt(1.0/(8.0*pi*rng))*exp(1i*kh.'*rng)./csrt(kh.');
   for j=1:NDEP
   % Here is the modal summation over MS
   cpr(j,n)=sum(cdfs(1:MS).*cdfr(j,1:MS).*crf(1:MS));
   end
end