function [H, Ht, Hs] = mld(p,st,varargin);
% function [H, Ht, Hs] = mld(p,st,T,S,meth,chk,bsz,sm5);
% 4 May 2004 - David S. Ciochetto (Dalhousie University - Physical Oceanography)
%
%	A function to calculate the mixed layer depth (H) from CTD profiles of pressure (p) and potential density (st).
% Multiple methods are implimented to estimate the mixed layer depth (mld). The option to use T is retained
% by asking for it to be sent to the function. 'chk' is an optional flag to make a plot to compare the data
% to the MLD estimate. If it is not desired to send T or to modify this to calculate T depths, then set it to
% the integer value of zero or to a surface value. T is used for method 4 and T and S are used for method 5.
% If T is not supplied for method 4, it can use a single T and S value at the surface to calculate the MLD.
% If T and S profiles are supplied to method 5, it will return depths to which these are uniform as well, Ht
% for T and Hs for S.
%
%	The estimate for the MLD (or other depths) will use the entire data vector sent here. If it is desired to
% exclude points at the top or bottom, you must trim them away before calling this routine.
%
% AVAILABLE METHODS:
%		(To choose the method, enter the number below for 'meth')
% 1. max N from a sigma_T profile smoothed by a Moving Average
% 2. CVP (Casult, Venzenia and Petrie, 2003 for the Scotian Shelf)
%	3. Modified CVP method where the bin size matches the MA length of method 1 and interpolate to the criteria.
% 4. KRH method (Kara et al. 2000) - specify delta T as a criteria for delta sigma_T
%	5. TF method (Thomson and Fine 2003) - specify an error norm EN to converge to. (calls their set of m-files)
%
% INPUTS:
%		p			pressure data for the profile [db or m]
%		st		sigma_T for the profile [kg/m^3]
%		T			Temperature profile [deg C]
%		S			Salinity profile 
% 	meth	desired method to estimate the MLD
%		chk		flag to create a plot of the result, 'y' if a check plot is desired.
%					It will add it to the current figure or make a new one. There are options to the plotting.
%					'y1' = just add the mld estimate, a colored line and symbol with a red dot in the center
%					'y2' = add the above plus text in the lower left hand corner with the methd
%					'y3' = add the above plus any fit to the profile
%					'y4' = add the above plus the orginal data profile
%
%					They can be called independently e.g. 'y13', 'y4', 'y1' or 'y1234' which is the same as 'y'
%		bsz		smoothing interval (bin size), if you wish to send the depth to average over, smooth, in units of [m]
%						or the units of the depth variable. Commonly used are 2 and 1. If 0 is sent, it will not filter
%						the data. Use this option for previously averaged data. If it is not choosen, a default of 2 is
%						set as this was found to produce the best results for method 1 and a MA low pass filter.
%
%		sm5		Controls the smoothing of method 5. Options = 'n','yma','yba' where 'n' is for no smoothing,
%					'yma' is to apply a moving average zero-phase filter like in method 1 of length
%
% OUTPUT:
%			H			estimate of the mixed layer depth [db or m based on units of p]
%			Ht		estimate of the isothermal layer depth (ILD) [db or m based on units of p]
%			Hs		estimate of the isohaline layer depth (MLD for the S profile) [db or m based on units of p]
%
%	NOTES:
%		1. Additional details of each method are in the comments in each respective section of the code.
%		2. There are settings that a user may wish to set within the code for each method. If you choose a MA filter
%				to smooth the data, you may want to verify the internal settings for instrument vertical velocity and
%				sampling rate. Defaults are 0.5 m/s at 8 Hz or 16 samples/m.
%		3. CVP is set for use with water 200m or less in depth. It automatically bin averages into 5m bins.
%		4. Method 5 m-files were supplied by K. Thomson. If distribuited, they are not to be included without permission.
%		5. KRH has no options to filter the data in this m-file. Send only the data here that you want to evaluate.
%				So prefilter or pre-smooth the data. Good results are obtained with 8Hz, 0.5 m/s CTD data bin averaged to 1 m.
%
%	REFERENCES:
%
%		Casault, B., A. F. Vezina and Petrie, B. (2003) "Atlas of surface mixed layer characteristics for the
%			Scotian Shelf and the Gulf of Maine," Bedford, N.S., Bedford Institute of Oceanography: 306.
%
%		Kara, A., Rochford, P. and Hurlburt, H. (2000) "An opitmal definition for ocean mixed layer depth,"
%			J. Geophys. Res., Vol. 105, No. C7 p. 16,803-21.
%
%		Thomson, R.E. and I. Fine, 2003 "Estimating Mixed Layer Depth from Oceanic Profile Data,"
%			Journal of Atmospheric and Oceanic Technology, Vol. 20, No. 2, pp. 319-329. 
%
% USES: pbin02.m, denan.m, samppm.m, kara_mld_1.m
%			(req permission for) by_s_m.m, get_xx_norm, s_m_p.m, r2b.m, spl.m, merge_.m
%     Corrected versions of these files may be supplied from David Ciochetto
%     upon confirmation of permission from Dr. Thomson.

