kill(all)$

load(lrats)$

single stage rocket,version 1.0 (gravity = -g)
peter.vlasschaert@gmail.com,17/07/2018.

depends([M,v,u],t);

depends(c,t);

p1:v-u=c;

u : velocity of the fuel 
v : velocity of the rocket 
relative to rocket  velocity : 'v-u' = c

p2:F=M*('diff(v,t,1))+('diff(M,t,1))*(v-u);

F : external force.

p3:subst([p1],p2);

Thrust term = T_h = ('diff(M,t,1))*c

M :total mass of the rocket
======================
P   : (p)ayload : satellite,nuclear warhead.
Mf : Mf(t),f  =fuel
Mc: Mc       = casing

depends(Mf,t);

p4:M=P+Mf++Mc;

rem : Mf only time dependent variable of M.

p5:Mf+Mc = Mo;

p6   :Mf(0)=ε*Mo;
p6a :Mc    = (1-ε)*Mo;

parameter  ε  is ratio = initial fuel supply/initialrocket mass
   =>  0 < ε < 1
pratical values : ε   ' 0.7 < ε < 0.9'

p7:1-ε;

1-ε : structural factor.

p8:lhs(p4)=subst([p6a],rhs(p4));

assume : fuel burnt and ejected = at a constant rate = k
=>  dM/dt = - k

p9:'diff(lhs(p8),t) = diff(rhs(p8),t);

p9a:lhs(p9) = (-1)*k;

p9b:subst([p9a],p9);

How to find , function : Mf=Mf(t)  ?
=============================

need to solve differential equation : 'p9b'
--------------------------------------------------------------
--------------------------------------------------------------

p10:ode2(p9b,Mf,t);

find %c: ?   :  t=0. use : initial fuel from 'p6'
----------------------------------------------------------------

p10a:subst([t=0,Mf=rhs(p6)],p10);

p10b:subst([rhs(p10a)=lhs(p10a)],p10);

after time      :   t= τ ,all fuel burn-out.  Mf = 0.
-----------------------------------------------------------------

p10c:subst([t=τ,Mf=0],p10b);

p10d:solve(p10c,τ)[1];

rocket mass  M(t) :  0 ≤ t ≤τ   ( = before time τ)
=====================================

p11:subst([p10b],p8);

p11a:expand(p11);

rem : no gravitational term (first approximation) 
           F = 0, into  'p3'

assume :  F = 0
=============
=============

p12:subst([F=0],p3);

1e)  dm/dt =-k

p12a:subst([p9a],p12);

2e) M=-k*t+P+Mo

p12b:subst([p11a],p12a);

How to solve : differential equation 'p12b',  v0≤v≤v
==========================================
rem : t=0  =>  v=v0 ( initial condition of  'ode ' )
----------------------------------------------------------------------------
time t:             0≤t≤τ
----------------------------------------------------------------------------

p13:ode2(p12b,v,t);

p13a:subst([t=0,v=v0],p13);

p13b:solve(p13a,%c);

p13c:subst(p13b,p13);

p13d:lhs(p13c)=logcontract(rhs(p13c));

increment of velocity :       Δv = v-v0,all full burn-out
===========================================

rem : payload ( = with payload )
************************************************

1e)  Δv
---------

p14:p13d-v0;

p14a:subst([v-v0=Δv],p14);

2e)  t=τ, see  'p10d'
----------------------------

p14b:subst([t=rhs(p10d)],p14a);

p14b1:ratsimp(p14b);

rem : payload ( = no payload = P=0)
*****************************************************

p14c:subst([P=0],p14b);

p14d:ratsimp(p14c);

assume :  F = - M*g ( = effect of gravity,downwards) ' negative against y direction'
====================================================================
====================================================================

p15:subst([F=-M*g],p3);

1e)  dM/dt = - k

p15a:subst([p9a],p15);

2e)  M=-k*t+P+Mo

p15b:subst([p11a],p15a);

g =constant acceleration per unit mass.

How to solve : differential equation 'p15b',  v0≤v≤v
==========================================
rem : t=0  =>  v=v0 ( initial condition of  'ode ' )
----------------------------------------------------------------------------
time t:             0≤t≤τ
----------------------------------------------------------------------------

p15c:ode2(p15b,v,t);

p15d:subst([t=0,v=v0],p15c);

p15e:solve(p15d,%c);

p15f:subst(p15e,p15c);

increment of velocity :       Δv = v-v0,all full burn-out
===========================================

p16:p15f-v0;

p16a:subst([lhs(p16)=Δv],p16);

p16b:logcontract(p16a);

p16c:ratsimp(p16b);

at full burn-out : t=τ

p16d:subst([t=rhs(p10d)],p16c);

moving one direction :  y direction
=============================

Δv = dy/dt, general t.

depends(y,t);

p16e:subst([Δv='diff(y,t)],p16c);

solve : ode 'p16e'
==============

p17:ode2(p16e,y,t);

?  How to find : t=0, y=y0

p17a:subst([t=0,y=y0],p17);

rem : solve  'p17a' for  %c
---------------------------------------

p17b:solve(p17a,%c);

p17c:subst(p17b,p17);

assume : y0 =0
-----------------------

p17d:subst([y0=0],p17c);

the height reached : h , when full -burnout  (y=h , t= τ =>'p10d')

p18:subst([y=h,t=rhs(p10d)],p17d);

p18aa:logcontract(p18);

max height =h  at burn-out  ( t = τ)
============================
============================

assume :  P <> 0 ( payload)
-----------------------------------------

p18a:ratsimp(p18);

p18b:logcontract(p18a);

assume :  P = 0 ( payload)
---------------------------------------

p18c:subst([P=0],p18b);

p18d:logcontract(p18c);

rem : magnitude of the thrust = c*k.
           -----------------------------------------------------------------
           I = specific impulse  = (c.k)/(g*k), 'seconden'
           -----------------------------------------------------------------

========================================
result :  I = c/g = 'specific impulse'
========================================

<==================================================================================================================>