function DY=top3d_elim_ode(t,Y,fa,c,cJd,invcJd,rGCdtilde,rGCd);
% function DY=top3d_elim_ode(t,Y,fa,c,cJd,invcJd,rGCdtilde,rGCd);
%
%  ■　(g) 微分方程式の定義
%
% 数値積分ルーチン ode45 を用いて、数値積分を行う微分方程式を定義する 
% 式(4.14)(4.31) 参照
%
% Copyright (C) 2007 清水 信行, 曽我部 潔, 竹原 昭一郎, 井上 剛志 All Rights Reserved.

% =================================================================================
% プログラムの使用条件
% 1. 本プログラムおよびデータを商用で，頒布してはならない．
% 2. 本プログラムおよびデータを著作者の許可なく営利目的で利用してはならない．
% 3. 本プログラムおよびデータの改変は，営利目的でない限り自由である．
% 4. 本プログラムおよびデータを使用することにより生じた損害などについては，著者，
%    プロクラム著作者，(社)日本機械学会，および(株)コロナ社は一切の責任を負わない．
% 5. 本プログラムの著作権は，本書編著者ならびにプログラム著作者に帰す．
% =================================================================================

               % オイラー角は3-1-3（Z'X'Z'）変換
thOG=Y(1:3);   % 状態変数から オイラー角(φ,θ,ψ)を抽出
omOGd=Y(4:6);  % 状態変数から 角速度(ωx',ωy',ωz')を抽出
C=cos(thOG);   % cos(φ),cos(θ),cos(ψ)
S=sin(thOG);   % sin(φ),sin(θ),sin(ψ)
nd=-c*omOGd;   % 回転に対する外トルク（この場合は減衰による外トルク）
%
AOG=[ C(1)*C(3)-S(1)*C(2)*S(3)  -C(1)*S(3)-S(1)*C(2)*C(3)   S(1)*S(2);
      S(1)*C(3)+C(1)*C(2)*S(3)  -S(1)*S(3)+C(1)*C(2)*C(3)  -C(1)*S(2);% 式(4.35)
      S(2)*S(3)                  S(2)*C(3)                  C(2)];% 座標変換マトリックス
%
invGed=1/S(2)*[S(3)        C(3)          0;    %オイラー角の時間微分を
               S(2)*C(3)  -S(2)*S(3)     0;    %求めるための角速度係数行列
              -C(2)*S(3)  -C(2)*C(3)     S(2)];%
%
omOGdtilde=tilde(omOGd);      % 角速度omOGdから作った歪み対称マトリックス
%
%　オイラー角の微分と角速度の関係（式(4.14)）
DthOG=invGed*omOGd;
%
%　角速度の時間微分について（式(4.31)）
DomOGd=invcJd*(nd-rGCdtilde*AOG'*fa-omOGdtilde*cJd*omOGd);
%
%　状態変数の時間微分
DY=[DthOG;DomOGd];            % ode45で数値積分する微分方程式
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%