# Best way to illustrate Keplers 2nd law with tikz

I am trying to visualize Keplers second law of planetary motion with tikz.

I managed to do the following:

\documentclass{article}

\usepackage{tikz}

\begin{document}

\begin{tikzpicture}
\draw [rotate around={0:(1.5,0)}] (1.5,0) ellipse (2.5cm and 2cm);
%
\fill (0,0) coordinate (O) circle (2pt) node[below =7pt] {sun};%
\coordinate (A1) at (1.62,2) ;%%
\coordinate (A2) at (2.15,1.93);%%
\coordinate (B1) at (-0.81,0.77);%%
\coordinate (B2) at (-0.9,0.56);%%
\coordinate (B3) at (-0.95,0.39);%%
\coordinate (B4) at (-0.99,0.14);%%
\coordinate (B5) at (-1,-0.05);%%
\coordinate (B6) at (-0.97,-0.31);%%
\coordinate (B7) at (-0.92,-0.52);%%
%
%
\coordinate (C1) at (3.25,-1.43);%%
\coordinate (C2) at (3.51,-1.18);%%
%
\coordinate (P) at (3.42,1.28) ;%%
\fill (P) circle (1pt) node[above right] {planet};%
%
%
\filldraw[fill=blue,opacity=0.5] (O) -- (A1) -- (A2) --cycle;%
\filldraw[fill=blue,opacity=0.5] (O) -- (B1) -- (B2) -- (B3) -- (B4) -- (B5) -- (B6) -- (B7)
--cycle;%
\filldraw[fill=blue,opacity=0.5] (O) -- (C1) -- (C2) --cycle;%
%

\draw (1.5,0) coordinate (M) --node[above]{\footnotesize $a$}
(4,0);%
\fill (M) circle (1pt);
\end{tikzpicture}

\end{document}


This gives the following output:

However there are multiple issues:

### How to fill the ellipse sector?

To fill the ellipse sectors I approximate them with polygons manually, in particular for the "larger" (but of course with equal area) part I used multiple manual "helper points". I did that by first drawing the image with geogebra, exporting it and cleaning up the code a bit. So is there a smart way to fill the true ellipse segment?

### How to get the equal areas automatically?

To make the areas of (approximately) same size I roughly measured the approximations with geogebra and adjusted it with the mouse (very painful :-)).

In an ideal solution the following would be possible:

• Specify two points of the ellipse (or the corresponding parameter of a parametrization) then the corresponding swept out area is filled automatically.

• Specify another point of the ellipse and a corresponding swept out area is drawn automatically (with the same area as the first one). You can do this multiple times as you like.

### It should be easy to modify the picture

For example when I change the sizes of the axes of the ellipse everything should change accordingly. This is not true in my solution because I specified the points manually.

Perhaps it is a good idea to externalize some calculations for this. That's why I added the luatex and python tags.

Edit:

-
If you wish to replicate the animation in the link, well that's a job for mozartstraße :) – Claudio Fiandrino Jul 2 '13 at 14:59
@ClaudioFiandrino I don't want to reproduce the animation in the link, I think the animation part would be another question or a kind of bonus to this question. I just added the link for those who don't know Keplers laws. The result should look just similar to mine but fix the coding issues... – student Jul 2 '13 at 15:00
How would you provide the coordinates optimally? Are they relative to the center of the ellipse or (as in your question and the answer) relative to one of the focus point? Or would you rather state two angles or one angle and a delta angle? – Qrrbrbirlbel Jul 2 '13 at 17:57
@Qrrbrbirlbel See my edit. – student Jul 2 '13 at 19:25

For the problem of filling the ellipse sectors, you can draw "bigger" triangles and then clip them to the elipse shape.

For this approach it is better to have the points A1, A2, B1, B7, C1 and C7 in polar coordinates. In fact, only the angle is important, since the radius will be made long enough to guarantee that the point is out of the ellipse. A radius of 5 is enough in this example.

