munk.m

function [val,z,zc,B]=munk(rxz,fld,arg,zc,B)
% [val,z,zc,B]=munk(rxz,fld,arg,zc,B)
%
% Makes the Munk (1979) sound speed profile (10.1121/1.1914492)
% also quoted by MIT Open Course Ware 2-S998 (Spring 2012, Lab 05)
%
% INPUT: 
%
% rxz         Single set of Cartesian position coordinates [x,z], in m, with
%             x the horizontal and z positive downward, OR, 
%  onez       a single depth in m, OR, a triplet
%  [zm dz zM] The minimum depth (positive m) [default: 0]
%             The depth discretization (positive m) [default: 10]
%             The maximum depth (positive m) [default: 4000]
%             and then a linearly spaced depth axis is made, OR
%  manyz      more than three depths that you care to specify
% fld         1 acoustic sound speed [default, only option]
% arg         1 returns the field
%             2 returns the field gradient d(fld)/dx
%             3 returns the field gradient d(fld)/dy
% zc          The depth of the minimum speed (m)
% B           The Munk (1979) scale depth parameter (m)
%
% OUTPUT:
%
% val         The sound speed or its derivatives at the requested depth(s)
% z           The requested depths (m)
% zc          The depth of the minimum speed (m)
% B           The Munk (1979) scale depth parameter (m)
%
% EXAMPLE:
%
% [cz,z,zc,B]=munk([0 10 -gebco(-149,-17)],[],1); plot(cz,z); axis ij; grid on
% xlabel('sound speed (m/s)'); ylabel('depth (m)'); title('Munk profile')
%
% [cz,z,zc,B]=munk([],[],3); plot(cz,z); axis ij; grid on
% xlabel('sound speed gradient (1/s)'); ylabel('depth (m)'); title('Munk profile')
%
% SEE ALSO:
%
% BULLEN, LINMOD
% 
% Last modified by fjsimons-at-alum.mit.edu, 10/26/2022

% Default values for the depth specification
defval('zm',0)
defval('dz',10)
defval('zM',4000)
% If empty makes a linear spacing with defaults
defval('rxz',zm:dz:zM);
% Interpret the input
if length(rxz)==1
  % If its is a single number, it's just a depth
  z=rxz;
elseif length(rxz)==2
  % If it is a single two-set, the second coordinate is depth used
  z=rxz(2);
elseif length(rxz)==3
  % If it is a triplet makes a linear spacing with the inputs
  z=rxz(1):rxz(2):rxz(3);
elseif length(rxz)>3
  % Any and all the depths, including those generated by default
  z=rxz;
end

% Default values
defval('fld',1)
defval('arg',1)

% Default values for the ocean model
defval('zc',1300)
defval('B',1300)

% The sound speed profile
epsilon=0.00737;
% The rescaled depth
zbar=2*(z-zc)/B;

% There is no distinction in the psord argument here
switch arg
  case 1
  % The sound speed
  val=1500*[1+epsilon*(zbar-1+exp(-zbar))];
 case 2
  val=zeros(size(zbar));
 case 3
  % The vertical derivative aka GRADIENT(cz,z)
  val=1500*epsilon*[1-exp(-zbar)]*2/B;
end