% 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 CROSS-SLOPE CASE

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

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.

% THIS SCRIPT IS FOR theta = 0

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

% track alignment

theta = 0;    
R_max = 25000;   % 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;  % P-wave speed in bottom
cs = 800;  % S-wave speed in bottom, cs = 0 for fluid bottom
beta_p = 0; % bottom attenuation for P-wave, dB/wavelength
beta_s = 0; % bottom attenuation for S-wave, dB/wavelength

eta = 1/(40*pi*log10(exp(1)));
cb = cb/(1+1i*beta_p*eta); % complex bottom sound speed
cs = cs/(1+1i*beta_s*eta);

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

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 = 50;      % number of terms in the Bessel series, Eq. (23)
pars.nb_max = 1;   % number of possible reflection from the bottom (total number of image sources is 4*nb_max + 2, see Eq. (25) )
pars.OmThre = 100;   % a threshold value of \Omega in Eq. (23). See comments in the function PenWedgeField( ) for details

pars.isave = 1;     % save the field of each image source (1) or not (0)

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);
xlabel('r (from the source), km');
ylabel('TL, dB re 1 m');
grid on;
axis ij;

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