Using a matrix of solution elements outside of ode45

1 view (last 30 days)
I am working on a problem for a falling parachute where the forces are proportional to velocity squared. I am using ode45 which calls a function which includes the differential equations of motion. ode45 outputs to the main program time, distance, and velocity. Within the function, I have built a matrix which includes time, distance, velocity, glide path angle, and acceleration using persistent variables. I would like to use this complete matrix after ode45 has finished in the main program but I only get the variables for the last time period and not the whole matrix. How do I make the whole matrix from the function inside the ode45 loop available outside ode45 in the main program?

Answers (2)

Mischa Kim
Mischa Kim on 10 Dec 2016
Duncan, how about simply re-building the matrix after the ode45-call? I assume all those variables are dependent on time, distance, and velocity.
For a more detailed analysis, please attach your code.
  2 Comments
Duncan McIntosh
Duncan McIntosh on 11 Dec 2016
Edited: Walter Roberson on 12 Dec 2016
function [T,T1,Y,Y1]=solution_module;
data_module;
global ACCEL r_o_enu;
v_0=0;a_0=0; alt=r_o_enu(3); % Initial Conditions for unsteady state ode. Parachute being deployed.
[T,Y]=ode45('usm_module',[0:1:5],[v_0;a_0]);
V=Y(6,2); A=ACCEL; % Initial Conditions (velocity, acceleration) for steady state ode. Parachute deployed.
[T1,Y1]=ode45('ssm_module',[5:1:450],[V;A])
end
$$$$$$$$$$$$$$$$$$$
function [F]=usm_module(t,x);
data_module;
global W_pl S_bar_pl S_bar_c b e r_o_enu g ACCEL;
alt=r_o_enu(3)+x(1); [rho] = rho_module(alt);
AR=(b.^2)./S_bar_c; m=W_pl./g;
A=.1.*rho.*S_bar_c.*t;
B=(rho./2).*(S_bar_c.*((.01+(.19./5).*t)+(.2.*t).^2./(pi.*e.*AR))+S_bar_pl.*2);
gama=-atan(B./A).*(180./pi);
k=sqrt(A.^2+B.^2)./(m.*((sind(gama)).^2));
F1=[x(2)];
F2=[k*(x(2)^2)-g];
F=[F1;F2];
ACCEL=[F2];
end
$$$$$$$$$$$$$$$$$$$$$
function [FF]=ssm_module(t,x);
data_module;
global W_pl S_bar_pl S_bar_c b e r_o_enu g;
F1=zeros(500,1); F2=zeros(500,1);
alt=r_o_enu(3)+x(1); [rho] = rho_module(alt);
AR=(b.^2)./S_bar_c; m=W_pl./g; t_d=5;
A=.1.*rho.*S_bar_c.*t_d;
B=(rho./2).*(S_bar_c.*((.01+(.19./5).*t_d)+(.2.*t_d).^2./(pi.*e.*AR))+S_bar_pl.*2);
gama=-atan(B./A).*(180./pi);
k=sqrt(A.^2+B.^2)./(m.*((sind(gama)).^2));
FF1=[x(2)];
FF2=[k*(x(2)^2)-g];
FF=[FF1;FF2];
end
$$$$$$$$$$$$$$$$
function data_module; % The Data Module contains all of the fixed parameters for all the other
modules. Sets sets these parameters to Global Variables.
global r W_pl S_bar_pl Cd_pl Cs_pl S_bar_c b Cd_c_o Cl_c_o Cd_c_i Cl_c_i e dt_infl lat_o long_o r_o_enu v_ac_atm g;
r=20925670; % (f), radius of the earth
W_pl=100; % (#f), weight of the payload
S_bar_pl=3; % (f^2), characteristic area of payload
Cd_pl=2; % Payload Drag Coefficient
Cs_pl=.5; % Payload Side Force Coefficient
S_bar_c=20; % (f^2), characteristic area of canopy
b=10; % (f), canopy wing span
Cd_c_o=.01; % Canopy Drag Coefficient, Initial
Cl_c_o=0; % Canopy Lift Coefficinet, Initial
Cd_c_i=.2; % Canopy Drag Coefficient, Inflated
Cl_c_i=1; % Canopy Lift Coefficient, Inflated
e=1; % Oswald Efficiency Factor
dt_infl=5; % (s), Inflation Time
lat_o=44.8831; % (degrees),Initial Latitude (N)
long_o=-68.6719; % (degrees), Initial Longitude (W)
r_o_enu=[0 0 15000]; % (f), Initial position vector (from aircraft)
v_ac_atm=[200 100 0]; % (f/s), Aircraft Initial velocity
g=32.174; % (f/s^2), earth gravitation constant
end
$$$$$$$$$$$$$$$$$$$$$$$$$$$$
function [rho] = rho_module(alt); % (slug/f^3) Converts geometric to geopotential altitude and then calculates density. Returns density and geometric altitude.
global r g;
temp_sl=518.69; % (degR), sea level temperature
rho_sl=.002377; % (slug/f^3), sea level density
a=-.00356; % (degR/f), standard temperature lapse rate
R=1716.49; % (f #f/slug degR)
% Converts geometric to geopotential altitude
alt_gp = (r./(r+alt)).*alt; % (ft), geopotential altitude
% Calculates density using geopotential altitude
temp=temp_sl+a.*alt_gp; % Calculates standard temperature
rho=rho_sl.*((temp./temp_sl).^(-((g./(a.*R))+1))); % Calculates standard density
end
Duncan McIntosh
Duncan McIntosh on 11 Dec 2016
May not be the best approach but I am looking at writing the matrix to Excel from within usm_module and ssm_module and then reading the Excel files from solution_module.

Sign in to comment.


Walter Roberson
Walter Roberson on 12 Dec 2016
Edited: Walter Roberson on 12 Dec 2016

Tags

Community Treasure Hunt

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

Start Hunting!