function sol=Newmark(M,C,KL,eq,F,init,period,h)
% function sol=Newmark(M,C,KL,eq,F,init,period,h)
%
% 方程式　M*xdd+C*xd)+KL*x+ｆint(x)＝F　を時間積分します
% t=0から刻み時間hでt=periodまで積分する。
% 初期値はinitで,ｆint(x)はeqで与える。　
% 
% Copyright (C) 2007 清水 信行, 曽我部 潔, 藤川 猛 All Rights Reserved.

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

Fint=char(eq(1));                   %非線形復元力
KN=char(eq(2));                     %接線係数（ばね定数）
sol=[];hh=h*h;
eps=1e-6;                           %ニュートン法収束許容地
  D=init(1);D1=D;                   %初期変位
  V=init(2);                        %初期速度
  A=inv(M)*(F-C*V-KL*D-eval(Fint)); %初期加速度
t=0 ; 
while t < period 
% 状態量の保存 
    D1=D; SF=KL*D1+eval(Fint);
    sol=[sol, [t;D;V;SF]];          % 状態量の保存 
% tn+hにおける近似解の設定
    A1=A;
    V1=V+h*A;
    D1=D+h*V+0.5*hh*A;
% ニュートンラプソン法による近似解の改善 
for inew=1:10;
    G=M*A1+C*V1+KL*D1+eval(Fint)-F; %誤差ベクトル
    if abs(G)<eps ;break ;end
    matJ=M+0.5*h*C+0.25*hh*(KL+eval(KN));  %ヤコビアンマトリックス
    dA=-inv(matJ)*G;                %solve dA  
    A1=A1+dA;                       %修正
    V1=V1+0.5*h*dA;
    D1=D1+0.25*hh*dA;
end
  A=A1;  V=V1;  D=D1;               %更新 
  t=t+h; 
end 
