diff --git a/Afshar_Distinguishibility.m b/Afshar_Distinguishibility.m
new file mode 100644
index 0000000000000000000000000000000000000000..76ddd0d8525819db3f0fab27759652afa853d6a1
--- /dev/null
+++ b/Afshar_Distinguishibility.m
@@ -0,0 +1,432 @@
+% 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