Hi @moinmoinnoob69, Hoping you are practic with ordinary differential equations, I think that it is appropriate to approach the script with the two body problem structure with a perturbation ,which in this case is the thrust. I will assume that you have the thrust curve of each of the phases of your project , from which you will compute the acceleration dividing the thrust by the mass, each instant. You are lucky because I am writing a thesis about a spacecraft which orbits around Mercury, and also if it seems very different, the structure of the script and the functions are not:
This is the main script
clc,clear
format long
mu_m = ...
s_0 = ...
%%%%%%%%%% Integrazione del de-orbiting
options = odeset('RelTol', 1.e-16, 'AbsTol', 1.e-16);
[t,S_1] = ode113(@ODE_1_LanderOrbit, t_span1 , s_0' , options) ;
%%%%%%%%%% Integrazione della fase propulsa a solido e coasting %%%%%%%%%%
[t,m_2] = ode113(@ODE_SRMSupport, 0:0.01:(t_f2-t_i2) , 3490 ) ;
t_span_SRM_int = 0:2.5:65 ;
F_SRM_int = 4.44822*[0 17500 18000 20000 22000 23500 26000 25000 26000 26500 27500 28500 30000 31000 31000 32000 32500 34000 34000 34000 34000 32500 32000 32000 31000 15000 0 ]; % La spinta in libbre , bha !
F_SRM_eva = interp1(t_span_SRM_int,F_SRM_int,0:0.01:(t_f2-t_i2),'linear','extrap') ;
for i2 = 1:length(0:0.01:(t_f2-t_i2))
a_SRM(i2)= F_SRM_eva(i2)/m_2(i2) ;
end
[t, S_2] = ode113(@(t, x) ODE_2_LanderSRM( t, x, a_SRM , t_span2 ), t_span2, S_1(end, 1:6));
[t, S_3] = ode113(@ODE_3_LanderCoasting, t_span3, [S_2(end, 1:6) ]);
%%%%%%%%%% Integrazione della fase propulsa a liquido %%%%%%%%%%
[t,m_3] = ode113(@ODE_LRESupport, 0:0.01:(t_f4-t_i4) , m_2(end)-560*0.453592 ) ;
t_span_LRE_int = [0 0.1 0.5 0.7 5 7.5 12 13 13.5 15 18 19 19.001] ;
F_LRE_int = 1.03*4.44822*[0 1400 1450 1400 2100 2400 2800 2825 2850 2850 2850 1000 0 ];;
F_LRE_eva = interp1(t_span_LRE_int,F_LRE_int,0:0.01:(t_f4-t_i4),'linear','extrap');
for i_2 = 1:length(0:0.01:(t_f4-t_i4))
a_LRE(i_2)= F_LRE_eva(i_2)/m_3(i_2);
end
[t, S_4] = ode113(@(t, x) ODE_4_LanderLRE(t, x, a_LRE,t_span4), t_span4, S_3(end, 1:6));
Here it is an example of the function that will you provide the mass each instant
function dA = ODE_SRMSupport(t,A1)
g_0 = 3.7;
m_p = 6205 * 0.453592 ;
I_t = 4.44822 * 1911070 ;
x= 0:2.5:65 ;
T = 4.44822*[0 17500 18000 20000 22000 23500 26000 25000 26000 26500 27500 28500 30000 31000 31000 32000 32500 34000 34000 34000 34000 32500 32000 32000 31000 15000 0 ]; % La spinta in libbre , bha !
i = interp1(x,T,t,'linear','extrap');
A5 = i*m_p/I_t ;
dA = [ -A5 ]' ;
end
And here an example of the SRM stage function, while the orbiting one can be easily provided removing the a_p. (the LRE is obviously very similar to the SRM one)
function dx = ODE_SRMSupport(t, x, ts,t_span_f)
istat = interp1(t_span_f,ts,t,'linear','extrap');
mu_m = ...
grav = - mu_m / (norm(x(1:3),2))^3;
versore = x(4:6)/norm(x(4:6)) ;
a_p = versore * istat / 1000;
A = [ 0,0,0, 1,0,0
0,0,0, 0,1,0
0,0,0, 0,0,1
grav,0,0, 0,0,0
0,grav,0, 0,0,0
0,0,grav, 0,0,0];
dx = A*x+[0 0 0 -a_p']';
Obviously it will not run because i rudimentally modificated it , deleting the lines that are not necessary for your problem, but I think it is a very good starting point (I also didn't delete my parameters, I am sorry XD). Hope it will help! I am sorry if i will be not very active to answer your questions this week but I am new to this forum and I only got registered because I had an issue about the thesis I am writing. Therefore, I will do my best to solve other problems you will propose. I am not familiar with some good literature about what you have to do, I am sorry about that.
I didn't add the perturbation of the drag (because Mercury doesen't have an atmosphere hahaha), but it can be easily computed as
a_drag = - velocity_versor * (module of the drag, function of the position of the rocket, provided by the integrator itself)
I am sorry for my bad language but I am italian :) .
Good luck!
