4.3  Lab 3, Simulation of crank-piston system

  4.3.1  Problem description
  4.3.2  Animation
  4.3.3  Part 1
  4.3.4  Part 2
  4.3.5  Part 3
  4.3.6  Appendix, source code

4.3.1  Problem description

pict

pict

4.3.2  Animation

The GUI based simulation was written in Matlab 2011a. A number of animated GIF files were made. In the table below, clicking on any image will start a gif animation in your browser.

These are relatively large animation GIF files and might take 1-2 seconds to load depending on network bandwidth. To run them locally, they can be downloaded as well.




pict

Using low spring stiffness to reach 5000 RPM in one second. (5 MB)

pict

Using low spring stiffness, disk is heavier than above, was not able to reach 5000 RPM in 1 second. (600 KB)

pict

Same as above, but the torque was removed at 500 RPM.




4.3.3  Part 1

The Lagrangian energy method was used to derive the equations of motion. The physical model used to represent the connection between the piston and the crank arm uses the method of Knarnopp-Margolis. This method uses a virtual spring and virtual damper for the connection between the piston and the crank arm instead of a rigid connection. This leads to a simpler mathematical model.

The equation of motion will now be derived based on this method.

Let \(f(\theta )\) be the distance which is shown in the diagram above as \(x^{\prime }\), The kinematic constraint is the following\[ f\left ( \theta \right ) =R\sin \theta +\sqrt{L^{2}-R^{2}\cos ^{2}\theta }\] The extension in the spring is \[ \Delta =x-f\left ( \theta \right ) \] Hence the spring force is\begin{equation} F_{s}=k\left [ x-f\theta )\right ] \tag{1} \end{equation} The damping force is\begin{align} F_{d} & =b\dot{x}\nonumber \\ & =b\frac{d}{dt}\left ( x-f\left ( \theta \right ) \right ) \nonumber \\ & =b\left ( \dot{x}-\frac{df\left ( \theta \right ) }{d\theta }\dot{\theta }\right ) \tag{2} \end{align}

The system is a 2 DOF, with generalized coordinates \(\theta ,x\).

Let \(I=\frac{MR^{2}}{2}\) be the moment of inertia of the disk (called \(J_{fw}\) in the above diagram) around an axis through the disk center, where \(M\) is the mass of the disk (called \(m_{fw}\) in the diagram). Let \(m\) be the mass of the cylinder (called \(m_{p}\) in the problem diagram).

The kinetic energy of the system is\[ T=\frac{1}{2}I\dot{\theta }^{2}+\frac{1}{2}m\dot{x}^{2}\] The potential energy is\begin{align*} V & =\frac{1}{2}k\Delta ^{2}\\ & =\frac{1}{2}k\left ( x-f\left ( \theta \right ) \right ) ^{2} \end{align*}

Therefore the Lagrangian is\begin{align*} L & =T-V\\ & =\frac{1}{2}I\dot{\theta }^{2}+\frac{1}{2}m\dot{x}^{2}-\frac{1}{2}k\left ( x-f\left ( \theta \right ) \right ) ^{2} \end{align*}

For \(\theta \) we find\begin{align*} \frac{\partial L}{\partial \dot{\theta }} & =I\dot{\theta }\\ \frac{d}{dt}\frac{\partial L}{\partial \dot{\theta }} & =I\ddot{\theta } \end{align*}

and\begin{align*} \frac{dL}{d\theta } & =-k\left ( x-f\left ( \theta \right ) \right ) \left ( -\frac{df\left ( \theta \right ) }{d\theta }\right ) \\ & =k\left ( x-f\left ( \theta \right ) \right ) \frac{df\left ( \theta \right ) }{d\theta } \end{align*}

Therefore the EQM for \(\theta \) is\begin{equation} \frac{d}{dt}\frac{\partial L}{\partial \dot{\theta }}-\frac{dL}{d\theta }=Q_{\theta }\tag{3} \end{equation} To find the generalized force \(Q_{\theta }\), we make a virtual \(\delta \theta \) displacement while keeping \(x\) fixed and then find the virtual work done. This results in\begin{align*} \delta W & =\tau \delta \theta +F_{d}\left ( \delta f\left ( \theta \right ) \right ) \\ & =\tau \delta \theta +F_{d}\frac{df\left ( \theta \right ) }{d\theta }\delta \theta \\ & =\left ( \tau +F_{d}\frac{df\left ( \theta \right ) }{d\theta }\right ) \delta \theta \end{align*}

Hence from the above we see that \[ Q_{\theta }=\tau +F_{d}\frac{df\left ( \theta \right ) }{d\theta }\] Eq. (3) becomes\begin{equation} I\ddot{\theta }-k\left ( x-f\left ( \theta \right ) \right ) \frac{df\left ( \theta \right ) }{d\theta }=\tau +F_{d}\frac{df\left ( \theta \right ) }{d\theta }\tag{4} \end{equation} Now we evaluate \(\frac{df\left ( \theta \right ) }{d\theta }\)\begin{align} f\left ( \theta \right ) & =R\sin \theta +\sqrt{L^{2}-R^{2}\cos ^{2}\theta }\nonumber \\ \frac{df\left ( \theta \right ) }{d\theta } & =R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\tag{5} \end{align}

