Can someone help me with this integration

In addition to getting an error in this, what I want to essentially do is replace f1(1, i) in the integ2 variable, but I keep getting f1(1, i) as it is. Its like I am telling MATLAB to evaluate f1(1,i) with the actual f1(1,i) from the Workspace but its just leaving it as it is.

clear
clc
beep off
theta_s = 0:0.1:pi/2;
syms theta_v phi 
f1 = ((((1/(2.*pi)).*((pi-phi).*cos(phi) + sin(phi)).*tan(theta_s).*tan(theta_v)) - (1/pi).*(tan(theta_s)+tan(theta_v)+sqrt(tan(theta_v).^2 + tan(theta_s).^2 - 2.*tan(theta_s).*tan(theta_v).*cos(phi)))));

for i = 1:length(theta_s)
    integ2 = @(theta_v, phi) f1(1,i).*cos(theta_v).*sin(theta_v);
    r2(i)= integral2(integ2, 0, pi/2, 0, pi);
end

Error using integral2Calc>integral2t/tensor (line 231)
Input function must return 'double' or 'single' values. Found 'sym'.

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 105)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

Answer 1

[Torsten]

THETA_S = 0:0.1:pi/2;
for i = 1:length(THETA_S)
    theta_s = THETA_S(i);
    f1 = @(theta_v,phi) ((((1/(2.*pi)).*((pi-phi).*cos(phi) + sin(phi)).*tan(theta_s).*tan(theta_v)) - (1/pi).*(tan(theta_s)+tan(theta_v)+sqrt(tan(theta_v).^2 + tan(theta_s).^2 - 2.*tan(theta_s).*tan(theta_v).*cos(phi)))));
    integ2 = @(theta_v, phi) f1(theta_v,phi).*cos(theta_v).*sin(theta_v);
    r2(i)= integral2(integ2, 0, pi/2, 0, pi);
end

[Ali Almakhmari]

https://youtu.be/pPM5NnkoTfI?t=14

Answer 2

[AndresVar]

you can do it without symbolic variables

clear
theta_s_vec = 0:0.1:pi/2;
for ii = 1:numel(theta_s_vec)
    theta_s = theta_s_vec(ii);
    integ2 = @(phi,theta_v) cos(theta_v).*sin(theta_v).*((((1/(2.*pi)).*((pi-phi).*cos(phi) + sin(phi)).*tan(theta_s).*tan(theta_v)) - (1/pi).*(tan(theta_s)+tan(theta_v)+sqrt(tan(theta_v).^2 + tan(theta_s).^2 - 2.*tan(theta_s).*tan(theta_v).*cos(phi)))));
    r2(ii) = integral2(integ2,0,pi,0,pi/2);
end
r2

r2 = 1×16

-1.5708 -1.5728 -1.5788 -1.5893 -1.6048 -1.6265 -1.6559 -1.6956 -1.7494 -1.8238 -1.9300 -2.0887 -2.3439 -2.8051 -3.8411 -7.9164

You can try symbolically and approximate, you get the same result. But can be slower sometimes.

clear
theta_s = 0:0.1:pi/2;
syms theta_v phi real positive
f1 = ((((1/(2.*pi)).*((pi-phi).*cos(phi) + sin(phi)).*tan(theta_s).*tan(theta_v)) - (1/pi).*(tan(theta_s)+tan(theta_v)+sqrt(tan(theta_v).^2 + tan(theta_s).^2 - 2.*tan(theta_s).*tan(theta_v).*cos(phi)))));
integ2 = f1*cos(theta_v)*sin(theta_v);

for ii = 1:numel(theta_s)
    r2(ii)=int(int(integ2(ii),phi,[0 pi]),theta_v,[0 pi/2]);
end
vpa(r2,5)

ans =