%==================================================================
% 
% TRACEO
% Faro, Sun Jul  3 18:26:46 WEST 2011
% Written by Tordar 
% 
%==================================================================

clear all, close all 

disp('Deep water test:') 

case_title = '''Munk profile & transmission loss''';

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

freq   =    50;
Rmaxkm =   101; Rmax = Rmaxkm*1000;
Dmax   =  5000; 

ray_step = Rmax/1000; 

zs = 1000; rs = 0;
np2 = 51; thetamax = 14; la = linspace(-thetamax,thetamax,np2);

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

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

altimetry(1,:) = [rs-2 Rmax+2];
altimetry(2,:) = [0         0];

surface_data.type  =   '''V'''; % 
surface_data.ptype =   '''H'''; % Homogeneous
surface_data.units =   '''W'''; % (Attenuation Units) Wavelenght
surface_data.itype =  '''FL''';
surface_data.x     = altimetry; % Surface coordinates
surface_data.properties = [0 0 0 0 0.0]; % Dummy parameters

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

c1 = 1500; z1 = 1300;

depths = linspace(0,Dmax,1001);

c = munk( depths, z1, c1 );
 
ssp_data.cdist = '''c(z,z)'''; % Sound speed profile
ssp_data.cclass = '''TABL''';
ssp_data.z    = depths(:);
ssp_data.r    = [];
ssp_data.c    = c(:);

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

object_data.nobjects = 0; % No objects

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

bathymetry(1,:) = [rs-2 Rmax+2];
bathymetry(2,:) = [Dmax   Dmax];

bottom_data.type   = '''E''' ;
bottom_data.ptype  = '''H''' ; % Homogeneous bottom
bottom_data.units  =  '''W'''; % (Attenuation Units) Wavelenght
bottom_data.itype  = '''FL'''; % Bottom interpolation type 
bottom_data.x      = bathymetry;     % Bottom coordinates 
bottom_data.properties = [1550.0 600.0 2.0 0.1 0.0]; % Bottom properties (speed, speed, density, absorption coefficient)

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

nra = 501; nza = 501;
rarray = linspace(0,100*1000,nra);
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('munk.in',case_title,source_data,surface_data,ssp_data,object_data,bottom_data,output_data);

disp('Calling TRACEO...')

!runtraceo munk

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

counter = 0;

load cpr

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

tej = flipud( jet( 256 ) );

counter = counter + 1; figure(counter), hold on
plot(rs,zs,'ko',rs,zs,'m*','MarkerSize',16) 
pcolor(rarray,zarray,tl), shading interp, colormap( tej ), caxis([50 120]), colorbar
axis([0 Rmax 0 Dmax])
hold off
view(0,-90)
xlabel('Range (m)')
ylabel('Depth (m)')
title('TRACEO - Munk profile')

disp('done.')
