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

polar equation : relation with ellipse,a,b,c,e,θ,r , 'version 1.0'
peter.vlasschaert@gmail.com,26/09/2018

polar equation :  r ?
=================

e = eccentricity = constant ratio  ( = definition)
======================================

p1:  r/k  = e;
p2: l/k1 = e;

solve for : k,k1
----------------------

p1a:solve(p1,k)[1];
p2a:solve(p2,k1)[1];

distance (GG1) = distance(FL1) + distance(L1 ?)
distance (GG1) = distance(FL1) + distance(H H1)
======================================== 
 k1                       = distance(FL1) + k = |FL1| + k
========================================

p3a: dist_FL1 = r*cos(θ);
p3b: k1 = rhs(p3a) + k;
p3c: subst([p1a,p2a],p3b);

find r : ?
=======

p4:solve(p3c,r);
p4a:p4[1];

e = eccentricity ( relation 'conic section')
=================================
circle          :  e =0
ellipse        : 0 < e < 1
parabola   :  e = 1
hyperbola : e > 1
=================================

circle  :
=====

p5:subst([e=0],p4a);

rem : r = 'radius' = l

ellipse :
======

elements of ellipse : a,b,c,e,θ
========================

nomenclature
---------------------

a = semi-major axis
b= semi-minor axis
c = focus distance 
θ = true anomaly
r  = radius
l = semilatus rectum = |FF1|

How to find l : ?
=============

vector 'v' , angle : θ = 0 (into 'p4a')
---------------------------------------------------

p5:subst([θ=0,r=r[p]],p4a);

rem : |FG| = r[p] = 'perigee'

vector 'u' , angle : θ = π (into 'p4a')
---------------------------------------------------

p5a:subst([θ=π,r=r[a]],p4a);

rem : |FG| = r[a] = 'apogee'

from definition 'ellipse' (see Figure  2)
===================

p5b: p5+p5a;

p5c:2*a = lhs(p5b);
p5c1:2*a=rhs(p5b);

result : l

p5d:(solve(p5c1,l));
p5e:p5d[1];

polar equation ' ellipse'  , use 'p4a'
====================

p6:subst(p5d,p4a);

How to find b : ?   => b = f1(a,e)
==========================

first method : pythagoras and definition 'ellipse'
------------------------------------------------------------------------

rem : see triangle : F'p F  = F'p 0 + 0 p F (both triangles are equal triangle = right triangles)
--------------------------------------------------------------------------------------------------------------------------------------
1e) definition         :  2*a = |F'p| + |pF| = g+h  
2e) right triangles : 2.a = 2*|pF| = 2*h =2*g
--------------------------------------------------------------------------------------------------------------------------------------

definition : right triangle
=====================
b^2 + c^2 = |pF|^2 = a^2

p7: b^2 +c^2 = a^2;
p7a:p7-part(lhs(p7),1);

find c ?  c=f2(e,a)
--------------------------
a=|0G| = r[p] + c

from 'p5' = r[p] use 'l'

p7b:subst([p5d],p5);

p7c:ratsimp(p7b);

p7d : a = rhs(p7c)+c;

p7e:solve(p7d,c);
p7f:p7e[1];

result : b (use 'p7a') and 'p7f'
--------------------------------------------

p8:subst(p7e,p7a);

p8a:solve(p8,b);

p8b:p8a[2];

second method : calculus , maximum
----------------------------------------------------------
----------------------------------------------------------

first step       :start from : 'p6'
=======================

p9:p4a;

p9a:x=r*cos(θ);
p9b:y=r*sin(θ);

second step : differentation with respect to θ
======================================
rem : 'b' y - direction , 'a' x- direction
---------------------------------------------------------------------

p9c:subst([p9],p9b);

thrid step : find max ' 'diff(y,θ) ' 
============================

p9d:'diff(lhs(p9c),θ) = diff(rhs(p9c),θ);
p9e:ratsimp(p9d);

'diff(y,θ) ' = 0

p9f:part(rhs(p9e),1) = 0;
p9g:trigsimp(p9f);
p9h:factor(p9g);

l ≠ 0, (cos(θ)+e) = 0

p10:part(lhs(p9h),2) = rhs(p9h);

p10a:solve(p10,cos(θ));
p10a1:p10a[1];

fourth step : plugin 'p10a' into 'p9c'
=============================

p10b: subst(p10a,p9c);

rem : sin(θ)  ,from  cos(θ)=-e  'p10a1'

identity  : geometric circle.
=>  (sin(α))^2+(cos(α))^2 = 1;

p10c: (sin(θ))^2+(cos(θ))^2 = 1;
p10d:subst(p10a,p10c);
p10e:solve(p10d,sin(θ));
p10f:p10e[2];

fifth step : max 'y' = b use l=a-a*e^2
==============================

final result b :

assume e < 1

assume(e<1 and e > 0,a>0.);

p11:subst([lhs(p9c)=b,p10f,p5e],p10b);

negdistrib;

p11a:1-e^2 =m;

p11b:lratsubst([p11a],p11);

p11c:subst([rhs(p11a) =lhs(p11a)],p11b);

method 1 = method 2  =>  b=a*sqrt(1-e^2)

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

rem : negdistrib  ?

negdistrib : not negdistrib;

p12:factor(p11);

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