% gillCoin.m
%*******************************************************
%         Rolling Coin  (GILL)    R-K-G法
%*******************************************************
% Copyright (C) 2007 清水 信行, 曽我部 潔 All Rights Reserved.

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

function [xa,cq4]=gillCoin(t, x0, r, m, g, rho, k, h, cq0);
%-------------------------------------------
% ルンゲ・クッタ・ギルによる数値積分計算
% 呼び出される関数
% subCoin.m : 転がり運動するコインのDAE(微分代数方程式)の定義
%-------------------------------------------

[hsdot]=subCoin(t, x0, r, m, g, rho, k );
for j=1:10;
  xk1=h*hsdot(j);
  r1=(xk1-2*cq0(j))/2;
  y1(j)=x0(j)+r1;
  cq1(j)=cq0(j)+3*r1-xk1/2;
end;

[hsdot]=subCoin(t+h/2, y1, r, m, g, rho, k );
for j=1:10;
  xk2=h*hsdot(j);
  r2=0.292893219*(xk2-cq1(j));
  y2(j)=y1(j)+r2;
  cq2(j)=cq1(j)+3*r2-0.292893219*xk2;
end;

[hsdot]=subCoin(t+h/2, y2,r, m, g, rho, k );
for j=1:10;
  xk3=h*hsdot(j);
  r3=1.707106781*(xk3-cq2(j));
  y3(j)=y2(j)+r3;
  cq3(j)=cq2(j)+3*r3-1.707106781*xk3;
end;

[hsdot]=subCoin(t+h, y3, r, m, g, rho, k );
for j=1:10;
  xk4=h*hsdot(j);
  r4=(xk4-2*cq3(j))/6;
  xa(j)=y3(j)+r4;
  cq4(j)=cq3(j)+3*r4-xk4/2;
end;