Calling a function inside an integration

조회 수: 10 (최근 30일)
Ali Almakhmari
Ali Almakhmari 2022년 3월 23일
댓글: Walter Roberson 2022년 3월 24일
So I posted a snipet of my code below. The issue I am having here is regarding how to syntax my code correctly in the integration line inside the for loop. I want to call the function for each "i" iteration and for each iteration of the integral to determine if "groundTemp(yi, theta, phi)" should be replaced by 90 or 250 depending on the condition stated in the function. I just don't know how to syntax this in a way to get rid of the error and do what I want. Any help would be appreciated.
t = 36:1:64;
v = 7.65;
h = 435;
r = tand(15)*h;
lambda = 0.21;
Aeff = 0.57;
for i = 1:length(t)
yi = r - v*(t(i));
T_AP(i) = ((lambda^2)/Aeff) .* integral2(@(theta, phi) groundTemp(yi, theta, phi) .* exp(-2.77.*(theta./(pi/12)).^2).*sin(theta), 0, pi/2, 0, 2*pi);
end
function [temp] = groundTemp(yi, theta, phi)
range = (435).*tan(theta).*cos(phi);
if(range > yi)
temp = 90;
elseif(range < yi)
temp = 250;
end
end
  댓글 수: 1
Walter Roberson
Walter Roberson 2022년 3월 24일
I had hoped to make progress using the Symbolic Toolbox. It turned out to get stuck, but the intermediate form is one I have never seen before:
Q = @(v) sym(v); %convert to symbolic
v = Q(7.65);
h = Q(435);
r = tand(Q(15))*h;
lambda = Q(0.21);
Aeff = Q(0.57);
syms phi theta T real
Pi = Q(pi);
yi = r - v*(T);
expr = groundTemp(yi, theta, phi) .* exp(-Q(2.77).*(theta./(Pi/12)).^2).*sin(theta)
tic, inner = int(expr, theta, 0, Pi/2), toc
tic, outer = int(inner, phi, 0, 2*pi), toc
term = outer;
T_AP = ((lambda^2)/Aeff) .* term
function [temp] = groundTemp(yi, theta, phi)
range = (435).*tan(theta).*cos(phi);
temp = piecewise(range > yi, 90, range < yi, 250, nan)
end
This generates
(147*int(intlib::intOverSet(250*exp(-(9972*theta^2)/(25*pi^2))*sin(theta), theta, ([0, pi/2] minus solve(-((51*T)/2900 + 3^(1/2) - 2)/cos(phi) < tan(theta), theta, NoSpuriousSolutions)) intersect solve(tan(theta) < -((51*T)/2900 + 3^(1/2) - 2)/cos(phi), theta, NoSpuriousSolutions)) + intlib::intOverSet(90*exp(-(9972*theta^2)/(25*pi^2))*sin(theta), theta, [0, pi/2] intersect solve(-((51*T)/2900 + 3^(1/2) - 2)/cos(phi) < tan(theta), theta, NoSpuriousSolutions)), phi, 0, pi/2))/1900 + (147*int(intlib::intOverSet(250*exp(-(9972*theta^2)/(25*pi^2))*sin(theta), theta, ([0, pi/2] minus solve(-((51*T)/2900 + 3^(1/2) - 2)/cos(phi) < tan(theta), theta, NoSpuriousSolutions)) intersect solve(tan(theta) < -((51*T)/2900 + 3^(1/2) - 2)/cos(phi), theta, NoSpuriousSolutions)) + intlib::intOverSet(90*exp(-(9972*theta^2)/(25*pi^2))*sin(theta), theta, [0, pi/2] intersect solve(-((51*T)/2900 + 3^(1/2) - 2)/cos(phi) < tan(theta), theta, NoSpuriousSolutions)), phi, (3*pi)/2, 2*pi))/1900 + (147*int(intlib::intOverSet(250*exp(-(9972*theta^2)/(25*pi^2))*sin(theta), theta, ([0, pi/2] minus solve(tan(theta) < -((51*T)/2900 + 3^(1/2) - 2)/cos(phi), theta, NoSpuriousSolutions)) intersect solve(-((51*T)/2900 + 3^(1/2) - 2)/cos(phi) < tan(theta), theta, NoSpuriousSolutions)) + intlib::intOverSet(90*exp(-(9972*theta^2)/(25*pi^2))*sin(theta), theta, [0, pi/2] intersect solve(tan(theta) < -((51*T)/2900 + 3^(1/2) - 2)/cos(phi), theta, NoSpuriousSolutions)), phi, pi/2, (3*pi)/2))/1900
The intlib::intOverSet is obviously something internal to MuPAD.
Unfortunately at the moment I can't seem to take it any further than that.

댓글을 달려면 로그인하십시오.

채택된 답변

Torsten
Torsten 2022년 3월 24일
편집: Torsten 2022년 3월 24일
t = 36:1:64;
v = 7.65;
h = 435;
r = tand(15)*h;
lambda = 0.21;
Aeff = 0.57;
for i = 1:length(t)
yi = r - v*(t(i));
T_AP(i) = lambda^2/Aeff*integral2(@(theta,phi)fun(yi,theta,phi), 0, pi/2, 0, 2*pi);
end
plot(t,T_AP)
function value = fun(yi,theta,phi)
for i=1:size(theta,1)
for j=1:size(theta,2)
temp = groundTemp(yi, theta(i,j), phi(i,j));
value(i,j) = temp*exp(-2.77.*(theta(i,j)/(pi/12))^2)*sin(theta(i,j));
end
end
end
function temp = groundTemp(yi, theta, phi)
range = (435).*tan(theta).*cos(phi);
if range > yi
temp = 90;
elseif range <= yi
temp = 250;
end
end

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by