%% Q1

clear;

t=3600;
D = 1e-5;
M = 1;
x = -2:0.01:2;

C = M/sqrt(4*pi*D*t)*exp(-x.^2/(4*D*t));

figure(1)
plot(x,C,'k-');
xlim([-1 1]);
xlabel('x (m)','fontsize',12);
ylabel('C(x,t)','fontsize',12);
title('1D diffusion, M=1, D=1x10^{-5}, C(x,t=3600)');

%% Q2

clear;

t=3600;
D_x = 1e-5;
D_y = 2e-5;
M = 1;
[x,y] = meshgrid([-2:0.1:2],[-2:0.1:2]);
 
C = M/(4*pi*sqrt(D_x*D_y)*t)*exp(-x.^2/(4*D_x*t)).*exp(-y.^2/(4*D_y*t));
 
figure(2)
surf(x,y,C);
xlim([-2 2]); ylim([-2 2]);
xlabel('x (m)','fontsize',12);
ylabel('y (m)','fontsize',12);
zlabel('C(x,y,t)','fontsize',12);
title('2D diffusion, M=1, D_{x} = 2x10^{-5}, D_{y} = 1x10^{-5}, C(x,y,t=3600)')


%% Q3

clear;

D = 1e-5;
M = 1;
x = -2:0.001:2;
t = [60:60:3600];

% % extra) 2*sigma = sqrt(8Dt)
% charac_d = sqrt(8*D*t);

% a) FWHM
FWHM = nan(length(t),1);
for i=1:length(t)
    C = M/sqrt(4*pi*D*t(i))*exp(-x.^2/(4*D*t(i)));
    FWHM(i) = 2*abs(x(find(C-(M/sqrt(4*pi*D*t(i))/2)>0,1,'first')));
end

% b) C_thresh = 0.3
C_thresh1 = nan(length(t),1);
for i=1:length(t)
    C = M/sqrt(4*pi*D*t(i))*exp(-x.^2/(4*D*t(i)));
    C_thresh1(i) = abs(2*x(find(C-0.3>0,1,'first')));
end

% c) C_thresh = 0.03
C_thresh2 = nan(length(t),1);
for i=1:length(t)
    C = M/sqrt(4*pi*D*t(i))*exp(-x.^2/(4*D*t(i)));
    C_thresh2(i) = abs(2*x(find(C-0.03>0,1,'first')));
end

figure(3); hold on; box on;
% plot(t,charac_d,'k-');
plot(t,FWHM,'r--');
plot(t,C_thresh1,'b:');
plot(t,C_thresh2,'g-.');
xlabel('t (s)','fontsize',12);
ylabel('Cloud width (m)');
% legend('2\sigma','FWHM','C_{thresh} = 0.3','C_{thresh} = 0.03','location','northwest');
legend('FWHM','C_{thresh} = 0.3','C_{thresh} = 0.03','location','northwest');