% Create a Z-Profile overturn in a linearly stratified water column
%
% 2005 Jan 10
% David S. Ciochetto
%
% USES: paper.m, 
%
% See Also: runlength.m, 
%

clear all; close all; paper

SR = 4; 						% Sampling Rate [Hz]
SS = 0.5;						% Sensor Speed [m/s]
dr = [-5 -40];			% Depth range of profile min max [m]
drho = [1023 1026];	% Density to match depth range [kg/m^3]
g = 9.81;						% Gravitational acceleration [m/s^2]
rho0 = 1026;				% Reference Density [kg/m^3]
Lz = -20;						% Location of center of overturn [m]
L = 3;							% Size of overturn [m]

% Generate the independent variable
z = flipud([dr(2):SS./SR:dr(1)]');

% Density gradient
drdz = diff(drho)./diff(dr);

% Buoyancy Frequency N = sqrt(-g/rho0*drho/dz) [1/s]
N = sqrt(-g./rho0.*drdz);
% Buoyancy Period [s]
T = 2*pi/N;

% Generate the linear density signal
r = drdz*(z-dr(1))+drho(1);


% Need some way to keep the noise realistic such as 10% of the signal or so? Do this with the amplitude.
% Possibly a percentage of the difference (not the reference value since that could be large). What is
% realistic to expect, smooth or small scale fluctuations.

% Add random noise in the signal to simulate small random natural fluctuations
a1 = 0.1;	% Random noise amplitude of the background, multiplies the random generator
t1 = 'n';	% Type of noise to add 'g' = Gaussian, 'u' = uniform and 'n' = no noise
switch t1
case 'g'		% Gaussian noise with mean zero, var = 1 and std = 1. Factor of 4 to be similar to rand
	r = r + a1/4*(randn(size(r))-0.5);
case 'u'		% Uniform noise that goes from 0 to 1 or -0.5 to 0.5
	r = r + a1*(rand(size(r))-0.5);
case 'n'		% Choosing this makes r the re-ordered linear profile for gravitational stability.
	r = r;
end

% Check
plot(r,z,'-k')
set(gca,'FontSize',12);
% NOTE: This is not the same as the re-ordered profile if noise is added to this step.


% Overturn the fluid
rho = r;											% rho simulates the measured signal
idx = z>=Lz-L/2 & z<=Lz+L/2;
ot = flipud(rho(idx));
rho(idx) = ot;

% Add Random noise to the signal (see test_runlength.m) to simulate random noise in the sampling methods.
a2 = 0.1;		% Amplitude of the random noise in the sampling.
t2 = 'u';	% Type of noise to add 'g' = Gaussian, 'u' = uniform and 'n' = no noise
switch t2
case 'g'
	rho = rho + a2/4*(randn(size(rho))-0.5);
case 'u'
	rho = rho + a2*(rand(size(rho))-0.5);
case 'n'
	rho = rho;
end

% Check
hold on; plot(rho,z,'r-'); hold off;
ylabel('z [m]'); xlabel('\rho [kg m^{-3}]');
title(['Z-overturn in a linearly stratified fluid with N = ' tsigfig(N,4) ' [s^{-1}]'])



% Perform Run-Length test (as in bottom of test_runlength.m)
% Use the sorted profile to do the differnece since there is more noise introduced to the original r.
figure
if t1 == 'n'		% I can use the method for testing run lengths since r is the reordered profile.
	res = runlength(rho-r); % Calculate the runlength of the signal minus the trend (re-ordered profile)
	% Plot a histogram of the resulting PDF
	N=[1:max(res)]; H=hist(res,N); y = H/sum(H);
	bar(N,y); set(gca,'FontSize',12);
	title('Runlength analysis of fluctuating density signal for linear stratification with noise plus overturn')
	xlabel('n'); ylabel('Probability of run of length n, P(n)  [0-1]')
	hold on; plot(N,2.^(-N),'y.-'); hold off;	% Add the theoretical distribution for noise in yellow
	pause
	% Make a semilog plot of this PDF including the threshold value.
	semilogy(N,y,'kx-',N,2.^(-N),'y.-',N,2*2.^(-N),'r:'); set(gca,'FontSize',12);
	title('Runlength analysis of fluctuating density signal for linear stratification with noise plus overturn.')
	ylabel('Probability of run of length n, P(n)  [0-1]')
	xlabel('n (Solid Line with dots is P(n) = 2^{-n}, dotted line is GK96 threshold value.)')
else
	% I need to gravitationally sort the profile and get the difference from that.
	
end

