hi ,i have this function that ode45 integrate, the numerical integration is correct, but in the main script I have to plot the vector 'alfa', and in order to do this i have to rewrite the same code, but i need to do a cicle 'for', how is this cicle ?

Question

Alessandro Vitiello 2018-10-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/422027-hi-i-have-this-function-that-ode45-integrate-the-numerical-integration-is-correct-but-in-the-main

关闭： MATLAB Answer Bot 2021-8-20

function [ derivate ] = eq_moto_vincolo_q_piatto(t, X )
format long
global mi m CQ T w A rho0 H R_earth beta a0 a1 b0 b1 b2 c0 c1 c2 c3  
global tD_ott alfaD_ott sigmaD_ott
global t_pre alfa_pre
%vettore incognita
r=X(1);
lambda=X(2);
lat=X(3);
gamma=X(4);
v=X(5);
zita=X(6);
calore_assorbito=X(7);
%%equazioni del moto
%equazioni cinematiche
if r<R_earth
  v=0;
end
dr=v*sin(gamma);
dlambda=v*cos(gamma)*cos(zita)/(r*cos(lat));
dlat=v*cos(gamma)*sin(zita)/r;
%densità
rho=rho0*exp(-(r-R_earth)/H);
%parametri aerodinamici
alfa_tilde=interp1(tD_ott,alfaD_ott,t,'spline');
flusso_tilde = ( c0 + c1.*(alfa_tilde) + c2.*(alfa_tilde).^2 + c3.*(alfa_tilde).^3 ) * ( 17700 .* sqrt(rho * 0.002970/1.75e9 ) .* (0.0001.*v * 3280.84 ).^(3.07) );    %BTU/ft^2 sec  
if flusso_tilde < 70
    alfa = alfa_tilde
else
      %c0 + c1*(alfa) + c2*(alfa)^2 + c3*(alfa)^3 = termine_noto ;
      termine_noto = 70 / ( 17700 * sqrt(rho * 0.002970/1.75e9 ) * (0.0001 * v * 3280.84 )^(3.07) ) ;
      coefficienti = [c3 c2 c1 c0-termine_noto];
alfa_sol = roots(coefficienti)
ind_alfa_sol=find(imag(alfa_sol));
alfa_sol_complessi=alfa_sol(ind_alfa_sol);
alfa_sol(ind_alfa_sol)=[]
ind = find ( min(abs(alfa_sol - alfa_pre)) );
alfa = alfa_sol(ind)
end
sigma_int=interp1(tD_ott,sigmaD_ott,t,'spline');
sigma = sigma_int * pi/180;
CL = a0 + a1*(alfa);
CD = b0 + b1*(alfa) + b2*(alfa)^2;
%forza aerodinamica
L=0.5*rho*CL*A*v^2;
D=0.5*rho*CD*A*v^2;
Q=0.5*rho*CQ*A*v^2;
%equazioni dinamiche
dgamma= ( -mi/(r^2*v) + v/r ) * cos(gamma) + ( T*cos(beta)*sin(alfa)+L*cos(sigma)+Q*sin(sigma) )/(m*v) + (w^2*r/v)* cos(lat)*( cos(lat)*cos(gamma) + sin(lat)*sin(gamma)*sin(zita) ) +2*w*cos(lat)*cos(zita); 
dv= -(mi/r^2)*sin(gamma) + ( T*cos(beta)*cos(alfa)-D )/(m) + w^2*r*cos(lat) * ( cos(lat)*sin(gamma) - sin(lat)*cos(gamma)*sin(zita) )  ; 
dzita= -v/r * tan(lat)*cos(gamma)*cos(zita) + ( T*sin(beta)-L*sin(sigma)+Q*cos(sigma) )/( m*v*cos(gamma) ) + 2*w*cos(lat)*tan(gamma)*sin(zita) - (w^2*r)/(v*cos(gamma))*sin(lat)*cos(lat)*cos(zita) - 2*w*sin(lat); 
%flusso di calore
flusso = ( c0 + c1.*(alfa) + c2.*(alfa).^2 + c3.*(alfa).^3 ) * ( 17700 .* sqrt(rho * 0.002970/1.75e9 ) .* (0.0001.*v * 3280.84 ).^(3.07) )   %BTU/ft^2 sec  
if t>t_pre
    t_pre = t;
    alfa_pre = alfa;
end
%funzione da integrare
derivate=[dr; dlambda; dlat; dgamma; dv; dzita; flusso];

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2018-10-3

You can use a global or a (better) a shared variable to keep a record of alfa, probably updating at your "if t>t_pre" . You should probably keep a record of the associated t value as well.

此问题已关闭。

hi ,i have this function that ode45 integrate, the numerical integration is correct, but in the main script I have to plot the vector 'alfa', and in order to do this i have to rewrite the same code, but i need to do a cicle 'for', how is this cicle ?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

回答（0 个）

另请参阅

标签

Community Treasure Hunt

hi ,i have this function that ode45 integrate, the numerical integration is correct, but in the main script I have to plot the vector 'alfa', and in order to do this i have to rewrite the same code, but i need to do a cicle 'for', how is this cicle ?

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

回答（0 个）

另请参阅

标签

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论