close all; clear; format compact; clc
% Matlab script
% Developed by Julio Antonio Lara-Hernandez (julioalh@atmosfera.unam.mx)
% and Jorge Zavala-Hidalgo (jzavala@atmosfera.unam.mx).
% Ocean-atmosphere interaction group, Center of Atmospheric Sciences,
% National Autonomous University of Mexico (UNAM), Mexico.
%
% This script shows that a Turbulent-Difussion velocity (v') can be
% equivalently constructed from random numbers of (1) uniform or (2) normal
% distributions. It also shows that by setting a = b/sqrt(dt), the effect of
% Turbulent-Difussion remains the same regardless the selected time-step (dt):
%
% (1) v' = Va * Ru
% Where: v' = Turbulent-Difussion velocity
%        Va = Advective velocity (e.g. velocity of an ocean current)
%        Ru = Random number of uniform distribution from -a to a
%        a  = b/sqrt(dt)
%        b  = User-defined parameter modulating the degree of Turbulen-Difussion
%        dt = Time-step defined by the user
%
% (2) v' = Rn
% Where: Rn = Random number of normal (Gaussian) distribution with a mean
%             equal to zero and a standard deviation std_Rn = Va*a/sqrt(3)
%
% In the resulting graphs after running the script, the central curve
% indicates the mean, the upper curve is the mean+std, and the lower curve
% is the mean-std. The overlapping of the curves from both approaches
% indicates that they are equivalent. Reproducing the same curves by
% considering different time-steps shows that setting a equal to b/sqrt(dt)
% works to achieve the same effect of Turbulent-Difussion regardless the 
% selected time-step.

%% User-defined parameters:
% Definig the number of particles to simulate
n_particles = 10000;
% Defining for how long the particles will be moving due to Turbulent-Difussion
simulationTime_h = 24; % Defined in hours
% Defining the value of parameter b
b = 4*1350^(1/2);
% Considering different time-steps (dt) to test
timeSteps = [1350/4, 1350/2, 1350, 1350*2]; % Defined in seconds
% Defining the value of the advective velocity (Va)
Va = 0.01; % Defined in m/s

%% Calculation section
% Asigning the advective velocity to each of the simulated particles
Va = zeros(1,n_particles) + Va; % Defined in m/s
% Initialize a cycle counter
counter = 0;
% Starting cycle for testing different time-steps
for dt = timeSteps
    % Update the cycle counter
    counter = counter + 1;
    % Calculating the number of time-steps needed to achieve the simulation time
    n_dt = (simulationTime_h*60*60)/dt;
    % To specify the x-axis in hours, considering that initial time is 0
    axisx_hrs = (0:dt:simulationTime_h*60*60)/60/60;
    % Calculating the value of parameter a
    a = b/dt^(1/2);
    % Calculating the standard deviation of a uniform distribution from -a to a
    std_Ru = a/sqrt(3);
    % Calculating the standard deviation for the random number of normal distribution
    std_Rn = std_Ru * Va(1);
    % Initializing the matrices that will contain the positions as computed from
    % approaches (1) v' = Va * Ru, and (2) v' = Rn 
    positions1 = zeros(n_dt+1,n_particles); % Defined in meters
    positions2 = zeros(n_dt+1,n_particles); % Defined in meters
    % Starting cycle for computing y displacements at each time-step
    for cicle_dt = 2:n_dt+1
        % Calculating the random velocity according to approach (1) v' = Va * Ru
        RandomVelocity1 = (-a+2*a.*rand(1,n_particles)) .* Va; % Defined in m/s
        % Calculating the random velocity according to approach (2) v' = Rn
        RandomVelocity2 = std_Rn.*randn(1,n_particles); % Defined in m/s
        % Total velocity for approach (1) is advective velocity plus random velocity
        TotalVelocity1 = Va + RandomVelocity1; % Defined in m/s
        % Total velocity for approach (2) is advective velocity plus random velocity
        TotalVelocity2 = Va + RandomVelocity2; % Defined in m/s
        % Calculating the new positions of particles considering approach (1)
        positions1(cicle_dt,:) = positions1(cicle_dt-1,:) + TotalVelocity1 .* dt; % Defined in meters
        % Calculating the new positions of particles considering approach (2)
        positions2(cicle_dt,:) = positions2(cicle_dt-1,:) + TotalVelocity2 .* dt; % Defined in meters
    end
    % Calculating the standard deviation of y positions at each time-step from approach (1)
    stde_appr1 = std(positions1,0,2);
    % Calculating the mean of y positions at each time-step from approach (1)
    mean_appr1 = mean(positions1,2);
    % Calculating the standard deviation of y positions at each time-step from approach (2)
    stde_appr2 = std(positions2,0,2);
    % Calculating the mean of y positions at each time-step from approach (2)
    mean_appr2 = mean(positions2,2);
    
    %% Graphing section
    % Plotting the mean+-standard deviations of the positions computed from both approaches
    subplot(2,2,counter)
    hold on
    title(['dt = ',num2str(dt),'s ; Va = ',num2str(Va(1)),'m/s'])
    % Plotting approach (1)
    plotAppr1=plot(axisx_hrs,[mean_appr1+stde_appr1,mean_appr1-stde_appr1,mean_appr1],'r','linewidth',3);
    % Plotting approach (2)
    plotAppr2=plot(axisx_hrs,[mean_appr2+stde_appr2,mean_appr2-stde_appr2,mean_appr2],'b','linewidth',1);
    % Add legends
    legend([plotAppr1(1),plotAppr2(1)],['v'' = Va*Ru ;  a = ',num2str(a)],...
        ['v'' = Rn ;  std = ',num2str(std_Rn)],'Location','northwest')
    % Setting limits in x axis
    xlim([0,simulationTime_h])
    % Add label to x axis
    xlabel('Time (h)')
    % Add label to y axis
    ylabel('Position in y (m)')
end
set(gcf,'position',[488,215,1043,555])