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

% 
% The following features are updated by Ph.D. candidate Jun Tang, College
% of Underwater Acoustic Engineering, Harbin Engineering University,
% Harbin, China, 2017.
% 1. This version is applicable to wedges with elastic bottom. To reduce to
%    fluid bottom, users can simply set cs to 0 in the launcher files.
% 2. Fields of each image source are solely saved in independent files for
%    possible future applications. 

%% 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 = [40 90];

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

isave = pars.isave;


%% auxiliary variables

% n receivers

nr = length(xr);
pr(1:nr) = 0;          % total field
prnl(1:nr) = 0;        % field of a certain image source defined by nb and ll

R = sqrt( (xr - xs).^2 + (yr-ys).^2  + (zr - zs).^2 );
prnl(1:nr) = exp(1i*k*R)./R;                        
ranges_plot = sqrt((xr-xr(1)).^2 + (yr-yr(1)).^2 + (zr-zr(1)).^2)/1000;

if isave > 0
dlmwrite(['nb=' num2str(0) ',' 'll=' num2str(-1) '.grid'],[ranges_plot' real(prnl)' imag(prnl)'],'\t');
end % save the field generated by S0,-1, which is the real source

pr = pr + prnl;

R = sqrt( (xr - xs).^2 + (yr-ys).^2  + (zr + zs).^2 ); 
prnl(1:nr) = -exp(1i*k*R)./R;      
ranges_plot = sqrt((xr-xr(1)).^2 + (yr-yr(1)).^2 + (zr-zr(1)).^2)/1000;

if isave > 0
dlmwrite(['nb=' num2str(0) ',' 'll=' num2str(2) '.grid'],[ranges_plot' real(prnl)' imag(prnl)'],'\t');
end % save the field generated by S0,2, which is the surface image of the real source

pr = pr + prnl;

% % % % % 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;

ll_vals = [-2 -1 1 2]; % ll = -2,-1,1,2

for nb = 1:nb_max
        
        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;    % different from that in Petrov's code, so is line 195.
        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;
    
    for j_ll = 1:4  
        
        ll = ll_vals(j_ll);
        % 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( (xrb - xis) + 1i*Rdir.*(yr-yis) ); % position of Rdir is different from that in Petrov's code
        
        
        
        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);  
        
        if isave > 0
        dlmwrite(['nb=' num2str(nb) ',' 'll=' num2str(ll) '.grid'],[ranges_plot' real(prnl)' imag(prnl)'],'\t');
        end; % save the field of the specific image source

        pr = pr + prnl;
        
    end;

end;   
    
disp('COMPLETED!');
toc;

    function V_phi_b = BotRefl
        
        c1 = MPars(1,1); 
        c2 = MPars(1,2); % P-wave speed
        c3 = MPars(1,3); % S-wave speed
        rho1 = MPars(2,1);
        rho2 = MPars(2,2);
        
        gam1 = k*cos_phi_b;
        gam2 = k*sqrt( (c1/c2)^2 -1 + cos_phi_b.^2); 
        gam3 = k*sqrt( (c1/c3)^2 -1 + cos_phi_b.^2);  
        
        DirectCond = ( real(gam1).*real(gam2) < 0);
        
        gam2(DirectCond) = - gam2(DirectCond);
        
        DirectCond = ( real(gam1).*real(gam3) < 0);
        
        gam3(DirectCond) = - gam3(DirectCond); 
        
        
        d = (1-2*((c3/c1)^2)*(1-cos_phi_b.^2));
        d = d.*d; % d = ( cos (2*phi_s) )^2, where phi_s is the trnamitted angle of S-wave in the bottom

        V_phi_b = ( rho2*d.*gam1.*gam3 + rho2*(1-d).*gam1.*gam2 - rho1*gam2.*gam3)./ (rho2*d.*gam1.*gam3 + rho2*(1-d).*gam1.*gam2 + rho1*gam2.*gam3); 
        % reflection coefficient of fluid-solid interface
        
        if abs(c3) < 10e-30 
            V_phi_b = (  (rho2/rho1)*gam1 - gam2  )./( (rho2/rho1)*gam1 + gam2 );
        end; % reflection coefficient of fluid-fluid interface is used when c3 = 0 
            
    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





