function[fig_num]=two_dof_frf_apply_initial(m,fn,omegan,omegad,damp,Q,fig_num,iu)
%
disp(' ');
disp(' Apply initial condition excitation? ');
disp(' 1=yes  2=no ');
init = input(' ');
%
if(init==1)
%
    iunit=iu;
    MST=Q';
%    
    clear d;
    clear v;
    disp(' ');
    if(iunit==1)
        d(1)=input(' Enter initial displacement of mass 1 (inch) ');
        d(2)=input(' Enter initial displacement of mass 2 (inch) ');
        disp(' ');
        v(1)=input(' Enter initial velocity of mass 1 (in/sec) ');
        v(2)=input(' Enter initial velocity of mass 2 (in/sec) ');
    else
        d(1)=input(' Enter initial displacement of mass 1 (m) ');
        d(2)=input(' Enter initial displacement of mass 2 (m) ');
        disp(' ');
        v(1)=input(' Enter initial velocity of mass 1 (m/sec) ');
        v(2)=input(' Enter initial velocity of mass 2 (m/sec) ');        
    end
%
    disp(' ');
    dur=input(' Enter duration (sec) ');
    disp(' ');
    disp(' Enter sample rate (samples/sec) '); 
    out1=sprintf(' (suggest > %9.5g)',20*fn(2,2));
    disp(out1);
    sr=input(' ');
    dt=1/sr;
%
    clear t;
    clear tt;
    clear nt;
    nt=round(dur/dt);
%
    clear DI;
    clear VI;
    clear nDI;
    clear nVI;
%
    DI=[d(1) d(2)];
    VI=[v(1) v(2)];
%
    sz=size(DI);
    if(sz(2)>sz(1))
       DI=DI';
    end
%
    sz=size(VI);
    if(sz(2)>sz(1))
      VI=VI';
    end
%
    nDI=MST*(m*DI);
    nVI=MST*(m*VI);
%
%
    clear x;
    clear v;
    clear a;
    clear A;
    clear B;
    clear om2;
    clear b1;
    clear b2;
    clear b3;
%
    om2=omegan.*omegan;
%
    A=nDI; 
    b1=omegan.*A;
    b2=damp.*b1;
    b3=b2+nVI;
%
    for(j=1:2)
        B(j)=b3(j)/omegad(j);
    end
%
    sz=size(A);
    if(sz(2)>sz(1))
        A=A';
    end
%
    sz=size(B);
    if(sz(2)>sz(1))
        B=B';
    end
    clear b1;
    clear b2;
    clear b3;
%
    om2=omegan.*omegan;
%
    A=nDI; 
    b1=omegan.*A;
    b2=damp.*b1;
    b3=b2+nVI;
%
    for(j=1:2)
        B(j)=b3(j)/omegad(j);
    end
%
    sz=size(A);
    if(sz(2)>sz(1))
        A=A';
    end
%
    sz=size(B);
    if(sz(2)>sz(1))
        B=B';
    end
%
    for(i=1:nt)
        t=dt*(i-1);
        tt(i)=t;
%
        arg=omegad*t;
%
        sz=size(arg);
        if(sz(2)>sz(1))
            arg=arg';
        end
%
        cos_arg=cos(arg);
        sin_arg=sin(arg);
        eee=exp(-damp.*omegan*t); 
%
        sz=size(cos_arg);
        if(sz(2)>sz(1))
            cos_arg=cos_arg';
        end    
%
        sz=size(sin_arg);
        if(sz(2)>sz(1))
            sin_arg=sin_arg';
        end  
%
        sz=size(eee);
        if(sz(2)>sz(1))
            eee=eee';
        end
%
        A_cos_arg=A.*cos_arg;
        A_sin_arg=A.*sin_arg;
        B_cos_arg=B.*cos_arg;    
        B_sin_arg=B.*sin_arg;   
%
        term1=A_cos_arg+B_sin_arg;
        n=eee.*term1;
%
        sz=size(n);
        if(sz(2)>sz(1))
            n=n';
        end      
%
        term2=-A_sin_arg+B_cos_arg;
%
        nv=-damp.*(omegan.*n) + omegad.*( eee.*term2 );
%
        na=-2*damp.*(omegan.*nv) - om2.*n;    
%
        x(i,:)=Q(:,1:2)*n(1:2);
        v(i,:)=Q(:,1:2)*nv(1:2);
        a(i,:)=Q(:,1:2)*na(1:2);  
%
    end
%
        figure(fig_num);
        fig_num=fig_num+1;
        plot(tt,x(:,1),tt,x(:,2));
        grid on;
        xlabel('Time(sec)');
        title('Displacement');
        if(iu==1)
            ylabel('Disp(in)');
        else
            ylabel('Disp(m)');    
        end
%
        figure(fig_num);
        fig_num=fig_num+1;
        plot(tt,v(:,1),tt,v(:,2));
        grid on;
        xlabel('Time(sec)');
        title('Velocity');
        if(iu==1)
            ylabel('Vel(in/sec)');
        else
            ylabel('Vel(m/sec)');    
        end
%
        figure(fig_num);
        fig_num=fig_num+1;
        if(iu==1)
            a(:,1)=a(:,1)/386;
            a(:,2)=a(:,2)/386;
        end
        plot(tt,a(:,1),tt,a(:,2));
        grid on;
        xlabel('Time(sec)');
        title('Acceleration');
        if(iu==1)
            ylabel('Accel(G)');
        else
            ylabel('Accel(m/sec^2)');    
        end
%
end