kill(all)$
load(lrats)$
load(vect)$

launch trajectories, (no rotation of the earth )'version 1.0a',corrected
peter.vlasschaert@gmail.com,12/09/2018

 forces acting on the rocket. (L,T,G,D)

s : N to M  :down track.
----------------------------------------------------------------------
T  =>    T[V] = T cos(α)  , T[L]  = T*sin(α) ,orange
G  =>    g[V] = -G*sin(γ) ,  g[L]  = -Cos(γ) , purple
----------------------------------------------------------------------
θ = theta
γ = gamma

see also : website 
----------------------------
Polar coordinates and Normal-Tangential coordinates ( = velocity),version 1.0

polar coordinates: (r,θ)
===================
θ      unit vectors       : e[θ]
 r     unit vectors       : e[r]
R : vector , r = radius

p1: 'diff(R,t) = 'diff(r,t)*e[r] + r*'diff(θ,t)*e[θ];

rem :  ref  x -axis : rotation trough angle of R : θ  , 'velocity'

Normal-Tangential coordinates 
==========================
Tangential unit vectors       : e[v]
 Normal      unit vectors       : e[L]

rem : ref horizontal : rotation trough angle of V:γ,'acceleration'

Why : replace : R-> V, θ -> γ,e[θ] ->e[V],e[r] -> e[L]  , 'p1'

p2: V = v*cos(γ)*e[θ] + v*sin(γ)*e[r];

rem :  V = d R/dt

p2a:rhs(p1) = rhs(p2);

coef : from  e[θ]
=============

p2b:coeff(p2a,e[θ]);

coef : from  e[r]
=============

p2c:coeff(p2a,e[r]);

see also : website 
----------------------------
rotating coordinates system,version 1.0

Why rotation of the system because : unit vectors (e[V],e[L]), change in position in space
--------------------------------------------------------------------------------------------------------------------------------------
(dV/dt)'fixed' = inertial part + noninertial part = (dV/dt)'moving'+ωxV
--------------------------------------------------------------------------------------------------------------------------------------
rem : x = 'cross product' between (ω,V)

1e) inertial part : acceleration
========================

p3:[R = V,r=v,θ = γ, e[θ] =e[V],e[r] = e[L]];
p3a:subst(p3,p1)$
p3b:subst([e[γ]=e[L]],p3a);

2e) non inertial part : acceleration      cross( ω *V) =vec[ω] x V
===================================================

figure 2: e[z] : above the purple plane ( e[z] > 0 )
                 e[V] perpendicular to e[L]
                 e[V],e[L] are in purple plane
=> e[V] perpendicular e[z] , e[L] perpendicular e[z].

a)   V , e[V]
--------

rem :  small v : direction  e[V]  : 

p4:  V= v*e[V];

b) vec[ω] ,e[z]
--------------

e[z], vector above the plane.(out the plane e[z] >0  )
                                                     (in    the plane  e[z] < 0)

p4a:vec[ω] = ω*e[z];

c) vec[ω] x V  :  e [z] x e[v]
---------------------------------------------------------
use : right hand rule for cross product
----------------------------------------------------------
result : e[z] x e[V] = - e[L]

p4b:p4*p4a;

p4c:lratsubst([e[z]*e[V] = -e[L]],p4b);

rem :  ω = 'diff(θ,t) , rotation  θ
---------------------------------------------

