
%   mdof_base_function.m

function[acc,ra,rd]=mdof_base_function(Q,omegan,dm,dt,dof,nt,yy,base)

    a=zeros(nt,dof);
    x=zeros(nt,dof);
    acc=zeros(nt,dof);   
    ra=zeros(nt,dof);
    rd=zeros(nt,dof);

    for j=1:dof

                omegad=omegan(j)*sqrt(1.-(dm(j)^2));
  
                cosd=cos(omegad*dt);
                sind=sin(omegad*dt);

                domegadt=dm(j)*omegan(j)*dt;

                eee1=exp(-domegadt);
                eee2=exp(-2.*domegadt);

                ecosd=eee1*cosd;
                esind=eee1*sind; 

                a1= 2.*ecosd;
                a2=-eee2;    

                omeganT=omegan(j)*dt;

                phi=2*(dm(j))^2-1;
                DD1=(omegan(j)/omegad)*phi;
                DD2=2*DD1;
    
                df1=2*dm(j)*(ecosd-1) +DD1*esind +omeganT;
                df2=-2*omeganT*ecosd +2*dm(j)*(1-eee2) -DD2*esind;
                df3=(2*dm(j)+omeganT)*eee2 +(DD1*esind-2*dm(j)*ecosd);

                MD=omegan(j)^3*dt;
                df1=df1/MD;
                df2=df2/MD;
                df3=df3/MD;

                af1=esind/(omegad*dt);
                af2=-2*af1;
                af3=af1;

                yj=yy(j,:);
  
                %   Relative Acceleration Response
 
                forward=[ af1,  af2,  af3 ];    
                back   =[     1, -a1, -a2 ];    
  
                a(:,j)=filter(forward,back,yj);

               %   Relative Displacement Response   
               
               rd_forward=[ df1,  df2,  df3 ];    
               rd_back   =[     1, -a1, -a2 ];     
  
               x(:,j)=filter(rd_forward,rd_back,yj);

    end

    for i=1:nt 
                rd(i,:)=Q*x(i,:)';
                ra(i,:)=Q*a(i,:)';  
                acc(i,:)=ra(i,:);
    end

    for i=1:dof
                acc(:,i)=acc(:,i)+base;
    end
