In non-impulsive orbital maneuvers the thrust acts over a significant period of time. In this case we can’t ignore the change in position vector. Given the initial state of the spacecraft we calculate the state after the burn
clc; clear all; global g0 mu T Isp M0; mu = 398600; %[km^3/s^2] Earth’s gravitational parameter M0 = 4000; %[kg] Gross mass before ignition R = [ 7200, 0, 0]; %[km] Position vector: V = [ 0, 7.6, 0]; %[km/s] Velocity vector g0 = 9.81; %[km/s^2] T = 20; %[kN] Thrust, const t_b = 200; %[s] Burn time m = 1000; %[kg] Mass after burn dm = (M0 - m)/t_b; Isp = T/(dm*g0); tspan = [0,t_b]; % Constructing initial state vector y0 = [R(1),R(2),R(3),V(1),V(2),V(3),M0]; % Using ODE45 to solve ordinary differential equations of motion % for the given initial values [time,f_y] = ode45(@f_yt,tspan,y0); fs = f_y(end,:); Rf = [fs(1),fs(2),fs(3)]; Vf = [fs(4),fs(5),fs(6)]; Mf = fs(7); fprintf('Position vector after burn R = [%4.2f %4.2f %4.2f] km\n',Rf); fprintf('Velocity vector after burn V = [%4.4f %4.2f %4.2f] km/s\n',Vf); fprintf('Mass after burn M = %4.2f kg\n',Mf);
Position vector after burn R = [7035.78 1651.64 0.00] km Velocity vector after burn V = [-1.7367 9.26 0.00] km/s Mass after burn M = 1000.00 kg
function dfdt = f_yt(t,y) global mu g0 T Isp; r = (y(1)^2+y(2)^2+y(3)^2)^0.5; v = (y(4)^2+y(5)^2+y(6)^2)^0.5; m = y(7); dfdt = [y(4),y(5),y(6),-mu*y(1)/r^3 + T/m*y(4)/v,... -mu*y(2)/r^3 + T/m*y(5)/v, -mu*y(3)/r^3 + T/m*y(6)/v,... -T/(Isp*g0)]'; return