% A launcher for the 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 :-)

% THIS IS THE SCRIPT FOR THE CASE OF PROPAGATION ALONG THE SLOPE

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

close all;
clear variables;
clc;

set(0, 'DefaultAxesFontSize', 18, 'DefaultAxesFontName', 'Arial');
set(0, 'DefaultTextFontSize', 18, 'DefaultTextFontName', 'Arial');

% track direction: 
% theta is the angle between the direction of Y axis and the track
% theta = 0  -- track is aligned across the slope (i.e. parallel to wedge apex)
% theta = -90 deg -- track is aligned along the slope (orthogonal to the
% apex, from x=x_s, y=0 to x=0, y=0). See pseudographics below.

%                 x=x_s            x=0 -- wedge apex
%                  |----------------|
%                  |                |
%   X direction <--*source----------|x=0,y=0
%                  |                |
%                  |----------------|
%                                   V
%                                   Y direction   

% track alignment

theta = -90;    
R_max = 4000;   % max range
nr = 500;       % number of receivers

% sourcre position: vector [1 x 3]

xs = 4000;  
ys =0;
zs = 100;
freq = 25; % source frequency, Hz

pars.src = [xs ys zs]; 

% receiver positions

xr = xs + sin(pi*theta/180)*linspace(0,R_max,nr);
yr = ys + cos(pi*theta/180)*linspace(0,R_max,nr);
zr = 30 + linspace(0,0,nr);

pars.rec = [xr; yr; zr];

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

cw = 1500;  % sound speed in water
cb = 1700;  % sound speed in bottom

beta = 0.5; % bottom attenuation, dB/wavelength
eta = 1/(40*pi*log10(exp(1)));
cb = cb/(1+1i*beta*eta); % complex bottom sound speed, including attenuation

pars.MPars = [[cw cb]; [1 1.5]];

% wavenumber = 2*Pi*f/cw;

pars.k = 2*pi*freq/1500;

% wedge angle = 2.8624 degree, ASA benchmark

pars.alpha_w = pi*2.8624/180;

% approximation parameters
% the parameters below are set for the propagation along the slope
% for cross-slope, see comments in PenWedgeField( )

pars.nnu = 2;      % number of terms in the Bessel series, Eq. (23)
pars.nb_max = 9;   % number of possible reflection from the bottom (total number of image sources is 4*nb_max + 2, see Eq. (25) )
pars.OmThre = 3;   % a threshold value of \Omega in Eq. (23). See comments in the function PenWedgeField( ) for details


pr = PenWedgeField( pars );

% plot the results

ranges_plot = sqrt((xr-xr(1)).^2 + (yr-yr(1)).^2 + (zr-zr(1)).^2)/1000;
TL_plot = 20*log10(abs(pr));

figure;
plot(ranges_plot, TL_plot,'linewidth',2);
ylim([-90 -30]);
xlabel('r (from the source), km');
ylabel('TL, dB re 1 m');
grid on;

% comparison against COUPLE results

TLCouple = dlmread('CoupleAtt.TL');
hold on;
plot(TLCouple(:,1),-TLCouple(:,2),'r--','linewidth',2);
legend('Image source method','COUPLE');

dlmwrite('ImageSourcesAlongSlopeAtt.TL',[ranges_plot' TL_plot'],'\t');
