function [c,cz,czz] = munk( depths , z1 , c1 ) ;  

%MUNK: Calculates Munk's deep water profile  
%
% Usage:   [c,cz,czz] = munk( depths , z1 , c1 ) 
%
%           where z1 is the channel axis and c1 is the sound speed at z1 
%           (default values: z1 = 1200 , c1 = 1480); "cz" and "czz" are 
%           the corresponding sound speed first and second derivatives.    
%       
%           Example:  
%           depths = [ 0:10:4500 ]  ; z1 = 1400 ; c1 = 1450 ;
%           C0     = munk( depths ) ; 
%           C1     = munk( depths , z1 ) ;  
%           C2     = munk( depths , z1 , c1 ) ;     
%           figure(1), plot( C0 , depths , C1 , depths , '-.' , C2 , depths , '--' ) 
%           grid on, view( 0, -90 )  
%           xlabel( 'Sound Velocity (in m/s)' ), ylabel( 'Depth (in m)' )          
 
%***************************************************************************************
% Second version: 11/09/2006
% 
% Contact: orodrig@ualg.pt
% 
% Any suggestions to improve the performance of this 
% code will be greatly appreciated. 
% 
%***************************************************************************************

c = []; cz = c; czz = c;

B  = 1.3e3 ; B2 = B*B;

epsilon = 7.37e-3 ;  

if nargin == 1 , z1 = 1200 ; c1 = 1480 ; end   

if nargin == 2 , c1 = 1480 ; end   

eta = 2*( depths - z1 )/B ; 

c   = 1 + epsilon*( eta + exp( -eta ) - 1 ); 
c   = c1*c ; 
cz  = ( 2*c1*epsilon/B  ) * ( 1 - exp( -eta ) );
czz = ( 4*c1*epsilon/B2 ) * exp( -eta ); 
