diff --git a/Afshar_Visibility.m b/Afshar_Visibility.m
new file mode 100644
index 0000000000000000000000000000000000000000..7792a64846143fc01f93c4f7f5a79dcb741e822e
--- /dev/null
+++ b/Afshar_Visibility.m
@@ -0,0 +1,956 @@
+% 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 = 30E-3; %physical size of screen where diffraction pattern will be observed
+slitnumber=2; %NEED TO CHANGE FIT VALUE SEPARATELY number of slits
+slitwidth=40E-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 need at least 100000
+pointdensity2=1000000; %lens iterating points/meter need at least 100000
+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
+
+%pre-defining matrices to contain zeroes, speeds up computation.
+meshindexmax=3;
+yNormlens1=zeros(1,adatanumber,meshindexmax);
+probout2=zeros(lensnumber,adatanumber+1,meshindexmax);
+probout4=zeros(screennumber,adatanumber,meshindexmax,21);
+totaladataX=zeros(adatanumber+1,meshindexmax,21);
+totaladataV1=zeros(adatanumber+1,meshindexmax,21);
+
+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;
+
+%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];
+
+adata=zeros(adatanumber,6, 21);
+fitcurve=zeros(adatanumber,1);
+epsillon=0; %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:2
+r = 0; %-epsillon*(2);%(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
+
+%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=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(k)=sum(phioutER4);
+end
+
+%probability slit to lens, amplitude^2
+
+for j=1:lensnumber
+   probout2(j,1,meshindex,phaser)=j-1;
+   probout2(j,a+1,meshindex,phaser)=(phiout2_A(j)+phiout2_B(j))*conj((phiout2_A(j)+phiout2_B(j)));
+end
+
+%Normalization slit to lens
+xsl0N = lenssize/(lensnumber-1)*probout2(:,1,meshindex);
+ysl0N = probout2(:,a+1,meshindex);
+yNormlens1(1,a+1,meshindex,phaser)=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(k))*conj((phiout4(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);
+
+adata(a,1)=(a-1)*0.05;
+%[maxY index]=max(probout4a(:,a+1,meshindex,1));
+index=round(screennumber/2);
+maxY=probout4(index,a+1,meshindex,1);
+minY=probout4(index,a+1,meshindex,2);
+adata(a,4)=((maxY-minY)/(maxY+minY))^2;
+
+end
+
+end
+
+% Total, add V^2 and D^2
+totaladataX(:,meshindex,:)=adata(:,1,:);
+totaladataV1(:,meshindex,:)=adata(:,4,:);
+
+end
+
+figure('Name','bars,phase');
+xT=lenssize/(lensnumber-1)*probout2(:,1,1,1);
+yT2=NormArea(1,3,1,1)*(1/yNormlens1(1,3,1,1))*probout2(:,3,1,1); % (j,a+1,meshindex,phaser)
+yT3=NormArea(1,3,2,1)*(1/yNormlens1(1,3,2,1))*probout2(:,3,2,1);
+yT4=NormArea(1,3,3,1)*(1/yNormlens1(1,3,3,1))*probout2(:,3,3,1); 
+yT5=NormArea(1,3,1,2)*(1/yNormlens1(1,3,1,2))*probout2(:,3,1,2); 
+yT6=NormArea(1,3,2,2)*(1/yNormlens1(1,3,2,2))*probout2(:,3,2,2);
+yT7=NormArea(1,3,3,2)*(1/yNormlens1(1,3,3,2))*probout2(:,3,3,2); 
+sp1=subplot(3,1,1);
+hold on;
+plot(xT,yT2,'k','LineWidth',1.5);
+plot(xT,yT5,'r','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([.35*screensize 0.65*screensize 0 inf]);
+xlabel('(meters)') 
+hold off;
+sp2=subplot(3,1,2);
+hold on;
+plot(xT,yT3,'k','LineWidth',1.5);
+plot(xT,yT6,'r','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([.35*screensize 0.65*screensize 0 inf]);
+xlabel('(meters)') 
+hold off;
+sp3=subplot(3,1,3);
+hold on;
+plot(xT,yT4,'k','LineWidth',1.5);
+plot(xT,yT7,'r','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([.35*screensize 0.65*screensize 0 inf]);
+xlabel('(meters)') 
+hold off;
+
+%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)
+yT3=NormArea(1,3,1,1)*(1/yNorm2(1,3,1,1))*probout4(:,3,1,1);
+yT4=NormArea(1,4,1,1)*(1/yNorm2(1,4,1,1))*probout4(:,4,1,1); 
+yT5=NormArea(1,5,1,1)*(1/yNorm2(1,5,1,1))*probout4(:,5,1,1);
+yT6=NormArea(1,6,1,1)*(1/yNorm2(1,6,1,1))*probout4(:,6,1,1); 
+sp1=subplot(5,1,1);
+plot(xT,yT2,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp2=subplot(5,1,2);
+plot(xT,yT3,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp3=subplot(5,1,3);
+plot(xT,yT4,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp4=subplot(5,1,4);
+plot(xT,yT5,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp5=subplot(5,1,5);
+plot(xT,yT6,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+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 127 MICRON MESH
+figure('Name','100 MICRON BARS');
+xT=screensize/(screennumber-1)*probout4(:,1,1,1);
+yT2=NormArea(1,2,2,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)
+yT3=NormArea(1,3,2,1)*(1/yNorm2(1,3,1,1))*probout4(:,3,2,1);
+yT4=NormArea(1,4,2,1)*(1/yNorm2(1,4,1,1))*probout4(:,4,2,1); 
+yT5=NormArea(1,5,2,1)*(1/yNorm2(1,5,1,1))*probout4(:,5,2,1);
+yT6=NormArea(1,6,2,1)*(1/yNorm2(1,6,1,1))*probout4(:,6,2,1); 
+yT7=NormArea(1,7,2,1)*(1/yNorm2(1,7,1,1))*probout4(:,7,2,1);
+yT8=NormArea(1,8,2,1)*(1/yNorm2(1,8,1,1))*probout4(:,8,2,1); 
+yT9=NormArea(1,9,2,1)*(1/yNorm2(1,9,1,1))*probout4(:,9,2,1);
+yT10=NormArea(1,10,2,1)*(1/yNorm2(1,10,1,1))*probout4(:,10,2,1);
+yT11=NormArea(1,11,2,1)*(1/yNorm2(1,11,1,1))*probout4(:,11,2,1);
+yT12=NormArea(1,12,2,1)*(1/yNorm2(1,12,1,1))*probout4(:,12,2,1);
+yT13=NormArea(1,13,2,1)*(1/yNorm2(1,13,1,1))*probout4(:,13,2,1);
+yT14=NormArea(1,14,2,1)*(1/yNorm2(1,14,1,1))*probout4(:,14,2,1); 
+yT15=NormArea(1,15,2,1)*(1/yNorm2(1,15,1,1))*probout4(:,15,2,1);
+yT16=NormArea(1,16,2,1)*(1/yNorm2(1,16,1,1))*probout4(:,16,2,1); 
+yT17=NormArea(1,17,2,1)*(1/yNorm2(1,17,1,1))*probout4(:,17,2,1);
+yT18=NormArea(1,18,2,1)*(1/yNorm2(1,18,1,1))*probout4(:,18,2,1);
+yT19=NormArea(1,19,2,1)*(1/yNorm2(1,19,1,1))*probout4(:,19,2,1);
+yT20=NormArea(1,20,2,1)*(1/yNorm2(1,20,1,1))*probout4(:,20,2,1);
+yT21=NormArea(1,21,2,1)*(1/yNorm2(1,21,1,1))*probout4(:,21,2,1);
+yT22=NormArea(1,22,2,1)*(1/yNorm2(1,22,1,1))*probout4(:,22,2,1);
+sp1=subplot(5,5,1);
+plot(xT,yT2,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp2=subplot(5,5,2);
+plot(xT,yT3,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp3=subplot(5,5,3);
+plot(xT,yT4,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp4=subplot(5,5,4);
+plot(xT,yT5,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp5=subplot(5,5,5);
+plot(xT,yT6,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp6=subplot(5,5,6);
+plot(xT,yT7,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp7=subplot(5,5,7);
+plot(xT,yT8,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp8=subplot(5,5,8);
+plot(xT,yT9,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp9=subplot(5,5,9);
+plot(xT,yT10,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp10=subplot(5,5,10);
+plot(xT,yT11,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp11=subplot(5,5,11);
+plot(xT,yT12,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp12=subplot(5,5,12);
+plot(xT,yT13,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp13=subplot(5,5,13);
+plot(xT,yT14,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp14=subplot(5,5,14);
+plot(xT,yT15,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp15=subplot(5,5,15);
+plot(xT,yT16,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp16=subplot(5,5,16);
+plot(xT,yT17,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp17=subplot(5,5,17);
+plot(xT,yT18,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp18=subplot(5,5,18);
+plot(xT,yT19,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp19=subplot(5,5,19);
+plot(xT,yT20,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp20=subplot(5,5,20);
+plot(xT,yT21,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp21=subplot(5,5,21);
+plot(xT,yT22,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)')
+sgtitle('100 MICRON BARS')
+
+%FIGURE 254 MICRON MESH
+figure('Name','254 MICRON BARS');
+xT=screensize/(screennumber-1)*probout4(:,1,1);
+yT2=NormArea(1,2,3,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)
+yT3=NormArea(1,3,3,1)*(1/yNorm2(1,3,1,1))*probout4(:,3,3,1);
+yT4=NormArea(1,4,3,1)*(1/yNorm2(1,4,1,1))*probout4(:,4,3,1); 
+yT5=NormArea(1,5,3,1)*(1/yNorm2(1,5,1,1))*probout4(:,5,3,1);
+yT6=NormArea(1,6,3,1)*(1/yNorm2(1,6,1,1))*probout4(:,6,3,1); 
+yT7=NormArea(1,7,3,1)*(1/yNorm2(1,7,1,1))*probout4(:,7,3,1);
+yT8=NormArea(1,8,3,1)*(1/yNorm2(1,8,1,1))*probout4(:,8,3,1); 
+yT9=NormArea(1,9,3,1)*(1/yNorm2(1,9,1,1))*probout4(:,9,3,1);
+yT10=NormArea(1,10,3,1)*(1/yNorm2(1,10,1,1))*probout4(:,10,3,1);
+yT11=NormArea(1,11,3,1)*(1/yNorm2(1,11,1,1))*probout4(:,11,3,1);
+yT12=NormArea(1,12,3,1)*(1/yNorm2(1,12,1,1))*probout4(:,12,3,1);
+yT13=NormArea(1,13,3,1)*(1/yNorm2(1,13,1,1))*probout4(:,13,3,1);
+yT14=NormArea(1,14,3,1)*(1/yNorm2(1,14,1,1))*probout4(:,14,3,1); 
+yT15=NormArea(1,15,3,1)*(1/yNorm2(1,15,1,1))*probout4(:,15,3,1);
+yT16=NormArea(1,16,3,1)*(1/yNorm2(1,16,1,1))*probout4(:,16,3,1); 
+yT17=NormArea(1,17,3,1)*(1/yNorm2(1,17,1,1))*probout4(:,17,3,1);
+yT18=NormArea(1,18,3,1)*(1/yNorm2(1,18,1,1))*probout4(:,18,3,1);
+yT19=NormArea(1,19,3,1)*(1/yNorm2(1,19,1,1))*probout4(:,19,3,1);
+yT20=NormArea(1,20,3,1)*(1/yNorm2(1,20,1,1))*probout4(:,20,3,1);
+yT21=NormArea(1,21,3,1)*(1/yNorm2(1,21,1,1))*probout4(:,21,3,1);
+yT22=NormArea(1,22,3,1)*(1/yNorm2(1,22,1,1))*probout4(:,22,3,1);
+sp1=subplot(5,5,1);
+plot(xT,yT2,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp2=subplot(5,5,2);
+plot(xT,yT3,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp3=subplot(5,5,3);
+plot(xT,yT4,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp4=subplot(5,5,4);
+plot(xT,yT5,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp5=subplot(5,5,5);
+plot(xT,yT6,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp6=subplot(5,5,6);
+plot(xT,yT7,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp7=subplot(5,5,7);
+plot(xT,yT8,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp8=subplot(5,5,8);
+plot(xT,yT9,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp9=subplot(5,5,9);
+plot(xT,yT10,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp10=subplot(5,5,10);
+plot(xT,yT11,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp11=subplot(5,5,11);
+plot(xT,yT12,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp12=subplot(5,5,12);
+plot(xT,yT13,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp13=subplot(5,5,13);
+plot(xT,yT14,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp14=subplot(5,5,14);
+plot(xT,yT15,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp15=subplot(5,5,15);
+plot(xT,yT16,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp16=subplot(5,5,16);
+plot(xT,yT17,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp17=subplot(5,5,17);
+plot(xT,yT18,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp18=subplot(5,5,18);
+plot(xT,yT19,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp19=subplot(5,5,19);
+plot(xT,yT20,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp20=subplot(5,5,20);
+plot(xT,yT21,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp21=subplot(5,5,21);
+plot(xT,yT22,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sgtitle('200 MICRON BARS')
+
+%FIGURE 381 MICRON MESH
+figure('Name','381 MICRON BARS');
+xT=screensize/(screennumber-1)*probout4(:,1,1);
+yT2=NormArea(1,2,4,1)*(1/yNorm2(1,2,1,1))*probout4(:,2,4,1); % probout4(k,a+1,meshindex,phaser) yNorm2(1,a+1,meshindex,phaser) NormArea(1,a+1,meshindex,phaser)
+yT3=NormArea(1,3,4,1)*(1/yNorm2(1,3,1,1))*probout4(:,3,4,1);
+yT4=NormArea(1,4,4,1)*(1/yNorm2(1,4,1,1))*probout4(:,4,4,1); 
+yT5=NormArea(1,5,4,1)*(1/yNorm2(1,5,1,1))*probout4(:,5,4,1);
+yT6=NormArea(1,6,4,1)*(1/yNorm2(1,6,1,1))*probout4(:,6,4,1); 
+yT7=NormArea(1,7,4,1)*(1/yNorm2(1,7,1,1))*probout4(:,7,4,1);
+yT8=NormArea(1,8,4,1)*(1/yNorm2(1,8,1,1))*probout4(:,8,4,1); 
+yT9=NormArea(1,9,4,1)*(1/yNorm2(1,9,1,1))*probout4(:,9,4,1);
+yT10=NormArea(1,10,4,1)*(1/yNorm2(1,10,1,1))*probout4(:,10,4,1);
+yT11=NormArea(1,11,4,1)*(1/yNorm2(1,11,1,1))*probout4(:,11,4,1);
+yT12=NormArea(1,12,4,1)*(1/yNorm2(1,12,1,1))*probout4(:,12,4,1);
+yT13=NormArea(1,13,4,1)*(1/yNorm2(1,13,1,1))*probout4(:,13,4,1);
+yT14=NormArea(1,14,4,1)*(1/yNorm2(1,14,1,1))*probout4(:,14,4,1); 
+yT15=NormArea(1,15,4,1)*(1/yNorm2(1,15,1,1))*probout4(:,15,4,1);
+yT16=NormArea(1,16,4,1)*(1/yNorm2(1,16,1,1))*probout4(:,16,4,1); 
+yT17=NormArea(1,17,4,1)*(1/yNorm2(1,17,1,1))*probout4(:,17,4,1);
+yT18=NormArea(1,18,4,1)*(1/yNorm2(1,18,1,1))*probout4(:,18,4,1);
+yT19=NormArea(1,19,4,1)*(1/yNorm2(1,19,1,1))*probout4(:,19,4,1);
+yT20=NormArea(1,20,4,1)*(1/yNorm2(1,20,1,1))*probout4(:,20,4,1);
+yT21=NormArea(1,21,4,1)*(1/yNorm2(1,21,1,1))*probout4(:,21,4,1);
+yT22=NormArea(1,22,4,1)*(1/yNorm2(1,22,1,1))*probout4(:,22,4,1);
+sp1=subplot(5,5,1);
+plot(xT,yT2,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp2=subplot(5,5,2);
+plot(xT,yT3,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp3=subplot(5,5,3);
+plot(xT,yT4,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp4=subplot(5,5,4);
+plot(xT,yT5,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp5=subplot(5,5,5);
+plot(xT,yT6,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp6=subplot(5,5,6);
+plot(xT,yT7,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp7=subplot(5,5,7);
+plot(xT,yT8,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp8=subplot(5,5,8);
+plot(xT,yT9,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp9=subplot(5,5,9);
+plot(xT,yT10,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp10=subplot(5,5,10);
+plot(xT,yT11,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp11=subplot(5,5,11);
+plot(xT,yT12,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp12=subplot(5,5,12);
+plot(xT,yT13,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp13=subplot(5,5,13);
+plot(xT,yT14,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp14=subplot(5,5,14);
+plot(xT,yT15,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp15=subplot(5,5,15);
+plot(xT,yT16,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp16=subplot(5,5,16);
+plot(xT,yT17,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp17=subplot(5,5,17);
+plot(xT,yT18,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp18=subplot(5,5,18);
+plot(xT,yT19,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp19=subplot(5,5,19);
+plot(xT,yT20,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp20=subplot(5,5,20);
+plot(xT,yT21,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp21=subplot(5,5,21);
+plot(xT,yT22,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sgtitle('300 MICRON BARS')
+
+%FIGURE 508 MICRON MESH
+figure('Name','400 MICRON BARS');
+xT=screensize/(screennumber-1)*probout4(:,1,1);
+yT2=NormArea(1,2,5,1)*(1/yNorm2(1,2,1,1))*probout4(:,2,5,1); % probout4(k,a+1,meshindex,phaser) yNorm2(1,a+1,meshindex,phaser) NormArea(1,a+1,meshindex,phaser)
+yT3=NormArea(1,3,5,1)*(1/yNorm2(1,3,1,1))*probout4(:,3,5,1);
+yT4=NormArea(1,4,5,1)*(1/yNorm2(1,4,1,1))*probout4(:,4,5,1); 
+yT5=NormArea(1,5,5,1)*(1/yNorm2(1,5,1,1))*probout4(:,5,5,1);
+yT6=NormArea(1,6,5,1)*(1/yNorm2(1,6,1,1))*probout4(:,6,5,1); 
+yT7=NormArea(1,7,5,1)*(1/yNorm2(1,7,1,1))*probout4(:,7,5,1);
+yT8=NormArea(1,8,5,1)*(1/yNorm2(1,8,1,1))*probout4(:,8,5,1); 
+yT9=NormArea(1,9,5,1)*(1/yNorm2(1,9,1,1))*probout4(:,9,5,1);
+yT10=NormArea(1,10,5,1)*(1/yNorm2(1,10,1,1))*probout4(:,10,5,1);
+yT11=NormArea(1,11,5,1)*(1/yNorm2(1,11,1,1))*probout4(:,11,5,1);
+yT12=NormArea(1,12,5,1)*(1/yNorm2(1,12,1,1))*probout4(:,12,5,1);
+yT13=NormArea(1,13,5,1)*(1/yNorm2(1,13,1,1))*probout4(:,13,5,1);
+yT14=NormArea(1,14,5,1)*(1/yNorm2(1,14,1,1))*probout4(:,14,5,1); 
+yT15=NormArea(1,15,5,1)*(1/yNorm2(1,15,1,1))*probout4(:,15,5,1);
+yT16=NormArea(1,16,5,1)*(1/yNorm2(1,16,1,1))*probout4(:,16,5,1); 
+yT17=NormArea(1,17,5,1)*(1/yNorm2(1,17,1,1))*probout4(:,17,5,1);
+yT18=NormArea(1,18,5,1)*(1/yNorm2(1,18,1,1))*probout4(:,18,5,1);
+yT19=NormArea(1,19,5,1)*(1/yNorm2(1,19,1,1))*probout4(:,19,5,1);
+yT20=NormArea(1,20,5,1)*(1/yNorm2(1,20,1,1))*probout4(:,20,5,1);
+yT21=NormArea(1,21,5,1)*(1/yNorm2(1,21,1,1))*probout4(:,21,5,1);
+yT22=NormArea(1,22,5,1)*(1/yNorm2(1,22,1,1))*probout4(:,22,5,1);
+sp1=subplot(5,5,1);
+plot(xT,yT2,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp2=subplot(5,5,2);
+plot(xT,yT3,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp3=subplot(5,5,3);
+plot(xT,yT4,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp4=subplot(5,5,4);
+plot(xT,yT5,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp5=subplot(5,5,5);
+plot(xT,yT6,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp6=subplot(5,5,6);
+plot(xT,yT7,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp7=subplot(5,5,7);
+plot(xT,yT8,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp8=subplot(5,5,8);
+plot(xT,yT9,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp9=subplot(5,5,9);
+plot(xT,yT10,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp10=subplot(5,5,10);
+plot(xT,yT11,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp11=subplot(5,5,11);
+plot(xT,yT12,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp12=subplot(5,5,12);
+plot(xT,yT13,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp13=subplot(5,5,13);
+plot(xT,yT14,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp14=subplot(5,5,14);
+plot(xT,yT15,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp15=subplot(5,5,15);
+plot(xT,yT16,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp16=subplot(5,5,16);
+plot(xT,yT17,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp17=subplot(5,5,17);
+plot(xT,yT18,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp18=subplot(5,5,18);
+plot(xT,yT19,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp19=subplot(5,5,19);
+plot(xT,yT20,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp20=subplot(5,5,20);
+plot(xT,yT21,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)') 
+sp21=subplot(5,5,21);
+plot(xT,yT22,'-k','LineWidth',1.5);
+xline(testline1,'--r');
+xline(testline2,'--r');
+axis([0 0.03 0 inf]);
+xlabel('Screen (meters)')
+sgtitle('400 MICRON BARS')
+
+
+%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 6 V^2 VS A FOR VARIOUS SLIT INTENSITY
+figure('Name','Visibility vs a');
+xadX1=totaladataX(:,1,1);
+yadV1=totaladataV1(:,1,1); %(a,meshindex,phaser)
+yadV2=totaladataV1(:,2,1);
+yadV3=totaladataV1(:,3,1);
+hold on;
+plot(xadX1, yadV1,'*-')
+plot(xadX1, yadV2,'*-')
+plot(xadX1, yadV3,'*-')
+hold off;
+xlabel('a-value') 
+ylabel('V^2') 
+title('Visibility')
+legend('0 micron-p1','127 micron-p1','508 micron-p1')
+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