%% $Id: pst-kepler.tex 1227 2026-01-01 07:39:19Z herbert $
%%
%% This is file `pst-kepler.tex',
%%
%% IMPORTANT NOTICE:
%%
%% Package `pst-kepler.tex'
%%
%% Authors  :   Jürgen Gilg (8.2.1966 - 6.5.2022)
%%              Manuel Luque <Manuel.Luque27@gmail.com>
%%              Thomas Ossig 
%%              Jean-François Burnol
%%              Herbert Voss <hvoss@tug.org>
%%
%% This program can be redistributed and/or modified under the terms
%% of the LaTeX Project Public License Distributed from CTAN archives
%% in directory macros/latex/base/lppl.txt.
%%
%% DESCRIPTION:
%%   `pst-kepler' is a PSTricks package to draw kepler's view of the world
%%
\csname PSTkeplerLoaded\endcsname
\let\PSTkeplerLoaded\endinput
\ifx\PSTricksLoaded\endinput\else\input pstricks.tex\fi
\ifx\PSTSOLIDESIIIDLoaded\endinput\else\input pst-solides3d.tex\fi
\ifx\PSTXKeyLoaded\endinput\else\input pst-xkey.tex \fi
%
\def\fileversion{0.01}
\def\filedate{2026/01/01}
\message{`PST-kepler' v\fileversion, \filedate\space (jg,hv,jfb,ml,to)}
%
\edef\PstAtCode{\the\catcode`\@} \catcode`\@=11\relax
\pst@addfams{pst-kepler}
%

%Icosahedron, dodecahedron, octahedron, hexahedron, tetrahedron

\newcount\hvSet@type
\define@choicekey*+[psset]{pst-kepler}{type}[\val\nr]{world,tetrahedron,hexahedron,octahedron,%
   dodekahedron,icosahedron}[world]{% 
  \global\hvSet@type=\nr
}
{\PackageWarning{pst-kepler}{erroneous input (#1) for type ignored. Using world.}%
  \global\hvSet@type=0} % preset to world


\define@boolkey[psset]{pst-kepler}[Pst@]{tetrahedron}[true]{}
\define@boolkey[psset]{pst-kepler}[Pst@]{hexahedron}[true]{}
\define@boolkey[psset]{pst-kepler}[Pst@]{octahedron}[true]{}
\define@boolkey[psset]{pst-kepler}[Pst@]{dodecahedron}[true]{}
\define@boolkey[psset]{pst-kepler}[Pst@]{icosahedron}[true]{}
\define@key[psset]{pst-kepler}{scaleTetra}[2]{\def\psk@kepler@scaleTetra{#1 }}
\define@key[psset]{pst-kepler}{scaleHexa}[4]{\def\psk@kepler@scaleHexa{#1 }}
\define@key[psset]{pst-kepler}{scaleOcta}[6]{\def\psk@kepler@scaleOcta{#1 }}
\define@key[psset]{pst-kepler}{scaleDodeca}[8]{\def\psk@kepler@scaleDodeca{#1 }}
\define@key[psset]{pst-kepler}{scaleIcosa}[10]{\def\psk@kepler@scaleIcosa{#1 }}
\psset{icosahedron, dodecahedron, octahedron, hexahedron, tetrahedron,
       scaleIcosa=10, scaleDodeca=8, scaleOcta=6, scaleHexa=4, scaleTetra=2}%

\def\psKeplerModel{\@ifnextchar[{\psKeplerModel@i}{\psKeplerModel@i[]}}
\def\psKeplerModel@i[#1]{%
  \@ifnextchar({\psKeplerModel@ii[#1]}{\psKeplerModel@ii[#1](0,0)}}
\def\psKeplerModel@ii[#1](#2){{%
  \psset{type=world,viewpoint=50 -20 30 rtp2xyz,Decran=50,#1}%
  \ifcase\hvSet@type
    \psKeplerWorld(#2)%
  \or
    \psKeplerTetrahedron(#2)%
  \or
    \psKeplerHexahedron(#2)%
  \or
    \psKeplerOctahedron(#2)%
  \or
    \psKeplerDodecahedron(#2)%
  \or
    \psKeplerIcosahedron(#2)%
  \fi
}}

\def\psKeplerTetrahedron(#1){%
  \rput(#1){%
  \codejps{
	/Gray4 {0.4 setgray} def
	/Gray7 {0.7 setgray} def
  	0.001 setlinewidth
	/tet1 {%
    	3 newtetraedre
    	dup 0.9 solidaffine
    	dup videsolid
    	dup (Gray4) outputcolors
    	dup (Gray7) inputcolors
	} def
	tet1 drawsolid**
  }}%
}

\def\psKeplerHexahedron(#1){%
  \rput(#1){%
  \codejps{
	/Gray5 {0.5 setgray} def
	/Gray7 {0.7 setgray} def
	0.001 setlinewidth
/cub1 {%
    3 newcube
    dup 0.9 solidaffine
    dup videsolid
    dup (Gray5) outputcolors
    dup (Gray7) inputcolors
} def
%    
cub1 drawsolid**
}}}

\def\psKeplerOctahedron(#1){%
  \rput(#1){%
  \codejps{%
	/Gray8 {0.8 setgray} def
	0.001 setlinewidth
/octa {%
    3 newoctaedre
    dup 0.9 solidaffine
    dup videsolid
    dup (Gray) outputcolors
    dup (Gray8) inputcolors
} def
octa drawsolid**
}}}


\def\psKeplerDodecahedron(#1){%
  \rput(#1){%
  \codejps{%
	/Gray6 {0.6 setgray} def
	/Gray8 {0.8 setgray} def
	0.001 setlinewidth
/dod1 {%
    3 newdodecaedre
    dup 0.9 solidaffine
    dup videsolid
    dup (Gray6) outputcolors
    dup (Gray8) inputcolors
} def
%
dod1    drawsolid**
}}}


\def\psKeplerIcosahedron(#1){%
  \rput(#1){%
  \codejps{%
	/Gray7 {0.7 setgray} def
	/Gray8 {0.8 setgray} def
	0.001 setlinewidth
/ico1 {%
    3 newicosaedre
    dup 0.9 solidaffine
    dup videsolid
    dup (Gray7) outputcolors
    dup (Gray8) inputcolors
} def
%
ico1    drawsolid**
}}}


\def\psKeplerWorld(#1){%
  \rput(#1){%
  \codejps{%
	/Paleblue {0.8 0.8 1 setrgbcolor} def
	/Gray1 {0.1 setgray} def
	/Gray2 {0.2 setgray} def
	/Gray3 {0.3 setgray} def
	/Gray4 {0.4 setgray} def
	/Gray5 {0.5 setgray} def
	/Gray6 {0.6 setgray} def
	/Gray7 {0.7 setgray} def
	/Gray8 {0.8 setgray} def
	/Gray9 {0.9 setgray} def
	0.001 setlinewidth
	/sph1 {%
    	6 -90 0 [12 36] newcalottespherecreuse
    	dup (Gray4) outputcolors
    	dup (Gray8) inputcolors
	} def
%
/cub1 {%
    6 newcube
    dup 0.9 solidaffine
    dup videsolid
    dup (Gray5) outputcolors
    dup (Gray7) inputcolors
} def
%
/sph2 {%
    3 -90 0 [12 36] newcalottespherecreuse
    dup (bleu) outputcolors
    dup (Paleblue) inputcolors
} def
%
/tet1 {%
    3 newtetraedre
    dup 0.9 solidaffine
    dup videsolid
    dup (Gray4) outputcolors
    dup (Gray7) inputcolors
} def
%
/sph3 {%
    1.5 -90 0 [6 18] newcalottespherecreuse
    dup (Gray5) outputcolors
    dup (Gray8) inputcolors
} def
%
/dod1 {%
    1.3 newdodecaedre
    dup 0.9 solidaffine
    dup videsolid
    dup (Gray6) outputcolors
    dup (Gray8) inputcolors
} def
%
/sph4 {%
    0.85 -90 0 [6 18] newcalottespherecreuse
    dup (Gray6) outputcolors
    dup (Gray9) inputcolors
} def
%
/ico1 {%
    0.7 newicosaedre
    dup 0.8 solidaffine
    dup videsolid
    dup (Gray7) outputcolors
    dup (Gray8) inputcolors
} def
%
/sph5 {%
    0.5 -90 0 [6 18] newcalottespherecreuse
    dup (Gray7) outputcolors
    dup (White) inputcolors
} def
%
/octa {%
    0.4 newoctaedre
    dup 0.8 solidaffine
    dup videsolid
    dup (Gray) outputcolors
    dup (Gray8) inputcolors
} def
%
/sph6 {%
    0.2 -90 0 [6 18] newcalottespherecreuse
    dup (Gray) outputcolors
    dup (White) inputcolors
} def
%
/sph7 {%
    0.125 [6 18] newsphere
    dup (jaune) outputcolors
    dup (White) inputcolors
} def
%
    sph1 cub1 solidfuz
    sph2 solidfuz
    tet1 solidfuz
    sph3 solidfuz
    dod1 solidfuz
    sph4 solidfuz
    ico1 solidfuz
    sph5 solidfuz
    octa solidfuz
    sph6 solidfuz
    sph7 solidfuz
    drawsolid**
}}}


\colorlet{Red}{red!30}
\colorlet{Green}{green!30}
\colorlet{Blue}{blue!30}
\colorlet{Magenta}{magenta!30}
\colorlet{Cyan}{cyan!30}



\def\psKeplerShow{\@ifnextchar[{\psKeplerShow@i}{\psKeplerShow@i[]}}
\def\psKeplerShow@i[#1]{%
  \@ifnextchar({\psKeplerShow@ii[#1]}{\psKeplerShow@ii[#1](0,0)}}
\def\psKeplerShow@ii[#1](#2){{%
  \psset{unit=1.5,viewpoint=50 -20 30 rtp2xyz,Decran=50,#1}%
  \rput(#2){%
    \begin{pspicture}(-6,-6)(8,6)
	\codejps{
		/tet {
		  \psk@kepler@scaleTetra newtetraedre
  		  dup 0.8 solidaffine
  		  dup videsolid
          dup (red) outputcolors
          dup (Red) inputcolors
        } def
		/hexa {%
  		  \psk@kepler@scaleHexa  newcube
          dup 0.8 solidaffine
  		  dup videsolid
          dup (green) outputcolors
          dup (Green) inputcolors
		} def
		/octa {%
  		  \psk@kepler@scaleOcta newoctaedre
		    dup 0.9 solidaffine
  			dup videsolid
  			dup (blue) outputcolors
  			dup (Blue) inputcolors
		} def
		/dodeca {%
	      \psk@kepler@scaleDodeca newdodecaedre
		    dup 0.9 solidaffine
    		dup videsolid
    		dup (magenta) outputcolors
    		dup (Magenta) inputcolors
		} def
		/icosa {%
    	  \psk@kepler@scaleIcosa newicosaedre
		  dup 0.9 solidaffine
    	  dup videsolid
    	  dup (cyan) outputcolors
	      dup (Cyan) inputcolors
	    } def
	  \ifPst@tetrahedron tet \fi
	  \ifPst@hexahedron hexa solidfuz \fi
	  \ifPst@octahedron octa solidfuz \fi
	  \ifPst@dodecahedron dodeca solidfuz \fi
	  \ifPst@icosahedron icosa solidfuz \fi
	   drawsolid**
 	}
	\end{pspicture}%
  }%
}}

% polexpr not needed as \PolDecToString is available as \xintDecToString
% \usepackage{xstring,siunitx}

\newcommand\FloatToFix[1]{\xintieval{[4]\xintfloatexpr\empty#1\relax}}
\newcommand\FloatEdef[2]{\edef#1{\xintfloatexpr\empty#2\relax}}

\def\FrameRate{36}%  36 frames per second
\def\Period{10}%    10 seconds for one revolution = 360 frames

\edef\TotalNumberOfFrames{\the\numexpr\FrameRate*\Period}%

\FloatEdef\DeltaEinitial{360/\TotalNumberOfFrames}%
\FloatEdef\Einitial{0}

\def\AreaFractionOfTotal{30}% Gibt an, wie groß der Flächenanteil x/360 an der Gesamtfläche ist

% Ellipse and Kepler law of areas

\xintdefvar Ox, Oy := 7, 2;                                     % Center coordinates of the ellipse
\xintdefvar a_el, b_el := 5, 3;                                 % demi axes
% ATTENTION THAT IN FRANCE ECCENTRICITY "eps_el" IS OFTEN NOTED WITH LETTER e
\xintdefvar e_el   := sqrt(a_el^2-b_el^2);                      % eccentricity
\xintdefvar F_1x, F_1y, F_2x, F_2y := Ox-e_el, Oy, Ox+e_el, Oy; % foci
\xintdefvar eps_el := e_el/a_el;                                % numeric eccentricity
% ATTENTION THAT IN FRANCE ECCENTRICITY "eps_el" IS OFTEN NOTED WITH LETTER e

\xintdeffloatfunc fK(E) := E - eps_el * oneRadian * sind(E);% oneRadian = 180/Pi in xinttrig.sty

% E <-- E - (f(E)- 360 t/T)/f'(E)
\xintNewFunction{EKepler}[1]{iter(#1;subs((abs(D) < 1e-6)?
                                          {break(@-D)}% stop now
                                          {@ - D} % @ is previous value, replace it with @ -D
                                         , D = (fK(@)- #1)/(1 - eps_el * cosd(@))
                                         )% end of substitution of D value
                                    , i = 1++)}

\xintdeffloatfunc Ellipse_E_x(E) := Ox+a_el*cosd(E);
\xintdeffloatfunc Ellipse_E_y(E) := Oy+b_el*sind(E);

% with angle phi: attention that "e_el" is not French excentricity but "c"
\xintdeffloatfunc r_Ellipse(f) := b_el^2/(a_el+e_el*cosd(f)) ;
%\xintdeffloatfunc Ellipse_x(f) := F_2x+cosd(f)*r_Ellipse(f);
\xintdeffloatfunc Ellipse_x(f) := F_2x+(b_el^2 - a_el*r_Ellipse(f))/e_el;
\xintdeffloatfunc Ellipse_y(f) := F_2y+sind(f)*r_Ellipse(f);
% Umrechnung des Winkels E in den Winkel phi
\xintdeffloatfunc E_to_phi_Pol(f) :=  pArgd(a_el*cosd(f)+Ox-F_2x, b_el*sind(f)+Oy-F_2y);

\xintdeffloatfunc phi2E_i(c, s, r) :=
      pArgd(b_el*(-Ox+F_2x+c*r),a_el*(-Oy+F_2y+s*r));
\xintdeffloatfunc phi2E(c, s) := phi2E_i(c, s, b_el^2/(a_el+e_el*c));
\xintdeffloatfunc phi_Pol_to_E(f) := phi2E(cosd(f),sind(f));

\edef\Fiix{\FloatToFix{F_2x}}
\edef\Fiiy{\FloatToFix{F_2y}}
\edef\Ox{\FloatToFix{Ox}}
\edef\Oy{\FloatToFix{Oy}}
\edef\ael{\FloatToFix{a_el}}
\edef\bel{\FloatToFix{b_el}}

\xintNewFunction {AreaEllipseS}[3]
  {0.5*(#2-#1<2*#3)?
       {(sqr(r_Ellipse(#1))+4*sqr(r_Ellipse((#1+#2)/2))+sqr(r_Ellipse(#2)))/6*(#2-#1)}
       {subs(subs(
                (sqr(r_Ellipse(#1))
                +4*add(sqr((r_Ellipse(#1+D*i))), i=1..[2]..(2*N))
                +2*add(sqr((r_Ellipse(#1+D*i))), i=2..[2]..(2*N-2))
                +sqr(r_Ellipse(#2))
                )*D/3 %% size of Simpson unit is of width 2D thus 2D/6 gives D/3
                , D=(#2-#1)/(2*N)
                )
             , N= ceil((#2-#1)/(2*#3))
             )
            }
      *oneDegree}

\FloatEdef\EKepler{EKepler(\Einitial)}
\FloatEdef\EKeplerI{EKepler(\Einitial+\AreaFractionOfTotal)}
\FloatEdef\EKeplerphi{E_to_phi_Pol(\EKepler)}
\FloatEdef\EKeplerIphi{E_to_phi_Pol(\EKeplerI)}

% In fact here we have increment of 30 degrees, which is 1/12 of a revolution.
% We thus know the area will be one twelfth of area of ellipse which is πab
% DO NOT USE \FloatToFix HERE IT CUTS DOWN DRASTICALLY THE PRECISION
\FloatEdef\AreaStart{
    %  spaces do not matter here...
    AreaEllipseS(
% NO !
%    \FloatToFix{E_to_phi_Pol(EKepler(\Einitial))},
    \EKeplerphi,
% No No !
%    \FloatToFix{E_to_phi_Pol(EKepler(\Einitial+30))},
% use E_to_phi_Pol(EKepler(\Einitial+30))
% or rather as we have already computed it
    \EKeplerIphi,
    %2
    1
    )
}%



\newcommand\CalculeEtAffichePoint[4][cyan]{%
%
  \FloatEdef\EKepler{EKepler(#2)}%
  \FloatEdef\EKeplerI{EKepler(#2+#3)}%
  %  \FloatEdef\EKeplerM{EKepler(#2+(#3/2))}%
% some precaution due to mod 360 things
  \FloatEdef\EKeplerphi{E_to_phi_Pol(\EKepler)}%
  \FloatEdef\EKeplerIphi{E_to_phi_Pol(\EKeplerI)}%
% WE ASSUME #3 POSITIVE
  \xintifboolexpr{\EKeplerIphi<\EKeplerphi}
      {\FloatEdef\EKeplerIphi{360+\EKeplerIphi}}
      {}%
  \FloatEdef\EKeplerMphi{(\EKeplerphi + \EKeplerIphi)/2}%
%%%%\immediate\write16{\xintfloateval{\EKepler, fK(\EKepler), \Einitial}}%
%
  \FloatEdef\PX{Ellipse_E_x(\EKepler)}%
  \FloatEdef\PY{Ellipse_E_y(\EKepler)}%
  \FloatEdef\PXI{Ellipse_E_x(\EKeplerI)}%
  \FloatEdef\PYI{Ellipse_E_y(\EKeplerI)}%
  \FloatEdef\PXM{Ellipse_x(\EKeplerMphi)}%
  \FloatEdef\PYM{Ellipse_y(\EKeplerMphi)}%
%
  \pscustom[fillstyle=solid,fillcolor=#1,opacity=0.25,linestyle=dashed]{%
    \pstEllipse(\Ox,\Oy)(\ael,\bel)[\FloatToFix{\EKepler}][\FloatToFix{\EKeplerI}]
    \lineto(\Fiix,\Fiiy)
    \closepath
    }%
%
% #4 is last parameter, it is to print some label
% But this is cheating, real thing would be to compute really
% the area !
  \uput{0.35}[\FloatToFix{\EKeplerMphi+180}]{0}(\FloatToFix{\PXM},\FloatToFix{\PYM}){#4}
% printing the mid-phi point would be counter-intuitive as it does
% not follow the celestial laws exactly
% I considered adding this but I remove finally
%  \psdot[linecolor=red,dotsize=6pt](\FloatToFix{\PXM}, \FloatToFix{\PYM})%
%
  \psline(\Fiix,\Fiiy)(\FloatToFix{\PX}, \FloatToFix{\PY})%
  \psline(\Fiix,\Fiiy)(\FloatToFix{\PXI}, \FloatToFix{\PYI})%
  \psdot[linecolor=magenta,dotsize=6pt](\FloatToFix{\PX}, \FloatToFix{\PY})%
  \psdot[linecolor=cyan,dotsize=6pt](\FloatToFix{\PXI}, \FloatToFix{\PYI})%
}%

% Computes by numerical integration of 0.5 r^2 dphi the area
% from phi angle #1 to phi angle #2 expressed in degrees
% with a step of #3 degrees between evaluation points
\xintNewFunction {AreaEllipseS}[3]
  {0.5*(#2-#1<2*#3)?
       {(sqr(r_Ellipse(#1))+4*sqr(r_Ellipse((#1+#2)/2))+sqr(r_Ellipse(#2)))/6*(#2-#1)}
       {subs(subs(
                (sqr(r_Ellipse(#1))
                +4*add(sqr((r_Ellipse(#1+D*i))), i=1..[2]..(2*N))
                +2*add(sqr((r_Ellipse(#1+D*i))), i=2..[2]..(2*N-2))
                +sqr(r_Ellipse(#2))
                )*D/3 %% size of Simpson unit is of width 2D thus 2D/6 gives D/3
                , D=(#2-#1)/(2*N)
                )
             , N= ceil((#2-#1)/(2*#3))
             )
            }
      *oneDegree}

\def\psKeplerAnimate{\@ifnextchar[{\psKeplerAnimate@i}{\psKeplerAnimate@i[]}}
\def\psKeplerAnimate@i[#1]{{%
  \if$#1$\else\psset{#1}\fi
  \begin{animateinline}[controls,loop,
    begin={\begin{pspicture}[showgrid](1,-2)(13,6)
        \psset{dotscale=0.8}\psset{PointSymbol=*,plotpoints=600,arrowscale=1.4,arrowinset=0.05,arrowlength=2}
            \pnode(\Ox,\Oy){O}
            \pnode(\Fiix,\Fiiy){F_2}
%            \uput{0.15}[-90]{0}(F_1){$F_1$}
            \pstEllipse[linecolor=black,linewidth=1pt](O)(\ael,\bel)
    },
    end={%
            \psdot[dotstyle=o,dotscale=1,fillcolor=orange](F_2)
            \uput{0.15}[-90]{0}(F_2){$S$}
            \pstEllipse[linecolor=purple,linestyle=none,arrows=->](O)(\ael,\bel)[260][270]
            \pstEllipse[linecolor=purple,linestyle=none,arrows=->](O)(\ael,\bel)[80][90]
    \end{pspicture}}]{20}% 20 frames/s (velocity of the animation)
\multiframe{\TotalNumberOfFrames}{rB=\FloatToFix{\Einitial}+\FloatToFix{\DeltaEinitial}}{% number of frames
\CalculeEtAffichePoint
   {\FloatToFix{\rB}}
   {\AreaFractionOfTotal}
% last parameter should be a box because it does again and again same expansion
   {$A_{1/2} \approx \num{\xintDecToString{\xintfloateval{[3]\AreaStart}}}$}%
%
% This is cheating... (but the area should be computed from the E values via E - ecc sin(E) up to a factor).
%
\CalculeEtAffichePoint[green]
    {\FloatToFix{\rB+120}}
    {\AreaFractionOfTotal}
    {$A_{3/4}\approx \num{\xintDecToString{\xintfloateval{[3]\AreaStart}}}$}%
%
}
\end{animateinline}%
}}

\catcode`\@=\PstAtCode\relax
%
\endinput 




