MATLAB: Plotting the product of two functions.

curve fitting

Hi,
I have plotted two offset Bessel functions of the first kind (blue and red dotted curves) – see figure below.
I would like to produce something similar to the black curve shown in the second attached figure below, but for the two curves I have plotted. Can someone please point me in the correct direction?
Thank you
lambda = 500e-6;                     % wavelength [m]
D = 0.0001 ;                         % diameter D of the primary aperture [m]
theta1 = -20:.01:20;                 % angular radius θ (as measured from the primary aperture) [Degrees]
u1 = (pi/lambda)*D*theta1;
theta2 = -20:.01:20;  
u2 = (pi/lambda)*D*theta2;
x1 = u1./pi +0.61;
x2 = u2./pi -0.61;
I01 = 1.0; 
I02 = 1.0;
i1 = real((besselj(1,u1)./(u1)).^2);
i2 = real((besselj(1,u2)./(u2)).^2);
I1 = I01*(i1./max(i1));
I2 = I02*(i2./max(i2));
colormap('default')
figure(1);
%clf;
hold on;
plot(x1,I1,'b:','LineWidth',1.5);
plot(x2,I2,'r:','LineWidth',1.5);
hold off;
grid on; 
set(gca,'FontSize',14);
%axis([-4 4 0 1.05])
xlabel('\theta (\lambda/D)', 'fontsize', 16);
ylabel('Normalized Intensity', 'fontsize', 16);
legend('Source 1','Source 2','fontsize', 10);
legend;
figure(2);
clf;
[X,Y]=meshgrid(-10:.1:10);
%[X,Y]=meshgrid(x);
R=sqrt(X.^2+Y.^2);
z=(besselj(1,R)./R).^2;
Z=z./max(max(z));
surf(X,Y,Z,'EdgeColor','none');
colorbar
camlight left;lighting phong;
set(gca,'XTick',[], 'YTick', [], 'ZTick', [])

Best Answer

As pointed out by John, this is the sum of these two curves. Check the following code. I had to do a bit of interpolating because of the way you defined x1 and x2. I have added comments on the modified line. The graph looks quite similar to the graph you want
lambda = 500e-6;                     % wavelength [m]
D = 0.0001 ;                         % diameter D of the primary aperture [m]
theta1 = -20:.01:20;                 % angular radius θ (as measured from the primary aperture) [Degrees]
u1 = (pi/lambda)*D*theta1;
theta2 = -20:.01:20;  
u2 = (pi/lambda)*D*theta2;
x1 = u1./pi +0.8;
x2 = u2./pi -0.8;
I01 = 1.0; 
I02 = 0.5;  % <<==== changed this parameter
i1 = real((besselj(1,u1)./(u1)).^2);
i2 = real((besselj(1,u2)./(u2)).^2);
I1 = I01*(i1./max(i1));
I2 = I02*(i2./max(i2));
% vvv Interpolated to make both vector equal
x = unique([x1 x2]);
I1_ = interp1(x1, I1, x);
I2_ = interp1(x2, I2, x);
colormap('default')
figure(1);
%clf;
hold on;
plot(x1,I1,'b:','LineWidth',1.5);
plot(x2,I2,'r:','LineWidth',1.5);
plot(x,I1_+I2_,'r:','LineWidth',1.5); % <<==== added this line
hold off;
grid on; 
set(gca,'FontSize',14);
%axis([-4 4 0 1.05])
xlabel('\theta (\lambda/D)', 'fontsize', 16);
ylabel('Normalized Intensity', 'fontsize', 16);
legend('Source 1','Source 2','fontsize', 10);
legend;
figure(2);
clf;
[X,Y]=meshgrid(-10:.1:10);
%[X,Y]=meshgrid(x);
R=sqrt(X.^2+Y.^2);
z=(besselj(1,R)./R).^2;
Z=z./max(max(z));
surf(X,Y,Z,'EdgeColor','none');
colorbar
camlight left;lighting phong;
set(gca,'XTick',[], 'YTick', [], 'ZTick', [])
Related Question