The following code implements this idea:

% We define the orbit as a macro because we will use it twice, first for clipping and then
% to actually draw the ellipse. This way we avoid inconsistencies.
\def\orbit{(1.5,0) ellipse(2.5cm and 2cm)}

\begin{tikzpicture}
\fill (0,0) coordinate (O) circle (2pt) node[below =7pt] {sun};%
\coordinate (A1) at (50.992527:5);
\coordinate (A2) at (41.913511:5);
\coordinate (B1) at (136.450216:5);
\coordinate (B7) at (-150.524111:5);
\coordinate (C1) at (-23.749494:5);
\coordinate (C2) at (-18.581735:5);
\coordinate (P) at (3.42,1.28) ;%%
\fill (P) circle (1pt) node[above right] {planet};%

\begin{scope}  % The blue shaded regions
\clip \orbit;
\filldraw[fill=blue,opacity=0.5] (O) -- (A1) -- (A2) -- cycle;
\filldraw[fill=blue,opacity=0.5] (O) -- (B1) -- (B7) -- cycle;%
\filldraw[fill=blue,opacity=0.5] (O) -- (C1) -- (C2) -- cycle;%
\end{scope}

% The ellipse
\draw \orbit;

\draw (1.5,0) coordinate (M) --node[above]{\footnotesize $a$}  (4,0);
\fill (M) circle (1pt);
end{tikzpicture}


This is the result:

## Update. Automate the finding of the sectors approximating them by tiangles or circular sectors

The following code implements some ideas, but the implementation is very hackish. These are the ideas:

1. Given a pair of initial and final angles (given indeed by two points located in the outside of the ellipse) the macro \ComputeArea calculates the area of the triangle formed by the sun and the two points in the orbit at those angles.
2. Given any other point in the orbit, the macro \ComputePointNextTo finds the next point in the orbit (in counterclockwise direction) which covers the same area computed before. In this case the sector is assumed to be a circular sector with center at the sun, instead of a triangle, to simplify the calculus.

In order to solve 1, I used the formula found here which gives the area of a triangle from the coordinates of its three vertex. In order to implement this in TikZ I had to first find the three points, which involves solving some intersections. The formula is implemented in a let...in path, and saved via \xdef in a macro called \area for later use.

In order to solve 2, I used the formula for the area of a circular sector of angle theta, which is area=(theta*r^2), given theta in radians. Finding for theta we have then: theta = 2*area/r^2. I implemented this formula again in a let...in path and from this value of theta I built a coordinate called (result) which lies in the appropiate angle at the outisde of the ellipse.

The complete code follows. In this case I kept the original figure with blue regions exactly as the one given by the OP, and added my computations.

The area of the "big" sector is computed, and the result is shown below the figure, for debugging purposes (the unit of length is pt, so the resulting area is in pt^2).

For each of the other blue sectors, I used the first point (A1) and (C1) as the "given point" and computed as described the other "next points" (A2) and (C2). The figure shows in two red lines over the blue sectors the directions in which the found points are.

As you can see, the approximation is good enough unless the figure has to be used to take precise measures on it.

### Code:

\def\orbit{(1.5,0) ellipse(2.5cm and 2cm)}

