Calling a function inside an integration
7 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Ali Almakhmari
il 23 Mar 2022
Commentato: Walter Roberson
il 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 Commento
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.
Risposta accettata
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
0 Commenti
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!