%   Runge-Kutta4 Method
%   NONLinear Case for Wen Model
% My God it is OK!!!!! 04.05.28
%**********************************************************************
%*                   RESPONSE CALCURATION PROGRAM                     *
%*                                 OF                                 *
%*   NONLINEAR RESTORING FORCE SYSTEM WITH MULTI DEGREES OF FREEDOM   *
%*                      MODEL : VERSATILE MODEL                       *  
%*                      METHOD:RINGE-KUTTA                            *      _
%*                                     FILE     : RungeKutta.FOR      *     |_|  3rd
%*                                     BASEDATA : ELCENT.DAT          *      |    
%*                                                                    *      _ 
%*                       CODE BY TANG HESHENG                         *     |_|  2nd
%*                             2004.1.29                              *      |
%**********************************************************************      _
%**********************************************************************     |_|  1st
%*                   MODEL AND PARAMETERS SETTING                     *      |                                                                         *
%*   DOT2{X(I)}+2*H(I)*W(I)*U(I)+Z(I)/M(I)-(1-DELTin)*M(I+1)/M(I)*    *   --------
%*      [2*H(I+1)*W(I+1)*DOT{U(I+1)}+Z(I+1)/M(I+1)]=-DOT2(Xg)         *   ////////
%*                                                                    *     Model    
%*                        U(I)=X(I)-X(I-1)                            *
%*                     DELTin=1,i=n,ELSE DELTin=0                     *
%*                                                                    *
%*        DOT{Z}=K*DOT{U}-ALPHA*ABS(DOT{U})*ABS{Z}**(N-1)             *
%*                       -BETA*DOT{U}*ABS{Z}**N                       *
%*                                                                    *
%*         M     :MASS                                                *
%*         N     :VERSATILE MODEL PARAMETER                           *
%*        ALPHA  :VERSATILE MODEL PARAMETER                           *
%*        BETA   :VERSATILE MODEL PARAMETER                           *
%*         H     :EQUALITY OF ATTENUATION PARAMETER                   *
%*                     H(I)=C(I)/2*SQRT(M(I)*K(I))                    *
%*         W     :EQUALITY OF NATURAL RADIAN FREQUENCY                *
%*                        W(I)=SQRT(K(I)/M(I))                        *
%*         G     :EARTHQUAKE INPUT  1 X NSTEP vector                  *
%**********************************************************************


function [DISP,VEL,ACC,Z]=KUTTA(M,K,C,NDOF,ALPHA,BETA,N,G,NSTEP,DT)
H=C./(2*sqrt(M.*K));
OMIGA=sqrt(K./M);
ACC=zeros(NDOF,NSTEP); % absolute value
VEL=zeros(NDOF,NSTEP);
DISP=zeros(NDOF,NSTEP);
ACC=zeros(NDOF,NSTEP);
Z=zeros(NDOF,NSTEP);