\def\ComputeArea#1#2{
\path[name path=orbit] \orbit;
\path[name path=line1] (O) -- (#1);
\path[name path=line2] (O) -- (#2);
\path[name intersections={of=orbit and line1,by=aux1}];
\path[name intersections={of=orbit and line2,by=aux2}];
\path let \p1=(O),
\p2=(aux1),
\p3=(aux2),
\n1 = {abs(\x1*(\y2-\y3)+\x2*(\y3-\y1)+\x3*(\y1-\y2))/2.0}
in  node[above]  {\pgfmathparse{\n1}\xdef\area{\pgfmathresult}};
}

\def\ComputePointNextTo#1{
\path[name path=line1] (O) -- (#1);
\path[name intersections={of=orbit and line1,by=aux1}];
\path let \p1=($(aux1)-(O)$),
\n1 = {veclen(\p1)},    % r
\n2 = {atan2(\x1,\y1)}, % initial angle
\n3 = {deg(2*\area/\n1/\n1)}    % angle to cover
in coordinate (result) at (\n2+\n3:5);
}

\usetikzlibrary{intersections,calc}

\begin{tikzpicture}

% Original figure (using the clipping technique)
\fill (0,0) coordinate (O) circle (2pt) node[below =7pt] {sun};%
\coordinate (A1) at (41.913511:5);
\coordinate (A2) at (50.992527:5);
\coordinate (B1) at (136.450216:5);
\coordinate (B7) at (-150.524111:5);
\coordinate (C1) at (-23.749494:5);
\coordinate (C2) at (-18.581735:5);

\coordinate (P) at (3.42,1.28) ;%%
\fill (P) circle (1pt) node[above right] {planet};%

\begin{scope}
\clip \orbit;
\filldraw[fill=blue,opacity=0.5] (O) -- (A1) -- (A2) -- cycle;
\filldraw[fill=blue,opacity=0.5] (O) -- (B1) -- (B7) -- cycle;%
\filldraw[fill=blue,opacity=0.5] (O) -- (C1) -- (C2) -- cycle;%
\end{scope}

\draw \orbit;
\draw (1.5,0) coordinate (M)
--node[above]{\footnotesize $a$} (4,0);
\fill (M) circle (1pt);

% Added. Trying to automatically find (A2) and (C2)
% from (A1) and (C1), such that the area is equal to the
% sector from (B1) to (B7)

\ComputeArea{B1}{B7}
\node[right] at (0,-2.3) {Area: \area};  % Show it, for debugging

\ComputePointNextTo{A1}
\draw[red] (O) -- (result);

\ComputePointNextTo{C1}
\draw[red] (O) -- (result);
\end{tikzpicture}


### Result:

-
Not that anything stands against saving the path in a macro (TikZ does support this) you can also save it as a style with the help of insert path: the saveuse path key. Might be helpful in some cases. :) – Qrrbrbirlbel Jul 2 '13 at 15:25
Is there a way to automate finding equal areas? So if someone specifies one sector, TikZ will find another sector of equal area? – dustin Jul 2 '13 at 15:31
Would the trajectory equation help? – dustin Jul 2 '13 at 16:16
In the last post here seems to be an answer: mathforum.org/library/drmath/view/53635.html – student Jul 2 '13 at 17:01
Very interesting. By the way, you can draw an ellipsoid arc just the same if you know the angles (from its center though): arc[start angle=<sa>, end angle=<ea>, x radius=<xr>, y radius=<yr>] (or (<sa>:<ea>:<xr> and <yr>)). I wonder if it would help to change TikZ’ coordinate system here with [x=2.5cm, y=2cm] then all ellipses and circles with radius 1 (and not 1cm!) would be stretched accordingly while polar coordinate are especially easy. Speaking of “unit circle” … – Qrrbrbirlbel Jul 2 '13 at 17:47

The Keplerequation (link to German wikipedia, which is unusually more informative than the English one on this topic) has no algebraic/closed solution. There are good approximations, but if one has to be approximative from the beginning, one can also simulate the physics instead of doing the math:

