kill(all)$

load("lrats")$

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

=>   rem : list of rational subst. =  lratsubst ' command'.

filename : Ballisticmiissiletrajectoriesneglectburnoutreentryversion1final.wxmx

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

Ballistic Missile Trajectories        'Version 1.0'
Neglecting (Burnout and Re-entry) ellipse with focus at center of the earth
peter.vlasschaert@gmail.com,                               13/06/2018

(u)nit vector : x-direction : i_u
                           y-direction :j_u
                           z-direction :k_u

unit vector : u

u:[i_u,j_u,k_u];

components : x,y,z  direction

pa1:u[1];
pa2:u[2];
pa3:u[3];

equation ellipse , polar coordinates.(focus at the center of the earth).

pb:r= p/(1+e*cos(θ));

initial conditions : fig 1

p1:θ_0+ψ /2= %pi;

v_0 : initial speed 'missile'
R_0  : radius of the earth (Sphere).

rem : β_0 : ' angle with local horizon ' (= perpendicular to R_0)
components V0  : x-direction : v_0*cos(β_0).
                                   y-direction : v_0.sin(β_0).

How to find : Ø_0  ? fig1 
====================
Because we need components vector V0 =V0(Ø_0)
============================================================
@: use (cos,sin), components ' angle measured with ref : ' pos x-axis'
      result : (cos(°)^2+sin(°)^2=1)   ,{°  = stands for any angle}
============================================================
1e)  θ_0  => θ_0-π/2        : see fig 1 (use @ )
2e)  β_0 : complement  :  π/2-β_0

p2:Ø_0 =( θ_0 - %pi/2)+(%pi/2-β_0);

Components V0 : speed missile
==============

V0:[v_0*cos(Ø_0),v_0*sin(Ø_0)];

x-direction , y-direction
-----------------------------------

p3: V0[1];
p4:V0[2];

use : 'p2' , components : x,y ?
---------------------------------------------

p3a:sublis([p2],p3);
p4a:sublis([p2],p4);

Components r0 :  free flight position
==============
rem : @  'see above':θ_0 => θ_0-π/2
rem : radius of the earth = R_0

x-direction,y-direction
---------------------------------

p5:r0:[R_0*cos(θ_0-%pi/2),R_0*sin(θ_0-%pi/2)];
p5a:r0[1];
p5b:r0[2];

calculation : h_vec  = r0 x V0
==============================================================
rem : {x = cross product  = vector product}
reason : vector : r0,vector : V0, h_vec =stands perpendicular both vectors
                 rem : r0  perpendicular h_vec
                 rem : V0   perpendicular h_vec
physical interpretation : r0 and V0 same plane ( = orbital plane)
==============================================================

How to find : component (k_u direction)

3x3 matrix : How to build
--------------------------------------
Menu : Algebra : Enter matrix
Rows       : 3
Columns : 3
general
Name        : h_vec
---------------------------------------------------------------------------------------------------------------
first        row : pa1,pa2,pa3
second row :p5a,p5b,0
tirhd       row :p3a,p4a,0

p6:h_vec=matrix(
 [pa1,pa2,pa3], 
 [p5a,p5b,0], 
 [p3a,p4a,0]
);

use ' maxima ' ,rhs = right hand side.
                            ,lhs = left    hand side.

p6a : rhs(p6);
p6b: lhs(p6);

menu : Algebra: Determinant.

p7:p6b=determinant(p6a);

How to do this : steps
                            1e)  high light expression after '(p7)' 
                            2e)  Menu : Simplify : Trignometric Simplification :Reduce Trignometric.

trigreduce(h_vec=k_u*(R_0*v_0*sin(θ_0)*sin(θ_0-β_0)+R_0*v_0*cos(θ_0)*cos(θ_0-β_0)));

rem : %   ,use equation from just before.

p8: %;

result : k_u component  from h_vec = h_k

