%==================================================================
%  
%  TRACEO: Pekeris block waveguide example.
%  Faro, Mon Jan  6 15:14:15 WET 2014
%  Written by Tordar 
%  
%==================================================================

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

clear all, close all

case_title = '''Pekeris waveguide & arrival calculations''';

freq       =  100;
Rmaxboxkm  = 5.1; Rmaxbox = Rmaxboxkm*1000;
Rmaxhrykm  = 5.0; Rmaxhry = Rmaxhrykm*1000;
Dmax   =  1000;
zs = 500; rs = 0; thetamax = 85; nthetas = 2001;
m = 1001; n =  501;
ranges = linspace(-Rmaxhry,Rmaxhry,m);
zr     = linspace(0,Dmax,n);

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

ray_step = Rmaxbox/1000; 

la = linspace(-thetamax,thetamax,nthetas);

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

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

altimetry = [-Rmaxbox-1 Rmaxbox+1;0 0];

surface_data.type  =   '''V''';
surface_data.ptype =   '''H'''; % Homogeneous
surface_data.units =   '''W'''; % (Attenuation Units) dB/Wavelenght
surface_data.itype =  '''FL'''; % Sea surface interpolation type
surface_data.x     = altimetry; % Surface coordinates
surface_data.properties = [0 0 0 0 0]; % [cp cs rho ap as]

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

c0 = 1500;

ssp_data.cdist = '''c(z,z)'''; % Sound speed profile

ssp_data.cclass = '''ISOV'''; % Isovelocity profile

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

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

object_data.nobjects = 0; % Number of objects

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

bathymetry = [-Rmaxbox-1 2000 2010 2990 3000 Rmaxbox+1;...
               Dmax      Dmax  500  500 Dmax Dmax];

bottom_data.type   = '''E''' ; 
bottom_data.ptype  = '''H''' ; % Homogeneous bottom
bottom_data.units  = '''W''' ; % Bottom attenuation units
bottom_data.itype  = '''2P'''; % Bottom interpolation type 
bottom_data.x      = bathymetry; % Bottom coordinates 
bottom_data.properties = [1700.0 0 1.7 0.7 0]; % [cp cs rho ap as]

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

output_data.array_shape = '''VRY''';
output_data.ctype       = '''ADP''';
output_data.r           = Rmaxhry  ;
output_data.z           = Dmax*0.75;

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

load aad

rayname = 'aad00000';

if neigrays >  0 

for i = 1:neigrays

    si = num2str(i); lsi = length(si); 
    eval( [ 'ray = ' rayname(1:8-lsi) si ';' ] ); 
    tau(i) = ray(3);
    amp(i) = ray(4) + 1i*ray(5);
    amp(i) = amp(i)*exp( 1i*ray(6) );
    
end
figure(1)
stem( tau, abs( amp ) )
box on, grid on
xlabel('Travel time (s)')
ylabel('Amplitude')
title('TRACEO - Pekeris waveguide, delay calculations')

else
disp('no arrivals found...')
end 

disp('done.')
