%==================================================================
%  
%  TRACEO: Pekeris flat waveguide example with 1 Object
%  Gambelas, Tue Jan  7 13:47:57 WET 2014
%  Written by Tordar 
%  
%==================================================================

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

clear all, close all 

case_title = '''Pekeris waveguide & coherent pressure''';

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

freq   = 70; % Correspont to propeller's rotation frequency 
% Setting range "box"
Rmaxkm = 0.1; 
Rmax = Rmaxkm*1000;
Dmax   =  100; % Depth of 100 m obtained from Walid & Sylvain, 2013, p.15

ray_step = Rmax/1000; 

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

source_data.ds       = ray_step;
source_data.position = [rs zs];
source_data.rbox     = [rs-1 Rmax+1];
source_data.f        = freq;
source_data.thetas   = la;

%==================================================================
%  
%  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]: 
% [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 = 50; 
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*R0*cos( linspace(pi,0,npo) ) + ro;
%zup   =         R0*sin( linspace(pi,0,npo) ) + Dmax/2;
%zdn   =        -R0*sin( linspace(pi,0,npo) ) + Dmax/2;

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];
zdn   = [zo zo-0.5 zo-0.5 zo];
zup   = [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 = [1700 0 1.7 0.7 0]; % [cp cs rho ap as]

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

nra = 201; rarray = linspace(0,Rmax,nra);
nza = 201; zarray = linspace(0,Dmax,nza);

output_data.ctype       = '''CPR'''; 
output_data.array_shape = '''RRY''';
output_data.r           = rarray;
output_data.z           = zarray;
output_data.miss        = 0.5;

%==================================================================
%  
%  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');

disp('Reading the output data...')

load cpr

tl = 20*log10( abs( pressure ) );

figure(1), hold on
plot(rs,zs,'ko',rs,zs,'m*','MarkerSize',16) 
pcolor(rarray,zarray,tl), shading interp, colorbar, caxis([-60 -10])
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, acoustic pressure')

disp('done.')
