
%% Problem 1
% Laser light, which can be modeled as a plane wave, is incident on a
% finite width double slit.  The wavelength of light is 633 nm.  The slit
% width is 20 microns wide and the center of the slits are 100 microns apart.
% Use Huygen's principle to find the intensity of the light
% on a screen 10 cm away. Compare to the expected intensity disribution
% from theory. Normalize the intensity so that the maximum is equal to 1.
% Hint: Only calculate the intensity at the screen, not
% everywhere in between.

%Use microns as units
Lambda=.633;
a=20;
d=100;
L=1e5;

%Setup an array of sources at opening of slits
Nsources=50;
YSourceVec1=linspace(d/2-a/2,d/2+a/2,Nsources)';
YSourceVec2=-linspace(d/2-a/2,d/2+a/2,Nsources)';
YSourceVec=cat(1,YSourceVec1,YSourceVec2);

%Screen coordinates
YScreen=(-1e4:.1:1e4)';

%Function to calculate amplitude
wave_amp=@(YSource,Phase) 1./sqrt(L.^2+(YScreen-YSource).^2) .* ...
    exp(i*(2*pi/Lambda*sqrt(L.^2+(YScreen-YSource).^2)+Phase));

WaveInt=0
for Phase=(.1:.1:2*pi); %Time Average
    WaveAmpAll=0;
    for ii=1:length(YSourceVec) %Sum over sources
        WaveAmpAll=WaveAmpAll+wave_amp(YSourceVec(ii),Phase);
    end
    WaveInt=WaveInt+abs(WaveAmpAll).^2;
end
%normalize
WaveInt=WaveInt/max(WaveInt);

%Plot
figure
plot(YScreen,WaveInt)

%Theory
SinTheta=YScreen./sqrt(L.^2+YScreen.^2);
doubleSlitDif=cos(pi/Lambda*SinTheta*d).^2.*sinc(SinTheta/Lambda*a).^2;
hold on
plot(YScreen,doubleSlitDif,'--','linewidth',2)

xlabel('Screen Position (cm)')
ylabel('Normlized Intensity')
legend('Huygens'' Principle','Fraunhofer Diffration Theory')


%% Problem 2
%
% The file 'HW12Data.mat' contains two variables.  'Tvec' indicates the measurement time.
% 'Volts' is the measured data.  Use a Fourier transform to find the 
% four frequency components that make up the signal.

clear all
clc
close all

load HW12Data.mat

F_Volts=fft(Volts);

%This is the frequency
Fvec=1/max(Tvec)*(0:length(Tvec)-1);

%Take everything below Nyquist Frequency (Positive frequences)
DT=Tvec(2)-Tvec(1);
mask=Fvec<1/DT/2;

%Plot 
FigH=figure
plot(Fvec(mask),abs(F_Volts(mask)))
xlabel('Frequency (Hz)')
ylabel('FFT(Volts)')

%Find where the peaks are located. Many ways to do this. Including just reading
% off the plot.  Here is one:

% Simple spike detector 
FVabs=abs(F_Volts(mask));
FVdiv=(FVabs(2:end-1))>FVabs(1:end-2)&(FVabs(2:end-1)>FVabs(3:end));
FVdiv=cat(2,0,FVdiv,0); % pad ends that were cut off above
figure
plot(FVdiv)

Peaks=FVdiv>0;
[IDs]=find(Peaks)

%Here are the frequencies in Hz.  
FoundFreq=Fvec(IDs)

figure(FigH)
for ii=1:length(IDs)
    s=sprintf('%g Hz',FoundFreq(ii));
    X=abs(Fvec(IDs(ii)));
    Y=abs(F_Volts(IDs(ii)));
    text(X,Y,s)   
end

























