kill(all)$
load(lrats)$

RE-ENTRY equations of motion 2D( assume flat earth,vertical motion).
peter.vlasschaert@gmail.com,20/08/2018.

see : RE-ENTRY equations of motion 2D( assume flat earth).
we start from : ODE1,ODE2

p1:'diff(Z,t,2)/g=(V^2*%e^(-Z/H)*sin(θ)*ρ0)/(2*B)-1;
p1a:'diff(X,t,2)/g=-(V^2*%e^(-Z/H)*cos(θ)*ρ0)/(2*B);

g : Z 'direction' ( 'Gravitation')
B = ballistic coef =  'B=(g*m)/(Cd*A[fr])', A[fr]=Af
V = velocity 
ρ = Density ( ρ0 @  Z = 0)
H = atmospheric scale height =' H=(R*T)/g ' = cte
θ = re-entry angle
X,Z = direction 

start : vertical motion
--------------------------------
================================
=>   V = dZ/dt , X =cte  take X=0, θ = π/2
================================

p2: V = 'diff(Z,t);
p2a: θ = π/2;

RE-ENTRY equations of motion 2D : 'p2,p2a', vertical
============================================

p3:subst([p2,p2a],p1);
p3a:subst([p2,p2a],p1a);

How to solve ,'p3'
===============

try : ode2 
========

p4:ode2(p3,Z,t);

p4a:p4[1];
p4b:p4[2];

rem : no solution, 'first order equation not linear in y'

see : mathematics books, maxima help  , 'gamma_incomplete'.
====================================================
------------------------------------------------------------------------------------------------
 => see wikipedia.
info : https://en.wikipedia.org/wiki/Incomplete_gamma_function.
-------------------------------------------------------------------------------------------------

try method 2 :
===========

1e)  use the change : Rule.
----------------------------------------

p5:'diff(V,t) ='diff(rhs(p2),t);

rem : math : 'diff(Z,t,2) = dV/dt = dV/dZ*dZ/dt  = V*dV/dZ
---------------------------------------------------------------------------------
=============================================
=>      Maxima ( find : V function of Z) = V(Z).
=============================================

p5a:rhs(p5) = V*'diff(V,Z);

