Lab 1, Simulate spring pendulum. GUI application using Matlab

We are asked to write a program to simulation the equation of motion of the following spring pendulum shown in the textbook on page 42

4.1.2 Code

Here is 45 MB AVI movie showing the program I wrote to solve this problem. pendulum.avi

To run the program, download this code.zip and extract it. It will create a new folder called code

Notice that the \(r\) generalized coordinate is measured from the corner where the spring is attached to and not from the static equilibrium position. This diagram below helps to illustrate this more.

Now the equations of motion are found and then they are converted to state form description and then ODE solver was used to integrate them.

4.1.3 Mathematical model

This system is 2 degree of freedom system, since two generalized coordinates are required to uniquely locate the mass \(m\) at any point of time. Hence, two EQM’s (equation of motions) are required.

\[ -T+mg\cos \theta =m(r^{\prime \prime }-r\left ( \theta ^{\prime }\right ) ^{2}) \]

\begin{equation} r^{\prime \prime }=r\left ( \theta ^{\prime }\right ) ^{2}+g\cos \theta -\frac{k}{m}\left ( r-l_{0}\right ) \tag{1} \end{equation}

Applying \(F=ma\) in the perpendicular direction to the spring results in the following EQM

\begin{align} -mg\sin \theta & =m\left ( r\theta ^{\prime \prime }+2r^{\prime }\theta ^{\prime }\right ) \nonumber \\ \theta ^{\prime \prime } & =-\frac{g}{r}\sin \theta -\frac{2r^{\prime }\theta ^{\prime }}{r}\tag{2} \end{align}

To use the ODE solver, the above second order ODE’s, Eqs., (1) and (2), need to be converted to ﬁrst order ODE. This is done by introducing state variables as follows

\[ \left . \begin{array} [c]{c}x_{1}=r\\ x_{2}=r^{\prime }\\ x_{3}=\theta \\ x_{4}=\theta ^{\prime }\end{array} \right \} \overset{\frac{d}{dt}}{\rightarrow }\left . \begin{array} [c]{l}x_{1}^{\prime }=x_{2}\\ x_{2}^{\prime }=x_{1}x_{4}^{2}+g\cos x_{3}-\frac{k}{m}\left ( x_{1}-l_{0}\right ) \\ x_{3}^{\prime }=x_{4}\\ x_{4}^{\prime }=-\frac{g}{r}\sin x_{3}-\frac{2x_{2}x_{4}}{x_{1}}\end{array} \right \} \]

\begin{equation} \begin{pmatrix} x_{1}^{\prime }\\ x_{2}^{\prime }\\ x_{3}^{\prime }\\ x_{4}^{\prime }\end{pmatrix} =\begin{pmatrix} x_{2}\\ x_{1}x_{4}^{2}+g\cos x_{3}-\frac{k}{m}\left ( x_{1}-l_{0}\right ) \\ x_{4}\\ \frac{g}{r}\sin x_{3}-\frac{2x_{2}x_{4}}{x_{1}}\end{pmatrix} \tag{3} \end{equation}

The above system is solved using ODE solver and the result are shown in the following section.

4.1.4 Result of simulation

Using the following parameters to generate ﬁgure (3.2) and ﬁgure 3.3 for the very stiﬀ spring

4.1.5 Discussion of results

When the spring is made very stiﬀ, then the change of the pendulum length can be seen from the ﬁgure 3.3 result to be very small. This is the same result as a normal pendulum will have. This is also veriﬁed by ﬁgure 3.2 which shows a simple harmonic motion as shown in ﬁgure 3.2 (pendulum angle vs. time).

When the spring stiﬀness is reduced (\(f=1Hz\)), then we can see that the motion is no longer a simple harmonic motion as can be seen in ﬁgure (3.3 and 3.4). Changing the length of the initial pendulum also resulted is completely diﬀerent proﬁle of the motion.

4.1.6 Appendix, Surce code listing

function varargout = eme_121_lab1(varargin)
% EME_121_LAB1 M-file for eme_121_lab1.fig
%      EME_121_LAB1, by itself, creates a new
%      EME_121_LAB1 or raises the existing singleton*.
%
%      H = EME_121_LAB1 returns the handle to a new
%          EME_121_LAB1 or the handle to the existing singleton*.
%
%      EME_121_LAB1('CALLBACK',hObject,eventData,handles,...)
%        calls the local function named CALLBACK in EME_121_LAB1.M
%        with the given input arguments.
%
%      EME_121_LAB1('Property','Value',...) creates a new
%           EME_121_LAB1 or raises the existing singleton*.
%           Starting from the left, property value pairs are applied
%           to the GUI before eme_121_lab1_OpeningFcn gets called.  An
%           unrecognized property name or invalid value makes
%           property application stop.  All inputs are passed
%           to eme_121_lab1_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.
%       Choose "GUI allows only one instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help eme_121_lab1