% Set the default values in case the arguments are not passed (T, S, meth, chk, bsz).
defaultValues = {0,0,1,'n',2,'yba'};
nonemptyIdx = ~cellfun('isempty',varargin);
defaultValues(nonemptyIdx) = varargin(nonemptyIdx);
[T S meth chk bsz sm5] = deal(defaultValues{:});

%% Setup, things not included in defaults above.
SR = 4;					% Instrument sampling rate [samples/s] - change if raw instrument data is used.
W = 0.5;				% Instrument translation speed [m/s]	 - change if raw instrument data is used.
deltaT = 0.8;		% Delta T for method of Kara et al. 2000
EN = 0.01;			% Error norm to converge the segments to in the Thomson and Fine method
% sm5 = 'n';	% Option to pre-smooth the Thompson method inputs at a size bsz with bin average/MA/none. ['yba'/'yma'/'n']
% 							% Note to smooth, both this must be 'y**' and bsz must be non-zero.
fail = 'n';			% Stop the routine if an estimate fails to work ('y'), otherwise if it fails return NaN ('n')
		% (%) of the screen along the x-axis to put text identifying the estimate type (for chk='y2')
px = 0.05;
		% Check plot common title if the data profile is plotted from this routine
tit = 'MLD estimate(s) for data (blk-.) based on smoothed data (see method color code)';
%%


%% USER SHOULD NOT HAVE TO MODIFY THE CODE BELOW %%

H=[]; Ht=[]; Hs=[];

if nargin < 4
	disp('Not enough inputs to choose the method: using default: method 1')
	meth = 1;
end

if nargin < 2
	error('ERROR: Not enough inputs to the function to continue.')
end

switch meth
case 1			% MLD = depth of max N
	% NOTES: 	% This method will allow for profile data with NaN values.
	
	% Smooth the data with a M point moving average (MA) where M is the number of samples averaged over.
	% Smooth to bsz m, how many samples (M) do I smooth over? It can differ for a model output or experimental data.
	sn = samppm(p);		% Find the number of samples per meter as indicated by the pressure data.
	
	if sn > 1 	% Most likely experimental data since I would make model output 1 m resolution
		%% % Calculate M to average over a depth = bsz of data with a/descent rate W and sample rate SR.
		M = SR/W*bsz;
	else				% Not experimental data, data every 1 m or more
		M = sn*bsz;
	end
	% M has to be greater than or equal to 1
	if M < 1 & bsz~=0
		error('The number of points for the MA is less than 1')
	end
	
	if bsz==0		% Indication to not filter the data.
		y=st;
	else 	% Apply a zero phase moving average (MA) filter of size M to the density data to smooth it (Low pass)
		B=ones(M,1)/M; y=filtfilt(B,1,denan(st));
	end
	% Find the location of the max of the gradient of y (sigma_T) by first difference
	ddy = diff(y); idx = ddy==max(ddy);
	H = p(idx);	% The MLD estimate
	% In case there are more than one point, take the first one.
	if length(H)>1
		H=H(1);
	end
	% Note no interpolation since this method assumes a high spatial sampling of the water column.
	% I may want to check if the value of p at idx is a NaN... but if data for st exists, data for p should too.
	% May need H = denan(p(idx));
	
	% CHECK by plotting if desired in call to program
	if exist('chk') & lower(chk(1))=='y'
		hold on;
		if (strcmpi(chk,'y') & isempty(findstr(chk,'3')>0)) | findstr(chk,'3')>0
			plot(y,denan(p),'-g+')
		end
		if (strcmpi(chk,'y') & isempty(findstr(chk,'4')>0)) | findstr(chk,'4')>0
			plot(st,p,'-k.'); grid on; 
			set(gca,'FontSize',12); xlabel('\sigma_T [kg/m^3]'); ylabel('p [m]');	title(tit);
		end
		xl = get(gca,'XLim'); yl = get(gca,'YLim'); 
		if (strcmpi(chk,'y') & isempty(findstr(chk,'1')>0)) | findstr(chk,'1')>0
			 line(xl,[H H],'Color','g'); plot(y(idx),H,'go',y(idx),H,'r.');
		end
		if (strcmpi(chk,'y') & isempty(findstr(chk,'2')>0)) | findstr(chk,'2')>0
			py = 0.24; 	% (%) along the y-axis to place the text identifying this method
			TxT = ['Method 1:  d\rho/dz_{max} (grn-+) ' num2str(bsz) 'm MA filter'];
			text(min(xl)+px*(max(xl)-min(xl)),min(yl)+py*(max(yl)-min(yl)),TxT,'FontWeight','bold','color','g');
		end
		hold off;
	end
	
	
