%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                           %%
%% This is file `pst-nutation.tex',                                          %%
%%                                                                           %%
%% IMPORTANT NOTICE:                                                         %%
%%                                                                           %%
%% Package `pst-nutation'                                                    %%
%%                                                                           %%
%% Manuel Luque, Herbert Voß                                                 %%
%%                                                                           %%
%% February 04, 2022                                                         %%
%%                                                                           %%
%% 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-nutation' is a PSTricks tools                                       %%
%%  Package to illustrate the concepts of rotation, precession and nutation. %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\csname PSTNutation\endcsname
\let\PSTNutation\endinput

% Require PSTricks
\ifx\PSTricksLoaded\endinput\else\input pstricks.tex\fi
\ifx\PSTSOLIDESIIIDLoaded\endinput\else\input pst-solides3d.tex\fi
\ifx\PSTXKeyLoaded\endinput\else\input pst-xkey\fi

\def\fileversion{0.01}
\def\filedate{2025/12/20}
\message{`PST-nutation' v\fileversion, \filedate\space (ml,hv)}

\edef\PstAtCode{\the\catcode`\@} \catcode`\@=11\relax

\pstheader{pst-nutation.pro}

\addtosolideslistobject{parallel,meridian}
%% la macro \psRotationIIID
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rotation autour d'un axe passant par O
\define@key[psset]{pst-solides3d}{angle}{\def\psk@solides@angle{#1 }}% rotation autour de l'axe en degres
\psset[pst-solides3d]{angle=0}
%
\define@boolkey[psset]{pst-solides3d}[Pst@]{quaternion}[true]{}
\psset{quaternion=false}%
%\psRotationIIID(x,y,z)
% les coordonnées du vecteur axe sont données par axe=ux uy uz
\def\psRotationIIID{\pst@object{psRotationIIID}}
\def\psRotationIIID@i(#1,#2,#3){{%
  \pst@killglue%
  \use@par%
\addto@pscode{
    \tx@optionssolides
    SolidesDict begin
    /Angle \psk@solides@angle def
    /cosT {Angle cos} bind def
    /sinT {Angle sin} bind def
\ifPst@quaternion
    /Angle_2 {Angle 2 div} bind def
    /p_ {Angle_2 cos} bind def
    /sinA_2 {Angle_2 sin} bind def
    axe unitaire3d
    /u3 exch sinA_2 mul def /u2 exch sinA_2 mul def /u1 exch sinA_2 mul def
    /t2  {p_ u1 mul} bind def
    /t3  {p_ u2 mul} bind def
    /t4  {p_ u3 mul} bind def
    /t5  {u1 dup mul neg} bind def
    /t6  {u1 u2 mul} bind def
    /t7  {u1 u3 mul} bind def
    /t8  {u2 dup mul neg} bind def
    /t9  {u2 u3 mul} bind def
    /t10 {u3 dup mul neg} bind def
/Rotation3D {
4 dict begin
   /M defpoint3d % on récupère les coordonnées
   M /z exch def
     /y exch def
     /x exch def
 t8 t10 add x mul t6 t4 sub y mul add t3 t7 add z mul add 2 mul x add % x'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 t4 t6 add x mul t5 t10 add y mul add t9 t2 sub z mul add 2 mul y add % y'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 t7 t3 sub x mul t2 t9 add y mul add t5 t8 add z mul add 2 mul z add  % z'
end} def
\else
   axe unitaire3d
    /u3 exch def /u2 exch def /u1 exch def
/Rotation3D {
5 dict begin
   /M defpoint3d % on récupère les coordonnées
   M /z exch def
     /y exch def
     /x exch def
     /S u1 x mul u2 y mul add u3 z mul add def
   u1 S mul 1 cosT sub mul x cosT mul add u2 z mul u3 y mul sub sinT mul add % x'
   u2 S mul 1 cosT sub mul y cosT mul add u3 x mul u1 z mul sub sinT mul add % y'
   u3 S mul 1 cosT sub mul z cosT mul add u1 y mul u2 x mul sub sinT mul add % z'
end} def
\fi
    end}%
\psSolid[transform=Rotation3D](#1,#2,#3)
}\ignorespaces%
}

\pst@addfams{pst-nutation}

\define@key[psset]{pst-nutation}{theta}{\edef\pst@nutation@theta{#1 }}
\psset[pst-nutation]{theta=-30}% inclination of the axis with reference to Oz
\define@key[psset]{pst-nutation}{dtheta}{\edef\pst@nutation@dtheta{#1 }}
\psset[pst-nutation]{dtheta=5}% amplitude the variation of the inclination
\define@key[psset]{pst-nutation}{alpha}{\edef\pst@nutation@alpha{#1 }}
\psset[pst-nutation]{alpha=0}% rotation angle about the axis Oz in degrees (see documentation)
\define@key[psset]{pst-nutation}{omega}{\edef\pst@nutation@omega{#1 }}
\psset[pst-nutation]{omega=0}% rotation angle about the axis in degrees (see documentation)
\define@key[psset]{pst-nutation}{frequency}{\edef\pst@nutation@frequency{#1 }}
\psset[pst-nutation]{frequency=6}% number of axis oscillations in one revolution (axis oscillation frequency)
\define@boolkey[psset]{pst-nutation}[Pst@nutation@]{LargeSphere}[true]{}
\psset[pst-nutation]{LargeSphere=true}

\def\psNutation{\def\pst@par{}\pst@object{psNutation}}
\def\psNutation@i{%
 \addbefore@par{solidmemory,lightsrc=,resolution=360,linecolor=red}%
    \begin@SpecialObj
\pstVerb{%
         /pi 3.14159 def
         /rd {180 pi div mul} def        %% transforme des rd en degres
         /deg {pi mul 180 div} def       %% transforme des degres en rd
         /inclinaison \pst@nutation@theta def %% par rapport à Oz (colatitude)
         /dinclinaison \pst@nutation@dtheta def % amplitude des ondulations
         /inclinaisonmax 90 inclinaison dinclinaison add sub def
         /inclinaisonmin 90 inclinaison dinclinaison sub sub def
         /nS \pst@nutation@frequency def % nombre d'ondulations sur un tour
         /rotAxis \pst@nutation@omega def
         /rotZ \pst@nutation@alpha deg def
         /inclinaisonrad inclinaison deg def
         /dinclinaisonrad dinclinaison deg def
}%
\pscircle[linecolor=\pslinecolor,linewidth=0.05]{2}%
%  M alpha_x alpha_y alpha_z rotateOpoint3d --> M'
\codejps{/AXE {0 0 1 inclinaison dinclinaison nS \pst@nutation@alpha mul sin mul add 0 \pst@nutation@alpha rotateOpoint3d } def
         AXE /zaxe exch 2 mul def /yaxe exch 2 mul def /xaxe exch 2 mul def}%
\psSolid[object=vecteur,
        definition={[.01 .01]},
        linecolor=orange,linewidth=0.05,
        args=xaxe neg yaxe neg  zaxe neg](xaxe neg,yaxe neg,zaxe neg)
\defFunction[algebraic]{nutationS}(t)%
{-4*sin(t)*sin(inclinaisonrad + dinclinaisonrad*sin(nS*t))}%
{4*cos(t)*sin(inclinaisonrad + dinclinaisonrad*sin(nS*t))}%
{-4*cos(inclinaisonrad + dinclinaisonrad*sin(nS*t))}%
\psSolid[object=courbe,r=0,range=0 rotZ,function=nutationS,fillcolor=orange,linewidth=0.05]%
\psSolid[object=sphere,r=2,ngrid=9 18,RotX=inclinaison dinclinaison nS \pst@nutation@alpha mul sin mul add,RotZ=\pst@nutation@alpha,name=S,action=none]%
\psRotationIIID[object=load,load=S,axe=AXE,angle=\pst@nutation@omega,linecolor=\pslinecolor](0,0,0)
\psSolid[object=vecteur,
        linecolor=orange,linewidth=0.05,
        args=xaxe yaxe zaxe](xaxe,yaxe,zaxe)
\ifPst@nutation@LargeSphere
\pscircle[linecolor=gray]{4}
\multido{\i=-90+30}{6}{\psSolid[object=parallel,r=4,phi=\i,linecolor=gray]}%
\multido{\i=0+30}{12}{\psSolid[object=meridian,r=4,theta=\i,linecolor=gray]}%
\psSolid[object=parallel,r=4,phi=inclinaisonmax,linecolor={[rgb]{0.5 0.5 1}}]
\psSolid[object=parallel,r=4,phi=inclinaisonmin,linecolor={[rgb]{0.5 0.5 1}}]
\fi
\defFunction[algebraic]{nutationN}(t)%
  {4*sin(t)*sin(inclinaisonrad + dinclinaisonrad*sin(nS*t))}%
  {-4*cos(t)*sin(inclinaisonrad + dinclinaisonrad*sin(nS*t))}%
  {4*cos(inclinaisonrad + dinclinaisonrad*sin(nS*t))}%
\psSolid[object=courbe,r=0,range=0 rotZ,function=nutationN,fillcolor=orange,linewidth=0.05]%
\end@SpecialObj
  \ignorespaces}%
\catcode`\@=\PstAtCode\relax
\endinput