kill(all)$

load(lrats)$

rem : load the 'packages' to use :     lratsubst([,,],equation)
l = list,rat =rational,subst = substitute.

simple model ,powered flight with Drag, 2D  model, 'version 1.0'
peter.vlasschaert@gmail.com,19/06/2018

filename : simplemodelofpoweredflightwithoutdragfinalversion1

equation : use also basic equation , regular flight.
first : approximation, (lift = 0 )
rem : are the forces on the plane ,see fig 1 .

?  D = drag :against the direction of motion .

vec(T)  = Thrust of powered flight,angle ' θ           ' with pos ' x axis'. (2 components)
vec(D)  = Drag                                    ,angle ' θ+%pi ' with pos 'x axis '. (2 components)
vec(G)  = gravity  ' - ', y direction                                                        (2 components)
vec(V)  =  angle 'γ' with pos ' x axis'                                                 (2 components)

rem : θ = pitch angle   (x : component  = T*cos(θ)=T_x,y:component = T*sin(θ)=T_y)
                                         (x : component  = -D*cos(θ)=D_x,y:component = -D*sin(θ)=D_y)
           V = 'velocity'  => (x: component  = V*cos(γ)=V_x,y:components=v*sin(γ)=V_y)
           G = gravity     => ( x:component  = 0,                           y:components = - g             )

rem : θ ≠ γ
v = ' velocity ' ,unit :  m/s
g = ' gravity  ' ,unit : m^2/s

Newton Second law :  F=m*a = m*dv/dt= sum of all external forces.
external forces  = {T,m*g} =>                                         , unit : kg*m^2/s

assume(m>0);

p1:m*'diff(vec(V),t) = vec(T)+vec(D)+m*vec(G);

p2:solve(p1,'diff(vec(V),t))[1];

vectorial : differential equation of motion ' with Thrust  '
=============================================

p3:ratexpand(p2);

differential equations 'p3' will result '   "two" differential equations: from the components: 2D'

depends([V_x,V_y],t);

p4:vec(V)=[V_x,V_y];
p5:vec(T)=[T_x,T_y];
p5a:vec(D)=[D_x,D_y];
p6:vec(G)=[G_x,G_y];

x: components  'differential equation'
----------------------------------------------------------

p7:[lhs(p4)=rhs(p4)[1],lhs(p5)=rhs(p5)[1],lhs(p6)=rhs(p6)[1],lhs(p5a)=rhs(p5a)[1]];

p8:lratsubst(p7,p3);

p9:ratexpand(p8);

trignometric :  'see above'

p10:T_x = T*cos(θ);
p10a:D_x = D*cos(θ+%pi);
p11:G_x = 0;

p12:lratsubst([p10,p10a,p11],p9);

y: components  'differential equation'
----------------------------------------------------------

p7a:[lhs(p4)=rhs(p4)[2],lhs(p5)=rhs(p5)[2],lhs(p6)=rhs(p6)[2],lhs(p5a)=rhs(p5a)[2]];

p8a:lratsubst(p7a,p3);

p9a:ratexpand(p8a);

g = constant  = 9.81 m^2/s 

assume(g>0);

trignometric :  'see above'

p10a:T_y = T*sin(θ);
p10b:D_y= D* sin(θ+%pi);
p11a:G_y = -g;

p12a:lratsubst([p10a,p10b,p11a],p9a);

p13a:ratexpand(p12a);

result  : two differential equations 'x,y' ,components.

p14:p12;
p14a:p13a;

assume :θ - γ, differences between two angles : very small value.
======================================================
x: component :cos(θ) = cos(γ)= V_x/V
y: component :sin(θ) = sin(γ)  = V_y/V

1e)  x-component :cos(θ) = cos(γ)
2e)  y-component :sin(θ) =sin(γ)

gamma = γ, equations.

p15:lratsubst([cos(θ) = cos(γ)],p14)$
p16:factor(expand(p15));
p15a:lratsubst([sin(θ) =sin(γ)],p14a)$
p16a:expand(p15a);

