%
%   PENDULUM_1.M  Evaluate and plot the dynamics of a simple linear pendulum
%
%   This example evaluates and plots the angular position and velocity of a
%   pendulum with a point mass at the end.  The mathematical model was
%   linearized about the equilibrium position (pendulum at rest).  The model
%   was developed by performing an appropriate force balance and the
%   resulting linear constant coefficient 2nd order IVP was solved
%   analytically.  The exact expressions for the angular position and
%   velocity versus time are evaluated and plotted in this file (see the
%   lecture notes for the formal development and analytical expressions).
%
%   The primary purpose of this file is to serve as a typical example for
%   function evaluation and plotting in Matlab.  This is a common need for
%   many engineering design and analysis problems.  Thus, the student should
%   master this capability ASAP to become an effective "problem-solver" in
%   school and on the job.  In this particular example, as a bonus, you may
%   also gain a little insight into the actual dynamics of pendulum motion.
%   This is my intention anyway...
%
%   File prepared by J. R. White, UMass-Lowell (April 2003)
%

%
%   getting started
      clear all;   close all;   nfig = 0;
%
%   define problem parameters
      m = 1;                    % pendulum mass (kg)
      c = 2;                    % friction coeff (kg/s)
      L = 1;                    % length of pendulum arm (m)
      g = 10;                   % gravitational acceleration (m/s^2)
      poso = 30*pi/180;         % initial position (degrees converted to radians)
      velo = 0;                 % initial velocity (radians/s)
%
%   evaluate roots of characteristic equation
      a = -c/(2*m);   b = sqrt(g/L - a^2);
%
%   compute equation coefficients
      c1 = poso;          c2 = (velo - a*c1)/b;
      d1 = b*c2 + a*c1;   d2 = a*c2 - b*c1;

%   define discrete time domain variable (with time in seconds)
      Nt = 101;   to = 0;   tf = 5;   t = linspace(to,tf,Nt);
%
%   evaluate analytical expressions for angular postion and velocity
%   (be careful with the "dot" operations -- only need one dot in each
%   equation)
      pos = exp(a*t).*(c1*cos(b*t) + c2*sin(b*t));
      vel = exp(a*t).*(d1*cos(b*t) + d2*sin(b*t));
%
%   let's also compute the energy components versus time
      Ep = m*g*L*(1 - cos(pos));    % potential energy
      Ek = 0.5*m*L^2*vel.^2;        % kinetic energy
      Etot = Ep + Ek;               % total energy
%
%   now let's plot the results in a variety of ways
%   position (in degrees) vs time
      nfig = nfig+1;   figure(nfig)
      plot(t,pos*180/pi,'r-','LineWidth',2),grid    
      title('PENDULUM\_1: Pendulum Dynamics  --  \theta(t) vs t')
      xlabel('Time (sec)'),ylabel('Angular Position (degrees)')
%   position and velocity vs time  (use subplot command)
      nfig = nfig+1;   figure(nfig)
      subplot(2,1,1),plot(t,pos*180/pi,'r-','LineWidth',2),grid    
      title('PENDULUM\_1: Pendulum Dynamics  --  \theta(t) and \omega(t) vs t')
      ylabel('Angular Position (degrees)')
      subplot(2,1,2),plot(t,vel*180/pi,'g--','LineWidth',2),grid 
      xlabel('Time (sec)'),ylabel('Angular Velocity (deg/sec)')
%   just for fun, let's look at the phase plane diagram (velocity vs
%   position)
      nfig = nfig+1;   figure(nfig)
      plot(pos*180/pi,vel*180/pi,'r-','LineWidth',2),grid    
      title('PENDULUM\_1: Phase Plane Plot  --  \omega(t) vs \theta(t)')
      xlabel('Angular Position (degrees)'),ylabel('Angular Velocity (deg/sec)')
%   finally, let's plot the energy components (units -> J = N-m)
      nfig = nfig+1;   figure(nfig)
      plot(t,Ep,'r-',t,Ek,'g--',t,Etot,'b-.','LineWidth',2),grid 
      rr = axis;  rr(2) = 3;  axis(rr);   
      title('PENDULUM\_1: Energy Components  --  E_{tot} = E_p + E_k')
      xlabel('Time (sec)'),ylabel('Energy (J)')
      legend('Potential Energy -> E_p','Kinetic Energy -> E_k',...
      'Total Energy -> E_{tot}')
%
%   end of problem