Substituting the above back into Eq. (2) gives\begin{equation} F_{d}=b\left ( \dot{x}-\left ( R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \dot{\theta }\right ) \tag{6} \end{equation} Substituting Eqs. (6,5) into Eq. (4) gives \begin{multline*} I\ddot{\theta }-k\left ( x-f\left ( \theta \right ) \right ) \left ( R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) =\\ \tau +\left ( b\left ( \dot{x}-\left ( R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \dot{\theta }\right ) \right ) \left ( R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \end{multline*} Simplifying the above gives the final EQM for \(\theta \)\begin{align} \ddot{\theta } & =\frac{\tau }{I}+\frac{kR}{I}\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \left ( \sqrt{L^{2}-R^{2}\cos ^{2}\theta }-R\sin \theta +x\right ) \nonumber \\ & +\frac{bR}{I}\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \left ( \dot{x}-R\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \dot{\theta }\right ) \tag{7} \end{align}

Now to find the EQM for \(x\)\begin{align*} \frac{\partial L}{\partial \dot{x}} & =m\dot{x}\\ \frac{d}{dt}\frac{\partial L}{\partial \dot{x}} & =m\ddot{x} \end{align*}

and\[ \frac{dL}{dx}=-k\left ( x-f\left ( \theta \right ) \right ) \] Therefore the EQM for \(x\) is\begin{equation} \frac{d}{dt}\frac{\partial L}{\partial \dot{x}}-\frac{dL}{dx}=Q_{x}\tag{8} \end{equation} To find \(Q_{x}\) we make a virtual \(\delta x\) displacement while keeping \(\theta \) fixed and find find the virtual work. Hence\[ \delta W=-F_{d}\delta x \] Therefore \(Q_{x}=F_{d}\) which was found in Eq. (2) above. Hence \[ Q_{x}=-b\left ( \dot{x}-\frac{df\left ( \theta \right ) }{d\theta }\dot{\theta }\right ) \] And Eq. (8) becomes\[ m\ddot{x}+k\left ( x-f\left ( \theta \right ) \right ) =-b\left ( \dot{x}-\frac{df\left ( \theta \right ) }{d\theta }\dot{\theta }\right ) \] Substituting the expressions for \(f\left ( \theta \right ) \) and \(\frac{df\left ( \theta \right ) }{d\theta }\) found earlier into the above results in

\[ m\ddot{x}+k\left ( x-\left ( R\sin \theta +\sqrt{L^{2}-R^{2}\cos ^{2}\theta }\right ) \right ) =-b\left ( \dot{x}-\left ( R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \dot{\theta }\right ) \] Therefore\begin{equation} \ddot{x}=\frac{k\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}{m}+\frac{kR\sin \theta }{m}-\frac{kx}{m}-\frac{b\dot{x}}{m}+\frac{bR\dot{\theta }\cos \theta }{m}+\frac{b\dot{\theta }R^{2}\cos \theta \sin \theta }{m\sqrt{L^{2}-R^{2}\cos ^{2}\theta }} \tag{9} \end{equation}

4.3.4  Part 2

Now Eqs. (7) and (9) above will be cast into first order form.\[\begin{pmatrix} x_{1}\\ x_{2}\\ x_{3}\\ x_{4}\end{pmatrix} =\begin{pmatrix} \theta \\ x\\ \theta ^{\prime }\\ x^{\prime }\end{pmatrix} \overset{\frac{d}{dt}}{\rightarrow }\begin{pmatrix} \dot{x}_{1}\\ \dot{x}_{2}\\ \dot{x}_{3}\\ \dot{x}_{4}\end{pmatrix} =\begin{pmatrix} x_{3}\\ x_{4}\\ \theta ^{\prime \prime }\\ x^{\prime \prime }\end{pmatrix} \] Hence\begin{align*} \begin{pmatrix} \dot{x}_{1}\\ \dot{x}_{2}\\ \dot{x}_{3}\\ \dot{x}_{4}\end{pmatrix} & =\begin{pmatrix} x_{3}\\ x_{4}\\ \frac{\tau }{I}+\frac{kR}{I}\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \left ( \sqrt{L^{2}-R^{2}\cos ^{2}\theta }-R\sin \theta +x\right ) \\ +\frac{bR}{I}\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \left ( \dot{x}-R\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \dot{\theta }\right ) \\ \frac{k\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}{m}+\frac{kR\sin \theta }{m}-\frac{kx}{m}-\frac{b\dot{x}}{m}+\frac{bR\dot{\theta }\cos \theta }{m}+\frac{b\dot{\theta }R^{2}\cos \theta \sin \theta }{m\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\end{pmatrix} \\ & =\begin{pmatrix} x_{3}\\ x_{4}\\ \frac{\tau }{I}+\frac{kR}{I}\cos x_{1}\left ( 1+\frac{R\sin x_{1}}{\sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}}\right ) \left ( \sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}-R\sin x_{1}+x_{2}\right ) \\ +\frac{bR}{I}\cos x_{1}\left ( 1+\frac{R\sin x_{1}}{\sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}}\right ) \left ( x_{4}-R\cos x_{1}\left ( 1+\frac{R\sin x_{1}}{\sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}}\right ) x_{3}\right ) \\ \frac{k\sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}}{m}+\frac{kR\sin x_{1}}{m}-\frac{kx_{2}}{m}-\frac{bx_{4}}{m}+\frac{bRx_{3}\cos x_{1}}{m}+\frac{bx_{3}R^{2}\cos x_{1}\sin x_{1}}{m\sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}}\end{pmatrix} \end{align*}

4.3.5  Part 3

The system was simulated on the computer. The first order equations above were solved using ode numerical solver. A GUI program was written to make it easier to do the simulation and observe the results. The following describes the simulation program user interface

pict

Crank to 5000 RPM

We are asked to find the torque needed to reach \(5000\) RPM in one second. By trial and error, a torque of \(\tau =7\) \(NM\) was found to meet this condition, with \(f=10Hz\), and damping ratio \(\xi =0.2\), and using the same other parameters shown in the problem table. Using this torque, the \(5000\) RPM was reached at about \(0.5\) seconds. The force on the piston started to become lower after the torque was removed, as well as the piston acceleration.

The following diagram shows the final plots at the end of the one second run.

pict

Simulation to keep spring deflection below fraction of mm

We are also asked to find the parameters which will cause the spring deflection to remain below fraction of a millimeter. It was found that the spring \(k\) must be made very large to simulate the effect of a rigid connection between the piston and the crank arm in order to keep the deflection very small. But when this was done, the simulation started to take very long time. Simulating \(1\) second on my computer would more than one hour.

