%==================================================================
%  
%  TRACEO: Pekeris flat waveguide example.
%  Gambelas, Mon Jan  6 12:38:36 WET 2014
%  Written by Tordar 
%  
%==================================================================

system('rm -rf rco.mat');

clear all, close all 

case_title = '''Pekeris waveguide & ray coordinates''';

%==================================================================
%  
%  Define source data:
%  
%==================================================================

freq   =   70;  % Correspont to propeller's rotation frequency 
% Setting range "box"
Rmaxkm =   1; 
Rmax = Rmaxkm*1000;
Dmax   =  100; 

ray_step = Rmax/1000; 

zs = 3;  % considering the source close to the sea surface
rs = 0; 
thetamax = 30; 
nthetas = 50; 
thetas = linspace(-thetamax,0,nthetas); % Zero is used to avoid reflections at the surface

%source_data.ds       = 0; % If set zero, model uses the ranges of the rbox to calculate a preliminary step.
% ?????????????????????? 1 de 2 »»» Não está funcionando com ds = 0
source_data.ds       = ray_step; % If set zero, model uses the ranges of the rbox to calculate a preliminary step.
source_data.position = [rs zs];
source_data.rbox     = [rs-1 Rmax+1];
source_data.f        = freq;
source_data.thetas   = thetas;

%==================================================================
%  
%  Define altimetry data:
%  
%==================================================================

altimetry = [rs-2 Rmax+2;0 0];

surface_data.type  =   '''V'''; % Surface type: V --» Vacuum over surface
surface_data.ptype =   '''H'''; % Surface Properties: H --» Homogeneous surface
surface_data.units =   '''W'''; % Attenuation Units: W --» dB/Wavelenght
surface_data.itype =  '''FL'''; % Sea surface interpolation type: FLat surface
surface_data.x     = altimetry; % Surface coordinates
surface_data.properties = [0 0 0 0 0]; % Dummy values
% Properties and coordinates of the surface [cp cs rho ap as] (for homogeneous surface): 
% [compressional velocity, shear velocity, density, compressional attenuation, shear attenuation]

%==================================================================
%  
%  Define sound speed data:
%  
%==================================================================

c0 = 1500;

ssp_data.cdist = '''c(z,z)'''; % Type of sound speed distribution: Sound speed profile c = c(z)

ssp_data.cclass = '''ISOV'''; % Class of sound speed: ISOVelocity profile c(z) = c0 = constant

ssp_data.z   = [0 Dmax]';
ssp_data.r   = [];
ssp_data.c   = [c0 c0]';

%==================================================================
%  
%  Define object data:
%  
%==================================================================

object_data.nobjects = 1; % Number of objects inside the waveguide

cobjp = 2000; % cp --» compressional velocity of the object
cobjs = 200; % cs --» shear velocity of the object
rhoo = 5;  % density of the object
alphao = 0;

npo = 4; 

R0 = 5; % number of coordinates of the object
factor = 1; 
ro = 10; % horizontal distance from boat's propeller to left edge of object
zo = 6; % depth of the left edge of object

robj  = factor*[ro ro+0.5 ro+2.5 ro+3];
zup   = [zo zo+0.5 zo+0.5 zo];
zdn   = [zo zo-0.5 zo-0.5 zo];

xobj = [robj(:)';zdn(:)';zup(:)']; 

object_data.itype     = '''2P'''; % Objects interpolation type: 2P --» Piecewise linear interpolation
object_data.npobjects = npo; % Number of points in each object
object_data.type      = '''R'''; % Object Type: R --» rigid surface (**same as Altimetry Block**)
object_data.units     = '''W'''; % Attenuation Units: W --» dB/Wavelenght (**same as Altimetry Block**)
object_data.properties = [cobjp cobjs rhoo alphao alphao];
% Properties and coordinates of the object [ocp ocs orho oap oas]: 
% [compressional velocity, shear velocity, density, compressional attenuation, shear attenuation]
object_data.x(1,1:3,:) = xobj;

%==================================================================
%  
%  Define bathymetry data:
%  
%==================================================================

bathymetry = [rs-2 Rmax+2;Dmax Dmax];

bottom_data.type   = '''E''' ; % Bottom type: E --» Elastic surface
bottom_data.ptype  = '''H''' ; % Bottom Properties: H --» Homogeneous bottom
bottom_data.units  = '''W''' ; % Bottom Attenuation Units: W --» dB/Wavelenght
bottom_data.itype  = '''FL'''; % Bottom interpolation type: FLat surface
bottom_data.x      = bathymetry; % Bottom coordinates 
bottom_data.properties = [1566.7 0 1.775 0.14 0]; % [cp cs rho ap as]
% Properties and coordinates of the bottom [cp cs rho ap as] (for homogeneous bottom): 
% [compressional velocity, shear velocity, density, compressional attenuation, shear attenuation]
% Compressional velocity of 1566.7 obtained from Walid & Sylvain, 2013, p.18 (1520 + 1580 + 1600)/3
% Bottom density of 1.775 g/cm3 obtained from Walid & Sylvain, 2013, p.18 (1.75 + 1.08)/2
% Compressional Attenuation of 0.14 dB/Wavelenght obtained from Walid & Sylvain, 2013, p.18 (0.13 + 0.15)/2

%==================================================================
%  
%  Define output data:
%  
%==================================================================

nra = 1; % Number of hydrophones along range
rarray = linspace(0,Rmax,nra);
% Rmax = 1km is used here just for simplicity. The real range from surce to array is about 3 Km (Figure 1 -Gingras, 1994)
nza = 48;  % Number of hydrophones along depth
zarray = linspace(18.7,112.7,nza);

output_data.ctype       = '''RCO'''; % Output: RCO --» Ray COordinates
output_data.array_shape = '''VRY'''; % Array Type: VRY --» Verical arRaY
output_data.r           = rarray;
output_data.z           = zarray;
output_data.miss        = 0.5; % used as a threshold to find eigenrays and to calculate arrivals
% ????????????????????????? 2 de 2 »»» O array é plotado na saída?

%==================================================================
%  
%  Call the function:
%  
%==================================================================

disp('Writing TRACEO waveguide input file...')

wtraceoinfil('pekeris.in',case_title,source_data,surface_data,ssp_data,object_data,bottom_data,output_data);

disp('Calling TRACEO...')

system('runtraceo pekeris');
%system('cp pekeris.in WAVFIL');
%system('/home/vbarroso/bin/traceo.exe');
%system('mv LOGFIL pekeris.log');
%system('rm WAVFIL');

figure(1), hold on
plot(rs,zs,'ko',rs,zs,'m*','MarkerSize',16) 
disp('Reading the output data...')

load rco

rayname = 'ray00000';

figure(1), hold on
plot(rs,zs,'ko',rs,zs,'m*','MarkerSize',16) 
for i = 1:nthetas

    si = num2str(i); lsi = length(si); 
    eval( [ 'ray = ' rayname(1:8-lsi) si ';' ] ); 
    r = ray(1,:);
    z = ray(2,:);
    plot(r,z)
    clear r z
    
end

grid on, box on
fill(xobj(1,:),xobj(2,:),'k')
fill(xobj(1,:),xobj(3,:),'k')
axis([0 Rmax 0 Dmax])
view(0,-90)
hold off 
xlabel('Range (m)')
ylabel('Depth (m)')
title('TRACEO - Pekeris waveguide, ray coordinates')

disp('done.')
