function pr = PenWedgeField( pars )

% PenWedgeField( ) function that implements the image source solution
% of G. Deane and M. Buckingham (see J. Acoust. Soc. Am. V. 93. No. 3. 1993, Pp. 1319-1328)
% The code is written by P.S. Petrov, V.I.Il'ichev Pacific Oceanological
% institute, Vladivostok, Russia, 2015.
% Take it and use it for any peaceful purposes :-)

% In the comments below I use the references to the equations from the
% paper of G. Deane and M. Buckingham


%% input
% receiver data: an array of the size [nr x 3]

xr = pars.rec(1,:);
yr = pars.rec(2,:);
zr = pars.rec(3,:);

% sourcre position: vector [1 x 3]

xs = pars.src(1);
ys = pars.src(2);
zs = pars.src(3);

yr = yr - ys;   % if ys <> 0, make shift the coordinate system to make ys = 0.

% media parameters (sound speeds, densities, etc -- are used for reflection
% coefficient computation)

MPars = pars.MPars;

% wavenumber = 2*Pi*f/cw;

k = pars.k;

% wedge angle

alpha_w = pars.alpha_w;

% parameters of the image source solution

nb_max = 10;     % number of possible reflections from the bottom 
                 % nb_max ~ 10 is OK to compute TL along the slope
                 % nb_max ~ 30 (full circle) is required for the
                 % cross-slope benchmark
                 
if isfield(pars,'nb_max')
    nb_max = pars.nb_max;
end;

nnu = 5;       % number of Bessel functions in the series; for a given value of \nu, the 2\nu + 2 Bessel terms will be used
                 % for the ASA benchmark along the slope (as described in
                 % the paper by Jensen and Ferla) nnu ~ 5 is OK, for the
                 % cross-slope full 3D computation nnu ~ 100 is required,
                 % this is very slow.

if isfield(pars,'nnu')
    nnu = pars.nnu;
end;
                 
                 
                 
OmThre = 5;    % for \Omega < OmThre, the Bessel functions will be computed by the call of Matlab function besselj()
                 % for \Omega >= OmThre the recurrence relations will be
                 % used to increase the efficiency
                 % nnu should be roughly from nnu to 2*nnu

if isfield(pars,'OmThre')
    OmThre = pars.OmThre;
end;                 


TLlims = [-90 -40];

if isfield(pars,'TLlims')
    TLlims = pars.TLlims;
end;                 




%% auxiliary variables

% n receivers

nr = length(xr);
pr(1:nr) = 0;


% source cylindrical coordinates

alpha_s = alpha_w - atan(zs/xs);
Rs = sqrt(xs^2 + zs^2);

% integration variables

nph = 1001;
phi = linspace(0,pi,nph);   % phi values
dph = phi(2) - phi(1);

nth1 = 3001;    % number of integration points WRT theta on [0, pi/2], see Eq. (23) from the paper
nth2 = 4001;    % number of integration points WRT theta on [pi/2, pi/2 - \infty], see Eq. (23) from the paper
nth = nth1 + nth2;
Mth = pi/2;    % theta "infinity" value
theta(1:nth) = [linspace(0,pi/2,nth1)  linspace(pi/2,pi/2 - 1i*Mth,nth2)];    % theta values
dth1 = theta(2) - theta(1);
dth2 = theta(end) - theta(end-1);


V_totalm(1:nth,1:nph) = 1;
V_totalp(1:nth,1:nph) = 1;

ns = 0;
prnl(1:nr) = 0;

tic;


