# farfield faurhofer diffraction power ratio mismatch

2 ビュー (過去 30 日間)
Sean 2023 年 5 月 14 日
コメント済み: Sean 2023 年 5 月 15 日
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)]);

サインインしてコメントする。

### 回答 (1 件)

Shaik 2023 年 5 月 14 日
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 件のコメント-1 件の古いコメントを表示-1 件の古いコメントを非表示
Sean 2023 年 5 月 15 日
Thank you, i tried with this code but power fraction is still very low for near field and 1 for farfield

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeBessel functions についてさらに検索

### Community Treasure Hunt

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

Start Hunting!

Translated by