\documentclass{standalone}
\usepackage{etoolbox}
\usepackage{tikz}
\gdef\myposx{10.0}
\gdef\myposy{0.0}
\gdef\vx{0.0}
\gdef\vy{4.6}
\gdef\forcefactor{150}
\gdef\deltat{0.01}
\gdef\smallmass{1}
\gdef\startone{100}
\gdef\endone{200}
\gdef\starttwo{1800}
\gdef\endtwo{1900}
\gdef\pathone{}
\gdef\pathtwo{}
\begin{document}
\begin{tikzpicture}[scale=0.2]
\filldraw(0,0)circle(0.1);
\foreach \n in {1,...,3625}
{
\pgfextra{%
\global\let\oldx\myposx
\global\let\oldy\myposy
\pgfmathsetmacro{\currentsquareddistance}{\myposx*\myposx+\myposy*\myposy}
\pgfmathsetmacro{\currentforce}{\forcefactor/\currentsquareddistance}
\pgfmathsetmacro{\currentangle}{atan2(\myposx,\myposy)}
\pgfmathsetmacro{\currentforcex}{-1*\currentforce*cos(\currentangle)}
\pgfmathsetmacro{\currentforcey}{-1*\currentforce*sin(\currentangle)}
\pgfmathsetmacro{\currentvx}{\vx+\deltat*\currentforcex/\smallmass}
\pgfmathsetmacro{\currentvy}{\vy+\deltat*\currentforcey/\smallmass}
\pgfmathsetmacro{\currentposx}{\myposx+\deltat*\currentvx}
\pgfmathsetmacro{\currentposy}{\myposy+\deltat*\currentvy}
\global\let\myposx\currentposx
\global\let\myposy\currentposy
\global\let\vx\currentvx
\global\let\vy\currentvy
\global\let\forcex\currentforcex
\global\let\forcey\currentforcey
\global\let\myangle\currentangle
\ifnumequal{\n}{\startone}{%
\global\let\startonex\oldx
\global\let\startoney\oldy
\xappto{\pathone}{(\oldx,\oldy)}
}{}
\ifnumequal{\n}{\starttwo}{%
\global\let\starttwox\oldx
\global\let\starttwoy\oldy
\xappto{\pathtwo}{(\oldx,\oldy)}
}{}
\ifnumequal{\n}{\endone}{%
\global\let\endonex\myposx
\global\let\endoney\myposy
\xappto{\pathone}{,(\myposx,\myposy)}
}{}
\ifnumequal{\n}{\endtwo}{%
\global\let\endtwox\myposx
\global\let\endtwoy\myposy
\xappto{\pathtwo}{,(\myposx,\myposy)}
}{}
}
% \draw[very thin,->](\oldx,\ol dy)--++(\forcex,\forcey);
\ifnumgreater{(\n-\startone)*(\endone-\n)}{-1}
{
\pgfextra{%
\xappto{\pathone}{,(\myposx,\myposy)}
}
}
{}
\ifnumgreater{(\n-\starttwo)*(\endtwo-\n)}{-1}
{
\pgfextra{%
\xappto{\pathtwo}{,(\myposx,\myposy)}
}
}
{}
\draw(\oldx,\oldy)--(\myposx,\myposy);
}
\begin{scope}[red]
\filldraw (0,0)%
\foreach \point in \pathone
{%
--\point
}--(0,0);
\filldraw (0,0)%
\foreach \point in \pathtwo
{%
--\point
}--(0,0);
\end{scope}
\end{tikzpicture}

\end{document}


The values at the beginning are in totally arbitrary units. Realism would do much better, because our planets have very inexcentric orbits, so there wouldn't be much to see.

What's still to do:

