Info

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 view (last 30 days)
Alessandro Vitiello on 3 Oct 2018
Closed: MATLAB Answer Bot on 20 Aug 2021
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 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
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 = 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];
Walter Roberson on 3 Oct 2018
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.