farfield faurhofer diffraction power ratio mismatch

2 visualizaciones (últimos 30 días)
Sean
Sean el 14 de Mayo de 2023
Comentada: Sean el 15 de Mayo de 2023
i calculated power ratio of diffraction but there seems to be some problem as power fraction is almost zero and power fraction farfield is almost one.
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;
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
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_farfield)]);

Respuestas (1)

Shaik
Shaik el 14 de Mayo de 2023
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)]);
Power fraction within cone: 5.9128e-05
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.
  1 comentario
Sean
Sean el 15 de Mayo de 2023
Thank you, i tried with this code but power fraction is still very low for near field and 1 for farfield

Iniciar sesión para comentar.

Categorías

Más información sobre Bessel functions en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by