for IT=1:NSTEP-1
 %----------   0   ----------------   
    U0(1)=VEL(1,IT);
    for I=2:NDOF
        U0(I)=VEL(I,IT)-VEL(I-1,IT);
    end
    for L=1:NDOF
        X0(L)=VEL(L,IT)*DT;
       if L==NDOF
           Y0(L)=(-2*H(L)*OMIGA(L)*U0(L)-Z(L,IT)/M(L)-G(IT))*DT;
       else
           Y0(L)=(-2*H(L)*OMIGA(L)*U0(L)-Z(L,IT)/M(L)...
            +M(L+1)/M(L)*(2*H(L+1)*OMIGA(L+1)*U0(L+1)+Z(L+1,IT)/M(L+1))-G(IT))*DT; 
       end 
        Z0(L)=(M(L)*OMIGA(L)^2*U0(L)-ALPHA(L)*abs(U0(L))*abs(Z(L,IT))^(N(L)-1)*Z(L,IT)...
            -BETA(L)*U0(L)*abs(Z(L,IT))^N(L))*DT;
    end
 %----------   1   ----------------   
    U1(1)=U0(1)+Y0(1)*0.5;
    for I=2:NDOF
        U1(I)=U0(I)+Y0(I)*0.5-Y0(I-1)*0.5;
    end
    for L=1:NDOF
        X1(L)=(VEL(L,IT)+Y0(L)*0.5)*DT;
       if L==NDOF
           Y1(L)=(-2*H(L)*OMIGA(L)*U1(L)-(Z(L,IT)+0.5*Z0(L))/M(L)-(G(IT)+G(IT+1))*0.5)*DT;
       else
           Y1(L)=(-2*H(L)*OMIGA(L)*U1(L)-(Z(L,IT)+0.5*Z0(L))/M(L)...
            +M(L+1)/M(L)*(2*H(L+1)*OMIGA(L+1)*U1(L+1)+(Z(L+1,IT)+0.5*Z0(L+1))/M(L+1))-(G(IT)+G(IT+1))*0.5)*DT; 
       end 
        Z1(L)=(M(L)*OMIGA(L)^2*U1(L)-ALPHA(L)*abs(U1(L))*abs((Z(L,IT)+0.5*Z0(L)))^(N(L)-1)*(Z(L,IT)+0.5*Z0(L))...
            -BETA(L)*U1(L)*abs((Z(L,IT)+0.5*Z0(L)))^N(L))*DT;
    end
 %----------   2   ----------------   
    U2(1)=U0(1)+Y1(1)*0.5;
    for I=2:NDOF
        U2(I)=U0(I)+Y1(I)*0.5-Y1(I-1)*0.5;
    end
    for L=1:NDOF
        X2(L)=(VEL(L,IT)+Y1(L)*0.5)*DT;
       if L==NDOF
           Y2(L)=(-2*H(L)*OMIGA(L)*U2(L)-(Z(L,IT)+0.5*Z1(L))/M(L)-(G(IT)+G(IT+1))*0.5)*DT;
       else
           Y2(L)=(-2*H(L)*OMIGA(L)*U2(L)-(Z(L,IT)+0.5*Z1(L))/M(L)...
            +M(L+1)/M(L)*(2*H(L+1)*OMIGA(L+1)*U2(L+1)+(Z(L+1,IT)+0.5*Z1(L+1))/M(L+1))-(G(IT)+G(IT+1))*0.5)*DT; 
       end 
        Z2(L)=(M(L)*OMIGA(L)^2*U2(L)-ALPHA(L)*abs(U2(L))*abs((Z(L,IT)+0.5*Z1(L)))^(N(L)-1)*(Z(L,IT)+0.5*Z1(L))...
            -BETA(L)*U2(L)*abs((Z(L,IT)+0.5*Z1(L)))^N(L))*DT;
    end
 %----------   3   ----------------   
    U3(1)=U0(1)+Y2(1);
    for I=2:NDOF
        U3(I)=U0(I)+Y2(I)-Y2(I-1);
    end
    for L=1:NDOF
        X3(L)=(VEL(L,IT)+Y2(L))*DT;
       if L==NDOF
           Y3(L)=(-2*H(L)*OMIGA(L)*U3(L)-(Z(L,IT)+Z2(L))/M(L)-G(IT+1))*DT;
       else
           Y3(L)=(-2*H(L)*OMIGA(L)*U3(L)-(Z(L,IT)+Z2(L))/M(L)...
            +M(L+1)/M(L)*(2*H(L+1)*OMIGA(L+1)*U3(L+1)+(Z(L+1,IT)+Z2(L+1))/M(L+1))-G(IT+1))*DT; 
       end 
        Z3(L)=(M(L)*OMIGA(L)^2*U3(L)-ALPHA(L)*abs(U3(L))*abs((Z(L,IT)+Z2(L)))^(N(L)-1)*(Z(L,IT)+Z2(L))...
            -BETA(L)*U3(L)*abs((Z(L,IT)+Z2(L)))^N(L))*DT;
    end
%-------------- result -------------
    for L=1:NDOF
        DISP(L,IT+1)=DISP(L,IT)+(X0(L)+2*X1(L)+2*X2(L)+X3(L))/6;
        VEL(L,IT+1)=VEL(L,IT)+(Y0(L)+2*Y1(L)+2*Y2(L)+Y3(L))/6;
        Z(L,IT+1)=Z(L,IT)+(Z0(L)+2*Z1(L)+2*Z2(L)+Z3(L))/6;
    end
%-------------- ABSOLUTE ACC -------------
    UU(1)=VEL(1,IT+1);
    for I=2:NDOF
        UU(I)=VEL(I,IT+1)-VEL(I-1,IT+1);
    end
    for L=1:NDOF
       if L==NDOF
           ACC(L,IT+1)=-2*H(L)*OMIGA(L)*UU(L)-(Z(L,IT+1))/M(L);
       else
           ACC(L,IT+1)=-2*H(L)*OMIGA(L)*UU(L)-(Z(L,IT+1))/M(L)...
            +M(L+1)/M(L)*(2*H(L+1)*OMIGA(L+1)*UU(L+1)+Z(L+1,IT+1)/M(L+1)); 
       end 
    end
end