case 2		% MLD estimated by Casult et al. 2003 method for the Scotian Shelf
	% They use a bin size of 5 m for data between 0-200 m depth (due to the use of old sparse data)
	bsz = 5;	% Set the bin size to 5 m
	[bz, br] = pbin02(p,st,bsz,'n');	% Bin the profile data
	
	% To extend the analysis with their bin structure, have to build  these vectors in chunks.
	%% (CVP) Table 1. Depth structure used for bin averaging of profile data.
	%% Depth range Depth increment Bin range
	%% 0 to 200 m			 5 m ± 2.5 m
	%% 200 to 300 m		10 m ± 5 m
	%% 300 to 500 m 	20 m ± 10 m
	
	% They filtered out profiles with less than three data points
	if length(bz) < 3
		if fail==lower('y')
			error('CVP Method requires more than three data points per profile')
		else
			H=NaN; Ht=NaN; Hs=NaN;
			return
		end
	end
	
	% They defined the surface mixed layer as the water between the surface and the pycnocline. The MLD
	% was determined as the first significant inflection point encountered traversing downward in the profile
	% coinciding with the top of the pycnocline.
	
	% Step 1. Calculate the density gradient at the midpoint between two successive data points
	% They use a first difference in sigma-t divided by a first difference in p [db]
	% Treat the pressures as positive values such that the derivative is (+)
	drdz = diff(br)./diff(abs(bz));
	% Need a new coordinate system for the midpoints to go with drdz
	zp = bz(1:end-1) + diff(bz)/2;
	% CHECK
	% subplot 121, plot(drdz,zp(1:end-1),'.-'); grid on; subplot 122, plot(br,bz,'s',st,p,'.-'); grid on
	% plot(st,p,'c.-',br,bz,'bs',2*drdz+22.5,zp(1:end-1),'r.-'); grid on;
	
	% Step 2. Estimate the MLD - MLD = point where the density gradient was >= 0.01 kg/m^4
	H = zp(drdz >= 0.01);
	if isempty(H)				% If the criteria is not met
		H=zp(end);				% Set it to the middle of the lowest bin depth calculated.
	else
		H = H(1);
	end
	
	% Step 3 - ommitted - They calculate a value of the mean mixed layer T, S and SigmaT by integrating
	%											each respective property from the mixed layer to the surface. (They do not state
	%											it but they must divide by the MLD then to get the correct average.
	
	% CHECK by plotting if desired in call to program
	if exist('chk') & lower(chk(1))=='y'
		hold on;
		if (strcmpi(chk,'y') & isempty(findstr(chk,'3')>0)) | findstr(chk,'3')>0
			plot(br,bz,'-bv')
		end
		if (strcmpi(chk,'y') & isempty(findstr(chk,'4')>0)) | findstr(chk,'4')>0
			plot(st,p,'-k.'); grid on; 
			set(gca,'FontSize',12); xlabel('\sigma_T [kg/m^3]'); ylabel('p [m]');	title(tit);
		end
		xl = get(gca,'XLim'); yl = get(gca,'YLim'); 
		if (strcmpi(chk,'y') & isempty(findstr(chk,'1')>0)) | findstr(chk,'1')>0
			idx=find(zp==H); rH=interp1([0 1],br([idx idx+1]),0.5);
			line(xl,[H H],'Color','b'); plot(rH,H,'bo',rH,H,'r.');
		end
		if (strcmpi(chk,'y') & isempty(findstr(chk,'2')>0)) | findstr(chk,'2')>0
			py = 0.20;
			TxT = ['Method 2: CVP (blu-v) ' num2str(bsz) 'm bin avg, d\rho/dz >= 0.1 kg/m^4'];
			text(min(xl)+px*(max(xl)-min(xl)),min(yl)+py*(max(yl)-min(yl)),TxT,'FontWeight','bold','color','b');
		end
		hold off;
	end
	
	
case 3		% Modified CVP method for higher data density (bin size to match the smoothing applied in meth = 1
					% Also do linear interpolation to match the criteria value better. (see notes for meth 2)
					% This may go horribly wrong as their value of drho/dz may be calibrated for the 5 m bin.
	if bsz ~= 0
		[bz, br] = pbin02(p,st,bsz,'n');	% Bin the profile data to match the size of method 1.
	else	% Do not filter the data if bsz = 0
		bz=p; br=st;
	end
	
	if length(bz) < 3									% Filter out profiles with less than three data points
		if fail==lower('y')
			error('CVP Method requires more than three data points per profile')
		else
			H=NaN; Ht=NaN; Hs=NaN;
			return
		end
	end
	% Step 1. Calculate the density gradient (Define the estimates location is at the midpoint between two successive data points)
	drdz = diff(br)./diff(abs(bz));
	zp = bz(1:end-1) + diff(bz)/2;
	% Step 2. Estimate the MLD - MLD = first point where the density gradient was >= 0.01 kg/m^4
	H = zp(drdz >= 0.01); 
	if isempty(H)				% If the criteria is not met
		H=zp(end);				% Set it to the middle of the lowest bin depth calculated.
		crit = 'n';
	else
		H = H(1);
		crit ='y';
	end

	% Linearly interpolate between the two bounding points to get H such that it matches the criteria exactly.
	idx=find(zp==H);
	if crit == 'y';
		% pbin02 may give NaN values in the surface if there is no data there, it takes all of the data in the profile!!!
		if sum(isnan(drdz([idx-1 idx])))==0	% If the point is not at the top of the profile, e.g. a value other than a NaN above it
			H=interp1(drdz([idx-1 idx]),zp([idx-1 idx]),0.01);	% Ok to interpolate
			rH=interp1(bz([idx-1 idx]),br([idx-1 idx]),H);
		else
			rH=interp1(bz([idx idx+1]),br([idx idx+1]),H);			% Otherwise take the top of the top of the profile and do not interpolate.
		end
	else
		rH=interp1(bz([idx idx+1]),br([idx idx+1]),H);
	end
	
	% CHECK by plotting if desired in call to program
	if exist('chk') & lower(chk(1))=='y'
		hold on;
		if (strcmpi(chk,'y') & isempty(findstr(chk,'3')>0)) | findstr(chk,'3')>0
			plot(br,bz,'-m*')
		end
		if (strcmpi(chk,'y') & isempty(findstr(chk,'4')>0)) | findstr(chk,'4')>0
			plot(st,p,'-k.'); grid on; 
			set(gca,'FontSize',12); xlabel('\sigma_T [kg/m^3]'); ylabel('p [m]');	title(tit);
		end
		xl = get(gca,'XLim'); yl = get(gca,'YLim'); 
		if (strcmpi(chk,'y') & isempty(findstr(chk,'1')>0)) | findstr(chk,'1')>0
			line(xl,[H H],'Color','m'); plot(rH,H,'mo',rH,H,'r.');
		end
		if (strcmpi(chk,'y') & isempty(findstr(chk,'2')>0)) | findstr(chk,'2')>0
			py = 0.15;
			TxT = ['Method 3: CVP^+ (mag-*) ' num2str(bsz) 'm bin avg, d\rho/dz >= 0.01 kg/m^4'];
			text(min(xl)+px*(max(xl)-min(xl)),min(yl)+py*(max(yl)-min(yl)),TxT,'FontWeight','bold','color','m');
		end
		hold off;
	end
	
	
case 4		% Method of Kara 2000
	% They determine the Isothermal Layer Depth (ILD) first from the T profile, then use the density criteria 
	% based on computing sigma_T from T and the equation of state.
	% If T is supplied calculate the ILD and use the T criteria to calculate the density.
	%deltaT = 0.8;		% Delta T criteria suggested in Kara et al. 2000
	minZ = 0;				% Min depth in abs val of data to process (0 - process all of the profile data)
	idx = abs(p) >= minZ;			% This takes all of the NaN data out of the array I send to the subroutine!!!
	if length(T) >= 2
		% Compute the ILD and the value of T at the ILD
		ttt=T(idx); sss=S(idx); rrr=st(idx);
		[Ht,Tt] = kara_mld_1(T(idx),p(idx),deltaT,ttt(1)-deltaT);
		% Compute the MLD and the value of st at the MLD
		% Not stated in the paper, they calculate the delsigmaT with the surface values of T and S
		%disp('MLD')
		deltaST = abs( (sw_dens(sss(1),ttt(1)+deltaT,0)-1000)-rrr(1) );
		[H,Rt] = kara_mld_1(st(idx),p(idx),deltaST,rrr(1)+deltaST);
	elseif T ~= 0 & S ~= 0	% Surface values are given to compute the MLD from density only.
		rrr=st(idx);
		deltaST = abs( (sw_dens(S,T+deltaT,0)-1000)-rrr(1) );
		[H,Rt] = kara_mld_1(st(idx),p(idx),deltaST,rrr(1)+deltaST);
	else
		disp('The method of Kara et al. 2000 is based on a temperature criteria.')
		disp('You must supply a surface T and S value for the method to run.')
		if fail==lower('y')
			error(' ')
		else
			H=NaN; Ht=NaN; Hs=NaN;
			return
		end
	end
	
	% CHECK by plotting if desired in call to program
	if exist('chk') & lower(chk(1))=='y'
		hold on;
		clr = [255 170 1]./255;
		if (strcmpi(chk,'y') & isempty(findstr(chk,'3')>0)) | findstr(chk,'3')>0
			% no fit method applied
		end
		if (strcmpi(chk,'y') & isempty(findstr(chk,'4')>0)) | findstr(chk,'4')>0
			plot(st,p,'-k.'); grid on; 
			set(gca,'FontSize',12); xlabel('\sigma_T [kg/m^3]'); ylabel('p [m]');	title(tit);
		end
		xl = get(gca,'XLim'); yl = get(gca,'YLim');
		if (strcmpi(chk,'y') & isempty(findstr(chk,'1')>0)) | findstr(chk,'1')>0
			 line(xl,[H H],'Color',clr); plot(Rt,H,'s','color',clr); plot(Rt,H,'r.');
		end
		if (strcmpi(chk,'y') & isempty(findstr(chk,'2')>0)) | findstr(chk,'2')>0
			py = 0.10;
			TxT = ['Method 4: KRH (org-s) \DeltaT = ' num2str(deltaT) ' ^oC'];
			text(min(xl)+px*(max(xl)-min(xl)),min(yl)+py*(max(yl)-min(yl)),TxT,'FontWeight','bold','color',clr);
		end
		hold off;
	end
	
	
case 5		% Method of Thomson 2003
	% Call their program or sub functions	(use it to determine a depth for S, T and sigma_T)
	% Use the full length of the data sent to the function.
	% May have trouble if the z-axis is not an increasing function of index...
	% The method does not like NaN values...

	if lower(sm5)=='yma'
		% Smooth the Data first with a MA filter
		sn = samppm(p);		% Find the number of samples per meter as indicated by the pressure data.
		
		if sn > 1 	% Most likely experimental data since I would make model output 1 m resolution
			%% % Calculate M to average over a depth = bsz of data with a/descent rate W and sample rate SR.
			M = SR/W*bsz;
		else				% Not experimental data, data every 1 m or more
			M = sn*bsz;
		end
		% M has to be greater than or equal to 1
		if M < 1 & bsz~=0
			error('The number of points for the MA is less than 1')
		end
		
		if bsz==0		% Indication to not filter the data - alternate to sm5 to match what I do above.
			y=denan(st);
		else 	% Apply a zero phase moving average (MA) filter of size M to the density data to smooth it (Low pass)
			B=ones(M,1)/M; y=filtfilt(B,1,denan(st));
		end
		% Replace the data sent to the function with the filtered data to be sent to the depth finding routine.
		st=y;
		
		if length(T) > 1
			clear y;
			y=filtfilt(B,1,denan(T));
			if bsz ~= 0
				T=y;
			else
				T=denan(T);
			end
		end
		
		if length(S) > 1
			clear y;
			y=filtfilt(B,1,denan(S));
			if bsz ~= 0
				S=y;
			else
				S=denan(S);
			end
		end
		
	elseif lower(sm5)=='yba'		% Bin average rather than moving average

		if bsz ~= 0
			[bz, br] = pbin02(p,st,bsz,'n');	% Bin the profile data to match the size of method 1.
		else	% Do not filter the data if bsz = 0
			bz=denan(p); br=denan(st);
		end
		p=bz; st=br; clear bz br;

		if length(T) > 1
			[bz, bT] = pbin02(p,T,bsz,'n');
			if bsz ~= 0
				T=bT;
			else
				T=denan(T);
			end
			clear bT bz
		end
		
		if length(S) > 1
			[bz, bS] = pbin02(p,S,bsz,'n');
			if bsz ~= 0
				S=bS;
			else
				S=denan(S);
			end
			clear bS bz
		end
	
	end
	
	% It does not like NaN values so deNaN the data.
	if sum(isnan(p))>0
		idx = ~isnan(st);
		p=p(idx);
		st=st(idx);
		if length(T)>1	% Assuming that all the data have the same data length.
			T=T(idx);
		else
			T=denan(T);
		end
		if length(S)>1
			S=S(idx);
		else
			S=denan(S);
		end
	end
	
	N=length(st);  % EN=0.01;	% Send all of the data that is sent here to the routines, so deNaN first.
	
	% Make sure that the depth is positive (Thomson and Fine use positive depths)
	pos = 'y'; idx = p<0;
	if sum(idx)>0.1*length(p)								% If over 10% of the values are (+) then invert it.
		pos = 'n'; p = -p;
	end
	
	if length(T) > 1
		[segsT, zht, Dht, Ht] = by_s_m(N,EN,min(p),max(p),p,T);
	end
	
	if length(S) > 1
		[segsS, zhs, Dhs, Hs] = by_s_m(N,EN,min(p),max(p),p,S);
	end
	
	[segs, zh, Dh, H] = by_s_m(N,EN,min(p),max(p),p,st);
	
	% Convert the depth back to a (+) number if (+) was fed to the routine.
	if pos == 'n'
		p = -p;		H = -H;		zh = -zh;
		if length(T)>1,	zht = -zht;		Ht = -Ht; end
		if length(S)>1, zhs = -zhs; 	Hs = -Hs; end
	end
	
	% CHECK by plotting if desired in call to program
	if exist('chk') & lower(chk(1))=='y'
		hold on;
		clr = [0 200 255]./255;
		if (strcmpi(chk,'y') & isempty(findstr(chk,'3')>0)) | findstr(chk,'3')>0
			plot(Dh, zh,'-d','color',clr);
		end
		if (strcmpi(chk,'y') & isempty(findstr(chk,'4')>0)) | findstr(chk,'4')>0
			plot(st,p,'-k.'); grid on; 
			set(gca,'FontSize',12); xlabel('\sigma_T [kg/m^3]'); ylabel('p [m]');	title(tit);
		end
		xl = get(gca,'XLim'); yl = get(gca,'YLim');
		if (strcmpi(chk,'y') & isempty(findstr(chk,'1')>0)) | findstr(chk,'1')>0
			idx=find(zh==H);
			line(xl,[H H],'Color',clr); plot(Dh(idx),H,'d','color',clr); plot(Dh(idx),H,'r.');
		end
		if (strcmpi(chk,'y') & isempty(findstr(chk,'2')>0)) | findstr(chk,'2')>0
			py = 0.05;
			TxT = ['Method 5: TF (lt. blu-d) Error Norm = ' num2str(EN)];
			text(min(xl)+px*(max(xl)-min(xl)),min(yl)+py*(max(yl)-min(yl)),TxT,'FontWeight','bold','color',clr);
		end
		hold off;
	end
	
otherwise		% If the wrong method is choosen
	disp(['The method choosen to calculate the MLD, ' num2str(meth) ' is not implimented'])
end