1e)  x-component : cos(γ) = V_x/V
2e)  y-component : sin(γ)  = V_y/V

p17:[cos(γ) = V_x/V];
p17a:[sin(γ) = V_y/V];

result  : two differential equations 'x,y' ,components.(assume: θ  = γ)

p18:lratsubst(p17,p16)$
p19:factor(ratexpand(p18));
p18a:lratsubst(p17a,p16a)$
p19a:ratexpand(p18a);

?  V 

a) 1e)   V_x  from   'p17'
    2e)   V_y  from   'p17a'

p20 :solve(p17[1],V_x)[1];
p20a:solve(p17a[1],V_y)[1];

b)  solve 'p20,p20a' for V

p21:(p20)^2;
p21a:(p20a)^2;

p22:factor(p21+p21a);

p23:trigreduce(p22);

p24:solve(p23,V)[2];

two differential equations 'x,y','p19,p19a' respect '*V_x,*V_y'

p25:factor(expand((V_x*p19)));
p25a:expand(V_y*p19a);

ℏ;

from : 'p26,p26a' or 'p27'  insert into 'p25,p25a'

p26:diff(V_x^2,t);
p26a:diff(V_y^2,t);

p27:'diff((V_x^2+V_y^2),t)=diff((V_x^2+V_y^2),t);

result  : two differential equations 'x,y' , V inside 'p25,p25a'.

p28:2*p25;
p28a:expand(2*p25a);

p29:p28+p28a;

p30:lratsubst([rhs(p27)=lhs(p27)],p29);

depends([V],t);

p31:lratsubst([p23],p30);

p32:ev(p31,diff);

p33:expand(p32);

p34:solve(p33,diff(V,t));

p35:expand(p34)[1];

trignometric : function : V_y/V = sin(γ)

1e)  insert  'sin(γ)' into  'p20a'

p36p:solve(p20a,sin(γ))[1];
p36q:solve(p20,cos(γ))[1];

2e) insert 'sin(γ)' into 'p35'

p36a:lratsubst([rhs(p36p)=lhs(p36p),rhs(p36q)=lhs(p36q)],p35);

We get : first differential equation .

p36b:expand(p36a);

We need to find : second differential equation. ?
    =>  diff(γ,t )  ?

1e) step :      p19*V_y,p19a*V_x

p37:factor(p19*V_y);
p37a:p19a*V_x$
p37b:expand(p37a);

2e) step :  subtract 'p37b-p37'

p38:p37b-p37;

3e) step : 2e/V

p39:expand(p38/V);

4e) step : V_x/V = cos(γ) subst in 'p39' into rhs

p40:lhs(p39)=factor(lratsubst([lhs(p20)/V=rhs(p20)/V],rhs(p39)));

p40a:lhs(p40) =expand(rhs(p40));

5e) find  lhs(p40)  =  ?  dγ/dt

a) tan(γ)  ? by definition

p41:sin(γ)/cos(γ) = tan(γ);

b) V_y,V_x  , tan(γ);

p42:p20a/p20;

c) combine a) and b)

p43:lhs(p42) =rhs(p41);

d) solve 'p43' for γ

p44:solve(p43,γ)[1];

e) differentiation 'p44' with respect  't'

p45:'diff(lhs(p44),t)=diff(rhs(p44),t);

f) simplification

a) radcan ' maxima'  see : menu

p46:lhs(p45)=radcan(rhs(p45));

b)  use 'p23' into 'p46'

p47:V*lratsubst([p23],p46);

g) observe  'p47' and 'p40'

p48:lhs(p47) = expand(rhs(p40a));

simplify ,two differential equations 
=============================

p49:p36b;
p50:expand(p48)/V;

rem: D can be function .  
   D= C_D*1/2*ρ*V^2*A_f
C_D =Drag Coef.
ρ       =Density
V       =Velocity
A_f   = frontal cross sectional area of vehicle

p52:cos(γ+%pi);
p52a:trigreduce(p52);

p53:sin(γ+%pi);
p53a:trigreduce(p53);

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