p9: h_k = coeff(rhs(p8),k_u);

rem : 'p7' , h_i = 0,h_j=0

We rename : h_k =h
=================

p10:sublis([h_k=h],p9);

recall  :p = h^2/μ
              μ = M_e *G.
              r = p/(1+e*cos(θ)).

Second Law Newton (Conservative field).
(rem : Attraction between M_e and m)
explain : Attraction direction of radius to center of the earth.
 rem : F =F(r) =F_r: force direction of 'r'

general law : 'Second Law of Newton'
================================

p11:F_r =-1*G*M_e*m/r^2 ;

with approximation : miissile mass = m , r=R_0
----------------------------------------------------------------------

p12a:sublis([r=R_0],p11);

Weight of missile .W_m ' force '
----------------------------------------------

p12b: W_m= -1*m*g0;

follows : F_r =W_m

assume : m>0,M_e>0,G>0

assume(m>0,M_e>0,G>0);

p12c : rhs(p12a) = rhs(p12b);

p12c1:solve(p12c,g0);

p12d:p12c/m;

rem : rational subst .: G*M_e by μ
μ = G.M_e

p12e:lratsubst([G*M_e= μ],p12d);

p12f :solve(p12e,μ);

 r=r(θ)

P : point
p  : distance
PF +FA = 2*a (ellipse )
a = semimajor axis 
b = semiminor axis
How to find : p  to put θ =π/2,polar equation

p22:sublis([θ=%pi/2],pb);

general : r  distance from focus F to any point on the ellipse.

calculation : p = h^2/μ
===================

1e)  μ = ' from above calculation '

p23:rhs(p22) =h^2/μ;

p23a:sublis(p12f,p23);

2e)  h = ' from above calculation '

p23b:sublis([p10],p23a);

From :  Energy per unit mass =E_m , orbit conservation energy .(closed orbit)
               ==========================================================
               E_m= Energy = kinetic Energy+Potential Energy = -μ/2*a , (constant)

general :  (E)nergy (K)inetic = EK  
rem : not per unit mass .

p24: E_k = 1/2*m*v^2;

general :  (E)nergy (P)inetic = EK  
rem : not per unit mass 
F = - grad (v_p)  = > conservative field : F= -d(v_p)/dr)
rem : E_p :potential energy

?  How to find : Potential Energy
============================

1e) Newton Second Law : 'p11' ,F  = F_r

p24a: rhs(p11) = -'diff(E_p,r,1);

How to find : E_p

assume(r>0,q>0);

p24b: E_p=sublis([q=r],integrate(lhs(p24a),r,q,inf));
           

ET = E_K +E_p = constant

p24c:rhs(p24)+rhs(p24b) =constant;

E_m : total energy per unit mass

p24d: p24c/m;

p24d1:expand(p24d);

closed orbit = constant/m = -μ/(2*a)  ,v=v_0,r = R_0

1e)  put in : -μ/(2*a) into 'p24d1'

p24e:lratsubst([rhs(p24d1)=-μ/(2*a)],p24d1);

p24e1:expand(p24e);

2e) put in : v=v_0,r=R_0

p24e2:sublis([v=v_0,r=R_0],p24e1);

3e)  replace : G.M_e by μ

p24e3:[G*M_e = μ];

p24e4:lratsubst(p24e3,p24e2);

p24e5:expand(p24e4);

p = a*(1-e^2)   ?  How to find

distance (FP)  plus distance(FA) = 2*a

1e) distance (FP)  = dist1 , θ = 0

p25:dist1=sublis([θ=0],rhs(pb));

2e) distance (FA)  = dist2 , θ = π

p25a:dist2=sublis([θ=%pi],rhs(pb));

3e) dist1 plus dist2  = 2*a (definition :ellipse)

p25b:rhs(p25)+rhs(p25a) = 2*a;

p25c:solve(p25b,p)[1];

p25d:lhs(p25c)=collectterms(rhs(p25c),a);

