function varargout = Mass_Spring_System(varargin)
% MASS_SPRING_SYSTEM M-file for Mass_Spring_System.fig
%      MASS_SPRING_SYSTEM, by itself, creates a new MASS_SPRING_SYSTEM or raises the existing
%      singleton*.
%
%      H = MASS_SPRING_SYSTEM returns the handle to a new MASS_SPRING_SYSTEM or the handle to
%      the existing singleton*.
%
% How to use the GUI
%    - Use the sliders to change the values of the parameters c/m,
%    k/m or N.
%    - Click on the button labeled "Infinite" to change the forcing
%    to a delta function. Click on "Finite" to come back to a square
%    wave forcing.
%    - Click on "Close up centered at (x,y) = (1,0)" to zoom in near
%    the origin.
%
% Developed by Benjamin Dyhr, Mathematics Department, University
% of Arizona, Spring 2007.
%

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

global forcing_function
global zoom

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

% Update handles structure
guidata(hObject, handles);

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

set(handles.figure1,'Name',' Mass-Spring System');
zoom = 1;
forcing_function = 2;
[x,cm,km,N,F,f2,x_min,x_max,y_min,y_max]= modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);

% --- Outputs from this function are returned to the command line.
function varargout = Mass_Spring_System_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;



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

% Hint: place code in OpeningFcn to populate Phi_m


% --- Executes on slider movement.
function slider1_Callback(hObject, eventdata, handles)
% hObject    handle to slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global forcing_function
% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

[x,cm,km,N,F,f2,x_min,x_max,y_min,y_max]= modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);

% *** set(handles.text1,'String',['Mode # ' int2str(sl1)]);

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

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end


% --- Executes on slider movement.
function slider2_Callback(hObject, eventdata, handles)
% hObject    handle to slider2 (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,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

cm = get(handles.slider2,'Value');
set(handles.slider2,'Value',cm);

[x,cm,km,N,F,f2,x_min,x_max,y_min,y_max]= modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);

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

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end


% --- Executes on slider movement.
function slider3_Callback(hObject, eventdata, handles)
% hObject    handle to slider3 (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,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

km = get(handles.slider3,'Value');
set(handles.slider3,'Value',km);

[x,cm,km,N,F,f2,x_min,x_max,y_min,y_max]= modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);

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

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end


% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global forcing_function

forcing_function = 1;

[x,cm,km,N,F,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);


% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global forcing_function
forcing_function = 2;

