%bessel3D.m
%For 3D plotting of diffraction from a cicular aperture
%Initial concept made by Alex Fritzler
%Extensively modified and adapted for 3D printing by Steven Sahyun
%10/25/2015

%Physical modeling
D = 69;%distance to screen in cm
rmax = 2;%scan distance in cm

[X,Y] = meshgrid(-rmax:0.01:rmax);
r = sqrt(X.^2 + Y.^2)+eps;

thetar = atan(r./D);

d = 60000; %Aperture size in nm
lambda = 655;%wavelength in nm
k = 2*pi/lambda;

I_0 = 200;

J_1r = besselj(1, k*d*sin(thetar)); %Bessel function of the first kind.
I_thetar = I_0*((2.*(J_1r))./(k*d*sin(thetar))).^2;% Eq. 6.18 p. 143

%phi = 0:0.1:2*pi;

figure(1)
plot(r, J_1r); %Plot of a bessel function
title('Bessel J_1(krsin\theta)')

figure(2)
plot(r, I_thetar); %Plot of the intensity
title('Intensity I(\theta)')

figure(3)
z = I_thetar; %3D graph of the intensity
meshz(z);
title('Intensity I(\theta)')

%Our eyes have a logrithmic response to intensity so we can see
%dim objects near bright ones. Plotting the square root of the 
%intensity better represents our perception of the pattern.
%It is similar to plotting out an "overexposed" graph to enhance
%the weaker diffraction pattern lines.

figure(4)
z = sqrt(I_thetar); %What we percieve with our eyes.
meshz(z);
title('Intensity sqrt(I(\theta))')

%convert to a 3D printable file; You need to have stlwrite.m and surf2solid.m by
%Sven Holcombe in the same directory as this file.
%See: http://jmumakerlab.blogspot.com/2013/11/exporting-stl-from-matlab.html
% http://www.mathworks.com/matlabcentral/fileexchange/42876-surf2solid-make-a-solid-volume-from-a-surface-for-3d-printing
% http://www.mathworks.com/matlabcentral/fileexchange/20922-stlwrite-filename--varargin-

CircIntens = surf2solid(10*X, 10*Y, z, 'elevation', -1);
stlwrite('CircIntens.stl',CircIntens) 


%To graph as a smooth surface
%figure(4)
%surf(z);
%surf(z,'FaceColor','interp',...
%   'EdgeColor','none',...
%   'FaceLighting','phong')
%daspect([5 5 0.1])
%axis tight
%view(-50,30)
camlight left