How to calculate a double integral with a subfunction?

I try to calculate a double integral for the subfunction "pattern_element", whose varibles are "theta“ and ”phi”. Both varibles are within [0, 2π]. But an error exists when I used the "integral2". I am wondering why the following code does not work. Thank you for your anttenion

%% basic parameters

f0 = 5e9; % frequency:Hz

W=16.40*10^(-3); % Width：m

L=13.52*10^(-3); % Length：m

obj=integral2(@(theta,phi)pattern_element(f0,theta,phi,L,W),0, 2*pi,0, 2*pi)

function [FieldThetaPhi]=pattern_element(f0,theta,phi,L,W)

c0 = 299792458; %velocity of light

lamda0 = c0/f0; %wavelength

k=2*pi/lamda0; %wavenumber

cW=0.5*k*W*sin(theta)*sin(phi);

cL=0.5*k*L*sin(theta)*cos(phi);

sincW=sin(cW)/cW;

if cW==0

sincW=1;

end

FieldTheta = sincW*cos(cL)*cos(phi);

FieldPhi = -sincW*cos(cL)*cos(theta)*sin(phi);

FieldThetaPhi= sqrt(FieldTheta.^2+FieldPhi.^2);

end

Answer 1

[Paul]

Hi Dewen Yu,

As the error message indicated, the function pattern_element needs to use element-wise math operations to compute the integrand for multiple values of the variables. The code already did that in the computation of FieldThetaPhi using the .^ operators, but nowhere else. I modifed the code to use the elementwise operators everywhere, and changed the check on the sinc calculation to a logical indexing that works over a vector rather than the if statemement that was more suitable for a scalar.

%% basic parameters
f0 = 5e9;            % frequency:Hz
W=16.40*10^(-3);     % Width：m
L=13.52*10^(-3);     % Length：m
obj=integral2(@(theta,phi)pattern_element(f0,theta,phi,L,W),0, 2*pi,0, 2*pi)

obj = 30.3904

function [FieldThetaPhi]=pattern_element(f0,theta,phi,L,W)  
c0 = 299792458;        %velocity of light
lamda0 = c0./f0;        %wavelength
k=2.*pi/lamda0;         %wavenumber
cW=0.5.*k.*W.*sin(theta).*sin(phi);
cL=0.5.*k.*L.*sin(theta).*cos(phi);
sincW=sin(cW)./cW;
%if cW==0
%    sincW=1;
%end
sincW(cW==0) = 1;
FieldTheta = sincW.*cos(cL).*cos(phi);
FieldPhi   = -sincW.*cos(cL).*cos(theta).*sin(phi);
FieldThetaPhi = sqrt(FieldTheta.^2+FieldPhi.^2);
end

[Dewen Yu]

Thank you very much for your answer, Paul. It works. I have learned from your code.

[Walter Roberson]

This is the way. integral2 passes in arrays so you need to use element-wise operations

Answer 2

[Torsten]

Use

obj = integral2(@(theta,phi)arrayfun(@(x,y)pattern_element(f0,x,y,L,W),theta,phi),0, 2*pi,0, 2*pi)

instead of

obj=integral2(@(theta,phi)pattern_element(f0,theta,phi,L,W),0, 2*pi,0, 2*pi)

[Dewen Yu]

Your great anwer is appreciated. Thank you very much for your valuable suggestion.