% HW8
%% Problem #1
% In this problem I take 10^-30 Kg = 1 (which is approximately the mass of electron),
% 10^-10 m = 1 (which is one Angstrom) and 1.5e-11 s = 1.

Hydrogen = HW8_wave_decay;
Hydrogen.N = 1000; %number of samples
Hydrogen.DeltaX = 200; %length of the spacial interval in Angstrom
Hydrogen.Sigma = 0.5; %Bohr's radius in Angstrom
Hydrogen.M = 1; %mass of electron which is approximately 10^-30
Hydrogen.T = 0.000052; %runnung time in 1.5e-11 s

Hydrogen = Hydrogen.gauss_wave;
Hydrogen.dec_wave;

% As you can see in the plot, it takes 8*10^-16 seconds for the wave-packet to
% spread from 0.5 A to 5 A~10*0.5 A. Supposing that it is linear then we can
% easily say that 1m = (2*10^10)*0.5 A and the total time that it takes to
% spread to 1m would be

figure
plot(Hydrogen.Tvec*1.5e-11,Hydrogen.SigmaX)
title('spread over time is linear')
ylabel('\sigma'); xlabel('time')

Time_to_spread_to_1m = ((Hydrogen.T)*1.5e-16/5)*2e10 % seconds

%% problem #2

Events = double(rand(1,1001)<0.05); %signal generated from a constant rate process
figure
plot(Events)
title('Signal generated from a constant rate process')

Tt = (0:1:1000);
Ins_Res = (1-exp(-Tt)).*exp(-0.11*Tt); %Photmultiplier tubes responce
figure
plot(Tt,Ins_Res)
title('Instrument Responce')
xlabel('t');ylabel('Responce')

FEvents = (fft(Events)); %Fourier transform
FIns_Res = (fft(Ins_Res)); %Fourier transform
Convu = ifft(FEvents.*FIns_Res); %Convolution in frequency domain reduces to a simple product
figure
plot(Tt,Convu)
title('Convolution')
xlabel('time')


%% Problem #3
% In theory Intensity is proportional to cos(pi*d*sin(theta)/lambda)^2*sinc(pi*b*sin(theta)/lambda)^2
% where d is the distance between the slits and b is the slits' width and
% sin(theta) = x/D. D is the distance between the slights and the screen

Slits = zeros(201,201); % An opeque obstacle
Slits(:,60:70) = 1; %carving the frist slit
Slits(:,130:140) = 1; % carving the second slit
F = fftshift(fft2(Slits)); % the diffraction pattern is simply the Fourier transform of the obstacle
II = abs(F).^2; %Intensity
D = 100000; %the distance of slits from the screen with the diffraction pattern (micron)
Lambda=0.633; %wavelength in micron

figure
imshow(II,[]);title('Intensity')
figure
imshow(Slits,[]);title('Slits')
figure
imshow(log(II),[]);title('log(I)')

I = II(101,:);
DX = 1/201;
X = (-100:100)*DX*Lambda*D; %distance from the middle of the screen 
figure; plot(X,I/max(I))
xlabel('Position on the screen(\mum)');ylabel('Intensity')
% The diffraction pattern of the two slits is given by the following function
% see 'https://en.wikipedia.org/wiki/Double-slit_experiment'
Theory_I = @(x,y)(cos(x).*(sin(y)./y)).^2;
b = 10; % the distance between the centers of the two slits
d = 70; % the slits' width
y = pi*b*X./(Lambda*sqrt(D^2+X.^2));
x = pi*d*X./(Lambda*sqrt(D^2+X.^2));
ITheory = Theory_I(x,y);
ITheory(101) = cos(x(101))^2;
ITheory_norm = ITheory/max(ITheory);
hold
plot(X,ITheory_norm,'r--')
legend('Numeric','Theory')

%% Problem #4

close all
clear 
clc

A = load('HW8_NoisyData.mat'); %loading the data
T = 2; %the whole time interval is two seconds
N = length(A.V); %Number of elements in data vector
Leng = length(A.V);
figure
subplot(1,2,1)
time = linspace(0,2,Leng); %time vector
plot(time,A.V)
title('noisy signal in time domain')
xlabel('Time (s)'); ylabel('Signal (Volts)')

FV = fft(A.V); % the Fourier transform of the noisy data
F = (0:N-1)*(1/T); %Frequencies

subplot(1,2,2)
plot(F,abs(FV))
title('Fourier Transform of V')
xlabel('Frequency (Hz)');ylabel('Amplitude')

% We can see directly from the plot that the peak, corresponing to the
% signal, is at 4 Hz.  

% But Let's get this using code...
F_ID=find(max(abs(FV))==abs(FV));
FSignal = F(F_ID);
% We just need the positive frequency:
fprintf('Signal Frequency: %g Hz\n',FSignal(1));

% Now let's find the amplitude.  
% From the definition of fft and remembering integral of sin^2 = L/2 (over
% full periods)

%Fourier coeffiecent
AF = abs(FV(F_ID(1)));
Amp = AF/(N/2);
fprintf('Signal Amplitude: %g Volts\n',Amp);


