kill(all)$

polar equation for ellipse (elements),version 1.0
peter.vlasschaert@gmail.com,11/07/2018

load(lrats)$

depends(u,θ);

differential equation : radial component       ' = radius  =>  u =  1/r '
======================================================
=>  ( ode 2 = ordinary differential equation ), two constants c1,c2

p1:'diff(u,θ,2)+u =μ/h^2;

rem : u=1/r

p2:ode2(p1,u,θ);

p3:lratsubst([%k1=c1,%k2=c2],p2);

p4:expand(p3);

rem :trick μ =0  to remove constant , because 'particular solution  = u_p'
rem : u = u_h+u_p
        ---------------------

find : u_p

p4a:lratsubst([u=u_p,c1=0,c2=0],p4);

find: u_h

p4b:lratsubst([u=u_h,μ=0],p4);

p4c1:coeff(rhs(p4b),c1);
p4c2:coeff(rhs(p4b),c2);

rem : c1,c2  => A,θ_0,μ/h^2 = constant
================================

p5: lhs(p4b)=A*cos(θ-θ_0) ;

p6:trigexpand(p5);

p7:expand(p6);

p8:rhs(p4b)=rhs(p7);

p9a:coeff(p8,p4c1);
p9b:coeff(p8,p4c2);

rem : c1,c2 into  'p4'

p10:subst([p9a,p9b],p4);

p11:trigreduce(u=μ/h^2+A*sin(θ)*sin(θ_0)+A*cos(θ)*cos(θ_0));

p12:subst([u=1/r,A=-A],p11);

p13:solve(p12,r)[1];

above and below divide by  'μ'

p14a:part(p13,2,1)/μ;
p14b:part(p13,2,2)/μ;
p14b1:expand(p14b);
p14b2:ratcoeff(p14b1,-cos(θ_0-θ))=e;
p14b3:lratsubst([p14b2,θ_0=ω],p14b1);

rem :(A*h^2)/μ = e  = constant.

p15:lhs(p13)= p14a/p14b3;

elliptical orbit.

vernal equinox : intersection (equator ≠ earth equator,ecliptic)
------------------------
path : great circle : called  ' ecliptic' which lies in the plane (earth orbit about the sun).
rem :  equator Celestical not the earth equator( ned: evenaar).

elements of the ellipse
===================

rem :  distance(A,H) = a
            distance(H,N) = b
            distance(F,Q) = p  (see 'p16a'),definition( 'see below')
            distance(F,B) = r
            distance(0,F) = distance(0,F') = a*e
angle        :  θ - ω  = β ,argument  of cos => | θ - ω | = β
angle        : red      :   ω
angle        : green :   θ
arc angle : black :    β
points  :  F,F' =focuses of ellipse
point    :  H = 0
Figure 1 : build (geogebra,open source).

see figure 1 : F'B +FB = 2*a = 7.14 +1.99  , a : semi-major axis.
                                                                           ,  b : semi-minor axis 
                         e : definition : e = FB/BB' =h/r= 1.99/3.08  < 1
rem : LM = directrix 

by definition :  p= h^2/μ
===================

p16:subst([p=p14a],p15);
p16a:lhs(p15) =p/p14b3;

rem :  cos(-α) =cos(α)

p17:cos(-α);

by definition:  p=a*(1-e^2) for ellipse.
==============================

p18:subst(p=a*(1-e^2),p16a);

rem : allowed choose ref other way

p19:  r=(a*(1-e^2))/(1+e*cos(ω-θ));

we see fig : |θ-ω|=β
================

p20:θ-ω =β;

p21:lratsubst([p20],p19);
p22:expand(p21);
p23:combine(p22);