To reduce vibration, the damping coefficient \(\xi \) was made larger than the \(0.2\) value in the table. But even when making \(\xi \) close to one, the simulation still took very long time to reach one second due to the large number of steps taken by the numerical solver. The spring stiffness was increased all the way to \(10^{7}\) N/m, and the force on the piston reached \(1.5MPa\) at one point. The value of \(f\) (spring frequency) to achieve this high spring stiffness was \(1000\) \(Hz.\) The damping ratio \(\xi \) was set to \(0.9\). The following is a plot of the simulation using the above parameters, showing that the deflection remained below half millimeter.

pict

Force and acceleration on piston

The acceleration on the piston was found numerically at each step by dividing the difference of the current speed of piston from the last time step speed, and dividing that over the difference in time as follows\[ a_{x}=\frac{\dot{x}\left ( t_{n}\right ) -\dot{x}\left ( t_{n-1}\right ) }{t_{n}-t_{n-1}}\] The force on the piston was found at each step as follows\[ F=-(F_{s}+F_{d}) \] Where \(F_{s}\) is the force in spring and \(F_{d}\) is the damping force. Expression for these forces were derived above.

The force on the piston depended on the torque used and the amount of spring stiffness and damping used. For example, in the first test case above when the spring was made very stiff, the largest force on the piston reached was \(1.5\times 10^{6}N\) and an acceleration of almost \(1000g\) was seen. After the torque was removed, the acceleration gradually dropped to \(500g\). The longer the simulation was run, these values became smaller. The simulation program updates both the acceleration and the force on the piston as the simulation is running.

General observation

The Knarnopp-Margolis model of the crank/piston connection was used to simulate the motion. This model is simpler mathematically than the rigid connection model. But it was hard to find the right parameters to approximate a rigid connection since making the spring stiffness large, causes the simulation to become very slow. When the spring stiffness was made very large to better approximate a rigid connection, and when the damping was made larger, the simulation would take about \(1\) hour for \(1\) second simulation time. Even when using a stiff ode solver such as ode15s the simulation was taking too long. Using a compiled language would eliminate this problem as it will be much faster, but this simulation was done using Matlab

The main advantage of the Knarnopp-Margolis model is that the implementation was simpler since the mathematical model was simpler. But more trial and error are needed to find optimal values for all the parameters of the problem which will cause the simulation time to become more reasonable, while approximating the rigid connection.

4.3.6  Appendix, source code

 
function varargout = nma_lab3_eme_121_fig(varargin) 
% NMA_LAB3_EME_121_FIG MATLAB code for nma_lab3_eme_121_fig.fig 
%      NMA_LAB3_EME_121_FIG, by itself, creates a new NMA_LAB3_EME_121_FIG or raises the existing 
%      singleton*. 
% 
%      H = NMA_LAB3_EME_121_FIG returns the handle to a new NMA_LAB3_EME_121_FIG or the handle to 
%      the existing singleton*. 
% 
%      NMA_LAB3_EME_121_FIG('CALLBACK',hObject,eventData,handles,...) calls the local 
%      function named CALLBACK in NMA_LAB3_EME_121_FIG.M with the given input arguments. 
% 
%      NMA_LAB3_EME_121_FIG('Property','Value',...) creates a new NMA_LAB3_EME_121_FIG or raises the 
%      existing singleton*.  Starting from the left, property value pairs are 
%      applied to the GUI before nma_lab3_eme_121_fig_OpeningFcn gets called.  An 
%      unrecognized property name or invalid value makes property application 
%      stop.  All inputs are passed to nma_lab3_eme_121_fig_OpeningFcn via varargin. 
% 
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one 
%      instance to run (singleton)". 
% 
% this is the GUI main line for EME 121 lab 3 
% by Nasser M. Abbasi 
% See also: GUIDE, GUIDATA, GUIHANDLES 
 
% Edit the above text to modify the response to help nma_lab3_eme_121_fig 
 
% Last Modified by GUIDE v2.5 07-May-2011 17:21:06 
 
% Begin initialization code - DO NOT EDIT 
gui_Singleton = 1; 
gui_State = struct('gui_Name',       mfilename, ... 
                   'gui_Singleton',  gui_Singleton, ... 
                   'gui_OpeningFcn', @nma_lab3_eme_121_fig_OpeningFcn, ... 
                   'gui_OutputFcn',  @nma_lab3_eme_121_fig_OutputFcn, ... 
                   'gui_LayoutFcn',  [] , ... 
                   'gui_Callback',   []); 
if nargin && ischar(varargin{1}) 
    gui_State.gui_Callback = str2func(varargin{1}); 
end 
 
if nargout 
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); 
else 
    gui_mainfcn(gui_State, varargin{:}); 
end 
% End initialization code - DO NOT EDIT 
 
 
% --- Executes just before nma_lab3_eme_121_fig is made visible. 
function nma_lab3_eme_121_fig_OpeningFcn(hObject, eventdata, handles, varargin) 
% This function has no output args, see OutputFcn. 
% hObject    handle to figure 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
% varargin   command line arguments to nma_lab3_eme_121_fig (see VARARGIN) 
 
% Choose default command line output for nma_lab3_eme_121_fig 
handles.output = hObject; 
 
set(handles.figure1, 'UserData',[]); 
set(handles.figure1,'Name','UC Davis, EME 121 lab#3, by Nasser M. Abbasi'); 
 
userData.stop = false; 
set(handles.figure1, 'UserData',userData); 
 
 
% Update handles structure 
guidata(hObject, handles); 
 
% UIWAIT makes nma_lab3_eme_121_fig wait for user response (see UIRESUME) 
% uiwait(handles.figure1); 
 
 
% --- Outputs from this function are returned to the command line. 
function varargout = nma_lab3_eme_121_fig_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT); 
% hObject    handle to figure 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Get default command line output from handles structure 
varargout{1} = handles.output; 
 
 
 