p4d:lhs(p4b) =subst([ω ='diff(θ,t)],rhs(p4c));

d) We want use , 'p2b' into 'p4d'
------------------------------------------------

p4e:solve(p2b,('diff(θ,t,1)))[1];

p4f:subst([p4e],p4d);

result :acceleration => Newton Law
-----------------------------------------------------
F = m*a = m*dV/dt

p5:a =(rhs(p3b)+rhs(p4f));

coef : from  e[v]
=============

p5a:coeff(rhs(p5),e[v]);

coef : from  e[L]
=============

p5b:coeff(rhs(p5),e[L]);
p5b1:part(p5b,1);
p5b2:part(p5b,2);

p5c:G =m*g;

Newton Law : ' v direction '
--------------------------------------------------------
from fig 1:  components : T[V],D, g[v] 
sign :   use : ref 
                     components ; T[V],-D,-g[v]

p6: m*p5a = T*cos(α)-D-G*sin(γ);
p6a:subst([p5c],p6);
p6b:p6a/m;
p6c:expand(p6b);

Newton Law : ' L direction '
--------------------------------------------------------
from fig 1:  components : T[L],L,g[L]
sign :   use : ref 
                     components ; T[L],L,-g[L]

p7: m*p5b = T*sin(α)+L-G*cos(γ);
p7a:subst([p5c],p7);
p7b:p7a/m;
p7c:expand(p7b);
p7d:p7c-p5b2;
p7e:p7d/v;
p7f:expand(p7e);

F = Gr*m1*m2/r^2;  Newton Law for attraction of two bodies.
------------------------------------------------------------------------------------------------------
?  find relation  g0 and g , m1 = mass (e)arth ,m2 = m ( = rocket,object)
g0 -> r[e] , g = r = r[e] +h , Gr = gravitation constant
---------------------------------------------------------------------------------------------------------

p8 :  Gr*m[e]/r[e]^2 *m = g0*m;
p8a: Gr*m[e]/r^2 *m = g*m;

p8b: p8a/p8;

p8c:solve(p8b,g)[1];

g0 : 9.81 m/s^2 down acceleration under gravitation

What about : D,L
================
ρ = density medium ( density of air)
C[L]  ,C[D] = constant 
v    = velocity
Afr = frontal area

p9: L = 1/2*ρ*v^2*Afr*C[L];
p9a:D = 1/2*ρ*v^2*Afr*C[D];

p10:[p8c,p9,p9a];

final equation : simulation ( by RK4 =Runge Kutta 4th order, ODE's)
========================================================

final result : dv/dt : ' p6c '
---------------------------
use : p8c,p9,p9a. = 'p10'

p11:subst(p10,p6c);

final result : dv/dt : ' p7f '
---------------------------
use : p8c,p9,p9a. = 'p10'

p12:subst(p10,p7f);

rem : γ = gamma

downtrack : s,along the earth,some correction: 
-------------------------------------------------------------------------
s = r*θ  => ds/dt = r*dθ/dt,                                    (+++)

p13: p2;
p13a:p2c;

p13b: s=r*θ;
p13b1: 'diff(s,t) =r*'diff(θ,t);

p13b2:lhs(p13b1)=r[e]/r*v*cos(γ);

see (+++) :  above :coef : from  e[θ], 'p2b'

change mass : dm/dt ( relation to T =thrust)

p14b:'diff(m,t) = -T/(g0*I[spc]);

negative sign (-) : mass burning during the flight
from : above , I[spc], units  ' second:s ',specefic for every propellant.
                           T        ,units   'kg*m/s^2 ','Newton'
                           g0     ,units   'm/s^2      ', acceleration nearby earth surface ( = sea level)
g0 : 9.81 m/s^2 down acceleration under gravitation

final ODEs , RK4,procedure
========================

p15:p11;
p15a:p12;
p15b:p13a;
p15c:p13b2;
p15d:p14b;

===============================================
 r =r[e]+h , h = height above the earth
 γ, angle between(local horizon,vector: V)
 α,angle between(vector          V,vector: T)
 r[e],earth radius
 ρ = density medium ( density of air)
 C[L]  ,C[D] = constant 
 v    = velocity
 Afr = frontal area
 m  = mass rocket +propellant (change during the flight)
 g0 : 9.81 m/s^2 down acceleration under gravitation
 s : along the earth surface
T :thrust
===============================================

corrected : 'p13b','p13b1','p13b2'  , version 1.0a 

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