Best way to illustrate Keplers 2nd law with tikz

Question

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:

Answer 1

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:

Answer 2

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]'");}
add_cus_dep("asy","eps",0,"asy");
add_cus_dep("asy","pdf",0,"asy");
add_cus_dep("asy","tex",0,"asy");

and run latexmk -pdf kepler.tex.

Answer 3

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)

Answer 4

Here's a metapost effort. The macro kep takes six arguments in two equal groups. In the first group there are (in order):

1. Left focus (i.e. position of sun)
2. Right focus
3. Length of semi-minor axis

In the second group there are positions of the planet given as metapost times where (in this case) 360 degrees is equal to 8 units of metapost time starting from the rightmost point of an axis aligned ellipse and travelling counterclockwise. The times given correspond to:

1. Initial position of planet.
2. Terminal position of planet.
3. Initial position of planet from which to calculate a terminal position resulting in equal swept area.

Compile with lualatex:

\documentclass{article}
\usepackage{luamplib}
\begin{document}
\mplibnumbersystem{double}
\mplibtextextlabel{enable}
\begin{mplibcode}
    color orange; orange:=(1,0.5,0);
    color yellow; yellow:=(.99,.93,0);
    color bluea; bluea:=(.74,.83,.9);
    color blueb; blueb:=(0,.47,.75);

    picture sun,planet;
    sun:=image(
    for i=1 upto 50:
        draw fullcircle scaled (i/50) withcolor (i/50)[yellow,orange] withpen pencircle scaled .1bp;
    endfor;
    );
    planet:=image(
    for i=1 upto 50:
        draw fullcircle scaled (i/50) withcolor (i/100)[bluea,blueb] withpen pencircle scaled .1bp;
    endfor;
    );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % https://tug.org/pipermail/metapost/2005-February/000222.html
    vardef area(expr p) = % p is a B\’ezier segment; result = \int y dx
        save xa, xb, xc, xd, ya, yb, yc, yd;
        (xa,20ya)=point 0 of p;
        (xb,20yb)=postcontrol 0 of p;
        (xc,20yc)=precontrol 1 of p;
        (xd,20yd)=point 1 of p;
        (xb-xa)*(10ya + 6yb + 3yc + yd)
        +(xc-xb)*( 4ya + 6yb + 6yc + 4yd)
        +(xd-xc)*( ya + 3yb + 6yc + 10yd)
    enddef;

    vardef Area(expr P) = % P is a cyclic path; result = area of the interior
        area(subpath (0,1) of P)
        for t=1 upto length(P)-1: + area(subpath (t,t+1) of P) endfor
    enddef;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    vardef ellipse(expr p,q,b)=
        save a_,c_; numeric a_,c_;
        c_=abs(p-q)/2;
        a_=c_++b;
        show(angle(q-p));
        fullcircle xscaled 2a_ yscaled 2b rotated angle(q-p) shifted .5[p,q]
    enddef;

    def kep(expr p,q,b)(expr ta,tb,tc)=
        save p_,x_,psize_; path p_[]; numeric x_[];
        psize_=.13u;
        x_1=ta; x_2=tb; x_3=tc;
        p_0=ellipse(p,q,b);
        p_1=p--subpath (x_1,x_2) of p_0--cycle;
        x_0=abs(Area(p_1));
        vardef eq_ar(expr t)=
            abs(Area(p--subpath (x_3,t) of p_0--cycle))<x_0
        enddef;
        x_4=solve eq_ar(x_3+eps,8+x_3-eps);
        p_2=p--subpath(x_3,x_4) of p_0--cycle;
        p_3=fullcircle scaled (psize_+2bp);
        p_4=subpath (x_1,x_2) of p_0 cutbefore p_3 shifted point x_1 of p_0 cutafter p_3 shifted point x_2 of p_0;
        p_5=subpath (x_3,x_4) of p_0 cutbefore p_3 shifted point x_3 of p_0 cutafter p_3 shifted point x_4 of p_0;
        for i=1,2: fill p_[i] withcolor .8; endfor;
        for i=1 upto 4: draw p--point x_[i] of p_0 withcolor .7; endfor;
        draw p_0 dashed evenly withcolor .6 withpen pencircle scaled .75bp;
        draw sun scaled psize_ shifted p;
        for i=4,5: drawarrow p_[i] withpen pencircle scaled 1bp withcolor blue; endfor;
        for i=1 upto 4: draw planet scaled psize_ shifted (point x_[i] of p_0); endfor;
    enddef;

    beginfig(0);
        u=3cm;
        z0=origin;
        z1=(2.2u,0);
        % kep(p,q,b)(t1,t2,t3)
        % p=left focus
        % q=right focus
        % b=length of semi-minor axis
        % t1=start time of known arc
        % t2=end time of known arc
        % t3=start time for second arc 
        kep(z0,z1,.75u)(1.5,3.8,6);
    endfig;
\end{mplibcode}
\end{document}