function mfwTag_Callback(hObject, eventdata, handles) 
% hObject    handle to mfwTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of mfwTag as text 
%        str2double(get(hObject,'String')) returns contents of mfwTag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function mfwTag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to mfwTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function fTag_Callback(hObject, eventdata, handles) 
% hObject    handle to fTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of fTag as text 
%        str2double(get(hObject,'String')) returns contents of fTag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function fTag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to fTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function zetaTag_Callback(hObject, eventdata, handles) 
% hObject    handle to zetaTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of zetaTag as text 
%        str2double(get(hObject,'String')) returns contents of zetaTag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function zetaTag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to zetaTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function torqueTag_Callback(hObject, eventdata, handles) 
% hObject    handle to torqueTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of torqueTag as text 
%        str2double(get(hObject,'String')) returns contents of torqueTag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function torqueTag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to torqueTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function rpmTag_Callback(hObject, eventdata, handles) 
% hObject    handle to rpmTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of rpmTag as text 
%        str2double(get(hObject,'String')) returns contents of rpmTag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function rpmTag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to rpmTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
% --- Executes on button press in runBtn. 
function runBtn_Callback(hObject, eventdata, handles) 
% hObject    handle to runBtn (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
userData = get(handles.figure1, 'UserData'); 
userData.stop = false; 
set(handles.figure1, 'UserData',userData); 
 
data.angleInitial=(pi/2)/180; 
data.angleSpeedInitial=0; 
 
data.M=str2num(get(handles.mfwTag,'String')); 
data.volumeOfCylinder=500*0.000001 ; %convery cc to m^3 
data.Dp=3*.0254; 
data.Ap=(pi/4)*data.Dp^2; 
data.stroke=data.volumeOfCylinder/data.Ap; 
data.R=data.stroke/2; 
data.L=3*data.R; 
data.I=data.M*data.R^2/2; 
set(handles.Itag,'String',sprintf('%3.3f',data.I)); 
 
data.m=0.25*data.I/data.R^2; 
set(handles.pistonMassTag,'String',sprintf('%3.3f',data.m)); 
 
data.f=str2num(get(handles.fTag,'String')); 
set(handles.pistonFreqTag,'String',sprintf('%3.3f',data.f)); 
 
data.omega=2*pi*data.f; 
set(handles.pistonWnTag,'String',sprintf('%3.3f',data.omega)); 
 
 
data.dampingRatio=str2num(get(handles.zetaTag,'String')); 
 
 
data.k=data.m*data.omega^2; 
set(handles.kTag,'String',sprintf('%3.3f',data.k)); 
 
data.b=2*data.dampingRatio*data.omega*data.m; 
set(handles.bTag,'String',sprintf('%3.3f',data.b)); 
 
data.torque = str2num(get(handles.torqueTag,'String')); 
data.maxRPM=str2num(get(handles.rpmTag,'String')); 
data.maxSimulationTime=str2num(get(handles.simTimeTag,'String')); 
data.handles=handles; 
data.currentAngle=data.angleInitial; 
data.minY=0; 
data.maxY=4; 
 
 
[g_msg,g_status]=nma_lab3_eme_121(data); 
 
 
% --- Executes on button press in stopTag. 
function stopTag_Callback(hObject, eventdata, handles) 
% hObject    handle to stopTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
userData = get(handles.figure1, 'UserData'); 
userData.stop = true; 
set(handles.figure1, 'UserData',userData); 
%enable_all(handles,'on'); 
 
 
 
function simTimeTag_Callback(hObject, eventdata, handles) 
% hObject    handle to simTimeTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of simTimeTag as text 
%        str2double(get(hObject,'String')) returns contents of simTimeTag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function simTimeTag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to simTimeTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function edit7_Callback(hObject, eventdata, handles) 
% hObject    handle to edit7 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of edit7 as text 
%        str2double(get(hObject,'String')) returns contents of edit7 as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function edit7_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to edit7 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function edit8_Callback(hObject, eventdata, handles) 
% hObject    handle to edit8 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of edit8 as text 
%        str2double(get(hObject,'String')) returns contents of edit8 as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function edit8_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to edit8 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function edit9_Callback(hObject, eventdata, handles) 
% hObject    handle to edit9 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of edit9 as text 
%        str2double(get(hObject,'String')) returns contents of edit9 as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function edit9_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to edit9 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function edit10_Callback(hObject, eventdata, handles) 
% hObject    handle to edit10 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of edit10 as text 
%        str2double(get(hObject,'String')) returns contents of edit10 as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function edit10_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to edit10 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function edit11_Callback(hObject, eventdata, handles) 
% hObject    handle to edit11 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of edit11 as text 
%        str2double(get(hObject,'String')) returns contents of edit11 as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function edit11_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to edit11 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function edit12_Callback(hObject, eventdata, handles) 
% hObject    handle to edit12 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of edit12 as text 
%        str2double(get(hObject,'String')) returns contents of edit12 as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function edit12_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to edit12 (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function Itag_Callback(hObject, eventdata, handles) 
% hObject    handle to Itag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of Itag as text 
%        str2double(get(hObject,'String')) returns contents of Itag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function Itag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to Itag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function pistonMassTag_Callback(hObject, eventdata, handles) 
% hObject    handle to pistonMassTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of pistonMassTag as text 
%        str2double(get(hObject,'String')) returns contents of pistonMassTag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function pistonMassTag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to pistonMassTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function pistonFreqTag_Callback(hObject, eventdata, handles) 
% hObject    handle to pistonFreqTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of pistonFreqTag as text 
%        str2double(get(hObject,'String')) returns contents of pistonFreqTag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function pistonFreqTag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to pistonFreqTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function pistonWnTag_Callback(hObject, eventdata, handles) 
% hObject    handle to pistonWnTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of pistonWnTag as text 
%        str2double(get(hObject,'String')) returns contents of pistonWnTag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function pistonWnTag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to pistonWnTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function kTag_Callback(hObject, eventdata, handles) 
% hObject    handle to kTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of kTag as text 
%        str2double(get(hObject,'String')) returns contents of kTag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function kTag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to kTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end 
 
 
 
function bTag_Callback(hObject, eventdata, handles) 
% hObject    handle to bTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    structure with handles and user data (see GUIDATA) 
 
% Hints: get(hObject,'String') returns contents of bTag as text 
%        str2double(get(hObject,'String')) returns contents of bTag as a double 
 
 
% --- Executes during object creation, after setting all properties. 
function bTag_CreateFcn(hObject, eventdata, handles) 
% hObject    handle to bTag (see GCBO) 
% eventdata  reserved - to be defined in a future version of MATLAB 
% handles    empty - handles not created until after all CreateFcns called 
 
% Hint: edit controls usually have a white background on Windows. 
%       See ISPC and COMPUTER. 
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) 
    set(hObject,'BackgroundColor','white'); 
end
 
function [g_msg,g_status]=nma_lab3_eme_121(data) 
% 
% process function called from GUI driver code to 
% plot simulation result 
% 
% EME 121, spring 2011 
% 
% by Nasser M. Abbasi, Matlab 2011 
% April 2011 
% 
 
[xInitial,~]=f(data.L,data.R,data.angleInitial); 
data.xInitial=xInitial; 
data.xSpeedInitial=0; 
 
data.timeStep=0; 
data.solver='ode15s'; 
 
g_status = true; 
g_msg=''; 
 
set(0,'CurrentFigure',data.handles.figure1); 
cla(data.handles.xAxesTag,'reset'); 
cla(data.handles.angleAxesTag,'reset'); 
cla(data.handles.rpmAxesTag,'reset'); 
cla(data.handles.mainAxes,'reset'); 
 
cla(data.handles.forceAxes,'reset'); 
cla(data.handles.accAxes,'reset'); 
 
IC = [data.angleInitial;data.xInitial;... 
    data.angleSpeedInitial;data.xSpeedInitial]; 
options = odeset('OutputFcn',@output); 
 
%Check if using auto step size of an explicit step size. see 
%GUI for input choice. 
if data.timeStep==0 
    timeSpan=[0 data.maxSimulationTime]; 
else 
    timeSpan=0:data.timeStep:data.maxSimulationTime; 
end 
 
%call  the solver 
 
SCREEN_CAPTURE=false;  %set to zero for animation capture 
 
currentIter=0; 
 
if SCREEN_CAPTURE 
    frameNumber = 0; 
    actualFrameNumber = 0; 
    theFrame=getframe(gcf); 
    [im,MAP]=rgb2ind(theFrame.cdata,32,'nodither'); 
    %MAP=colormap(gray); 
end 
[t,x] = feval(data.solver,@rhs,timeSpan,IC,options); 
 
if SCREEN_CAPTURE 
    imwrite(im,MAP,'demo.gif','DelayTime',0,'LoopCount',inf,... 
        'WriteMode','overwrite'); 
end 
 
%----------------------------------- 
% ODE solver function 
%----------------------------------- 
    function dx=rhs(t,x) 
        %the ode45 function 
 
        k=data.k; 
        R=data.R; 
        L=data.L; 
        I=data.I; 
        b=data.b; 
        m=data.m; 
        tao=data.torque; 
 
        [d,ext]=f(L,R,x(1)); 
        y=sqrt(L^2-R^2*(cos(x(1)))^2); 
        if x(3)/(2*pi)*60>data.maxRPM 
            data.torque=0; 
        end 
 
        %fprintf('angle=%3.3f\t x(2)=%3.3f\t f(theta)=%3.3f\n',mod(x(1)*180/pi,360),x(2),d); 
 
        dx=[x(3); 
            x(4); 
 
            tao/I + (k*R/I)*cos(x(1)) * ( 1+R*sin(x(1))/y )* ... 
            (-y-R*sin(x(1)) + x(2) ) + ... 
            b*R/I * cos(x(1)) * ( 1 + (R*sin(x(1))/y) ) *... 
            ( x(4) - R*cos(x(1)) * (1 + (R*sin(x(1))/y) ) * x(3) ); 
 
            k*y/m + k*R*sin(x(1))/m - k*x(2)/m - b*x(4)/m ... 
            + b*R*x(3)*cos(x(1))/m + b*x(3)*R^2*cos(x(1))*sin(x(1))/(m*y) 
            ]; 
 
    end 
%----------------------------------- 
%  ODE output function 
%----------------------------------- 
    function status= output(t,x,~) 
        %called by ode45 after each step. Plot the current 
        %pendulum position for simulation 
 
        userData = get(data.handles.figure1, 'UserData'); 
        if userData.stop 
            status=true; 
            g_status =true; 
            return; 
        end 
 
        if isempty(t) 
            status = true; 
        else 
 
            currentIter=currentIter+1; 
            if currentIter==1 
                status = false; 
                return; 
            end 
 
            [D,ext]=f(data.L,data.R,x(1,1)) ; 
 
            data.currentAngle=x(1); 
            % plot angle vs time 
 
            if currentIter>2 
                % position vs. time 
                set(0,'CurrentFigure',data.handles.figure1); 
                set(data.handles.figure1, 'CurrentAxes',... 
                    data.handles.xAxesTag); 
                plot([data.lastTime t], [data.lastPosition x(2)],'-'); 
                xlim([0 data.maxSimulationTime]); 
                %                 if abs(data.minY-data.maxY)>10*eps 
                %                     ylim([data.minY data.maxY]); 
                %                 end 
                xlabel('time (sec)'); 
                title({'Position Vs. time';sprintf('x = %3.7f meter',x(2))}); 
                hold on; 
                set(gca,'FontSize',7); 
                drawnow; 
 
                %plot vibration 
                set(0,'CurrentFigure',data.handles.figure1); 
                set(data.handles.figure1, 'CurrentAxes',... 
                    data.handles.angleAxesTag); 
                currentSpringExtensionIn_mm = 1000*(x(2)-D); 
 
                plot([data.lastTime t],... 
                [data.lastSpringExtension currentSpringExtensionIn_mm],'-'); 
 
                %ylabel('mm'); 
                xlabel('time (sec)'); 
                title({'virtual spring extension in mm';... 
                    sprintf('%3.7f',currentSpringExtensionIn_mm)}); 
                hold on; 
                set(gca,'FontSize',7); 
                %xlim([0 data.maxSimulationTime]); 
                drawnow; 
 
                % RPMangular speed vs. time plot 
                set(data.handles.figure1, 'CurrentAxes',... 
                    data.handles.rpmAxesTag); 
 
                plot([data.lastTime t],[data.lastRPM x(3)/(2*pi)*60],'-'); 
                title({'current RPM';sprintf('%3.3f ',x(3)/(2*pi)*60)}); 
                xlim([0 data.maxSimulationTime]); 
                xlabel('time (sec)'); 
                hold on; 
                set(gca,'FontSize',7); 
                drawnow; 
 
                %force on piston plot 
                set(0,'CurrentFigure',data.handles.figure1); 
                set(data.handles.figure1,'CurrentAxes',data.handles.forceAxes); 
                currentFd = currentDampingForce(x,data.L,data.R,data.b); 
                currentFs =  data.k*(x(2)-D); 
 
                % if sign(currentFd) ~= sign(currentFs) 
                %     fprintf('spring force=%3.6f, damping force=%3.6f\n',... 
                %     currentFs,currentFd); 
                % end 
 
                forceOnPiston = - (currentFs +currentFd); 
 
                plot([data.lastTime t], [data.lastTotalForceOnPiston... 
                    forceOnPiston],'-'); 
 
                data.lastTotalForceOnPiston=forceOnPiston; 
 
                set(gca,'FontSize',7); 
                title({'force on piston (Newton)';... 
                    sprintf('%3.7f',forceOnPiston)}); 
                hold on; 
                xlim([0 data.maxSimulationTime]); 
                drawnow; 
 
                % acc in g 
                set(data.handles.figure1, 'CurrentAxes',... 
                    data.handles.accAxes); 
 
                currentAcc=(x(4)-data.lastXspeed)/(t-data.lastTime); 
                currentAcc=currentAcc/9.8; 
                plot([data.lastTime t],[data.lastAcc currentAcc],'-'); 
                title({'piston acceleration (in g)';... 
                    sprintf('%3.7f',currentAcc)}); 
 
                xlim([0 data.maxSimulationTime]); 
                xlabel('time (sec)'); 
 
                if currentIter>3 
                    hold on; 
                end 
                set(gca,'FontSize',7); 
                drawnow; 
 
                data.lastTime=t; 
                data.lastPosition=x(2); 
                data.lastRPM=x(3)/(2*pi)*60; 
                data.lastAngle=mod(x(1)*180/pi,360); 
                data.lastSringExtension=1000*(x(2)-D); 
                data.lastXspeed=x(4); 
                data.lastAcc=currentAcc; 
            else 
                data.lastAngle=x(1); 
                data.lastPosition=x(2); 
                data.lastTime=t; 
                data.lastRPM=x(3)/(2*pi)*60; 
                data.lastTotalForceOnPiston = - ( data.k*(x(2)-D)+... 
                    currentDampingForce(x,data.L,data.R,data.b) ); 
                data.lastSpringExtension=1000*(x(2)-D); 
                data.lastXspeed=x(4); 
                data.lastAcc=0; 
            end 
 
            %---------------- central simulation now added 
            set(data.handles.figure1, 'CurrentAxes',data.handles.mainAxes); 
 
            % 
            %             %plot circle 
            z=-1:.01:1; 
            plot(data.R*sin(2*pi*z),data.R*cos(2*pi*z),'LineWidth',3); 
 
            axis equal; 
            hold on; 
            theta = pi/2-x(1); 
            % 
            % 
            %plot R line 
 
            plot([0 data.R*cos( theta) D],[0 data.R*sin( theta) 0],'k-'); 
            axis equal; 
            plot([data.R*cos(theta) D],[data.R*sin(theta) 0],'r-',... 
                'LineWidth',2); 
            plot([0 D],[0 0],'k-');   axis equal; 
 
            %plot piston housing 
            plot([0.8*data.L 1.37*data.L],... 
                [data.R/3 data.R/3 ],'-','LineWidth',4,'color','k'); 
 
            %plot piston housing 
            plot([1.4*data.L 2*data.L 2*data.L 1.4*data.L],... 
                [data.R/3 data.R/3 -data.R/3 -data.R/3],'-',... 
                'LineWidth',4,'color','k'); 
 
            %plot piston housing 
            plot([0.8*data.L 1.37*data.L],... 
                [-data.R/3 -data.R/3 ],'-','LineWidth',4,'color','k'); 
 
 
            plot([D x(2,1)],[0 0],'r');   axis equal; 
            %lengthOfSpring = abs(x(2,1)-D); 
            %if (x(2,1)-D)<0 
            %    fprintf('WARNING, solution x is smaller than f(theta)\n'); 
            %end 
 
            plot(data.R*cos( theta),data.R*sin( theta),'o','MarkerSize',6,... 
                'MarkerFaceColor','r');   axis equal; 
            plot(0,0,'o','MarkerSize',4,... 
                'MarkerFaceColor','k');   axis equal; 
 
            %spring removed for now. 
            %plot spring 
            %[xSpring,ySpring]=nma_spring.makeV2(lengthOfSpring,5,data.R/4); 
            % 
            %A = [cos(-pi/2) -sin(-pi/2);sin(-pi/2) cos(-pi/2)]; 
            %newSpringCoordinates=arrayfun(@(i) A*[xSpring(i);ySpring(i)],... 
            %        1:length(xSpring),'UniformOutput',false); 
            %        newSpringCoordinates= [newSpringCoordinates{:}]'; 
            % 
            %        xSpring=newSpringCoordinates(:,1); 
            %        ySpring=newSpringCoordinates(:,2); 
            %        xSpring=xSpring+(x(2,1)-lengthOfSpring); 
            %        line(xSpring,ySpring); 
 
            %make box 
            %plot the box 
            plot(x(2),0,'s','MarkerSize',28,... 
                'MarkerFaceColor','m'); 
 
            set(gca,'XTickLabel',[]); 
            set(gca,'YTickLabel',[]); 
            if data.torque>0 
                textToShow='torque is on'; 
            else 
                textToShow='torque is off'; 
            end 
 
            title({'Crank-Piston simulation using Karnopp-Margolis virtual piston mathematical model';... 
                'spring/damper attached to piston but not displayed'; 
                sprintf('Current time = [%4.3f sec], %s',t(1),textToShow)},... 
                'FontSize',8); 
            axis equal; 
 
            xlim([-1.3*data.R data.R+2*data.L]); 
 
            %draw torque if needed 
            %plot torque as arc around center. 
            if abs(data.torque)>10*eps 
                n = 100; % The number of points in the arc 
                x0=0;y0=0; x1=data.R*.6; y1=0; x2=-data.R*.6; y2=.1; 
 
                P0 = [x0;y0]; 
                P1 = [x1;y1]; 
                P2 = [x2;y2]; 
                v1 = P1-P0; 
                v2 = P2-P0; 
                c  = det([v1,v2]); 
                a  = linspace(0,atan2(abs(c),dot(v1,v2)),n); % Angle range 
                v3 = [0,-c;c,0]*v1; 
                v  = v1*cos(a)+((norm(v1)/norm(v3))*v3)*sin(a); 
                A = [cos(theta) -sin(theta);sin(theta) cos(theta)]; 
                xy_1=A*[v(1,:)+x0;v(2,:)+y0]; 
                line(xy_1(1,6:end),xy_1(2,6:end),... 
                    'LineWidth',1.5,'LineStyle','-',... 
                    'color','black'); 
 
                line(xy_1(1,1:5),xy_1(2,1:5),'LineWidth',2,... 
                    'LineStyle','none',... 
                    'Marker','>','MarkerFaceColor','red'); 
                %text('Interpreter','latex','String','$$\tau$$',... 
                %    'Position',[xy_1(1,end),xy_1(2,end)],... 
                %    'Fontsize',16); 
 
                hold on; 
                axis equal 
            end 
 
            axis equal; 
 
            drawnow; 
            hold off; 
            status = false; 
            if SCREEN_CAPTURE 
                frameNumber = frameNumber+1; 
                if mod(frameNumber,5)==0 
                    theFrame = getframe(gcf); 
                    actualFrameNumber = actualFrameNumber +1; 
                    im(:,:,1,actualFrameNumber ) = ... 
                        rgb2ind(theFrame.cdata,MAP,'nodither'); 
                    %im=rgb2ind(theFrame.cdata,MAP,'nodither'); 
                    %imwrite(im,MAP,[num2str(actualFrameNumber) '.gif'],'gif'); 
                end 
            end 
            %pause(.001); 
        end 
    end 
end 
%-------------------------------- 
% This is the f(theta) function in the textbook 
%------------------------------- 
function [d,ext]=f(L,R,theta) 
 
ext=R*sin(theta); 
d=sqrt(L^2-R^2*cos(theta)^2); 
d=d+ext; 
end 
%------------------ 
function fd=currentDampingForce(x,L,R,b) 
%see equation 6 in report 
a=x(1); 
 
fd=b*(x(4)-(R*cos(a)+ R^2*cos(a)*sin(a)/sqrt(L^2-R^2*cos(a)^2))*x(3)); 
 
end
 
classdef nma_spring 
    % 
    %static class to make spring for plotting animations 
    %by Nasser M. Abbasi 
 
 
    properties 
    end 
 
    methods(Static) 
 
        %--------------------------- 
        % 
        %--------------------------- 
        function [x,y] = make(r,theta,N) 
            %r: total length of spring 
            %theta: in radians, anticlock wise is positive, 
            %       theta zero is position x-axis 
            %N: number of twists in spring 
            % 
            %OUTPUT: 
            % x,y coordinates of line to use to plot spring 
            len = (4/6)*r; 
            p = zeros(N,2); 
            delr = len/N; 
 
            r0     = (1/6)*r; 
            p(2,1) = r0; 
            p(2,2) = theta; 
 
            for n=3:N-2 
                p(n,1)=r0+delr*n; 
                z=atan(2*delr/p(n,1)); 
                p(n,2)=theta+(-1)^n*z; 
            end 
            p(end-1,1)=(5/6)*r; p(end-1,2)=theta; 
            p(end,1)=r; p(end,2)=theta; 
 
            [x,y] = pol2cart(p(:,2),p(:,1)); 
        end 
 
        %---------------------------------------- 
        function [x,y] = makeBox(w,h) 
            x = zeros(5,1); 
            y = zeros(5,1); 
            x(1)= w/2; y(1)=0; 
            x(2)= w/2; y(2)=h/2; 
            x(3)= -w/2; y(3)=h/2; 
            x(4)=-w/2;  y(4)=0; 
            x(5)=x(1); y(5)=y(1); 
 
        end 
        %--------------------------- 
        % 
        %--------------------------- 
        function [x,y] = makeV2(len,N,w) 
            %len: total length of spring 
            %N: number of twists in spring 
            % 
            %OUTPUT: 
            % x,y coordinates of line to use to plot spring 
 
            nCoordinates = 2*N+4; 
            p = zeros(nCoordinates,2); 
            h = (1/6)*len; 
            del=(len-2*h)/(N+1); 
            done = false; 
            n = 1; 
            while not(done) 
 
                if n==1 
                   p(1,1)=0; 
                   p(1,2)=0; 
                   n = n + 1; 
                elseif n==2 
                   p(2,1)=0; 
                   p(2,2)=h; 
                   n = n + 1; 
                elseif n<nCoordinates-2 
                   p(n,1)=w; 
                   p(n,2)=p(n-1,2)+del; 
                   n = n  + 1; 
                   p(n,1)=-w; 
                   p(n,2)=p(n-1,2); 
                   n = n + 1; 
                elseif n==nCoordinates-1 
                   p(n,1) = 0; 
                   p(n,2)=p(n-1,2)+del; 
                   n = n + 1; 
                else 
                   p(n,1) = 0; 
                   p(n,2)=p(n-1,2)+h; 
                   done = true; 
                end 
            end 
            x=p(:,1); 
            y=p(:,2); 
 
        end 
 
 
 
        %------------------------ 
        % 
        %------------------------ 
        function test(theta) 
 
            close all; 
 
            theta = theta*pi/180; 
            figure; 
            set(gcf,'Units','normalized'); 
            delr=0.3; 
            for r=1:.1:3 
                cla; 
                [x,y] = nma_spring.make(r,theta,30); 
                line(x,y); 
                hold on; 
 
                %baseX = [2*delr*cos((pi/2)-theta) -2*delr*sin(theta)]; 
                %baseY = [-2*delr*sin((pi/2)-theta) 2*delr*cos(theta)]; 
 
                baseX = [-2*delr*sin(abs(theta)) 2*delr*sin(abs(theta))]; 
                baseY = [-2*delr*cos(theta)      2*delr*cos(theta)]; 
                line(baseX,baseY); 
 
                triX=[baseX(1) ... 
                    baseX(1) ... 
                    30*delr*cos(theta)-abs(baseX(1))... 
                    baseX(1)]; 
 
                triY=[baseY(1) ... 
                    -(30*delr*sin(abs(theta))+abs(baseY(1)))... 
                    -(30*delr*sin(abs(theta))+abs(baseY(1)))... 
                    baseY(1)]; 
 
                line(triX,triY); 
 
                L=sqrt(x(end)^2+y(end)^2); 
                xx0=L*cos(theta); 
                yy0=-L*sin(abs(theta)); 
                massX=[xx0 ... 
                    xx0+2*delr*sin(abs(theta)) ... 
                    xx0+2*delr*sin(abs(theta))+4*delr*cos(theta) ... 
                    xx0+2*delr*sin(abs(theta))+4*delr*cos(theta)-4*delr*sin(abs(theta))... 
                    xx0+2*delr*sin(abs(theta))+4*delr*cos(theta)-4*delr*sin(abs(theta))-4*delr*cos(theta)... 
                    xx0]; 
 
                massY=[ yy0 ... 
                    yy0+2*delr*cos(theta) ... 
                    yy0+2*delr*cos(theta)-4*delr*sin(abs(theta)) ... 
                    yy0+2*delr*cos(theta)-4*delr*sin(abs(theta))-4*delr*cos(theta)... 
                    yy0+2*delr*cos(theta)-4*delr*sin(abs(theta))-4*delr*cos(theta)+4*delr*sin(abs(theta))... 
                    yy0]; 
 
                line(massX,massY); 
 
                xlim([-3,10]); ylim([-10,3]); 
 
                set(gca,'DataAspectRatioMode','manual','DataAspectRatio',[1 1 1]); 
                drawnow; 
                pause(.1); 
            end 
 
        end 
 
        %------------------------------- 
        % 
        %------------------------------- 
        function testV2(len) 
 
            close all; 
 
            theta = 0; 
            figure; 
            set(gcf,'Units','normalized'); 
            delr=0.3; 
            for r=1:len/3:len 
                cla; 
                [x,y] = nma_spring.makeV2(r,10,1); 
                line(x,y); 
                hold on; 
 
                %baseX = [2*delr*cos((pi/2)-theta) -2*delr*sin(theta)]; 
                %baseY = [-2*delr*sin((pi/2)-theta) 2*delr*cos(theta)]; 
 
                baseX = [-2*delr*sin(abs(theta)) 2*delr*sin(abs(theta))]; 
                baseY = [-2*delr*cos(theta)      2*delr*cos(theta)]; 
                line(baseX,baseY); 
 
                triX=[baseX(1) ... 
                    baseX(1) ... 
                    30*delr*cos(theta)-abs(baseX(1))... 
                    baseX(1)]; 
 
                triY=[baseY(1) ... 
                    -(30*delr*sin(abs(theta))+abs(baseY(1)))... 
                    -(30*delr*sin(abs(theta))+abs(baseY(1)))... 
                    baseY(1)]; 
 
                line(triX,triY); 
 
                L=sqrt(x(end)^2+y(end)^2); 
                xx0=L*cos(theta); 
                yy0=-L*sin(abs(theta)); 
                massX=[xx0 ... 
                    xx0+2*delr*sin(abs(theta)) ... 
                    xx0+2*delr*sin(abs(theta))+4*delr*cos(theta) ... 
                    xx0+2*delr*sin(abs(theta))+4*delr*cos(theta)-4*delr*sin(abs(theta))... 
                    xx0+2*delr*sin(abs(theta))+4*delr*cos(theta)-4*delr*sin(abs(theta))-4*delr*cos(theta)... 
                    xx0]; 
 
                massY=[ yy0 ... 
                    yy0+2*delr*cos(theta) ... 
                    yy0+2*delr*cos(theta)-4*delr*sin(abs(theta)) ... 
                    yy0+2*delr*cos(theta)-4*delr*sin(abs(theta))-4*delr*cos(theta)... 
                    yy0+2*delr*cos(theta)-4*delr*sin(abs(theta))-4*delr*cos(theta)+4*delr*sin(abs(theta))... 
                    yy0]; 
 
                line(massX,massY); 
 
                xlim([-3,10]); ylim([-10,3]); 
 
                set(gca,'DataAspectRatioMode','manual','DataAspectRatio',[1 1 1]); 
                drawnow; 
                pause(.3); 
            end 
 
        end 
 
    end 
end