kill(all)$
load(lrats)$

sounding rocket (=vertical path) ,version 1.0
peter.vlasschaert@gmail.com, 13/08/2018

sounding rocket : ?
=================
 =>            https://en.wikipedia.org/wiki/Sounding_rocket

Show forces on a sounding rocket.

Differential equations : sounding rocket
------------------------------------------------------------
m[0] : initial mass  (t=0)
m      :            mass  (t =t)
m,m[0] >0, m[0] > m.
------------------------------------------------------------

--------------------------------------------------------------------------------------------------------------------

p1:'diff(y,t) = v;
p2:m*'diff(v,t) = -m*g-'diff(m,t)*V[e];

rem  :  minus sign before V[e] :  mass decrease. (lost of mass during the flight)

---------------------------------------------------------------------------------------------------------------------

find v: ( = vertical speed)
=====

p2a:expand(p2/m);

math notation :      dv/dt         = -V[e]*(1/dt)*dm/m -g
                                   dv/dt*dt    = -V[e]*dt/dt*dm/m -g*dt
                                   dv              = -V[e]*dm/m -g*dt

see : two ODEs
-----------------------

assume(m[0] > 0,m > 0,m[0] > m);

p2b: 'integrate(1,v,0,v) = -V[e]*'integrate(1/m,m,m[0],m)-'integrate(g,t,0,t);

p2c:ev(p2b,nouns);

p2d:logcontract(p2c);

thrust = constant    = >  mass flow rate : dm/dt = k  also constant.
==============                                            ====================

p3:'diff(m,t) = k;

find m: (from 'p3')
======

p4:ode2(p3,m,t);

from : differential equations 'see above '

find : %c
--------------

t=0  =>  m=m[0]  ( = initial conditions differential equation 'p3')
-----------------------

p4a:subst([t=0,m=m[0]],p4);

p4b:subst([rhs(p4a)=lhs(p4a)],p4);

solve : 'p4b' to 't'
--------------------------

p4c:solve(p4b,t);

find y: ( =distance)
=====

use :  'p1' and 'p2d'

p5:lhs(p1) = rhs(p2d);

math notation :      dy/dt         = -V[e]*ln(m[0]/m) -g*t               (1e)
                                   dy/dt*dt    = -V[e]*ln(m[0]/m)*dt -g*t*dt     (2e)
                                   dv              = -V[e]*ln(m[0]/m)*dt -g*t*dt     (3e)
=>    dt = dm/k  , math = ln  => maxima = log

p5a:'integrate(1,y) = -V[e]*'integrate(log(m[0]/m)/k,m,a1,a2)-'integrate(g*t,t,0,t);

condition : m[0] > m , m>0,m[0] > 0

assume( a2< a1,a2>0);

p5b:ev(p5a,nouns);

p5c:expand(p5b);

rem : manuel integration ' integration by parts '
=======================================
------------------------------------------------------------------------------------
rule : d(f*g)=g*df+f*dg => d(f*g)-f*dg =g*df =>  ' integrate '
------------------------------------------------------------------------------------

p5d:'integrate(log(x),x) =integrate(log(x),x)+C;

------------------------------------------------------------------------------------

input :   a1 = m[0] ,lowest boundary of integration
               a2  = m     ,highest boundary of integration.

p5e:subst([a1=m[0],a2=m],p5c);

rem : by notation : (m[0])[0] =m[0] )
-----------------------------------------------------

p5f:[  (m[0])[0] =m[0] ];

p5g:subst(p5f,p5e);

y = y(m,m[0]) ,with  use ' p4c ' = t
-------------------------------------------------

p5h:subst(p4c,p5g);

find v[bo]  = velocity at burnout ,  m=m[f] , f = fuel.
=========================================

before burnout:
=============

use : 'p2d', use of notation ' see p5f ', 'p4c'

p6:[p4c[1],v=v[bo],m=m[f],(m[f])[0] =m[0] ];
p6a:subst(p6,p2d);

after burnout :the mass of rocket is constant  (m=m[f]), ( = free flight)
=========================================================

free flight  ? 
==========

without  D = Drag  and T = Thrust term  => Energy is constant.= E[i] , i = 1,2.
----------------------------------------------------------------------------------------------------------------
suppose : begin state  , i=1.
                     end    state  , i=2
----------------------------------------------------------------------------------------------------------------

definition E[i]
===========

p7:E[i]  := 1/2*m[i]*v[i]^2 + m[i]*g *y[i];

p7a: E[1] =E[2];

---------------------------------------------------------------------------------------------------------------

after burnout  
===========

p7b: Energy = subst([v[1]=v[bo],y[1]=y[bo],m[1]=m[f]],E[1]);

free flight : g=0 , into ' p2d'  and  ' p5h '.
================================

p8:y[bo] = subst([m=m[f],(m[f])[0]=m[0]],rhs(p5h))$
p8a :subst([g=0],p8);

p9:subst([g=0],p6a);

result : Energy
============

E[i], i=1  ( = begin state)
----------

p10:lhs(p7b) = subst([p8a,p9],rhs(p7b));

find  y[max ] :
===========

E[i] , i=2   ( = end state)
------------

y[2] = y[max],v[2] = 0,m[2] = m[f].

p11:lhs(p7b) = subst([y[2] = y[max],v[2]=0,m[2]=m[f]],E[2]);

 result  y[max] ?
==============

p12:rhs(p10) = rhs(p11);

p12a:solve(p12,y[max]);

p12b:expand(p12a)[1];

p12c:logcontract(p12b);

rem :maximum altitude (= y[max]) occurs when speed  goes to zero ( i= 2)
          ==========================================================
          k     = mass flow  rate  ( = dm/dt) = constant
         V[e] = exhaust speed
         m[f]  = constant 
         m[0] = constant
         g      = 9.81 m/s^2
        ===========================================================

-----------------------------------------------------------------------------------------------------------------------------------------