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', [])