disp(' ')
disp(' sineforce_ode45.m ')
disp(' ver 1.0  December 16, 2007 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the response of a  ')
disp(' single-degree-of-freedom system to a sinusoidal ')
disp(' force excitation using the ode45 method. ')
disp(' ')
disp(' The ode45 method is a numerical method.')
disp(' ')
%
clear omegan;
clear omegad;
clear fn;
clear damp;
clear N;
clear t;
clear x;
clear w;
%
disp(' select units')
iunits=input(' 1=English  2=metric ');
%
if(iunits==1)
    disp(' ')
    mass=input(' input mass (lbm) ');
    mass=mass/386.;
    disp(' ')
    stiffness=input(' input stiffness (lbf/in) ');    
    disp(' ')
    v0 = input(' input initial velocity (in/sec) ');
    disp(' ')
    d0 = input(' input initial displacement (in) ');
else
    disp(' ') 
    mass=input(' input mass (kg) ');
    disp(' ')
    stiffness=input(' input stiffness (N/m) ');    
    disp(' ')
    v0 = input(' input initial velocity (m/sec) ');
    disp(' ')
    d0 = input(' input initial displacement (m) ');    
end
%
%% disp(' ')
%% fn=input(' Enter the natural frequency (Hz)  ');
%
disp(' ')
idamp=input(' Enter damping format:  1= damping ratio   2= Q  ');	
%
disp(' ')
if(idamp==1)
    damp=input(' Enter damping ratio (typically 0.05) ');
else
    Q=input(' Enter the amplification factor (typically Q=10) ');
    damp=1./(2.*Q);
end
%
disp(' ')
f=input(' Enter the forcing frequency (Hz)  ');
%
if(iunits==1)
    disp(' ')
    A = input(' input the force amplitude (lbf) ');
else
    disp(' ')
    A = input(' input the force amplitude (N) ');   
end
%
omegan=sqrt(stiffness/mass);
fn=omegan/(2.*pi);
disp(' ');
out1=sprintf(' fn = %8.4g Hz \n',fn);
disp(out1);
%
omega=2.*pi*f;
%
omegad=omegan*sqrt(1.-damp^2);
%
%   State Space Solution
%
if(fn>f)
   period=1/fn;
else
   period=1./f;
end
%
dt=period/64;
sr=1/dt;
dur=30*period;
N=dur/dt;
%
t=linspace(0,dur,N);
t=t';
%
w(1)=d0;
w(2)=v0;
%
clear b;
clear c;
%
A=A/mass;
%
b=-omegan^2;
c=-2.*damp*omegan;
domegan=damp*omegan;
%
[t,w] = ode45(@sineforce_function,t,w,[],[],b,c,A,omega);
%
disp(' ');
disp(' plot displacement?');
k=input(' 1=yes  2=no ');
if(k==1)
    figure(1);
    plot(t,w(:,1));
    if(iunits==1)
        ylabel('Displacement (in)');
    else
        ylabel('Displacement (m)');
    end
    xlabel('Time (sec)');
end
%        
disp(' ');
disp(' plot velocity?');
k=input(' 1=yes  2=no ');
if(k==1)
    figure(2);
    plot(t,w(:,2));
    if(iunits==1)
        ylabel('Velocity (in/sec)');
    else
        ylabel('Velocity (m/sec)');
    end
    xlabel('Time (sec)');
end
%        
disp(' ');
disp(' plot acceleration?');
k=input(' 1=yes  2=no ');
if(k==1)
    clear acc;
%
    for(i=1:length(t))
        acc(i)=A*sin(omega*t(i))+c*w(i,2)+b*w(i,1);
    end
%    
    figure(3);
    if(iunits==1)
        acc=acc/386.;
    end    
    plot(t,acc);
    if(iunits==1)
        ylabel('Acceleration (G)');
    else
        ylabel('Acceleration (in/sec^2)');
    end
    xlabel('Time (sec)');
end
%
disp(' ');
disp(' plot steady state transfer function magnitude (disp/force) ?');
k=input(' 1=yes  2=no ');
if(k==1)
    clear omega;
    clear mag;
    clear phase;
    clear trans;
    omega=linspace(0,  2*pi*fn*10, 1000 );
    figure(4);
%
    n1=0; 
    n2=0; 
    n3=1;
    d1=1;
    d2=2.*domegan; 
    d3=omegan^2;
%
    [mag,phase]=bode_transfer_function(n1,n2,n3,d1,d2,d3,omega);
%
    omega=omega/(2.*pi);
    mag=mag/mass;
    plot( omega,mag);
    grid;
%    
    title('Transfer Function Magnitude');
    if(iunits==1)
        ylabel('Disp/Force (in/lbf)');     
    else
        ylabel('Disp/Force (m/N)');       
    end   
    xlabel('Frequency (Hz)');    
%
    disp(' plot steady state transfer function phase angle?');
    kk=input(' 1=yes  2=no ');
    if(kk==1)
        figure(5);
        plot( omega,phase);
        grid;
%    
        title(' ');
        ylabel('Phase Angle (deg)');     
        xlabel('Frequency (Hz)');    
    end
end