Calling a function inside an integration

12 visualizzazioni (ultimi 30 giorni)
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 Commento
Walter Roberson
Walter Roberson il 24 Mar 2022
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.

Accedi per commentare.

Risposta accettata

Torsten
Torsten il 24 Mar 2022
Modificato: Torsten il 24 Mar 2022
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

Più risposte (0)

Categorie

Scopri di più su Programming in Help Center e File Exchange

Prodotti


Release

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by