• compute the full period OR replace the \foreach-loop by a while-construct that terminates, when the ellipse is full
• compute the start and end points out of the time phase
• replace all the global definitions by appropriate ones (I'm not good in TeX scoping)
-
+1 for the physical approach! – bersanri Jul 3 '13 at 8:05
I get TeX capacity exceeded, sorry [pool size=1179524]. Perhaps it would be better to externalize the calculations with luatex, phyton or something like that. – student Jul 3 '13 at 10:38
@student Yes, that can happen. Of course, you can externalize the computations. You can also reduce the exactness (given by \deltat. one needs to adapt the upper limit of the loop then). – Toscho Jul 3 '13 at 15:13

The basic struct PlanetaryMotion which handles the elliptic sector area calculations is defined in asydef environment, and two examples of the illustration are shown in two asy pictures.

kepler.tex:

\documentclass{article}
\usepackage{lmodern}
\usepackage{subcaption}
\usepackage[inline]{asymptote}
\usepackage[left=2cm,right=2cm]{geometry}

\begin{asydef}
import graph;
import gsl; // for newton() solver
size(200);

struct PlanetaryMotion{
real a,b,e;
real planetTime,sunR,planetR;
pair F0,F1;

guide orbit;

transform tr=scale(-1,-1); // to put the Sun in the left focus

pair ellipse(real t){
return (a*cos(t),b*sin(t));
}

real Area(real t){ // area 0..t
return a*b/2*(t-e*sin(t));
}

real calcArea(real t0,real t1){
return Area(t1)-Area(t0);
}

real AreaPrime(real t){
return 1/2*a*b*(1-e*cos(t));
}

real findTime(real areaToFit, real tstart){ // find time tend to fit areaToFit
real tend=newton(new real(real t){return calcArea(tstart,t)-areaToFit;}
,new real(real t){return AreaPrime(t);},tstart,tstart+2pi);
return tend;
}

void drawBG(){
draw(tr*orbit,darkblue);
filldraw(tr*shift(F0)*scale(sunR)*unitcircle,yellow,orange);
filldraw(tr*shift(ellipse(planetTime))*scale(planetR)*unitcircle,blue,lightblue);

//dot(tr*F1,UnFill);
label("$F_0$",tr*F0,3N);
//label("$F_1$",tr*F1,3N);
label("Sun",tr*F0,3S);
label("planet",tr*ellipse(planetTime),SW);

draw(((0,0)--(a,0)));
label("$a$",(a/2,0),N);
dot((0,0),UnFill);

}

void drawSector(real t0, real t1,pen p=blue+opacity(0.3)){
fill(tr*(F0--graph(ellipse,t0,t1)--cycle),p);
}

void operator init(real a, real b
,real planetTime
,real sunR=0.05a, real planetR=0.3sunR
){
this.a=a;
this.b=b;
this.planetTime=planetTime;
this.sunR=sunR;
this.planetR=planetR;
this.e=sqrt(a^2-b^2)/a;
this.F0=(a*e,0);
this.F1=(-a*e,0);
this.orbit=graph(ellipse,0,2pi);
}
}

\end{asydef}
\begin{document}
\begin{figure}
\captionsetup[subfigure]{justification=centering}
\centering
\begin{subfigure}{0.49\textwidth}
\begin{asy}
PlanetaryMotion pm=PlanetaryMotion(1,0.618,1.2pi);

pm.drawBG();

real t0,t1,t2,t3,t4,t5;

t0=-0.1pi;
t1= 0.1pi;

pm.drawSector(t0,t1);

real area0=pm.calcArea(t0,t1);

t2=0.7pi;
t3=pm.findTime(area0,t2);
pm.drawSector(t2,t3);

t4=1.5pi;
t5=pm.findTime(area0,t4);
pm.drawSector(t4,t5);
\end{asy}
\caption{}
\label{fig:1a}
\end{subfigure}
%
\begin{subfigure}{0.49\textwidth}
\begin{asy}
PlanetaryMotion pm=PlanetaryMotion(1,0.8,1.35pi,sunR=0.09);

pm.drawBG();

real t0,t1,t2,t3,t4,t5;

t0=-0.05pi;
t1= 0.17pi;

pm.drawSector(t0,t1);

real area0=pm.calcArea(t0,t1);

t2=0.4pi;
t3=pm.findTime(area0,t2);
pm.drawSector(t2,t3);

t4=1.7pi;
t5=pm.findTime(area0,t4);
pm.drawSector(t4,t5);
\end{asy}
\caption{}
\label{fig:1b}
\end{subfigure}
\caption{Illustration of Keplers 2nd law with \texttt{Asymptote}.}
\end{figure}

\end{document}


To process it with latexmk, create file latexmkrc:

sub asy {return system("asy '\$_[0]'");}

and run latexmk -pdf kepler.tex.