How to convert my 2d trajectory code into a 3d by banking

Question

Jose Camorlinga Medina 2018-11-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/429928-how-to-convert-my-2d-trajectory-code-into-a-3d-by-banking

I used the equations of motion to model a trajecroy of a spaceplan, the first 2 phases of the plane were kept in 2d for simplicity. In the third phase it is supposed to do a bank back maneuver and glide back to the takeoff site. I am not sure how to model the bank maneuver.

global Re ht go deg rhoo hs msr mpir mpL mor FTr Ispr dr Ar CDr mdr mfr nr me mpi mo FTv Ispv md d A CD mf n tb tbr k CL Aw
%% Constant
% General constants
Re = 6378000;                   % Radius of earth (m) 
ht = 160;                       % Height where gravity turn begins (m)
go = 9.81;                      % Sea-level acceleration due to gravity (m/s^2)
deg = pi/180;                   % Degree to radian convertion (rad)
rhoo = 1.225;                   % Sea-level pressure (Kg/m^3)                          
hs = 7.5e3;
r = 1.4;
R_o= 287.053; % (J/Kg/K)
k = 398600*1000^3;                     % Gravitational constant (km^3/s^2)
% Rocket Parameters
msr = 2400;                     % Mass of emtpty rocket structure (kg)
mpir = 3000;                     % Initial mass of rocket propellent (kg)
mpL = 1361;                     % Payload mass (kg)
mor = msr+mpir+mpL;             % Inital rocket mass (kg)
FTr = 50000;                    % Vacume Thrust rocket (N)
Ispr = 260;                     % Vacume specific impulse rocket (s)
mdr = FTr/(Ispr*go);            % Mass flow rate rocket (kg/s)
dr = 1;                         % Frontal Diameter rocket (m)
Ar = pi/4*dr^2;                 % Frontal area of rocket (m^2)
mfr = mor-mpir;                 % Burnout mass (kg)
nr = mor/mfr;                   % Mass ratio rocket ()
CDr = 0.1;                      % Drag coefficient ()
% Spaceplane Parameters
me = 63000;                    % Mass of empty structure (kg)
mpi = 78000;                    % Initial mass of propellent (kg)
mo = me+mpi+mor;             % Initial mass of launch vehicle
FTv = 2*10^6;                  % Vacume Thrust (N)
Ispv = 455;                     % Vacume specific impulse (s)
md = FTv/(Ispv*go);             % Mass flow rate (kg/s)
d = 5;                          % Frontal area of (m^2)
A = pi/4*d^2;                   % Frontal area of spaceplane (m^2)
mf = mo-mpi;                    % Mass of fuel (kg)
n = mo/mf;                      % Initial to final mass ratio ()
CD=0.5;                         % Drag coefficient rocket()
CL=1;
Aw=145;
%% Phase 1 (Ascent-Separation)
% time
to = 0;                         % Initial time (s)
tb = mpi/md;                    % Burnout time (s)
tspan = [to,tb];                % Time interval vector
% Initial condition
vo = 0;                         % velocity (m/s)
gamo = 89.85*deg;               % flight path angle (rad)
xo = 0;                         % downrange distance (km)
zo = 0;                         % altitude (km)
yo = [vo; gamo; xo; zo];                    % Initial condition vector
[t,y] = ode45(@EOMP1,tspan,yo);                       % First phase ode45 command
% Numerically solved solution tspan
v = y(:,1);               % Velocity (km/s)
gam = y(:,2);             % Flight path angle (degrees)
x = y(:,3);               % Downrange distance (km)
z = y(:,4);               % Altitude (km)
Z = z/1e3;
%% Phase 2 (Separation-Payload orbital injection) 
% Time 
tbr = tb+mpir/mdr;                 % Burn time (s)
tor = tb;                        % Initial time for the numerical integration
tspanr = [tor,tbr];             % Range of integration
% Initial conditions
vor = v(end);                   % velocity (m/s)
gamor = gam(end);               % flight path angle (rad)
xor = x(end);                   % downrange distance (m)
zor = z(end);                   % altitude (m)
yor = [vor; gamor; xor; zor];
[tr,yr] = ode45(@EOMP2, tspanr, yor);
vr = yr(:,1);                    % Velocity (m/s)
gamr = yr(:,2);                  % Flight path angle (rad)
xr = yr(:,3);                    % Downrange distance (m)
zr = yr(:,4);                    % Altitude (m)