Piston problem (main program source for Unix workstation)

%   This program is for 3.175 in "An Introduction to Dynamics, 3rd ed."
%   By David J. McGill and Wilton W. King.
%
%   part 1. constant torque (moment) on bar1
%

%
%   necessary input constant 
%

m = 1;                      %  mass of the piston
l1 = 1;                     %  length of the bar1 (left bar)
l2 = 2;                     %  length of the bar2 (right bar)
M = 5;                      %  constant moment at bar1
I = 3;                      %  inertia of bar1

tspan = [0,5];              %  time duration for calculation

%
%   Intitial conditions for theta and theta_dot
%

theta = 0;                  %  initial displacement of the bar
theta_dot = 0;              %  initial angular velocity of the bar1
theta_ddot = M/I;           %  initial angular acceleration of the bar1

%
%   Setup integration solver (ODE45)
%

x0 = [theta,theta_dot,theta_ddot];        

%
%  Perform ODE45 solver
%

options = odeset('Refine',10);  %  options 'Refine' was used to produce
                                %  good-looking plots
                                %  default value of Refine was 4

[t,x] = ode45('pist_sol',tspan,x0,options);
[n,m] = size(x);

Theta = x(:,1);
Phi = asin((l1/l2)*sin(Theta));
X = l1*cos(Theta)+l2*cos(Phi);

%
%  Plot the results
%

figure(1); clf; orient tall;
subplot(2,1,1),plot(t,Theta);
title('Problem 3.175'); 
xlabel('Time (sec.)');    
ylabel('Theta');

subplot(2,1,2),plot(t,Phi);
ylabel('Phi');
xlabel('Time (sec.)');      

figure(2); clf; orient tall;
subplot(2,1,1),plot(t,X);
title('Problem 3.175');
xlabel('Time (sec.)');
ylabel('X (m)');

subplot(2,1,2),plot(Theta,X);  
xlabel('Theta (rad.)');
ylabel('X (m)');

mxx = max(X);       mnx = min(X);
mxt = max(Theta);   mnt = min(Theta);
axis([mnt,mxt,mnx,mxx])