for nb = 0:nb_max
    
    ll_vals = [-1 2];
    if nb > 0
        ll_vals = [-2 -1 1 2];
        
        Theta_b = 2*alpha_w*(nb-1);
        cos_phi_b(1:nth,1:nph) = sin(theta.')*cos(phi)*sin(Theta_b) + cos(theta.')*ones(size(phi))*cos(Theta_b);
        V_totalm(1:nth,1:nph) = V_totalm.*BotRefl;
        
        Theta_b = -2*alpha_w*nb;
        cos_phi_b(1:nth,1:nph) = sin(theta.')*cos(phi)*sin(Theta_b) + cos(theta.')*ones(size(phi))*cos(Theta_b);
        V_totalp(1:nth,1:nph) = V_totalp.*BotRefl;
        
        
    end;
    
    
    
    for ll = ll_vals
        
        % the angular positions of the image sources, Eq. (25) 
        % please note that the expressions given in the paper are slightly
        % incorrect
        
        if ll == -2
            Phi_nl = -alpha_s - 2*alpha_w*(nb-1);
            ns = nb - 1;
        elseif ll == -1
            Phi_nl = alpha_s - 2*alpha_w*nb;
            ns = nb;
        elseif ll == 1
            Phi_nl = alpha_s + 2*alpha_w*nb; % in the paper this one is incorrect!
            ns = nb;
        elseif ll == 2
            Phi_nl = -alpha_s + 2*alpha_w*(nb+1); % in the paper this one is incorrect!
            ns = nb + 1;
        end;
                
        xis = Rs*cos(Phi_nl); % in the Cartesian coordinates associated with surface!
        yis = 0;
        zis = Rs*sin(Phi_nl);
        
        Rmatr = [[cos(alpha_w)  -sin(alpha_w)]; [sin(alpha_w)  cos(alpha_w)]];
        CTrans = Rmatr*[  xr; -zr ];
        
        xrb = CTrans(1,:);
        zrb = CTrans(2,:);
        
        
        Rdir = sign(zrb - zis);       % the formulae work for cos(zeta) > 0, thus we should choose the Cartesian coordinates so that projection of \bar(R) onto z axis is positive.
        
        R = sqrt( (xrb - xis).^2 + (yr-yis).^2  + (zrb - zis).^2 );
        zeta = acos( Rdir.*(zrb - zis)./R );
        xi = angle( Rdir.*(xrb - xis) + 1i*(yr-yis) );
        
        
        
        Omega_z(1:nth,1:nr) = (cos(theta).')*( k*R(1:nr).*cos( zeta(1:nr) ) );
        Omega(1:nth,1:nr) = (sin(theta).')*( k*R(1:nr).*sin( zeta(1:nr) ) ); %???
        
        J2nu(1:nth,1:nr) = besselj(0,Omega);
        J2nup1(1:nth,1:nr) = besselj(1,Omega);
        
        
        absOmSmall = (abs(Omega)<OmThre);
        
        WPnl(1:nth,1:nr) = 0;
        
        for nu = 0:nnu
            
            disp(['Computing solution: nb=' int2str(nb) '; l=' int2str(ll) '; nu=' int2str(nu)]);
            toc;
            
            if ll<0
                %W(1:nth,1) = dph*trapz(V_total(1:nth,1:nph).*repmat(cos(2*nu*phi),nth,1),2);
                W(1:nth,1) = dph*int_simpson(V_totalm(1:nth,1:nph).*repmat(cos(2*nu*phi),nth,1),2);
                a2nu = 2*((-1)^ns)*W(1:nth,1)*cos(  2*nu*xi(1:nr) );
                %W(1:nth,1) = dph*trapz(V_total(1:nth,1:nph).*repmat(cos((2*nu+1)*phi),nth,1),2);
                W(1:nth,1) = dph*int_simpson(V_totalm(1:nth,1:nph).*repmat(cos((2*nu+1)*phi),nth,1),2);
                a2nup1 = 2*((-1)^ns)*W(1:nth,1)*cos(  (2*nu+1)*xi(1:nr) );
            else
                
                W(1:nth,1) = dph*int_simpson(V_totalp(1:nth,1:nph).*repmat(cos(2*nu*phi),nth,1),2);
                a2nu = 2*((-1)^ns)*W(1:nth,1)*cos(  2*nu*xi(1:nr) );
                %W(1:nth,1) = dph*trapz(V_total(1:nth,1:nph).*repmat(cos((2*nu+1)*phi),nth,1),2);
                W(1:nth,1) = dph*int_simpson(V_totalp(1:nth,1:nph).*repmat(cos((2*nu+1)*phi),nth,1),2);
                a2nup1 = 2*((-1)^ns)*W(1:nth,1)*cos(  (2*nu+1)*xi(1:nr) );
                
                
            end;
            
            
            
            if nu == 0
                a2nu = a2nu/2;
            end;
            
            WPnl(1:nth,1:nr) = WPnl(1:nth,1:nr) + ((-1)^nu)*(  a2nu.*J2nu + 1i*a2nup1.*J2nup1  );  % expression under integral wrt theta
            
            %aJ2nu(1:nth,1:nr) = besselj(2*(nu+1),Omega);
            %aJ2nup1(1:nth,1:nr) = besselj(2*(nu+1)+1,Omega);
            J2nu(1:nth,1:nr) = (-J2nu(1:nth,1:nr).*Omega + 2*(2*nu+1)*J2nup1(1:nth,1:nr) )./Omega;
            J2nup1(1:nth,1:nr) = -J2nup1(1:nth,1:nr) + 2*(2*nu+2)*J2nu(1:nth,1:nr)./Omega;
            %J2nu((Omega==0)) = 0;   % recurrence relations do not work for zero argument! The J should be put to zero where Omega = 0.
            %J2nup1((Omega==0)) = 0; % when Omega is close to zero, rcurrence relations also produce errors!!
            J2nu(absOmSmall) = besselj(2*(nu+1),Omega(absOmSmall));
            J2nup1(absOmSmall) = besselj(2*(nu+1)+1,Omega(absOmSmall));
            
            
        end;
        
        WPnl(1:nth,1:nr) = (1i*k/pi)*WPnl.*exp( 1i*Omega_z ).*repmat(sin( theta.'),1,nr);
        
        %prnl(1:nr) = dth1*trapz(WPnl(1:nth1,1:nr),1) + dth2*trapz(WPnl(nth1+1:nth,1:nr),1);
        prnl(1:nr) = dth1*int_simpson(WPnl(1:nth1,1:nr),1) + dth2*int_simpson(WPnl(nth1+1:nth,1:nr),1);
        
        pr = pr + prnl;
        
    end;

figure;
plot( sqrt((xr-xr(1)).^2 + (yr-yr(1)).^2 + (zr-zr(1)).^2)/1000 ,20*log10(abs(pr)),'linewidth',2);
ylim(TLlims);
xlabel('r (from the source), km');
ylabel('TL, dB re 1 m');
grid on;
title(['The image source solution for ' int2str(nb) ' bottom reflections']);

end;

    
disp('COMPLETED!');
toc;

    function V_phi_b = BotRefl
        
        c1 = MPars(1,1);
        c2 = MPars(1,2);
        rho1 = MPars(2,1);
        rho2 = MPars(2,2);
        
        gam1 = k*cos_phi_b;
        gam2 = k*sqrt( (c1/c2)^2 -1 + cos_phi_b.^2); 
        
        DirectCond = ( real(gam1).*real(gam2) + imag(gam1).*imag(gam2) < 0);
        
        gam2(DirectCond) = - gam2(DirectCond);
        
        V_phi_b = (  (rho2/rho1)*gam1 - gam2  )./( (rho2/rho1)*gam1 + gam2 );
        
    end

    function Fint = int_simpson(ffun,fdim)
       
        [nx, ny] = size(ffun);
        
        if fdim == 1
            inds = 1:fix( (nx-1)/2);
            Fint(1,1:ny) = sum( ffun(2*inds-1,1:ny) + 4*ffun(2*inds,1:ny) + ffun(2*inds+1,1:ny),1 )/3;  
            
            if rem(nx-1,2)
                Fint(1,1:ny) = Fint + 0.5*(ffun(end-1,1:ny) + ffun(end,1:ny));
            end;
            
        elseif fdim == 2
            inds = 1:fix( (ny-1)/2);
            Fint(1:nx,1) = sum( ffun(1:nx,2*inds-1) + 4*ffun(1:nx,2*inds) + ffun(1:nx,2*inds+1),2 )/3;  
            
            if rem(ny-1,2)
                Fint(1:nx,1) = Fint + 0.5*(ffun(1:nx,end-1) + ffun(1:nx,end-1));
            end;
        else
            error('In function int_simpson argument fdim must be either 1 or 2!');
        end;
            
            
    end

end





