%HW_4
%% Problem # 1
 clear 
 close all
% mass of nitrogen molecules
M = 2*2.32*10^-26; %kg

%Boltzmann constant
K = 1.38*10^-23; %J/K

% Temperature at STP
T = 273; %K

%standard deviation
Sig = sqrt(K*T/M);

% Number of Particles
N=10e5

% drawing random velocities from normal dist. with sigma = Sig
Vel = Sig*randn(3,N);

%the speed (magnitude of the velocity) of the particles
Speed = sqrt(sum(Vel.^2));

%histogram of the speed 
NHistBins=100;
figure
hist(Speed,NHistBins)
xlabel('Velocity')
ylabel('Count')
hold on

%finding the area of histogrm
[HeightV,BinV] = hist(Speed,NHistBins);
BinVWidth = BinV(2)-BinV(1); %BinWidth
AreaV = sum(BinVWidth*HeightV);

%the given velocity distribution is normlized. to make the histogram and
%the normalized dist. the same they we can multiply by N
Dv=BinVWidth;
SpeedDist=@(v)(M/(2*pi*K*T))^(3/2)*4*pi.*v.^2.*exp(-M*v.^2/(2*K*T));
P = N*SpeedDist(BinV)*Dv;
fig = plot(BinV,P);
fig.Color = [1,0.2,0];
fig.LineWidth = 1.5;
title('Velocity (m/s)')

%the energy of the particles
Energy = 1/2*M*Speed.^2;

%histogram of energy
figure
hist(Energy,NHistBins)
xlabel('Energy')
ylabel('Counts')
hold on

%finding the area of the energy histogram
[HeightE,BinE] = hist(Energy,NHistBins);
BinEWidth = BinE(7)-BinE(6);
AreaE = sum(BinEWidth*HeightE);

%plotting the given energy dist.

DE=BinEWidth;
EDist=@(E)2*sqrt(E/pi).*(1/(K*T))^(3/2).*exp(-E/(K*T))
P = N*EDist(BinE)*DE;
plot(BinE,P,'r','LineWidth',1.5)
title('Energy (J)')

%% problem # 2

% CDF(X) = (1/sqrt(2*pi))*integral(exp(-t^2/2),[-inf  X])
% erf(X) = (2/sqrt(pi))*integral(exp(-t^2),[0   X])
% CDF(X) = Y = 0.5*(1+erf(X/sqrt(2)))
% X = sqrt(2)*erfinv(2*Y-1), where erfinv() is the inverse of the erf()

% Y is the integral of the normalized distribution and has values in the interval [0  1]
Y = rand(1,10000);
X = sqrt(2)*erfinv(2*Y-1);

NHistBins=100;
hist(X,20)

[Counts,BinCenter] = hist(X,NHistBins);
BinWidth=BinCenter(2)-BinCenter(1);
NormCounts=Counts/sum(Counts);
figure
bar(BinCenter,NormCounts,'barwidth',1)
xlabel('Value')
ylabel('Counts')

hold on
plot(BinCenter,normpdf(BinCenter)*BinWidth,'r')
legend('Our Distribution','Standard Normal')

fprintf('Mean: %g \n',mean(X));
fprintf('Standard Deviation: %g \n',std(X));


%% problem # 3

N=100000;
NHistBins=100;
Xrand = sum(rand(12,N))-6;
figure
hist(Xrand,20)
%We are summing distributions and according to central limit theorem this
%results in a result that gets closer to the normal distribution. 

%The mean is obviously zero

%The variance of the uniform distribution is:
syms x
VarUnif=int((x-1/2)^2,x,0,1)
StdUnif=sqrt(VarUnif)

%The summed distibution tends toward normal with standard deviation
%sqrt(n)*sqrt(sigma^2) where sigma^2 is variance of the original
%distribution

%In this case:
sqrt(12)*sqrt(VarUnif)

%Compare to the standard normal:
[Counts,BinCenter] = hist(Xrand,NHistBins);
BinWidth=BinCenter(2)-BinCenter(1);
NormCounts=Counts/sum(Counts);
figure
bar(BinCenter,NormCounts,'barwidth',1)
xlabel('Value')
ylabel('Counts')

hold on
plot(BinCenter,normpdf(BinCenter)*BinWidth,'r')
legend('Our Distribution','Standard Normal')

fprintf('Mean: %g \n',mean(Xrand));
fprintf('Standard Deviation: %g \n',std(Xrand));
%% Problem #4
% The number of events recorded over a time period will be poisson
% distributed.  From the PMF, we can get the probability of getting k, with
% an expected (average) count of lambda

P4 = poisspdf(4,4)

%% problem 5

RSpace = 2*(rand(5000,3)>.5)-1;     %producing random integer numbers in the interval [-1  1]
RanWalk = cumsum(RSpace);     %cumsum sums up the elements of an array in the columns
XWalk = RanWalk(:,1);                %X-position of the walker
YWalk = RanWalk(:,2);                %Y-position of the walker
ZWalk = RanWalk(:,3);                %Z-position of the walker
figure
plot3(XWalk,YWalk,ZWalk)
xlabel('X')
ylabel('Y')
zlabel('Z')
title('Random Walk')