[x,cm,km,N,F,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global zoom;

zoom = 1; 

[x,cm,km,N,F,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);


% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global zoom;

zoom = 2;

[x,cm,km,N,F,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);


function [x,cm,km,N,F,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles)

global forcing_function
global zoom

N = floor(get(handles.slider1,'Value'));
if N == 0
    N = 1;
end
set(handles.text10,'String',['N = ' int2str(N) ]);
set(handles.slider1,'Value',N);

cm = get(handles.slider2,'Value');
set(handles.text11,'String',['c/m = ' num2str(cm)]);

km = get(handles.slider3,'Value');
set(handles.text13,'String',['k/m = ' num2str(km)]);


%hammerblow (normal view)


if zoom == 1 & forcing_function == 1

    
    x_min = 0;
    x_max = 15;
    y_min = -1.1;
    y_max = 3.1;

    x=[x_min:(x_max-x_min)/1000:x_max];

    siz=size(x);
    F=ones(1,siz(2));
    f2=ones(1,siz(2));
    
       %delta function
 
    F = 0*F;
    F(68)=100;
    
        %Definition of Solution
    r1=1;
    r2=1;
         
    if cm^2-4*km < 0
     
    r1= .5*( sqrt(cm^2-4*km) - cm);
    r2= .5*(-sqrt(cm^2-4*km) - cm);
    
    for j=1:siz(2)
        
         if x(j) < 1
            
           f2(j)=0;
            
         else 
             
            f2(j)=    exp(r1*(x(j)-1) )/(-(r2-r1)) - exp(r2*(x(j)-1)  )/(-(r2-r1)) ;
            
         end
         end
    
    elseif cm^2-4*km >= 0 
        
    r1= .5*( sqrt(-cm^2+4*km) - cm);
    r2= .5*(-sqrt(-cm^2+4*km) - cm);     
    
     for j=1:siz(2)
        
        if      x(j) < 1
            
           f2(j)=0;
            
        else
                
              f2(j)= real( exp(r1*(x(j)-1) )/(-(r2-r1))   - exp(r2*(x(j)-1) )/(-(r2-r1))) ;
                                                                                                           
        end
        end
    
    end
    
  
   
%square impulse (normal view)

elseif  zoom == 1 & forcing_function == 2

    x_min = 0;
    x_max = 15;
    y_min = -1.1;
    y_max = 3.1;
            
    x=[x_min:(x_max-x_min)/(100*(x_max-x_min)):x_max];

    siz=size(x);
    F=ones(1,siz(2));
    f2=ones(1,siz(2));
    
    %Definition of Forcing Function
    
    for j=1:siz(2)
        
        if      x(j) <= 1
            
            F(j)=0;
            
        elseif x(j) <= 1+1/N
                
            F(j) = N;
                
        elseif x(j) >  1+1/N

            F(j) = 0;
            
        end
    end
    
    %Definition of Solution
    r1=1;
    r2=1;
    
   if cm^2-4*km < 0
     
    r1= .5*( sqrt(cm^2-4*km) - cm);
    r2= .5*(-sqrt(cm^2-4*km) - cm);
    
         for j=1:siz(2)
        
         if x(j) < 1
            
           f2(j)=0;
            
         elseif x(j) <= 1+1/N
                
         f2(j)=N*(1/(r1*r2) + exp(r1*(x(j)-1))/(r1*(r1-r2)) + exp(r2*(x(j)-1))/(r2*(r2-r1)));
                
         elseif x(j) >  1+1/N

            f2(j)= N*(   exp(r1*(x(j)-1)      )/(r1*(r1-r2)) + exp(r2*(x(j)-1)      )/(r2*(r2-r1)) ...
                      - (exp(r1*(x(j)-(1+1/N)))/(r1*(r1-r2)) + exp(r2*(x(j)-(1+1/N)))/(r2*(r2-r1))) );
            
         end
         end
    
    elseif cm^2-4*km>=0 
        
    r1= .5*( sqrt(-cm^2+4*km) - cm);
    r2= .5*(-sqrt(-cm^2+4*km) - cm);     
    
     for j=1:siz(2)
        
        if      x(j) < 1
            
           f2(j)=0;
            
        elseif x(j) <= 1+1/N
                
         f2(j)=N*real(...
             1/(r1*r2) + exp(r1*(x(j)-1))/(r1*(r1-r2)) + exp(r2*(x(j)-1))/(r2*(r2-r1))...
                                                                                           ) ;
                
        elseif x(j) >  1+1/N

            f2(j)= N*real(...
                   exp(r1*(x(j)-1)      )/(r1*(r1-r2)) + exp(r2*(x(j)-1)      )/(r2*(r2-r1)) ...
                      - (exp(r1*(x(j)-(1+1/N)))/(r1*(r1-r2)) + exp(r2*(x(j)-(1+1/N)))/(r2*(r2-r1))) ...
                                                                                                       );
            
        end
        end
    
    end
    
 





%hammerblow (zoomed)
elseif forcing_function == 1 & zoom==2
   
    x_min = 0.9;
    x_max = 1.1;
    y_min = -0.1;
    y_max = 0.1;

    x=[x_min:(x_max-x_min)/1000:x_max];

    siz=size(x);
    F=ones(1,siz(2));
    f2=ones(1,siz(2));
    
        %delta function
    
    F = 0*F;
    F(501)=100;
    
        %Definition of Solution
    r1=1;
    r2=1;
         
        if cm^2-4*km < 0
     
    r1= .5*( sqrt(cm^2-4*km) - cm);
    r2= .5*(-sqrt(cm^2-4*km) - cm);
    
    for j=1:siz(2)
        
         if x(j) < 1
            
           f2(j)=0;
            
         else 
             
            f2(j)=    exp(r1*(x(j)-1) )/(-(r2-r1)) - exp(r2*(x(j)-1)  )/(-(r2-r1)) ;
            
         end
         end
    
    elseif cm^2-4*km >= 0 
        
    r1= .5*( sqrt(-cm^2+4*km) - cm);
    r2= .5*(-sqrt(-cm^2+4*km) - cm);     
    
     for j=1:siz(2)
        
        if      x(j) < 1
            
           f2(j)=0;
            
        else
                
              f2(j)= real( exp(r1*(x(j)-1) )/(-(r2-r1))   - exp(r2*(x(j)-1) )/(-(r2-r1))) ;
                                                                                                           
        end
        end
    
    end
  
   
%square impulse (zoomed)

elseif forcing_function == 2 & zoom == 2

    x_min = 0.9;
    x_max = 1.1;
    y_min = -0.1;
    y_max = 0.1;
            
    x=[x_min:(x_max-x_min)/(100*(x_max-x_min)):x_max];

    siz=size(x);
    F=ones(1,siz(2));
    f2=ones(1,siz(2));
    
    %Definition of Forcing Function
    
    for j=1:siz(2)
        
        if      x(j) <= 1
            
            F(j)=0;
            
        elseif x(j) <= 1+1/N
                
            F(j) = N;
                
        elseif x(j) >  1+1/N

            F(j) = 0;
            
        end
    end
    
    %Definition of Solution
    r1=1;
    r2=1;
    
   if cm^2-4*km < 0
     
    r1= .5*( sqrt(cm^2-4*km) - cm);
    r2= .5*(-sqrt(cm^2-4*km) - cm);
    
         for j=1:siz(2)
        
         if x(j) < 1
            
           f2(j)=0;
            
         elseif x(j) <= 1+1/N
                
         f2(j)=N*(1/(r1*r2) + exp(r1*(x(j)-1))/(r1*(r1-r2)) + exp(r2*(x(j)-1))/(r2*(r2-r1)));
                
         elseif x(j) >  1+1/N

            f2(j)= N*(   exp(r1*(x(j)-1)      )/(r1*(r1-r2)) + exp(r2*(x(j)-1)      )/(r2*(r2-r1)) ...
                      - (exp(r1*(x(j)-(1+1/N)))/(r1*(r1-r2)) + exp(r2*(x(j)-(1+1/N)))/(r2*(r2-r1))) );
            
         end
         end
    
    elseif cm^2-4*km>=0 
        
    r1= .5*( sqrt(-cm^2+4*km) - cm);
    r2= .5*(-sqrt(-cm^2+4*km) - cm);     
    
     for j=1:siz(2)
        
        if      x(j) < 1
            
           f2(j)=0;
            
        elseif x(j) <= 1+1/N
                
         f2(j)=N*real(...
             1/(r1*r2) + exp(r1*(x(j)-1))/(r1*(r1-r2)) + exp(r2*(x(j)-1))/(r2*(r2-r1))...
                                                                                           ) ;
                
        elseif x(j) >  1+1/N

            f2(j)= N*real(...
                   exp(r1*(x(j)-1)      )/(r1*(r1-r2)) + exp(r2*(x(j)-1)      )/(r2*(r2-r1)) ...
                      - (exp(r1*(x(j)-(1+1/N)))/(r1*(r1-r2)) + exp(r2*(x(j)-(1+1/N)))/(r2*(r2-r1))) ...
                                                                                                       );
            
        end
        end
    
    end 
    
end


function plotmodes(hObject, eventdata, handles)

[x,cm,km,N,F,f2,x_min,x_max,y_min,y_max]=modes(hObject, eventdata, handles);

axes(handles.graph1)
plot(x,F,'-r',x,f2,'-b')
axis([x_min x_max y_min y_max])
