Integrating a Multi Variable Function

I currently have a working simulink model of a rocket which I'd now like to turn into MATlab code. The simulink uses a function which calculates atmospheric conditions such as drag, density and gravity, as well as the ODEs of velocity magnitude (v), downrange distance (x), flightpath angle (gam) and altitude (y), as well as some other things like the changing mass of the rocket, and drag + gravity loss (VG,VD).

Here's the function:

function [m0_d,y_d,gam_d,x_d,v_d,Vx,Vy,q,VG,VD] = fcn(m0,y,gam,v,t)

%mass_payload = 100 ; %mass of payload (kg)
mass_fuel_struc = 2500; %mass of the rocket structure (kg)
%m0 = mass_payload + mass_fuel_struc; %total rocket mass (kg)
R_E = 6.3781*10^6; %radius of the Earth (m)
M_E = 5.9722*10^24; %mass of the Eath (kg)
G = 6.674*10^-11; %universal gravitational constant (m^3 kg^-1 s^-2)
thrust = 30200 ; %initial engine thrust (N)
isp = 292; %ISP of fuel (s)

%area, drag coefficient, gravity
Area = pi*0.35^2; cd = 0.75; g_s = 9.81;

burntime = (mass_fuel_struc*g_s*isp)/thrust; %exact burntime (s)
mdot = mass_fuel_struc/burntime; %mass fuel consumption rate (kg/s)

% rocket mass
if t <= burntime
    m0_d = -mdot;
    T = thrust;
else %rocket has run out of fuel so mass is now constant
    m0_d = 0;
    T = 0;
end

%density calculator
if y <= 11000
    Temp = 15.04-0.00649*y;
    p = 101.29*((Temp+273.1)/288.08)^5.258;
elseif y <= 25000
    Temp = -56.46;
    p = 22.65*exp(1.73-0.000157*y);
else
    Temp = -131.21 + 0.00299*y;
    p = 2.488*((Temp+273.1)/216.6)^(-11.388);
end

rho = p/(0.2869*(Temp+273.1));
% gravity model
g = (G*M_E)/(R_E+y)^2;
%drag force D
D = 0.5*rho*v^2*Area*cd;

%speed magnitude v
v_d = T/m0 - D/m0 - g*sin(gam);

%path angle gam (from 90 to 0 degrees)
if t <= 30
    gam_d = 0; %vertical flight
elseif t <= 32
    gam_d = - 0.0048*gam; %kickover manouvre
else
    gam_d = - ((g/v - v/(R_E+y))*cos(gam)); %gravity turn equation
end

%x distance downrange
x_d = R_E/(R_E+y)*v*cos(gam);

%altitude y
y_d = v*sin(gam);
VG = g*sin(gam); %gravity loss
VD = D/m0; %drag loss
%dynamic pressure
q = (rho * v^2)/2;

%Vertical Velocity
Vy = v*sin(gam);
%Horizontal Velocity
Vx = v*cos(gam);

end

How could I then use this function in a MATLAB script to integrate the same variables as in the Simulink model? I'd ideally like to use the ode45 function, running the simulation for 900s. Initial conditions of the integration blocks are:

m0_d = 2600, y_d = 0, x_d = 0, gam_d = pi/2, v_d = 0, VG = 0, VD = 0.

I appreciate any help given.

Answer 1

[Abhijeet]

Hi Hamish,

I understand that you want to integrate multiple variables mentioned in your Simulink model into MATLAB code. You can use ‘integral3’ function to numerically evaluate triple integral. For example:-

%Define the anonymous function f(x,y,z)=ysinx+zcosx.
fun = @(x,y,z) y.*sin(x)+z.*cos(x)
% Integrate over the region 0≤x≤π, 0≤y≤1, and −1≤z≤1. 
q = integral3(fun,0,pi,0,1,-1,1)

Reference link

• https://in.mathworks.com/help/symbolic/index.html?s_tid=CRUX_lftnav ( Documentation for using Symbolic Variables)
• https://www.mathworks.com/help/matlab/ref/integral3.html (Documentation on using the ‘integral3’ function)
• https://in.mathworks.com/help/symbolic/integration.html (Documentation on using integration on single variable)