disp(' ')
disp(' pendulum_free_ode45.m ')
disp(' ver 1.1  July 28, 2006 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the free vibration of a  ')
disp(' simple pendulum using the ode45 method.  ')
disp(' The pendulum may have large angular displacement. ')
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 angle units')
aunits=input(' 1=degrees  2=radians ');
%
if(aunits==1)
    disp(' ')
    v0 = input(' input initial angular velocity (deg/sec) ');
    disp(' ')
    d0 = input(' input initial angular displacement (deg) ');
    v0=v0*pi/180;
    d0=d0*pi/180;
else
    disp(' ')
    v0 = input(' input initial angular velocity (rad/sec) ');
    disp(' ')
    d0 = input(' input initial angular displacement (rad) ');
end
%
disp(' ');
disp(' Select length units')
iunits=input(' 1=English  2=metric ');
%
if(iunits==1)
    disp(' ')
    G = 386.;
    L = input(' input pendulum length (in) ');
else
    disp(' ')
    G = 9.81;
    L = input(' input pendulum length (m) ');    
end
%
omegan=sqrt(G/L);
fn=omegan/(2.*pi);
%
%   State Space Solution
%
period=1/fn;
%
out4 = sprintf('\n fn = %f Hz   period = %f sec  ',fn,period);
disp(out4)
%
dt=period/64;
sr=1/dt;
dur=8*period;
N=dur/dt;
%
out6 = sprintf('\n dt = %f sec    dur = %f sec ',dt,dur);
disp(out6)
%
t=linspace(0,dur,N);
t=t';
%
w(1)=d0;
w(2)=v0;
%
clear accel;
clear b;
clear c;
%
omegan2=omegan^2;
[t,w] = ode45(@pendulum_free_function,t,w,[],[],omegan2);
%
for(i=1:N)
    accel(i)=-omegan2*sin(w(i,1));
end
%
if(aunits==1)
    w(:,1)=w(:,1)*180/pi;  
    w(:,2)=w(:,2)*180/pi;   
    accel=accel*180/pi; 
end
%
dmax=max(w(:,1));
dmin=min(w(:,1));
vmax=max(w(:,2));
vmin=min(w(:,2));
%
if(aunits==1)
    out6 = sprintf('\n max disp = %11.4g deg   min disp = %11.4g deg',dmax,dmin);
    disp(out6)    
    out6 = sprintf('\n max vel = %11.4g deg/sec   min vel = %11.4g deg/sec',vmax,vmin);
    disp(out6)  
else
    out6 = sprintf('\n max disp = %11.4g rad   min disp = %11.4g rad',dmax,dmin);
    disp(out6)      
    out6 = sprintf('\n max vel = %11.4g rad/sec   min vel = %11.4g rad/sec',vmax,vmin);
    disp(out6)     
end    
%
disp(' ');
disp(' plot angular displacement?');
k=input(' 1=yes  2=no ');
%
if(k==1)
    figure(1);
    plot(t,w(:,1));
    xlabel('Time (sec)'); 
    if(aunits==1)
        ylabel('Angular Disp (deg)');
    else
        ylabel('Angular Disp (rad)');      
    end    
    grid on;    
end
%
disp(' ');
disp(' plot angular velocity?');
k=input(' 1=yes  2=no ');
%
if(k==1)
    figure(2);
    plot(t,w(:,2));
    xlabel('Time (sec)'); 
    if(aunits==1)
        ylabel('Angular Vel (deg/sec)');
    else
        ylabel('Angular Vel (rad/sec)');      
    end    
    grid on;   
end
%
disp(' ');
disp(' plot angular acceleration?');
k=input(' 1=yes  2=no ');
%
if(k==1)
    figure(3);
    plot(t,accel);
    xlabel('Time (sec)'); 
    if(aunits==1)
        ylabel('Angular Accel (deg/sec^2)');
    else
        ylabel('Angular Accel (rad/sec^2)');      
    end    
    grid on;   
end
disp(' ');
disp(' plot phase portrait?');
k=input(' 1=yes  2=no ');
if(k==1)
    figure(4);
    plot(w(:,1),w(:,2));
    if(aunits==1)
        ylabel('Angular Vel (deg/sec)');
        xlabel('Angular Disp (deg)');     
    else
        ylabel('Angular Vel (rad/sec)');
        xlabel('Angular Disp (rad)');         
    end
    grid on;   
end