% File: examples.m
% Date: 22-Oct-07
% Author: Scott Sewell

% This matlab script is a container for various example plots
% designed to inspire ideas for physical modeling and analysis

% SECTION 1: simpoisson.m
% This  script demonstrates how to generate and 
% compare sample populations from the Parent PDF's from
% which they are drawn.

mu=1.67;
xmax=mu+8*sqrt(mu);
numsamples=1000;
x=[0:xmax];
sample=poissrnd(mu,1,numsamples);
x1=[0:xmax];
parent=poisspdf(x1,mu)*numsamples;
y=hist(sample,x);
subplot(2,2,1);
bar(x,y,.8,'g');
hold on;
plot(x1,parent,'ko:')
legend('Sample Population','Parent Distribution');
title('Sample and Parent Distributions - POISSON - \mu=1.67');
xlabel('Value')
ylabel('Number of Samples (of 1000 trials)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% singleslit(b,l,f) plots the single slit diffraction pattern for 
% a slit of width 'b' in mm
% a wavelength of 'l' in nm
% a slit-to-screen distance of 'f' in m.
% Produces a plot of intensity versus screen position in mm.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b1=.1; 	
b2=.05;
l=632; 	
f=0.5;  	

theta=[-.02:.0001:.02];
beta1=pi*(b1/1000)/(l/1e9).*sin(theta);
beta2=pi*(b2/1000)/(l/1e9).*sin(theta);
y1=(sin(beta1)./beta1).^2;
y2=(sin(beta2)./beta2).^2;
subplot(2,2,2);
plot(theta*f*1000,y1);
hold on;
plot(theta*f*1000,y2,'--');
legend('100\mum slit','50 \mum slit');
title('Modelling Diffraction from a Single Slit');
xlabel('Screen Position [mm]')
ylabel('Intensity');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% From Hecht, Optics, Edition 4, Equation 10.24
% doubleslit(b,a,l,f) plots the two slit diffraction pattern for 
% a slit of width 'b' in mm
% a slit separation of 'a' in mm
% a wavelength of 'l' in nm
% a slit-to-screen distance of 'f' in m.
% Produces a plot of intensity versus screen position in mm.
% Try doubleslit(0.1,0.4,632,.5)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=0.4;
theta=[-.02:.0001:.02];
beta=pi*(b1/1000)/(l/1e9).*sin(theta);
alpha=pi*(a/1000)/(l/1e9).*sin(theta);
y=4.*(sin(beta).^2./beta.^2).*(cos(alpha).^2);
subplot(2,2,3);
plot(theta*f*1000,y);
title('Modelling Two-Slit Diffraction');
xlabel('Screen Position [mm]')
ylabel('Intensity');

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Planck Radiation Law
% See Hecht 'Optics' for Details
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h = 6.626e-34;   	% Planck's Constant = 4.135 x 10^-15 eV s
c = 3e8;  		% speed of light (MKS)
T= 6000; 		% kelvin
k= 1.38066e-23; 	% Boltzmann constant in J/K
l=0:20e-9:2000e-9;
p=2*3.14*h*c*c./(l.^5);	% 
b6000=p./(exp(h*c./(l*k*T)-1));
b5000=p./(exp(h*c./(l*k*5000)-1));
b4000=p./(exp(h*c./(l*k*4000)-1));
subplot(2,2,4);	 
plot(l,b6000,'.');
title('Planck Radiation Law');
xlabel('Wavelength [m]')
ylabel('Spectral Irradiance [W m^{-2} sr^{-1} nmm^{-1}]');
hold on;
plot(l,b5000,'--');
plot(l,b4000,'+');
legend('6000K','5000K','4000K');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plots a maxwell-boltzmann distribution of neutron velocities
% See Kittel and Kroemer 'Thermal Physics' for details
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m= 29.0;   		% amu air = 29
M= 0.029;  		% kg/mol air = 0.029
T_neutrons= 273; 		% kelvin
k= 1.38066E-23; 		% Boltzmann constant in J/K
R= 8.3145;		% universal gas constant = kN(A) in J/mol K
vprob_273=sqrt(2*R*T_neutrons/M);	% m/s
v=[0:10:5*vprob_273];	% create a molecular speed variable of decent range
Nv_273=4*3.14*(M/(2*3.14*R*T_neutrons))^1.5*v.^2.*exp(-M*v.^2/(2*R*T_neutrons));
Nv_323=4*3.14*(M/(2*3.14*R*323))^1.5*v.^2.*exp(-M*v.^2/(2*R*323));	

figure(2);
subplot(2,2,1); 
plot(v,Nv_273,'.');
title('Maxwell-Boltzmann Distribution');
xlabel('Velocity [m s^{-1}]');
ylabel('N(v)dv');
hold on;
plot(v,Nv_323,'*');
legend('T=273K','T=323K');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulating an exponential decay curve
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=[0:20:2000];
y=100*exp(-.005*x);
noisyy=fix(y+sqrt(y).*randn(1,length(y)));
sig=sqrt(y);
subplot(2,2,2);
errorbar(x,y,sig,'.');
title('Simulated Exponential Decay with Poisson Uncertainties');
xlabel('Bin');
ylabel('Counts');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% File:   sampleerrorbars.m
% Date:   02-Nov-04
% Author: I. Chuang <ichuang@mit.edu>
%
% Demonstrate matlab plotting with vertical and horizontal error bars.
% Uses herrorbar() from www.mathworks.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
npts = 32;
xdat = linspace(0,2*pi,npts);
zdat = cos(xdat);
noise = (rand(1,npts)-0.5)/2;
ydat = zdat+noise;
subplot(2,2,3);
errorbar(xdat,ydat,0.2*ones(1,npts),'ko');
hold on;
herrorbar(xdat,ydat,0.1*ones(1,npts),'ko');
plot(xdat,zdat,xdat,ydat,'go');
grid on;
title('Sample matlab file demonstrating vertical and horozintal error bars');
