%Homework #7
%% Problem #1
clear all
clc

DesiredDigits=50;

%See HW7_NDigits.m 
A = HW7_NDigits(1,DesiredDigits); 
B = HW7_NDigits(3,DesiredDigits); 
C = A/B
D=A+C
E=D*B
F=A-B

%problem here:
X=HW7_NDigits(4,DesiredDigits);
Y=HW7_NDigits(6,DesiredDigits);
Z=HW7_NDigits(2,DesiredDigits);

Z+Y
W = Z/Y
W/X  % This is incorrect

V = HW7_NDigits(24,DesiredDigits);

%This won't work because the 0.16666... can only be represented by a
%double.
K = HW7_NDigits(0.16666666666666666666666666666666666666666666666666666,100)
F = HW7_NDigits(0.15555555555555555555555555555555555555555555555555555,100)

Z/V
Z/(X*Y) 

%%  Calculate the series expansion of e^1
%% finding e
clear 
clc
DesiredDigits=100

A=HW7_NDigits(1,DesiredDigits*2);
E=HW7_NDigits(1,DesiredDigits*2);
N=1;
while A.OoM>-DesiredDigits
   A=A/HW7_NDigits(N,DesiredDigits*2)
   E=E+A;
   E
   N=N+1;
end
E.N=DesiredDigits;
E.NArray=E.NArray(1:DesiredDigits);
E

% Compare to symbolic:
digits(DesiredDigits)
%Do the series expansion of e^1
syms k
MyE2 = symsum(1/factorial(k), k, 0, DesiredDigits);
%Make into vector
Y=MyE2;
MyEDouble=zeros(DesiredDigits,1);
for ii=1:DesiredDigits
    MyEDouble(ii)=double(floor(Y));
    Y=(Y-floor(Y))*10;
end
Sym = HW7_NDigits(1,DesiredDigits);
Sym.NArray=MyEDouble

%% Problem #2

%Dimensions selected to keep array size small

Lambda=.633;            %microns
D=Lambda*50;            %microns 
GridSize=D*10;          %microns 
GridSpacing=Lambda/2;   %microns
Z=D*0.8;               %ScreenDist 

%Create aperture function
[X_A Y_A]=meshgrid((-GridSize/2:GridSpacing:GridSize/2),(-GridSize/2:GridSpacing:GridSize/2));
R_A=sqrt(X_A.^2+Y_A.^2);
A = R_A>D/2;     %Un-obstructed area
figure;imshow(A);title('Blocking Disk')

% Calculate the intensity across a line centered on disk shadow
X=(-D:GridSpacing:D);
Y=0;

%Calculate the numerical integral
NX=length(X);
clear Huygens Huygens_Fresnel
h=waitbar(0,'calculating integral...');
for nnx=1:NX
    waitbar(nnx/NX)
    %Position at screen:
    XS=X(nnx);
    YS=Y;
    ZS=Z;
    %Distance from each point on aperture to screen:
    R=sqrt((X_A-XS).^2+(Y_A-YS).^2+ZS.^2);
    WaveContibutions=A.*exp(i*2*pi/Lambda*R)./R;
    Huygens(nnx)=abs(sum(WaveContibutions(:))).^2;
    
    CosAlpha=SZ./R;
    %Using Huygens_Fresnel principle
    WaveContibutions_F=WaveContibutions*(1+CosAlpha)*i/Lambda;
    %Using Huygens principle
    Huygens_Fresnel(nnx)=abs(sum(WaveContibutions_F(:))).^2;
end
close(h)
%Normalize for easy comparison
Huygens_Norm=Huygens/max(Huygens);
Huygens_Fresnel_Norm=Huygens_Fresnel/max(Huygens_Fresnel);

figure
plot(X,Huygens_Norm)
hold on
plot(X,Huygens_Fresnel_Norm)
xlabel('X (microns)')
ylabel('Intensity')
legend('Huygens','Huygens-Fresnel')


%% Problem #3
clc
% an arbitrary vector
A = [1.2,2,-1,0,5.3,-0.05];
%comparing the results of our function with fft
DFT_A = HW7_dft(A)
FFT_A = fft(A)