rem : V*('diff(V,Z,1)) = 1/2*d(V^2)/dZ.

p5b:let(V*('diff(V,Z,1)),1/2*'diff(V^2,Z));

use : ' p3 ' and 'p2'

p5c:subst([p5a],p3)$
p5d:subst([rhs(p2)=lhs(p2)],p5c);

use rule : 'p5b'

p5e:letsimp(p5d)*2*g;

p5f:expand(p5e);

p5g:lhs(p5f)-part(rhs(p5f),1) = part(rhs(p5f),2) ;

++ =rem : coef V^2 , 'lhs',p5g

p5h:ratcoeff(lhs(p5g),V^2);

rem : use ode2, no proper differential equation.

depends([G,F,V],Z);

p6:G(Z)*'diff(V^2*F(Z),Z) = rhs(p5g);

p7:ev(p6,diff);

p8:V*('diff(V,Z,1)) = 1/2*'diff(V^2,Z,1);

p9:lratsubst([p8],p7);

correlation between  'p5g' and 'p9'

p10:F(Z)*G(Z) =1;

p10a:solve(p10,G(Z));

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

rem : coef V^2 , 'lhs',p11 ,see above : ++

p12:ratcoeff(lhs(p11),V^2);

relation ,coef  V^2

p13:p12=p5h;

use relation : 'p10a' into 'p13'

assume(B>0,g>0,ρ0>0,H>0);

p14:subst(p10a,p13);

rem : math :d F(Z)/F(Z) = -(%e^(-Z/H)*g*ρ0)/B *dZ;

p14a:[F(Z)=F];

p14b:subst(p14a,p14)*F;

because : subst , p14a , we can use ,ode2 ( Maxima).
=>  Find  F ?  , from 'p7' .

p14c:ode2(p14b,F,Z);
p14d:subst([%c = c2],p14c);

F : rem F=F(Z)
---------------------

    find c2 : ?
===========

find coef : 'c2'

p15:1/ratcoeff(rhs(p14d),c2);

p15a:1/F = 1/c2*p15;

modification ,ODE, 'p6'   => F(Z)=F,G(Z) =G

p15b:[F(Z) = F,G(Z) = G,G=1/F];
p15c:subst(p15b,p6);

p15d:p15c*F;

assume(c2>0);

p15e:subst([p14d],p15d);

analysis : lhs, 'p15e ', to do integration  ' dZ'
====================================
=> How to use 'pickapart ' , see : help  'maxima'
=>  implement new constant : c3.

p15f:pickapart(p15e,2);

p15g : %t44;

solve ,p15e ,  by integration .
========================

p16:p15g = c3-'integrate(%t45,Z,Z1,Z2);

p16a:expand(p16/c2);
p16b:subst([c3/c2=c4],p16a);
p16c:solve(p16b,V^2)[1];
p16d:expand(p16c);

We need to find : C4 ?
=>  Z = infinity  ,V=V0

p17:factor(p16d);

p17a:subst([V=V0],lhs(p17)) ='limit(rhs(p17),Z,inf);

 I find pickapart : by trail ?
=====================

p17b:pickapart(p17,3);

 evaluation  ,  'p17b' 
------------------------------------------

p17c,limit  (Z=infinity) , %t54

p17c:limit(%t54,Z,inf);

V^2=-%t54*%t55: copy,past, into  'p18'.
=================================

p18:V^2=-%t54*%t55;

p18a:subst([%t54=p17c,Z1=Z4,Z2=Z3],p18);

p18a1:part(rhs(p18a),1);
p18a2:part(rhs(p18a),2);
p18a3:part(rhs(p18a),2,1);
p18a3:part(rhs(p18a),2,1,1);
p18a4:part(rhs(p18a),2,1,3);

p18b:sqrt(p17)$
p18c:(subst([abs(V) =-1*V,Z1=0,Z2=Z5],p18b))*(-1);

analysis : upper bound :  Z1= inf , into integral.
=======================================
rem : Z1= inf , V=V0; => C4
rem : abs(V) = -V

p19:subst([V=V0],p18a);
p19a:solve(p19,c4);
p19b:subst(p19a,p18c);

rem :integrate(M(Z),Z,0,inf) = integrate(M(Z),Z,0,Z) +  integrate(M(Z),Z,Z,inf) .
=============================================================
How to do above :
=>    Z4 = Z , Z3= inf ;  result :second integral  = 0 , by put  Z5 = 0.
-------------------------------------------------------------------------------------------------------------

p20:[Z4 = Z, Z3 = inf , Z5=0];
p20a:subst(p20,p19b);

p20b:I = subst(p20,p18a4);

rem : change variable 
==================
Menu : Calculus => change variable
------------------------------------------------------
input : integral /sum : rhs(p20b)
              old variable   : Z
              new variable : Y
                    Equation  : Y= Z/H
------------------------------------------------------

p20c:changevar(rhs(p20b), Y=Z/H, Y, Z);
p20c1: Y = Z/H;
p20d:I1= subst([rhs(p20c1)=lhs(p20c1)],p20c);

def  K :(H*g*ρ0)/B;
===============

p21:(H*g*ρ0)/B = K;

assume(K>0);

p21a:lratsubst([p21],p20d);

p21a1:l2=subst([K=1],rhs(p21a));

p21aa:l2a=subst([H=1,K=1],rhs(p21a));

use : taylor expantion  : exp(x)

use : menu =>  get series ?
input : variable : x
             point       : 0
             enable : power series.

p21b:x=K*%e^(-Y);

p22:niceindices(powerseries(exp(x), x, 0));

p23:exp(x) = p22;
p23a:subst([p21b],p23);

p23b:subst([p23a],p21a);

rem : don't know how to interchange position , 'sum,integral'
===================================================

p24:lhs(p23b) = H*sum(K^i/i!*' integrate(%e^(-Y*i),Y,Y,inf),i,0,inf);

assume(i>0,Y>0);

p24a:w1 = 'integrate(%e^(-Y*i),Y,Y,inf);

p25:changevar(integrate(%e^(-Y*i),Y,Y,inf), y=Q*i, Q, Y);

p26:[rhs(p24a) = p25];
p26a:subst(p26,p24);
p26b:subst([p20c1],p26a);

rem :  ? 'p21a' and 'p26b'
=====================

p27:rhs(p21a)= rhs(p26b);

final result  V : ?
==============

p28:p20a;

subst  K : 'V'.

p29:lratsubst([p21],p28);

final result  V : ?
==============

p30:lhs(p29)=part(rhs(p29),1);