Calling a function inside an integration
10 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Ali Almakhmari
le 23 Mar 2022
Commenté : Walter Roberson
le 24 Mar 2022
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 commentaire
Walter Roberson
le 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.
Réponse acceptée
Torsten
le 24 Mar 2022
Modifié(e) : Torsten
le 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
0 commentaires
Plus de réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!