How to find : e ,above parameters ={v_0,R_0,g0,β_0}

1e)   e from p25d

p26:solve(p25d,e);

take : pos one.

p26a:p26[2];

2e) put  'a' from  relation 'Energy per unit mass' = p24e5

a1)    solve : 'p24e5' for a

p26b:solve(p24e5,a)[1];

a2)    μ=R_0^2*g0  into  a

p26c:sublis(p12f,p26b);

3e)  p ?   'p23b'

p26d:p23b;

result : e  follows from 'p26c and p26d'
               in to e=-sqrt(1-p/a)
 ===============================

p26e:sublis([p26c,p26d],p26a);

p/a   = 'p26f'

p26e1:p26d/p26c;

p26e2:expand(p26e1);

find : coef  'cos(β_0)^2' from rhs of  'p26e2', Maxima : 'coeff'

p26e3:coeff(rhs(p26e2),cos(β_0)^2);

p26e4:lhs(p26e)=sqrt(1-p26e3*cos(β_0)^2);

   For circular orbit : dv/dt = 0 : radius =constant => dr/dt = 0: d/dt(dr/dt) = 0.
================================================================
acceleration   vec(ac)                   : (tangential component,normaal component)
unit vector acceleration               : (a_t,a_n)
a_n  : direction of r

p27:ac=diff(v,t)*a_t+v^2/r*a_n;

p27a:F_c= -1*coeff(rhs(p27),a_n)*m;

use : Newton second Law .

p27b:p11;

rem : F_c = F_r , missile in to orbit.

initial position : use , r = R_0,v=v_c;
rem : c : symbol into 'v_c'  for 'central'

p27c:sublis([r=R_0,v=v_c],rhs(p27a) = rhs(p27b));

p27d:lratsubst(p24e3,p27c);

find :v_c  ?
symbol  :Δ = v_0/v_c

p27e:solve(p27d,v_c)[2];

p27f:sublis(p12f,p27e);

def  :  Δ  = initial velocity/(circular speed) = v_0/sqrt(R_0*g0)
==================================================

p28: Δ  = v_0/rhs(p27f);

e : e(Δ,β_0)

p28a: lratsubst([rhs(p28) =lhs(p28)],p26e4);

see : above  e= sqrt(1-p/a)

p28b:lratsubst([rhs(p28) =lhs(p28)],p26e2);

p29b:ratsimp(p28b);

p30b:-1*factor(p29b);

p31b:p30b+1;

p31b1:lhs(p28a) = sqrt(rhs(p31b));

rem : trignometric identity : (cos(β_0))^2+(sin(β_0))^2=1

p31b2:1-p/a = Δ^2*(Δ^2-2)*cos(β_0)^2+(cos(β_0))^2+(sin(β_0))^2;

p31b3:lhs(p31b2)=ratsimp(rhs(p31b2));

p31b4:sqrt(p31b3);

p31b5:lhs(p28a)=rhs(p31b4);

use : trajectory equation   r=R_0, θ = θ_0

p32:[ r=R_0, θ = θ_0];

insert 'p32' into 'pb'

p32a:sublis(p32,pb);

ψ/2+θ_0=%pi =π  , from p1

p32b:expand(solve(p1,θ_0));

p32c:cos(p32b);

p32d:lratsubst(p32c,p32a);

?  find p from  'p32b' into  ' polar equation = R_0'

p32e:sublis([p26d],p32d);

Find :ψ  ,free flight angle.

p32f:solve(p32e,cos(ψ/2));

p32g: lratsubst([rhs(p28) =lhs(p28)],p32f);

p32e:ratsimp(p32g)[1];

p32f:expand(p32e);

find : ψ= free flight angle

p32g:solve(p32f,ψ);

