Skip to content
Snippets Groups Projects
Commit 3ad29f0a authored by bgergely's avatar bgergely
Browse files

Add new file

parent 6bc8b447
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Afshar_Distinguishibility.m 0 → 100644
+ 432
− 0
View file @ 3ad29f0a
% Afshar number 638nm, 40 micron slits d, center to center D 250 micron spacing,
% lens diameter 3cm, lens distance 0.55m, 100cm focal length, wire gauge 127 microns, detector
% seperation .6mm, first order should be ~x*lambda/d~1.4mm from center
% this takes the better part of a day to complete on a i7-3770 with 16gb of ram
clear %clears workspace
clc %clears command window
close all %deletes all open figures
format longE %Long scientific notation with 15 digits after the decimal point for double values
tic %times the program
fprintf('Program Start %s\n', datestr(now,'HH:MM:SS PM')) %print starting time of code
%Define physical constants and parameters
lambda = 638E-9; %wavelength of light
lenssize = 30E-3; %physical size of lens
screensize = 20E-4; %physical size of screen where diffraction pattern will be observed
slitnumber=2; %NEED TO CHANGE FIT VALUE SEPARATELY number of slits
slitwidth=50E-6; %NEED TO CHANGE FIT VALUE SEPARATELYindividual slit width
slitseparation=250E-6; %NEED TO CHANGE FIT VALUE SEPARATELY center to center slit separation
Distance1=0.55; %distance from slits to lens
Distance2=0.55; %Distance1*4; %distance from lens to screen
%Define computational parameters
pointdensity1=1000000; %source iterating points/meter
pointdensity2=1000000; %lens iterating points/meter
pointdensity3=1000000; %screen iterating points/meter
sourcesize=slitnumber*slitwidth+((slitseparation-slitwidth)*(slitnumber-1)); %total width of source including the center between slits
sourcenumber=round(pointdensity1*sourcesize+1); %calc number of points over entire source including center section
lensnumber=round(pointdensity2*lenssize+1); %calc number of iterating points at the lens
screennumber=round(pointdensity3*screensize+1); %calc number of iterating points at the screen
slit1end=round(sourcenumber*slitwidth/sourcesize); %calc number position of first slit
slit2start=round((sourcenumber*(slitwidth+(slitseparation-slitwidth))/sourcesize)); %calc position of 2nd slit
testline1=screensize/2+(lambda*(Distance1+Distance2)/slitseparation); %location of young's double slit
testline2=screensize/2-(lambda*(Distance1+Distance2)/slitseparation); %location of young's double slit
slitline1=screensize/2-slitwidth/2-slitseparation/2; %slit location at screen
slitline2=screensize/2-slitseparation/2+slitwidth/2; %slit location at screen
slitline3=screensize/2+slitseparation/2-slitwidth/2; %slit location at screen
slitline4=screensize/2+slitseparation/2+slitwidth/2; %slit location at screen
slit1center=slitline1+slitwidth/2;
slit2center=slitline3+slitwidth/2;
%define slit intensity bias
adatanumber=20; %number of a values in 0.05 increments starting from 0
meshindexmax=3; %number of mesh size steps
%pre-defining matrices to contain zeroes, speeds up computation.
adata=zeros(adatanumber+1,meshindexmax);
newGrating=zeros(1,lensnumber);
yNormlens1=zeros(1,adatanumber+1,meshindexmax);
NormArea=zeros(1,adatanumber+1,meshindexmax);
probout2=zeros(lensnumber,adatanumber+1,meshindexmax);
probout4=zeros(screennumber,adatanumber+1,meshindexmax,2);
phasefinalout=zeros(screennumber,adatanumber+1,meshindexmax,2);
totaladataX=zeros(adatanumber+1,meshindexmax);
% totaladataV=zeros(screennumber,adatanumber,meshindexmax);
totaladataD=zeros(adatanumber+1,meshindexmax);
% totaladataT=zeros(screennumber,adatanumber,meshindexmax);
for meshindex=1:meshindexmax
meshwidth=(meshindex-1)*127E-6; %width of wires
meshgap=lambda*Distance1/(slitseparation)-meshwidth; %spacing between wires placed at nodes of diffraction pattern
meshwnumber=round(pointdensity2*meshwidth); %number of computaional points that the wires represent
meshgnumber=round(pointdensity2*meshgap); %number of computational points that the gaps represent
meshnumber=fix(lensnumber/(meshwnumber+meshgnumber));
if 1 == rem(meshnumber,2)
else
meshnumber=meshnumber-1;
end
leftoverlensnumber=lensnumber-((meshwnumber+meshgnumber)*meshnumber);
meshcellwidth=(meshwnumber+meshgnumber);
totalmeshnumber=(meshcellwidth)*meshnumber;
epsillon=1; %epsillon is used to determine the lens strength factor
for a=1:adatanumber+1 %number of iterations for the slit intensity bias
for phaser=1:1
r = -epsillon*(100/55);%(epsillon*.1)*1/(2*Distance2); %Lens strength factor
%Setup intial double slit wave vector for slit 1
phiin_A = zeros(round(sourcenumber),1);
for n=1:slit1end+1
phiin_A(n)=2*((a-1)*0.05)*sqrt(1/(slit1end*slitnumber));
end
phiin_B = zeros(round(sourcenumber),1);
%Setup intial double slit wave vector for slit 1
for n=slit2start:round(sourcenumber)
phiin_B(n)=2*(1-((a-1)*0.05))*sqrt(1/(slit1end*slitnumber))*exp((phaser-1)*(pi)*(sqrt(-1)));
end
%Calcuate propagation kernel2 from source to lens
Kernel2=zeros(sourcenumber,lensnumber);
for j=1:lensnumber
for n=1:sourcenumber
Kernel2(n,j)=exp(2*pi*(sqrt(-1))*(1/lambda)*sqrt(((sourcesize/2-(sourcesize/(sourcenumber-1))*(n-1))-(lenssize/2-(lenssize/(lensnumber-1))*(j-1)))^2+Distance1^2));
end
end
%GRATING
vec=zeros(1,meshcellwidth);
positions=(meshwnumber+1:meshcellwidth);
vec(positions)=1;
Grating=repmat(vec, 1, meshnumber);
Grating(totalmeshnumber+meshwnumber) = 0;
startpoint=round((lensnumber-numel(Grating))/2);
newGrating=[ones(1, max(0, round((lensnumber-numel(Grating))/2))), Grating];
leftover=ones(1, lensnumber-numel(newGrating));
newGrating=[newGrating leftover];
%Calcuate propagation kernel3 from lens to screen
Kernel3=zeros(lensnumber,screennumber);
for j=1:lensnumber
for k=1:screennumber
Kernel3(j,k)=exp(2*pi*sqrt(-1)/lambda*(sqrt(((lenssize/2-lenssize/(lensnumber-1)*(j-1))-(screensize/2-screensize/(screennumber-1)*(k-1)))^2+Distance2^2)+r*(lenssize/2-lenssize/(lensnumber-1)*(j-1))^2));
end
end
%calc phi for slit(n) to lens(j) with wire grid
phiout2_A=zeros(lensnumber,1);
for j=1:lensnumber
phioutER2 = zeros(lensnumber,1);
for n=1:sourcenumber
phioutER2(n)=phiin_A(n)*Kernel2(n,j);
end
phiout2_A(j)=newGrating(j)*sum(phioutER2); %Grating
end
%calc phi for slit(n) to lens(j) with wire grid
phiout2_B=zeros(lensnumber,1);
for j=1:lensnumber
phioutER2 = zeros(lensnumber,1);
for n=1:sourcenumber
phioutER2(n)=phiin_B(n)*Kernel2(n,j);
end
phiout2_B(j)=newGrating(j)*sum(phioutER2); %Grating
end
%calc phi for lens(j) to screen (k) ONLY
phiout3=zeros(screennumber,2);
for k=1:screennumber
phioutER3 = zeros(screennumber,1);
for j=1:lensnumber
phioutER3(j)=Kernel3(j,k);
end
phiout3(k)=sum(phioutER3);
end
%calc phi slit to screen TOTAL
phiout4_A=zeros(screennumber,2);
for k=1:screennumber
phioutER4 = zeros(screennumber,1);
for j=1:lensnumber
phioutER4(j)=(phiout2_A(j)+phiout2_B(j))*Kernel3(j,k);
end
phiout4_A(k)=sum(phioutER4);
end
%probability slit to lens, amplitude^2
for j=1:lensnumber
probout2(j,a+1,meshindex)=(phiout2_A(j)+phiout2_B(j))*conj((phiout2_A(j)+phiout2_B(j)));
probout2(j,1,meshindex)=j-1;
end
%Normalization slit to lens
xsl0N = lenssize/(lensnumber-1)*probout2(:,1,meshindex);
ysl0N = probout2(:,a+1,meshindex);
yNormlens1(1,a+1,meshindex)=trapz(xsl0N,ysl0N);
NormArea(1,a+1,meshindex,phaser)=yNormlens1(1,a+1,meshindex)/yNormlens1(1,a+1,1);
%probability TOTAL at screen
for k=1:screennumber
probout4(k,1,meshindex,phaser)=k-1;
probout4(k,a+1,meshindex,phaser)=(phiout4_A(k))*conj((phiout4_A(k)));
phasefinalout(k,a+1,meshindex,phaser)=(phiout4_A(k));
end
%Normalization lens to screen
xT0N=screensize/(screennumber-1)*probout4(:,1,1,1);
yT0N=probout4(:,a+1,meshindex,phaser);
yNorm2(1,a+1,meshindex,phaser)=trapz(xT0N,yT0N);
end
adata(a,1)=(a-1)*0.05;
% refine this section
%generate detector function
slitnumber1=round(slitline1*(screennumber-1)/screensize)-5;
slitnumber2=round(slitline2*(screennumber-1)/screensize)+5;
slitnumber3=round(slitline3*(screennumber-1)/screensize)-5;
slitnumber4=round(slitline4*(screennumber-1)/screensize)+5;
Detector = zeros(round(screennumber),1);
for k=1:slitnumber1
Detector(k)=0;
end
for k=slitnumber1:slitnumber2
Detector(k)=1;
end
for k=slitnumber2:slitnumber3
Detector(k)=0;
end
for k=slitnumber3:slitnumber4
Detector(k)=1;
end
for k=slitnumber4:screennumber
Detector(k)=0;
end
for k=1:round(screennumber/2)
xSD1(k,1)=screensize/(screennumber-1)*(k-1);
ySAD1temp(k,1)=NormArea(1,a+1,meshindex,phaser)*(1/yNorm2(1,a+1,meshindex,phaser))*phiout4_A(k)*Detector(k); %
ySBD1temp(k,1)=NormArea(1,a+1,meshindex,phaser)*(1/yNorm2(1,a+1,meshindex,phaser))*phiout4_A(k)*Detector(k); %
end
for k=round(screennumber/2):screennumber
xSD2(k,1)=screensize/(screennumber-1)*(k-1);
ySAD2temp(k,1)=NormArea(1,a+1,meshindex,phaser)*(1/yNorm2(1,a+1,meshindex,phaser))*phiout4_A(k)*Detector(k);
ySBD2temp(k,1)=NormArea(1,a+1,meshindex,phaser)*(1/yNorm2(1,a+1,meshindex,phaser))*phiout4_A(k)*Detector(k);
end
ySAD1=trapz(xSD1,(ySAD1temp.*conj(ySAD1temp)));
ySBD1=trapz(xSD1,(ySBD1temp.*conj(ySBD1temp)));
ySAD2=trapz(xSD2,(ySAD2temp.*conj(ySAD2temp)));
ySBD2=trapz(xSD2,(ySBD2temp.*conj(ySBD2temp)));
adata(a,4)=(((ySAD1+ySBD1)-(ySAD2+ySBD2))/((ySAD1+ySBD1)+(ySAD2+ySBD2)))^2; %D^2 D = 0.5*abs{|<S1|D1>|^2-|<S1|D2>|^2}+0.5*abs{|<S2|D1>|^2-|<S2|D2>|^2}
end
% Total, add V^2 and D^2
totaladataX(:,meshindex)=adata(:,1);
totaladataD(:,meshindex)=adata(:,4);
%totaladataT(:,meshindex)=adata(:,5);
end
%FIGURE 0 MESH
figure('Name','0 MICRON BARS');
xT=screensize/(screennumber-1)*probout4(:,1,1,1);
yT2=NormArea(1,2,1,1)*(1/yNorm2(1,2,1,1))*probout4(:,2,1,1); % probout4(k,a+1,meshindex,phaser) yNorm2(1,a+1,meshindex,phaser) NormArea(1,a+1,meshindex,phaser)
yT12=NormArea(1,12,1,1)*(1/yNorm2(1,12,1,1))*probout4(:,12,1,1);
yT22=NormArea(1,22,1,1)*(1/yNorm2(1,22,1,1))*probout4(:,22,1,1);
y2T2=NormArea(1,2,1,1)*(1/yNorm2(1,2,1,1))*probout4(:,2,2,1); % probout4(k,a+1,meshindex,phaser) yNorm2(1,a+1,meshindex,phaser) NormArea(1,a+1,meshindex,phaser)
y2T12=NormArea(1,12,1,1)*(1/yNorm2(1,12,1,1))*probout4(:,12,2,1);
y2T22=NormArea(1,22,1,1)*(1/yNorm2(1,22,1,1))*probout4(:,22,2,1);
y3T2=NormArea(1,2,1,1)*(1/yNorm2(1,2,1,1))*probout4(:,2,3,1); % probout4(k,a+1,meshindex,phaser) yNorm2(1,a+1,meshindex,phaser) NormArea(1,a+1,meshindex,phaser)
y3T12=NormArea(1,12,1,1)*(1/yNorm2(1,12,1,1))*probout4(:,12,3,1);
y3T22=NormArea(1,22,1,1)*(1/yNorm2(1,22,1,1))*probout4(:,22,3,1);
sp1=subplot(3,3,1);
plot(xT,yT2,'-k','LineWidth',1.5);
xline(slitline1,'--r');
xline(slitline2,'--r');
xline(slitline3,'--r');
xline(slitline4,'--r');
axis([0 0.002 0 inf]);
xlabel('Screen (meters)')
sp2=subplot(3,3,4);
plot(xT,yT12,'-k','LineWidth',1.5);
xline(slitline1,'--r');
xline(slitline2,'--r');
xline(slitline3,'--r');
xline(slitline4,'--r');
axis([0 0.002 0 inf]);
xlabel('Screen (meters)')
sp3=subplot(3,3,7);
plot(xT,yT22,'-k','LineWidth',1.5);
xline(slitline1,'--r');
xline(slitline2,'--r');
xline(slitline3,'--r');
xline(slitline4,'--r');
axis([0 0.002 0 inf]);
sp4=subplot(3,3,2);
plot(xT,y2T2,'-k','LineWidth',1.5);
xline(slitline1,'--r');
xline(slitline2,'--r');
xline(slitline3,'--r');
xline(slitline4,'--r');
axis([0 0.002 0 inf]);
sp5=subplot(3,3,5);
plot(xT,y2T12,'-k','LineWidth',1.5);
xline(slitline1,'--r');
xline(slitline2,'--r');
xline(slitline3,'--r');
xline(slitline4,'--r');
axis([0 0.002 0 inf]);
sp6=subplot(3,3,8);
plot(xT,y2T22,'-k','LineWidth',1.5);
xline(slitline1,'--r');
xline(slitline2,'--r');
xline(slitline3,'--r');
xline(slitline4,'--r');
axis([0 0.002 0 inf]);
xlabel('Screen (meters)')
sp7=subplot(3,3,3);
plot(xT,y3T2,'-k','LineWidth',1.5);
xline(slitline1,'--r');
xline(slitline2,'--r');
xline(slitline3,'--r');
xline(slitline4,'--r');
axis([0 0.002 0 inf]);
xlabel('Screen (meters)')
sp8=subplot(3,3,6);
plot(xT,y3T12,'-k','LineWidth',1.5);
xline(slitline1,'--r');
xline(slitline2,'--r');
xline(slitline3,'--r');
xline(slitline4,'--r');
axis([0 0.002 0 inf]);
xlabel('Screen (meters)')
sp9=subplot(3,3,9);
plot(xT,y3T22,'-k','LineWidth',1.5);
xline(slitline1,'--r');
xline(slitline2,'--r');
xline(slitline3,'--r');
xline(slitline4,'--r');
axis([0 0.002 0 inf]);
xlabel('Screen (meters)')
sgtitle('0 MICRON BARS')
%The code below is commented out, after code completion highlight individual section below to plot what is needed
%{
%FIGURE GRATING
figure('Name','GRATING');
xsl = lenssize/(lensnumber-1)*probout2(:,1);
scatter(xsl,newGrating)
ylim([0 0.01]);
%run this after running both the Distinguishility and Visibility code
%FIGURE 4 a vs D^2 and V^2
figure('Name','Visibility');
xadX1=totaladataX(:,1);
yadV1=totaladataV(:,1);
yadD1=totaladataD(:,1);
yadT1=totaladataV(:,1)+totaladataD(:,1);
yadV2=totaladataV(:,2);
yadD2=totaladataD(:,2);
yadT2=totaladataV(:,2)+totaladataD(:,2);
yadV3=totaladataV(:,3);
yadD3=totaladataD(:,3);
yadT3=totaladataV(:,3)+totaladataD(:,3);
yadV4=totaladataV(:,4);
yadD4=totaladataD(:,4);
yadT4=totaladataV(:,4)+totaladataD(:,4);
yadV5=totaladataV(:,5);
yadD5=totaladataD(:,5);
yadT5=totaladataV(:,5)+totaladataD(:,5);
yadV6=totaladataV(:,6);
yadD6=totaladataD(:,6);
yadT6=totaladataV(:,6)+totaladataD(:,6);
hold on;
plot(xadX1, yadV1,'*-r')
plot(xadX1, yadD1,'*-b')
plot(xadX1, yadT1,'*-k')
plot(xadX1, yadV2,'*-r')
plot(xadX1, yadD2,'*-b')
plot(xadX1, yadT2,'*-k')
plot(xadX1, yadV3,'*-r')
plot(xadX1, yadD3,'*-b')
plot(xadX1, yadT3,'*-k')
plot(xadX1, yadV4,'*-r')
plot(xadX1, yadD4,'*-b')
plot(xadX1, yadT4,'*-k')
plot(xadX1, yadV5,'*-r')
plot(xadX1, yadD5,'*-b')
plot(xadX1, yadT5,'*-k')
plot(xadX1, yadV6,'*-r')
plot(xadX1, yadD6,'*-b')
plot(xadX1, yadT6,'*-k')
hold off;
xlabel('a-value')
title('Distinguishibility vs Visibility')
legend('V^2','D^2','Total')
xlim([0 1]);
ylim([0 1.05]);
%FIGURE 5 D^2 VS A FOR VARIOUS SLIT SIZE
figure('Name','Distinguishibility');
xadX1=totaladataX(:,1);
yadD1=totaladataD(:,1);
yadD2=totaladataD(:,2);
yadD3=totaladataD(:,3);
hold on;
plot(xadX1, yadD1,'*-')
plot(xadX1, yadD2,'*-')
plot(xadX1, yadD3,'*-')
hold off;
xlabel('a-value')
ylabel('D^2')
title('Distinguishibility')
legend('0 micron','127 micron','508 micron')f
xlim([0 1]);
ylim([0 1.05]);
%FIGURE 6 V^2 VS A FOR VARIOUS SLIT SIZE
figure('Name','Visibility');
xadX1=totaladataX(:,1);
yadV1=totaladataV(:,1);
yadV2=totaladataV(:,2);
yadV3=totaladataV(:,3);
yadV4=totaladataV(:,4);
yadV5=totaladataV(:,5);
yadV6=totaladataV(:,6);
hold on;
plot(xadX1, yadV1,'*-')
plot(xadX1, yadV2,'*-')
plot(xadX1, yadV3,'*-')
plot(xadX1, yadV4,'*-')
plot(xadX1, yadV5,'*-')
plot(xadX1, yadV6,'*-')
hold off;
xlabel('a-value')
ylabel('V^2')
title('Visibility')
legend('0 micron','100 micron','200 micron','300 micron','400 micron','500 micron')
xlim([0 1]);
ylim([0 1.05]);
%}
fprintf('Program End %s\n', datestr(now,'HH:MM:SS PM')) %print finishing time of code
toc %print total time to run
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment