Hi Sean,
The issue in your code lies in the calculation of the variable beta inside the second loop. Currently, you calculate beta outside of the loop, and its value remains constant throughout the loop iterations. To fix this, you need to calculate beta inside the loop so that it corresponds to the current value of q. Here's the modified code:
clc
close all
clear all
lambda = 332.8e-9; % wavelength of light in vacuum
k = 2*pi/lambda;
a = 3e-7; % radius of diffracting circular aperture
Io = 1.0; % relative intensity
R = 1; % distance of screen from aperture
half_angle = 13.3; % half angle of light cone
center_angle = 0; % center angle of light cone
Y = (-0.55e-0:1e-3:0.9e-0);
Z = Y; % coordinates of the screen
I(1:length(Y),1:length(Z))=0;
for i=1:length(Y)
for j=1:length(Z)
q = (Y(i).^2+Z(j).^2).^0.5;
beta = k*a*q/R; % Calculate beta inside the loop
I(i,j) = Io.*((besselj(1,(beta))/(beta+1e-12)).^2);
end
end
% Calculate fraction of power in cone
theta = atan2(Y,Z);
cone_indices = (theta > center_angle-half_angle) & (theta < center_angle+half_angle);
power_fraction = sum(sum(I(cone_indices)))/sum(sum(I));
% Calculate far field using Fraunhofer approximation
N = 1000; % number of points in far field
theta = linspace(-pi/2,pi/2,N); % angular coordinates in far field
E_theta = zeros(1,N);
for i=1:length(theta)
q = k*a*sin(theta(i)); % distance from center of aperture in far field
beta = k*a*q/R; % Calculate beta inside the loop
E_theta(i) = (Io/(2*pi*R*q))*((2*besselj(1,beta)/(beta+1e-12))*exp(-1i*k*q))/(q);
end
% Calculate power in cone in far field
theta = abs(theta)*180/pi;
cone_indices = (theta > center_angle-half_angle) & (theta < center_angle+half_angle);
power_fraction_farfield = sum(abs(E_theta(cone_indices)).^2)/sum(abs(E_theta).^2);
disp(['Power fraction within cone: ' num2str(power_fraction)]);
By calculating beta inside the loop, you ensure that it corresponds to the current value of q for each iteration. This should provide the correct power fraction within the cone and far field.
