kill(all)$
load(lrats)$

equation motion rocket 2D (flight,θ,ψ,simple),version 1.0
peter.vlasschaert@gmail.com, 19/08/2018.

 Forces on the rocket

design : Geogebra 5. ( open source).

assume(g>0,m>0,ρ>0);

forces   =>   direction of flight                  : v
-------------------------------------------------------------
Thrust  : F*cos(ψ-θ)  ,  'w1'
Drag     : -D             = -1/2 C[d]*ρ*v^2*A[fr].
Weight : -G*sin(θ) = - m*g*sin(θ)
rem :  D opposit direction of  'direction of flight'.
rem : (-) , neg direction       of  'direction of flight'.

A[fr] = frontal area
ρ       = density of the medium .
v       = velocity

p1:w1=F*cos(ψ-θ);
p2: D = 1/2*C[d]*ρ*v^2*A[fr];
p3: G  = m.g;
p4:g1=-1*m*g*sin(θ);

differential equations : 'v'
--------------------------------------

p5:m*'diff(v,t) = rhs(p1)-lhs(p2)+rhs(p4);

p6: expand(p5/m);

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

forces   =>   normal 'direction of flight  ': L
--------------------------------------------------------------
Thrust  : -F*sin(ψ-θ),                               'w2'
Lift         : L               = 1/2 C[l]*ρ*v^2*A[fr].
Weight : -G*cos(θ) = - m*g*cos(θ)
rem : (-) , θ change pos direction of ω 'Thrust'

p1a :w2=F*sin(ψ-θ);;
p2a :L= 1/2 *C[l]*ρ*v^2*A[fr];
p4a: g2 = -1*m*g*cos(θ);

differential equations                              : 'L'
-------------------------------------------------------------

acceleration : v^2/R ( pos direction of  L).
R = radius of curvature .( circle :the radius of the circle )

p5a : v= ω*R;
p6a :  ω = 'diff(θ,t);
p7a : subst([p6a],p5a);

find : x (use 'p7a' )
===============

p8a: v^2/R = v*x;
p9a:solve(p8a,x)[1];

p10a:lhs(p9a)=rhs(p7a)/R;

p11a : subst([p10a],p8a);

forces : m*v^2/R = m*v*('diff(θ,t,1))
============================

p12a : m*rhs(p11a) = -rhs(p1a)+L+rhs(p4a);

p13a:expand(p12a/m);

------------------------------------------------------------------------------------------------
two ode's : 'v' ,'L' ( use  'p2',' p2a' )
=============================
-------------------------------------------------------------------------------------------------

p14:subst([p2],p6);
p14a:subst([p2a],p13a);

assume constant =  thrust . ( F ) ,    'p14,p14a' . (one stage rocket)
======================================================
--------------------------------------------------------------------------------------------------------------------------
rem :  F =constant  >  0,  because = v[e] : exhaust gasses velocity (direction neg)
--------------------------------------------------------------------------------------------------------------------------

p15: F = - 'diff(m,t)*v[e];

p151: solve(p15,('diff(m,t,1)))[1];

p152:[rhs(p151) =lhs(p151)];

solve : ode2, ' p15'

p15a:ode2(p15,m,t);

find : %c  ?
==============
t = 0  =>  m = mi

p15b:subst([t=0,m=mi ],p15a);

p15c: subst([rhs(p15b)=lhs(p15b)],p15a);

p15d:lratsubst(p152,p15c);

We can set : ('diff(m,t,1)) = -k , k = constant
===================================
why : initial mass = m[i] goes down after burning fuel . (reason : neg sign)

p15e:'diff(m,t,1)=-k;
p15f:subst([p15e],p15d);

what value from :  k =  ?
====================
k=burning constant = m[fuell]/t[fb].
fb : t[fb],    when all fuel is burning (lost mass = m[fuel])

p16:k = m[fuel]/t[fb];

p16a:subst([p16],p15f);

find : F = F(m), forces function of the mass.
====================================

p17:subst([p16],p15e);

p17a:subst([p17],p15);

result : 2 ODES . ' p14,p14a' use of 'p16a,p17a'
==============

p18:subst([p16a,p17a],p14);
p18a:subst([p16a,p17a],p14a);

If  θ = ψ then  ?  2 ODE's
====================

p19:[ψ = θ];
p19a:subst(p19,p18);
p19b:subst(p19,p18a);

rem : specefic fuel:  ? v[e] = I[spc]*g[0].
================================
=>     I[spc]  = 'time'. (different for every fuel)

p20:v[e] =I[spc]*g[0];
p20a:subst([p20],p19a);
p20b:subst([p20],p19b);

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