Time  of Flight : t_ff  :  θ_0   =>   θ_0 + ψ ( freeflight time)
==============================================
t=t_0     at   θ  = θ_0   
t= T/2    at   θ =  π     
=>  halve of free flight time. = 1/2 time t_ff
t_ff   = 2*(T/2-t_0)
----------------------------------------------------------------------------------
rem :  t = T  =>  θ = 2*π
----------------------------------------------------------------------------------
==============================================
r =a*(1-e*cos(ψ)) =a*(1-e^2)/(1+e*cos(θ))
==============================================
sqrt(R_0*g0^2/a^3)*t= ψ - e*sin(ψ)
==============================================

assume(R_0>0,g0>0);

p1f:t_ff =2*(T/2-t_0) ;

p2f:expand(p1f);

p3f:sqrt(R_0*g0^2/a^3)*t= ψ - e*sin(ψ);

  t = T  => ψ = 2*π
===============

p4f: [t=T,ψ=2*%pi];
p5f:sublis(p4f,p3f);

 solve for  'T'
-------------------

p6f:solve(p5f,T)[1];

a) we want to remove : ' a ' from ' with Δ' => a=a(Δ)
===========================================

p7f:p26c;

a=a(Δ)

p8f1:lratsubst(rhs(p28) = lhs(p28),p7f);

p8f2:fullratsimp(p8f1);

b)  'p8f2 ' into 'p6f'  ( a  into  T )
=========================

p8f3: lratsubst([p8f2],p6f);

Find : T -> T_c ,circular orbit  :  a=R_0
==============================

p9f:[T=T_c,a= R_0];
p10f:sublis(p9f,p6f);

p11f:p8f3/p10f;

rem :  T/T_c =1/(2-Δ^2) ?  same as p11f

p12f:  xx = -1*sqrt(-1)/(-(2-Δ^2))^(3/2);

p12f1:xxx =-1*sqrt(-1)/((-1)^(3/2)*(2-Δ^2)^(3/2));

p13f:lratsubst([xxx = lhs(p11f)],p12f1);

 find t_ff  => t_ff=T-2*t_0
? T,t_0

1e)   find T.

p14f:solve(p13f,T)[1];

p15f:sublis([p14f],p2f);

2e)   a*(1-e*cos(ψ)) =a*(1-e^2)/(1+e*cos(θ))
         θ = θ_0  let   ψ =ψ_0

use relation : f(ψ,...) = g(θ,....)

p16f:[θ=θ_0,ψ=ψ_0];
p17f:a*(1-e*cos(ψ)) =a*(1-e^2)/(1+e*cos(θ));
p18f:sublis(p16f,p17f)/a;

p19f :solve(p18f,cos(ψ_0))[1];

remove : θ_0  with 'p1' => only function ψ_0 ('p19f')

p20f: p1;
p20f1:expand(solve(p1,θ_0));

p20f2:sublis(p20f1,p19f);

3e) find t =t_0 => ψ =ψ_0

p21f:[t=t_0,ψ=ψ_0];
p21f1:sublis(p21f,p3f);

4e)  insert a into  'p21f1' use 'p8f2'

p21f2:sublis([p8f2],p21f1);

5e)   T_c=(2*%pi*R_0)/g0 'p10f ', relation lhs(p21f2)

p22f:p10f;

p23f:solve(p22f,g0);

p24f:sublis(p23f,p21f2);

p25f:solve(p24f,t_0)[1];

p26f:factor(p25f);

rem :  2*pi*(Δ^2-2)^(3/2)   ? same pq1 ,see t_0 above.

pq1:yy = 2*%pi*((-1)^(3/2)*(2-Δ^2)^(3/2));

pq2:part(p26f,2,1);
pq3:rhs(pq1);

p26f : lhs(p26f)=pq2/pq3;

p27f:expand(p26f);

find : total result : t_ff  ?

p28f:p15f;

p29f:sublis([p27f],p28f);

p30f:lhs(p28f) = factor(rhs(p29f));

p31f:radcan(p30f);

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