% Last Modified by GUIDE v2.5 10-Apr-2011 23:32:56

% By Nasser M. Abbasi

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @eme_121_lab1_OpeningFcn, ...
    'gui_OutputFcn',  @eme_121_lab1_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 eme_121_lab1 is made visible.
function eme_121_lab1_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 eme_121_lab1 (see VARARGIN)

% Choose default command line output for eme_121_lab1
handles.output = hObject;

set(handles.figure1, 'UserData',[]);
set(handles.figure1,'Name','EME 121 first lab by Nasser M. Abbasi');

data.stop = false;
set(handles.figure1, 'UserData',data);

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes eme_121_lab1 wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = eme_121_lab1_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 mass_tag_Callback(hObject, eventdata, handles)
% hObject    handle to mass_tag (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 mass_tag as text
% str2double(get(hObject,'String')) returns contents of mass_tag as a double

% --- Executes during object creation, after setting all properties.
function mass_tag_CreateFcn(hObject, eventdata, handles)
% hObject    handle to mass_tag (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 f_tag_Callback(hObject, eventdata, handles)
% hObject    handle to f_tag (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 f_tag as text
% str2double(get(hObject,'String')) returns contents of f_tag as a double

% --- Executes during object creation, after setting all properties.
function f_tag_CreateFcn(hObject, eventdata, handles)
% hObject    handle to f_tag (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 length_tag_Callback(hObject, eventdata, handles)
% hObject    handle to length_tag (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 length_tag as text
% str2double(get(hObject,'String')) returns contents of length_tag as a double

% --- Executes during object creation, after setting all properties.
function length_tag_CreateFcn(hObject, eventdata, handles)
% hObject    handle to length_tag (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 angle_zero_tag_Callback(hObject, eventdata, handles)
% hObject    handle to angle_zero_tag (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 angle_zero_tag as text
% str2double(get(hObject,'String')) returns contents of angle_zero_tag as a double

% --- Executes during object creation, after setting all properties.
function angle_zero_tag_CreateFcn(hObject, eventdata, handles)
% hObject    handle to angle_zero_tag (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 angle_speed_zero_Callback(hObject, eventdata, handles)
% hObject    handle to angle_speed_zero (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 angle_speed_zero as text
%        str2double(get(hObject,'String')) returns contents of angle_speed_zero as a double

% --- Executes during object creation, after setting all properties.
function angle_speed_zero_CreateFcn(hObject, eventdata, handles)
% hObject    handle to angle_speed_zero (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 r_zero_tag_Callback(hObject, eventdata, handles)
% hObject    handle to r_zero_tag (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 r_zero_tag as text
%        str2double(get(hObject,'String')) returns contents of r_zero_tag as a double

% --- Executes during object creation, after setting all properties.
function r_zero_tag_CreateFcn(hObject, eventdata, handles)
% hObject    handle to r_zero_tag (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 r_speed_zero_Callback(hObject, eventdata, handles)
% hObject    handle to r_speed_zero (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 r_speed_zero as text
%        str2double(get(hObject,'String')) returns contents of r_speed_zero as a double

% --- Executes during object creation, after setting all properties.
function r_speed_zero_CreateFcn(hObject, eventdata, handles)
% hObject    handle to r_speed_zero (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 run_btn.
function run_btn_Callback(hObject, eventdata, handles)
% hObject    handle to run_btn (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[data,status] = parse_input(handles);
if not(status)
    return;
end

data.handles = handles;
data.g = 9.8;

enable_all(handles,'off');
userData = get(handles.figure1, 'UserData');
userData.stop = false;
set(handles.figure1, 'UserData',userData);

[g_msg,g_status]=process_eme_121_lab1(data);
if not(g_status)
    uiwait(errordlg(sprintf('Processing terminated: %s',g_msg),...
        'Bad Input', 'modal'));
    uicontrol(handles.f_tag);
end
enable_all(handles,'on');

%----------------------------
function [data,status]= parse_input(handles)

[data.mass,status]=verify_valid_positive_numeric...
    (get(handles.mass_tag,'String'),handles.mass_tag,...
    'Mass must be positive number');
if not(status)
    return;
end

[data.f,status]=verify_valid_positive_numeric...
    (get(handles.f_tag,'String'),handles.f_tag,...
    'natural frequency must be positive number');
if not(status)
    return;
end

[data.L,status]=verify_valid_positive_numeric...
    (get(handles.length_tag,'String'),handles.length_tag,...
    'pendulum length must be positive number');
if not(status)
    return;
end

[data.angle_zero,status]=verify_valid_numeric...
    (get(handles.angle_zero_tag,'String'),handles.angle_zero_tag,...
    'Initial angle must be numerical value');
if not(status)
    return;
end

if abs(data.angle_zero)>180
    uiwait(errordlg('Initial angle must be between 0 and 180 degrees only',...
        'Bad Input', 'modal'));
    status = 0;
    uicontrol(handles.angle_zero);
    return
else
    data.angle_zero = data.angle_zero*pi/180;
end

[data.r_zero,status]=verify_valid_positive_numeric...
    (get(handles.r_zero_tag,'String'),handles.r_zero_tag,...
    'initial r must be  numerical value');
if not(status)
    return;
end

[data.r_speed_zero,status]=verify_valid_numeric...
    (get(handles.r_speed_zero,'String'),handles.r_speed_zero,...
    'initial r speed must be numerical value');
if not(status)
    return;
end

[data.time_step,status]=verify_valid_numeric...
    (get(handles.time_Step_tag,'String'),handles.time_Step_tag,...
    'time step must be zero or larger');
if not(status)
    return;
end

[data.angle_speed_zero,status]=verify_valid_numeric...
    (get(handles.angle_speed_zero,'String'),handles.angle_speed_zero,...
    'initial angle speed must be numerical value');
if not(status)
    return;
end

[data.max_t,status]=verify_valid_positive_numeric...
    (get(handles.max_t_tag,'String'),handles.max_t_tag,...
    'maximum simulation time must be positive number');
if not(status)
    return;
end

contents = cellstr(get(handles.ode_solver_tag,'String'));
data.solver = contents{get(handles.ode_solver_tag,'Value')};

data.k = data.mass*(2*pi*data.f)^2;

%-------------------------------
function enable_all(handles,to)
set(handles.mass_tag,'Enable',to);
set(handles.f_tag,'Enable',to);
set(handles.length_tag,'Enable',to);
set(handles.angle_zero_tag,'Enable',to);
set(handles.r_speed_zero,'Enable',to);
set(handles.angle_speed_zero,'Enable',to);
set(handles.r_zero_tag,'Enable',to);
set(handles.max_t_tag,'Enable',to);
set(handles.ode_solver_tag,'Enable',to);
set(handles.time_Step_tag,'Enable',to);

function max_t_tag_Callback(hObject, eventdata, handles)
% hObject    handle to max_t_tag (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 max_t_tag as text
%        str2double(get(hObject,'String')) returns contents of max_t_tag as a double

% --- Executes during object creation, after setting all properties.
function max_t_tag_CreateFcn(hObject, eventdata, handles)
% hObject    handle to max_t_tag (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 reset_tag.
function reset_tag_Callback(hObject, eventdata, handles)
% hObject    handle to reset_tag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

data = get(handles.figure1, 'UserData');
data.stop = true;
set(handles.figure1, 'UserData',data);
enable_all(handles,'on');

% --- Executes on selection change in ode_solver_tag.
function ode_solver_tag_Callback(hObject, eventdata, handles)
% hObject    handle to ode_solver_tag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns ode_solver_tag contents as cell array
%        contents{get(hObject,'Value')} returns selected item from ode_solver_tag

% --- Executes during object creation, after setting all properties.
function ode_solver_tag_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ode_solver_tag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu 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 time_Step_tag_Callback(hObject, eventdata, handles)
% hObject    handle to time_Step_tag (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 time_Step_tag as text
%        str2double(get(hObject,'String')) returns contents of time_Step_tag as a double

% --- Executes during object creation, after setting all properties.
function time_Step_tag_CreateFcn(hObject, eventdata, handles)
% hObject    handle to time_Step_tag (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]=process_eme_121_lab1(data)
% This function is called by the Matlab GUI to solve
% the equation of motion for Lab 1 problem, EME 121,
% spring 2011, UC Davis.
%
% INPUT: data, a record which contains the problem parameters
%
% The function will call ode45 solver and plot the solution
%
% by Nasser M. Abbasi
% 4/12/2011

% Reset axes for plotting
cla(data.handles.axes,'reset');
cla(data.handles.time_angle_axes,'reset');
cla(data.handles.time_r_axes,'reset');
cla(data.handles.angle_r_axes,'reset');

%Check if using auto step size of an explicit step size. see
%GUI for input choice.
if data.time_step==0
    time_span=[0 data.max_t];
else
    time_span=0:data.time_step/1000:data.max_t;
end

%set the y-axis for plotting
max_y = data.L;

%set the initial conditions for ode45
IC = [data.r_zero;data.r_speed_zero;data.angle_zero;data.angle_speed_zero];
options = odeset('OutputFcn',@output);

g_status = true;
g_msg='';

%call  the solver
[t,x] = feval(data.solver,@rhs,time_span,IC,options,data);

%plot the result
N=length(t);
picks=1:N;

set(data.handles.figure1, 'CurrentAxes',data.handles.time_angle_axes);
angle=x(:,3);
plot(t(picks),angle(picks)*180/pi,'-');
title('pendulum angle vs time');
xlabel('Time (sec)');
ylabel('angle (degree)');
grid;
drawnow;

set(data.handles.figure1, 'CurrentAxes',data.handles.time_r_axes);
r=x(:,1);
stairs(t(picks),r(picks),'-');
title('pendulum length vs time');
xlabel('Time (sec)'); ylabel('length (m)');
grid;
drawnow;

set(data.handles.figure1, 'CurrentAxes',data.handles.angle_r_axes);

plot(angle(picks)*180/pi,r(picks),'-');
title('pendulum angle vs length');
xlabel('Angle (degree)'); ylabel('length (m)');
grid;
drawnow;

%-----------------------------------
%
%-----------------------------------
    function dx=rhs(~,x,data)
        %the ode45 function
        temp = x(1)-data.L;

        if abs(temp<2*eps)
            term = 0;
        else
            term = data.k/data.mass*temp;
        end

        dx=[x(2);
            x(1)*x(4)^2+data.g*cos(x(3))-term;
            x(4);
            -2*x(2)*x(4)/x(1)- data.g/x(1)*sin(x(3))
            ];
        if x(1)>max_y
            max_y = x(1)+0.2*x(1);
        end
    end

%-----------------------------------
%
%-----------------------------------
    function status= output(t,x,~,data)
        %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;
        else
            status = false;
            if not(isempty(t))
                plot_solution(x(1,end),x(3,end),t(1),data);
                %pause(0.1);
                if x(1,end)>10*data.L  %spring stretched too far
                    status = true;
                    g_status = false;
                    g_msg ='spring over stretched';
                end
            end
        end
    end

%-----------------------------------
%
%-----------------------------------
    function plot_solution(r,theta,t,data)
        %plot the final solution
        set(0,'CurrentFigure',data.handles.figure1);
        set(data.handles.figure1, 'CurrentAxes',data.handles.axes);

        [x_local,y] = pol2cart(theta-pi/2,r);
        plot(x_local,y,'o','MarkerSize',20,'MarkerEdgeColor','r',...
            'MarkerFaceColor','r','LineWidth',2);
        [x_local,y] = nma_spring.make(r,theta-(pi/2),30);
        hold on;
        line(x_local,y);
        plot(0,0,'o','MarkerSize',5,'MarkerEdgeColor','k',...
            'MarkerFaceColor','k','LineWidth',1);
        title(sprintf(...
            'time = %3.4f, angle(degree) = %3.3f, r = %3.3f, k=%3.4f',...
            t,theta*180/pi,r,data.k));

        xlim([-max_y max_y]);
        ylim([-max_y max_y]);
        %refreshdata
        grid;
        drawnow;
        hold off;
    end

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

l.198 — TeX4ht warning — “SaveEverypar’s: 2 at “begindocument and 3 “enddocument —


\(m\)	\(f\)	\(l_{0}\)	\(r\left ( 0\right ) \)	\(r^{\prime }\left ( 0\right ) \)	\(\theta \left ( 0\right ) \)	\(\theta ^{\prime }\left ( 0\right ) \)

\(5kg\)	\(1000\ Hz\)	\(1\ m\)	\(1\ m\)	\(0\)	\(45^{0}\)	\(0\)


\(m\)	\(f\)	\(l_{0}\)	\(r\left ( 0\right ) \)	\(r^{\prime }\left ( 0\right ) \)	\(\theta \left ( 0\right ) \)	\(\theta ^{\prime }\left ( 0\right ) \)

\(5kg\)	\(1\ Hz\)	\(1\ m\)	\(1\ m\)	\(0\)	\(45^{0}\)	\(0\)


\(m\)	\(f\)	\(l_{0}\)	\(r\left ( 0\right ) \)	\(r^{\prime }\left ( 0\right ) \)	\(\theta \left ( 0\right ) \)	\(\theta ^{\prime }\left ( 0\right ) \)

\(5kg\)	\(1\ Hz\)	\(1\ m\)	\(2\ m\)	\(0\)	\(55^{0}\)	\(0\)

4.1 Lab 1, Simulate spring pendulum. GUI application using Matlab

4.1.1 Problem description

4.1.2 Code

4.1.3 Mathematical model

4.1.4 Result of simulation

4.1.5 Discussion of results

4.1.6 Appendix, Surce code listing