function varargout = SB2013(varargin)
% SB2013 MATLAB code for SB2013.fig
%      SB2013, by itself, creates a new SB2013 or raises the existing
%      singleton*.
%
%      H = SB2013 returns the handle to a new SB2013 or the handle to
%      the existing singleton*.
%
%      SB2013('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SB2013.M with the given input arguments.
%
%      SB2013('Property','Value',...) creates a new SB2013 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before SB2013_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to SB2013_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 SB2013

% Last Modified by GUIDE v2.5 19-Oct-2015 17:07:52

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @SB2013_OpeningFcn, ...
    'gui_OutputFcn',  @SB2013_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 SB2013 is made visible.
function SB2013_OpeningFcn(hObject, eventdata, handles, varargin) %#ok<*INUSL>
% This function has no outputfig 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 SB2013 (see VARARGIN)

handles.output = hObject;

handles.K = 0;
handles.M = 0;
handles.T = 0;
handles.Disp = 0;
handles.Vel = 0;
handles.ACC = 0;
handles.Strain = 0;
handles.Stress = 0;

handles.GUILOCATION = pwd;

handles.FileNameTEMP = 'Default.mat';
c = strrep(handles.FileNameTEMP, '.mat', '');
set(handles.FileSaveName, 'string', c)
load([handles.GUILOCATION, '\Saved Inputs\', handles.FileNameTEMP])

handles.LocationControl_valTEMP = LocationControl_val;

handles.EarthquakeFile = EarthquakeFile;
set(handles.Location, 'string', handles.EarthquakeFile)

handles.EQ = EQ;
handles.TTEMP = TTEMP;
handles.gTEMP = g;
if g == 1
    set(handles.g, 'value', 0)
else
    set(handles.g, 'value', 1)
end

handles.unitTEMP = unit;
if unit == 1
    set(handles.meter, 'value', 1)
else
    set(handles.foot, 'value', 1)
end

handles.dtTEMP = dt;
set(handles.dt, 'string', num2str(handles.dtTEMP))

handles.ScaleTEMP = Scale;
set(handles.Scale, 'string', num2str(handles.ScaleTEMP))

handles.analysis = analysis;
if analysis == 1
    set(handles.Linear, 'value', 1)
    set(handles.su, 'enable', 'off')
    set(handles.GammaMax, 'enable', 'off')
    set(handles.MatProfile, 'enable', 'off')
    set(handles.suMulti, 'enable', 'off')
    set(handles.GammaMaxMulti, 'enable', 'off')
    set(handles.NYS, 'enable', 'off')
else
    set(handles.Nonlinear, 'value', 1)
    set(handles.su, 'enable', 'on')
    set(handles.GammaMax, 'enable', 'on')
    set(handles.MatProfile, 'enable', 'on')
    set(handles.suMulti, 'enable', 'on')
    set(handles.GammaMaxMulti, 'enable', 'on')
    set(handles.NYS, 'enable', 'on')
end

handles.NYSTEMP = NYS;
set(handles.NYS, 'string', num2str(handles.NYSTEMP))

handles.TolTEMP = Tol;
set(handles.Tol, 'string', num2str(handles.TolTEMP))

handles.BaseTEMP = Base;
if Base == 0
    set(handles.RigidBase, 'value', 1)
    set(handles.pr, 'enable', 'off')
    set(handles.Vsr, 'enable', 'off')
    set(handles.pr, 'string', '')
    set(handles.Vsr, 'string', '')
    set(handles.Incident, 'enable', 'off')
    
    set(handles.LucoDamping, 'enable', 'on')
    handles.DampType = DampType;
    
    if DampType == 1
        set(handles.RayleighDamping, 'Value', 1)
    else
        set(handles.LucoDamping, 'Value', 1)
    end
else
    
    set(handles.FlexibleBase, 'value', 1)
    set(handles.pr, 'enable', 'on')
    set(handles.Vsr, 'enable', 'on')
    set(handles.pr, 'string', num2str(pr))
    set(handles.Vsr, 'string', num2str(Vsr))
    set(handles.Incident, 'enable', 'on')
    
    set(handles.RayleighDamping, 'Value', 1)
    set(handles.LucoDamping, 'enable', 'off')
    
    handles.DampType = 1;
end

handles.IncidentTEMP =  Incident;
set(handles.Incident, 'value', handles.IncidentTEMP)

handles.hTEMP = h;
set(handles.h, 'string', num2str(handles.hTEMP))

handles.NTEMP = N;
set(handles.N, 'string', num2str(handles.NTEMP))

handles.prTEMP = pr;
handles.VsrTEMP = Vsr;

handles.mode_i_frequencyTEMP = mode_i_frequency;
handles.mode_j_frequencyTEMP = mode_j_frequency;
handles.mode_i_dampingTEMP = mode_i_damping;
handles.mode_j_dampingTEMP = mode_j_damping;
handles.damping = damping;

handles.MatNumTEMP = 1;
set(handles.MatNum, 'String', handles.MatNumTEMP)

handles.MatTEMP = MatList(1,2);
set(handles.Mat, 'String', handles.MatTEMP)

handles.MatList = MatList;
set(handles.MatTable,'Data',MatList)

handles.LayersTEMP = Layers;
set(handles.Layers, 'string', num2str(handles.LayersTEMP))

handles.VsTEMP = Vs;
set(handles.Vs, 'string', num2str(handles.VsTEMP))

handles.suTEMP = su;
set(handles.su, 'string', num2str(handles.suTEMP))

handles.GammaMaxTEMP = GammaMax;
set(handles.GammaMax, 'string', num2str(GammaMax))

handles.denTEMP = den;
set(handles.den, 'string', num2str(den))

handles.MatProfileTEMP = MatProfile;
set(handles.MatProfile, 'string', num2str(handles.MatProfileTEMP))

handles.LayersMultiTEMP = LayersMulti;
set(handles.LayersMulti, 'string', num2str(LayersMulti))

handles.VsMultiTEMP = VsMulti;
set(handles.VsMulti, 'string', num2str(VsMulti))

handles.suMultiTEMP = suMulti;
set(handles.suMulti, 'string', num2str(handles.suMultiTEMP))

handles.GammaMaxMultiTEMP = GammaMaxMulti;
set(handles.GammaMaxMulti, 'string', num2str(GammaMaxMulti))

handles.denMultiTEMP = denMulti;
set(handles.denMulti, 'string', num2str(denMulti))

handles = OutputOff(handles);
PlotModel(N,h,Base,handles)

guidata(hObject, handles);

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

%  -  -  -  Outputs from this function are returned to the command line.
function varargout = SB2013_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning outputfig 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 outputfig from handles structure
varargout{1} = handles.output;


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

handles.FileNameTEMP = FileSaveName;
guidata(hObject, handles);


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

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


function LocationControl_Callback(hObject, eventdata, handles)
% hObject    handle to LocationControl (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 LocationControl contents as cell array
%        contents{get(hObject, 'Value')} returns selected item from LocationControl
handles.output = hObject;
LocationControl_val = get(handles.LocationControl, 'Value');
LocationControl_str = get(handles.LocationControl, 'String');
switch LocationControl_str{LocationControl_val}
    case '1. User Defined'
        try
            [a b] = uigetfile('.txt', 'Select Earthquake File', [handles.GUILOCATION, '\Earthquake Data']);
            EQ = load(strcat(b, a));
            [m n] = size(EQ);
            if min(m, n) == 1
                set(handles.dt, 'enable', 'on')
                handles.EQ = EQ;
                handles.TTEMP = 0;
            elseif min(m, n) == 2
                handles.EQ = EQ(:, 2);
                handles.TTEMP = EQ(:, 1);
                handles.dtTEMP = handles.TTEMP(2) - handles.TTEMP(1);
                set(handles.dt, 'string', handles.dtTEMP)
                set(handles.dt, 'enable', 'off')
            else
                errordlg('There should be at most 2 columns. If there is only one column, then it is the Earthquake motion. If there are 2 columns, then the first column is time and second column is the Earthquake motion.',...
                    'Please check format of Earthquake Data.')
                return
            end
            
            set(handles.Location, 'string', a)
            handles.EarthquakeFile = a;
        catch %#ok<*CTCH>
            errordlg('Error loading Earthquake Data.',...
                'Error loading Earthquake Data')
            return
        end
        
    case '2. El Centro'
        set(handles.Location, 'string', 'elcentro.txt')
        handles.EQ = load([handles.GUILOCATION, '\Earthquake Data\elcentro.txt']);
        handles.TTEMP = handles.EQ(:,1);
        handles.EQ = handles.EQ(:,2);
        handles.dtTEMP = handles.TTEMP(2) - handles.TTEMP(1);
        set(handles.dt, 'string', handles.dtTEMP)
        set(handles.dt, 'enable', 'off')
end

handles = OutputOff(handles);

handles.LocationControl_valTEMP = LocationControl_val;
guidata(hObject, handles);

%  -  -  -  Executes during object creation, after setting all properties.
function LocationControl_CreateFcn(hObject, eventdata, handles)
% hObject    handle to LocationControl (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

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

% Hint: get(hObject, 'Value') returns toggle state of g
g = get(hObject, 'Value');

if g == 0
    handles.output = hObject;
    if get(handles.meter, 'value') == 1
        unit = 1;
    else
        unit = 2;
    end
    g = 1;
else
    if get(handles.meter, 'value') == 1
        unit = 1;
        g = 9.81;
    else
        unit = 2;
        g = 32.2*12;
    end
end

handles = OutputOff(handles);

handles.unitTEMP = unit;
handles.gTEMP = g;
guidata(hObject, handles);

%  -  -  -  Executes when selected object is changed in AccelerationControl.
function AccelerationControl_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in AccelerationControl
% eventdata  structure with the following fields (see UIBUTTONGROUP)
%	EventName: string 'SelectionChanged' (read only)
%	OldValue: handle of the previously selected object or empty if none was selected
%	NewValue: handle of the currently selected object
% handles    structure with handles and user data (see GUIDATA)
handles.output = hObject;
g = handles.gTEMP;
switch get(eventdata.NewValue, 'Tag')
    case 'meter'
        unit = 1;
        if g ~= 1
            g = 9.81;
        end
        set(handles.hstr, 'string', 'Height of ele (m)')
        set(handles.RockDenstr, 'string', 'Density of Rock (kg/m^3)')
        set(handles.Vsrstr, 'string', 'Vs of Rock (m/s)')
        set(handles.Vsstr, 'string', 'Shear Wave Velocity (m/s)')
        set(handles.sustr, 'string', 'Shear Strength (kPa)')
        set(handles.GammaMaxstr, 'string', 'Max Shear Strain (%)')
        set(handles.denstr, 'string', 'Mass Density (kg/m^3)')
    case 'foot'
        unit = 2;
        if g ~= 1
            g = 32.2*12;
        end
        set(handles.hstr, 'string', 'Height of ele (in)')
        set(handles.RockDenstr, 'string', 'Density of Rock (slug/in^3)')
        set(handles.Vsrstr, 'string', 'Vs of Rock (in/s)')
        set(handles.Vsstr, 'string', 'Shear Wave Velocity (in/s)')
        set(handles.sustr, 'string', 'Shear Strength (ksi)')
        set(handles.GammaMaxstr, 'string', 'Max Shear Strain (%)')
        set(handles.denstr, 'string', 'Mass Density (slug/in^3)')
end

handles = OutputOff(handles);

h = handles.hTEMP;
N = handles.NTEMP;
Base = handles.BaseTEMP;

plot(zeros(N, 1), 0:h:(N - 1)*h, 'o')
hold on
plot(zeros(N + 1, 1), 0:h:N*h, '-')
axis([ - 0.05 0.05 0 N*h]);
set(gca, 'YDir', 'reverse')
set(gca, 'xtick', []);

if unit == 1
    ylabel('depth (m)')
else
    ylabel('depth (in.)')
end

if Base == 1
    plot(0, (N*h), 'o')
end
hold off

handles.unitTEMP = unit;
handles.gTEMP = g;
guidata(hObject, handles);


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

if isempty(dt) == 1
    errordlg('Time step for EQ data should be a number.',...
        'Invalid Entry for dt (of EQ data)')
    set(handles.dt, 'String', handles.dtTEMP)
    return
end

if dt ~= handles.dtTEMP
    
    handles = OutputOff(handles);
    
end

handles.dtTEMP = dt;

guidata(hObject, handles)

%  -  -  -  Executes during object creation, after setting all properties.
function dt_CreateFcn(hObject, eventdata, handles)
% hObject    handle to dt (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 Scale_Callback(hObject, eventdata, handles)
% hObject    handle to Scale (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 Scale as text
%        str2double(get(hObject, 'String')) returns contents of Scale as a double
Scale = str2num(get(hObject, 'String'));

if isempty(Scale) == 1
    set(handles.Scale, 'string', '1')
    Scale = 1;
end

if Scale ~= handles.ScaleTEMP
    
    handles = OutputOff(handles);
    
end

handles.ScaleTEMP = Scale;
guidata(hObject, handles)

%  -  -  -  Executes during object creation, after setting all properties.
function Scale_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Scale (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 when selected object is changed in AnalysisControl.
function AnalysisControl_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in AnalysisControl
% eventdata  structure with the following fields (see UIBUTTONGROUP)
%	EventName: string 'SelectionChanged' (read only)
%	OldValue: handle of the previously selected object or empty if none was selected
%	NewValue: handle of the currently selected object
% handles    structure with handles and user data (see GUIDATA)
handles.output = hObject;
str = handles.LayersTEMP;
A = strread(str, '%s', 'delimiter', ' '); %#ok<REMFF1>

Vs = handles.VsTEMP;
su = handles.suTEMP;
suMulti = handles.suMultiTEMP;
GammaMax = handles.GammaMaxTEMP;
GammaMaxMulti = handles.GammaMaxMultiTEMP;
MatProfile = handles.MatProfileTEMP;
NYS = handles.NYSTEMP;

switch get(eventdata.NewValue, 'Tag')
    case 'Linear'
        analysis = 1;
        
        set(handles.su, 'enable', 'off')
        set(handles.su, 'string', num2str(ones(1, length(handles.VsTEMP))))
        handles.suTEMP = ones(1, length(handles.VsTEMP));
        
        set(handles.GammaMax, 'enable', 'off')
        set(handles.GammaMax, 'string', num2str(zeros(1, length(handles.VsTEMP))))
        handles.GammaMaxTEMP = zeros(1, length(handles.VsTEMP));
        
        set(handles.MatProfile, 'enable', 'off')
        set(handles.MatProfile, 'string', num2str(zeros(1, length(handles.VsTEMP))))
        handles.MatProfileTEMP = zeros(1, length(handles.VsTEMP));
        
        set(handles.suMulti, 'enable', 'off')
        set(handles.suMulti, 'string', '15')
        handles.suMultiTEMP = 15;
        
        set(handles.GammaMaxMulti, 'enable', 'off')
        set(handles.GammaMaxMulti, 'string', '0')
        handles.GammaMaxMultiTEMP = 0;
        
        set(handles.NYS, 'enable', 'off')
        set(handles.NYS, 'string', num2str(1))
        handles.NYSTEMP = 1;
        
    case 'Nonlinear'
        analysis = 2;
        
        set(handles.su, 'enable', 'on')
        set(handles.su, 'string', num2str(75*ones(1, length(handles.VsTEMP))))
        handles.suTEMP = 75*ones(1, length(handles.VsTEMP));
        
        set(handles.GammaMax, 'enable', 'on')
        set(handles.GammaMax, 'string', num2str(3*ones(1, length(handles.VsTEMP))))
        handles.GammaMaxTEMP = 3*ones(1, length(handles.VsTEMP));
        
        set(handles.MatProfile, 'enable', 'on')
        set(handles.MatProfile, 'string', num2str(MatProfile))
        handles.MatProfileTEMP = MatProfile;
        
        set(handles.suMulti, 'enable', 'on')
        set(handles.suMulti, 'string', '0')
        handles.suMultiTEMP = 0;
        
        set(handles.GammaMaxMulti, 'enable', 'on')
        set(handles.GammaMaxMulti, 'string', '0')
        handles.GammaMaxMultiTEMP = 0;
        
        set(handles.NYS, 'enable', 'on')
        set(handles.NYS, 'string', num2str(20))
        handles.NYSTEMP = 20;
end

handles = OutputOff(handles);

handles.analysis = analysis;
guidata(hObject, handles);


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

if isempty(NYS) == 1
    errordlg('NYS should be a number.', 'Invalid Entry for NYS')
    set(handles.NYS, 'String', handles.NYSTEMP)
    return
end

if rem(NYS,1) ~= 0
    errordlg('NYS should be an integer.', 'Invalid Entry for #')
    set(handles.NYS, 'String', handles.NYSTEMP)
    return
end

if NYS < 1
    errordlg('NYS should be 1 or greater. For linear analysis the shear strength will be set extremely high so it will not be reached.')
    set(handles.NYS, 'String', handles.NYSTEMP)
    return
end


if NYS ~= handles.NYSTEMP
    
    handles = OutputOff(handles);
    
end

handles.NYSTEMP = NYS;
guidata(hObject, handles)


%  -  -  -  Executes during object creation, after setting all properties.
function NYS_CreateFcn(hObject, eventdata, handles)
% hObject    handle to NYS (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 Tol_Callback(hObject, eventdata, handles)
% hObject    handle to Tol (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 Tol as text
%        str2double(get(hObject, 'String')) returns contents of Tol as a double
Tol = str2num(get(hObject, 'String'));

if isempty(Tol) == 1
    set(handles.Tol, 'string', '1e - 005')
    Tol = 1*10^ - 5;
end

if Tol ~= handles.TolTEMP
    
    handles = OutputOff(handles);
    
end

handles.TolTEMP = Tol;
guidata(hObject, handles)


%  -  -  -  Executes during object creation, after setting all properties.
function Tol_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Tol (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 when selected object is changed in Base.
function Base_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in Base
% eventdata  structure with the following fields (see UIBUTTONGROUP)
%	EventName: string 'SelectionChanged' (read only)
%	OldValue: handle of the previously selected object or empty if none was selected
%	NewValue: handle of the currently selected object
% handles    structure with handles and user data (see GUIDATA)
handles.output = hObject;
switch get(eventdata.NewValue, 'Tag')
    case 'RigidBase'
        
        handles.BaseTEMP = 0;
        set(handles.pr, 'enable', 'off')
        set(handles.Vsr, 'enable', 'off')
        set(handles.pr, 'string', '')
        set(handles.Vsr, 'string', '')
        set(handles.Incident, 'enable', 'off')
        set(handles.Incident, 'value', 0)
        handles.IncidentTEMP = 0;
        
        set(handles.LucoDamping, 'Enable', 'on')
        
    case 'FlexibleBase'
        
        handles.BaseTEMP = 1;
        set(handles.pr, 'enable', 'on')
        set(handles.Vsr, 'enable', 'on')
        set(handles.pr, 'string', '2000')
        set(handles.Vsr, 'string', '1000')
        handles.prTEMP = 2000;
        handles.VsrTEMP = 1000;
        %set(handles.Incident, 'enable', 'on')
        %set(handles.Incident, 'value', handles.IncidentTEMP)
        
        handles.DampType = 1;
        set(handles.RayleighDamping, 'Value', 1)
        set(handles.LucoDamping, 'Enable', 'off')
        
end

handles = OutputOff(handles);

h = handles.hTEMP;
N = handles.NTEMP;
Base = handles.BaseTEMP;

PlotModel(N,h,Base,handles)

guidata(hObject, handles);


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

% Hint: get(hObject, 'Value') returns toggle state of Incident
Incident = get(hObject, 'Value');

handles = OutputOff(handles);

handles.IncidentTEMP = Incident;
guidata(hObject, handles);


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

if isempty(h) == 1
    errordlg('Height of element should be a number.',...
        'Invalid Entry for Height of ele')
    set(handles.h, 'String', handles.hTEMP)
    return
end

if h ~= handles.hTEMP
    
    handles = OutputOff(handles);
    
end

N = handles.NTEMP;
Base = handles.BaseTEMP;

PlotModel(N,h,Base,handles)

handles.hTEMP = h;
guidata(hObject, handles);


%  -  -  -  Executes during object creation, after setting all properties.
function h_CreateFcn(hObject, eventdata, handles)
% hObject    handle to h (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 N_Callback(hObject, eventdata, handles)
% hObject    handle to N (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 N as text
%        str2double(get(hObject, 'String')) returns contents of N as a double
N = str2num(get(hObject, 'String'));

if isempty(N) == 1
    errordlg('Number of elements should be a number.',...
        'Invalid Entry for # of Ele')
    set(handles.N, 'String', handles.NTEMP)
    return
end

if rem(N,1) ~= 0
    errordlg('# of Ele should be an integer.', 'Invalid Entry for #')
    set(handles.N, 'String', handles.NTEMP)
    return
end

if N < 1
    errordlg('# of Ele should be 1 or greater.')
    set(handles.N, 'String', handles.NTEMP)
    return
end

if N ~= handles.NTEMP
    
    handles = OutputOff(handles);
    
end

h = handles.hTEMP;
Base = handles.BaseTEMP;

PlotModel(N,h,Base,handles)

handles.NTEMP = N;
guidata(hObject, handles);

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

pr = str2num(get(hObject, 'String'));

if isempty(pr) == 1
    errordlg('Density of Rock should be a number.',...
        'Invalid Entry for Density of Rock')
    set(handles.pr, 'String', handles.prTEMP)
    return
end

if pr ~= handles.prTEMP
    
    handles = OutputOff(handles);
    
end

handles.prTEMP = pr;
guidata(hObject, handles)


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

Vsr = str2num(get(hObject, 'String'));

if isempty(Vsr) == 1
    errordlg('Shear Velocity of Rock should be a number.',...
        'Invalid Entry for Velocity (for Rock)')
    set(handles.Vsr, 'String', handles.VsrTEMP)
    return
end

if Vsr ~= handles.VsrTEMP
    
    handles = OutputOff(handles);
    
end

handles.VsrTEMP = Vsr;
guidata(hObject, handles)

%  -  -  -  Executes during object creation, after setting all properties.
function Vsr_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Vsr (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 when selected object is changed in Base.
function Damping_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in Base
% eventdata  structure with the following fields (see UIBUTTONGROUP)
%	EventName: string 'SelectionChanged' (read only)
%	OldValue: handle of the previously selected object or empty if none was selected
%	NewValue: handle of the currently selected object
% handles    structure with handles and user data (see GUIDATA)
handles.output = hObject;
switch get(eventdata.NewValue, 'Tag')
    case 'RayleighDamping'
        handles.DampType = 1;
    case 'LucoDamping'
        handles.DampType = 2;
end

handles = OutputOff(handles);

guidata(hObject, handles);


% --- Executes on button press in CDC.
function CDC_Callback(hObject, eventdata, handles)

if handles.DampType == 1
    
    options.Resize='on';
    options.WindowStyle='normal';
    
    prompt = {'Frequency #1:','Frequency #1:','Damping at Frequency #1:',...
        'Damping at Frequency #2:'};
    dlg_title = 'Enter Rayleigh Damping Coefficients';
    num_lines = [1 50; 1 50; 1 50; 1 50];
    def = {num2str(handles.mode_i_frequencyTEMP), num2str(handles.mode_j_frequencyTEMP),...
        num2str(handles.mode_i_dampingTEMP),num2str(handles.mode_j_dampingTEMP)};
    answer = inputdlg(prompt,dlg_title,num_lines,def,options);
    
    if isempty(answer)
        return
    end
    
    mode_i_frequency = str2num(cell2mat(answer(1)));
    mode_j_frequency = str2num(cell2mat(answer(2)));
    mode_i_damping = str2num(cell2mat(answer(3)));
    mode_j_damping = str2num(cell2mat(answer(4)));
    
    if isempty(mode_i_frequency)
        errordlg('Frequency #1 should be a number.')
        return
    elseif length(mode_i_frequency) > 1
        errordlg('There should only be a single value for Frequency #1.')
        return
    elseif mode_i_frequency < 0
        errordlg('Frequency #1 should be non negative.')
        return
    elseif isempty(mode_j_frequency)
        errordlg('Frequency #2 should be a number.')
        return
    elseif length(mode_j_frequency) > 1
        errordlg('There should only be a single value for Frequency #2.')
        return
    elseif mode_j_frequency < 0
        errordlg('Frequency #2 should be non negative.')
        return
    elseif isempty(mode_i_damping)
        errordlg('Damping at Frequency #1 should be a number.')
        return
    elseif length(mode_i_damping) > 1
        errordlg('There should only be a single value for Damping at Frequency #1.')
        return
    elseif mode_i_damping < 0
        errordlg('Damping at Frequency #1 should be non negative.')
        return
    elseif isempty(mode_j_damping)
        errordlg('Damping at Frequency #2 should be a number.')
        return
    elseif length(mode_j_damping) > 1
        errordlg('There should only be a single value for Damping at Frequency #2.')
        return
    elseif mode_j_damping < 0
        errordlg('Damping at Frequency #2 should be non negative.')
        return
    end
    
    if handles.mode_i_frequencyTEMP ~= mode_i_frequency ||...
            handles.mode_j_frequencyTEMP ~= mode_j_frequency ||...
            handles.mode_i_dampingTEMP ~= mode_i_damping ||...
            handles.mode_j_dampingTEMP ~= mode_j_damping
        handles = OutputOff(handles);
    end
    
    handles.mode_i_frequencyTEMP = mode_i_frequency;
    handles.mode_j_frequencyTEMP = mode_j_frequency;
    handles.mode_i_dampingTEMP = mode_i_damping;
    handles.mode_j_dampingTEMP = mode_j_damping;
    
else
    
    options.Resize='on';
    options.WindowStyle='normal';
    
    prompt = {'Enter Damping (%):'};
    dlg_title = 'Enter Luco Damping Coefficient';
    num_lines = [1 50];
    def = {num2str(handles.damping)};
    answer = inputdlg(prompt,dlg_title,num_lines,def,options);
    
    if isempty(answer)
        return
    end
    
    damping = str2num(cell2mat(answer));
    
    if isempty(damping)
        errordlg('Luco Damping should be a number.')
        return
    elseif length(damping) > 1
        errordlg('There should only be one value for Luco Damping.')
        return
    elseif damping < 0
        errordlg('Luco Damping should be non negative.')
        return
    end
    
    if handles.damping ~= damping
        handles = OutputOff(handles);
    end
    
    handles.damping = damping;
    
end

guidata(hObject, handles);


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

MatNumTEMP = str2num(get(hObject, 'String'));

if isempty(MatNumTEMP) ~= 0
    errordlg('Material Number should be a number.', 'Invalid Entry for #')
    set(handles.MatNum, 'String', handles.MatNumTEMP)
    set(handles.MatNum, 'String', handles.MatNumTEMP)
    return
end

if MatNumTEMP < 1
    errordlg('Material Number should be greater than 0.', 'Invalid Entry for #')
    set(handles.MatNum, 'String', handles.MatNumTEMP)
    return
end

if rem(MatNumTEMP,1) ~= 0
    errordlg('Material Number should be an integer.', 'Invalid Entry for #')
    set(handles.MatNum, 'String', handles.MatNumTEMP)
    return
end

handles.MatNumTEMP = MatNumTEMP;
guidata(hObject,handles);


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

d = dir([pwd, '\Material\']);
str = {d.name};
str = str(3:end);
str(end+1) = {'Clear'};

[s1,v1] = listdlg('PromptString','Select a file:',...
    'SelectionMode','single',...
    'ListString',str);

if v1 ~= 0
    
    if s1 == length(str)
        File = '';
    else
        [Path, File, Ext]=fileparts(str{s1});
    end
    
    Mat = File;
    set(handles.Mat, 'String', File)
    MatNum = handles.MatNumTEMP;
    MatList = handles.MatList;
    
    MatList{MatNum,1} = MatNum;
    MatList{MatNum,2} = Mat;
    set(handles.MatTable,'Data',MatList)
    
    if MatNum < 13
        MatNum = MatNum + 1;
        set(handles.MatNum,'String',num2str(MatNum))
        handles.MatNumTEMP = MatNum;
    end
    
    handles.MatList = MatList;
    guidata(hObject,handles);
    
end


function Layers_Callback(hObject, eventdata, handles)
% hObject    handle to Layers (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 Layers as text
%        str2double(get(hObject, 'String')) returns contents of Layers as a double
str = get(hObject, 'String');
A = strread(str, '%s', 'delimiter', ' '); %#ok<REMFF1> % Layers is initially stored as a string.

for i = 1:length(A)
    
    Layers{i} = str2num(cell2mat(A(i))); %#ok<AGROW> % Changes Layers to numbers.
    
    if isempty(Layers{i}) == 1
        errordlg('Layers should be a vector of numbers. (Example: 1:3 4)',...
            'Invalid Entry for Layers.')
        set(handles.Layers, 'String', handles.LayersTEMP)
        return
    end
    
end

tot = 0;
for i = 1:length(A)
    tot = tot + length(Layers{i}); % Used to check length of layer
end

check = zeros(length(A), 1);
for i = 1:length(A)
    check(i) = max(Layers{i}); % Used to check if there is N number of elements is defined
end

n = 0;
for i = 1:length(A)
    test = Layers{i};
    for j = 1:length(test)
        n = n + 1;
        check2(n) = test(j); %#ok<AGROW> % Used to check for repeated values in Layers
    end
end

if length(check2) ~= length(unique(check2))
    errordlg('Please check Element vector for duplicates', ...
        'Duplicate Entries in Element Vector')
    return
end

if sum(mod(check2, 1)) ~= 0
    errordlg('Please check that all entries in element vector are integers.', ...
        'Noninteger Entries in Element Vector')
    return
end

if sum(str) ~= sum(handles.LayersTEMP)
    
    handles = OutputOff(handles);
    
    if handles.analysis == 1
        
        set(handles.su, 'enable', 'off')
        set(handles.su, 'string', num2str(ones(1, length(handles.VsTEMP))))
        handles.suTEMP = ones(1, length(handles.VsTEMP));
        set(handles.GammaMax, 'enable', 'off')
        set(handles.GammaMax, 'string', num2str(zeros(1, length(handles.VsTEMP))))
        handles.GammaMaxTEMP = zeros(1, length(handles.VsTEMP));
        set(handles.suMulti, 'enable', 'off')
        set(handles.suMulti, 'string', '15')
        handles.suMultiTEMP = 15;
        set(handles.GammaMaxMulti, 'enable', 'off')
        set(handles.GammaMaxMulti, 'string', '0')
        handles.GammaMaxMultiTEMP = 0;
        set(handles.NYS, 'enable', 'off')
        set(handles.NYS, 'string', num2str(1))
        handles.NYSTEMP = 1;
        
    end
    
end


handles.LayersTEMP = str;
guidata(hObject, handles);

%---Executes during object creation, after setting all properties.
function Layers_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Layers (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 Vs_Callback(hObject, eventdata, handles)
% hObject    handle to Vs (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 Vs as text
%        str2double(get(hObject, 'String')) returns contents of Vs as a double
handles.output = hObject;
Vs = str2num(get(hObject, 'String'));

if isempty(Vs) == 1
    errordlg('Shear Velocity be a number.',...
        'Invalid Entry for Shear Velocity')
    set(handles.Vs, 'String', num2str(handles.VsTEMP))
    return
end

if length(Vs) ~= length(handles.VsTEMP) || sum(Vs == handles.VsTEMP)/...
        length(Vs) ~= 1
    
    handles = OutputOff(handles);
    
end

if handles.analysis == 1
    
    set(handles.su, 'enable', 'off')
    set(handles.su, 'string', num2str(ones(1, length(Vs))))
    handles.suTEMP = ones(1, length(Vs));
    set(handles.GammaMax, 'enable', 'off')
    set(handles.GammaMax, 'string', num2str(zeros(1, length(Vs))))
    handles.GammaMaxTEMP = zeros(1, length(Vs));
    set(handles.suMulti, 'enable', 'off')
    set(handles.suMulti, 'string', '15')
    handles.suMultiTEMP = 15;
    set(handles.GammaMaxMulti, 'enable', 'off')
    set(handles.GammaMaxMulti, 'string', '0')
    handles.GammaMaxMultiTEMP = 0;
    set(handles.NYS, 'enable', 'off')
    set(handles.NYS, 'string', num2str(1))
    handles.NYSTEMP = 1;
    
end

handles.VsTEMP = Vs;
guidata(hObject, handles)


%---Executes during object creation, after setting all properties.
function Vs_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Vs (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 su_Callback(hObject, eventdata, handles)
% hObject    handle to su (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 su as text
%        str2double(get(hObject, 'String')) returns contents of su as a double
su = str2num(get(hObject, 'String'));

if isempty(su) == 1
    errordlg('Shear Strength should be a number.',...
        'Invalid Entry for Shear Strength')
    set(handles.su, 'String', num2str(handles.suTEMP))
    return
end

if length(su) ~= length(handles.suTEMP) || sum(su == handles.suTEMP)/...
        length(su) ~= 1
    
    handles = OutputOff(handles);
    
end

handles.suTEMP = su;
guidata(hObject, handles);


%---Executes during object creation, after setting all properties.
function su_CreateFcn(hObject, eventdata, handles)
% hObject    handle to su (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 GammaMax_Callback(hObject, eventdata, handles)
% hObject    handle to GammaMax (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 GammaMax as text
%        str2double(get(hObject, 'String')) returns contents of GammaMax as a double
GammaMax = str2num(get(hObject, 'String'));

if isempty(GammaMax) == 1
    errordlg('Shear Strain should be a number.',...
        'Invalid Entry for Shear Strain')
    set(handles.GammaMax, 'String', num2str(handles.GammaMaxTEMP))
    return
end

if length(GammaMax) ~= length(handles.GammaMaxTEMP) || ...
        sum(GammaMax == handles.GammaMaxTEMP)/length(GammaMax) ~= 1
    
    handles = OutputOff(handles);
    
end

handles.GammaMaxTEMP = GammaMax;
guidata(hObject, handles);


%---Executes during object creation, after setting all properties.
function GammaMax_CreateFcn(hObject, eventdata, handles)
% hObject    handle to GammaMax (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 den_Callback(hObject, eventdata, handles)
% hObject    handle to den (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 den as text
%        str2double(get(hObject, 'String')) returns contents of den as a double
den = str2num(get(hObject, 'String'));

if isempty(den) == 1
    errordlg('Density should be a number.', 'Invalid Entry for Density')
    set(handles.den, 'String', num2str(handles.denTEMP))
    return
end

if length(den) ~= length(handles.denTEMP) || ...
        sum(den == handles.denTEMP)/length(den) ~= 1
    
    handles = OutputOff(handles);
    
    if handles.analysis == 1
        
        set(handles.su, 'enable', 'off')
        set(handles.su, 'string', num2str(ones(1, length(handles.VsTEMP))))
        handles.suTEMP = ones(1, length(handles.VsTEMP));
        set(handles.GammaMax, 'enable', 'off')
        set(handles.GammaMax, 'string', num2str(zeros(1, length(handles.VsTEMP))))
        handles.GammaMaxTEMP = zeros(1, length(handles.VsTEMP));
        set(handles.suMulti, 'enable', 'off')
        set(handles.suMulti, 'string', '15')
        handles.suMultiTEMP = 15;
        set(handles.GammaMaxMulti, 'enable', 'off')
        set(handles.GammaMaxMulti, 'string', '0')
        handles.GammaMaxMultiTEMP = 0;
        set(handles.NYS, 'enable', 'off')
        set(handles.NYS, 'string', num2str(1))
        handles.NYSTEMP = 1;
        
    end
    
end


handles.denTEMP = den;
guidata(hObject, handles);

%---Executes during object creation, after setting all properties.
function den_CreateFcn(hObject, eventdata, handles)
% hObject    handle to den (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 MatProfile_Callback(hObject, eventdata, handles)
% hObject    handle to MatProfile (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 MatProfile as text
%        str2double(get(hObject,'String')) returns contents of MatProfile as a double
MatProfile = str2num(get(hObject, 'String'));

if isempty(MatProfile) == 1
    errordlg('Material should be a number.', 'Invalid Entry for Material.')
    set(handles.MatProfile, 'String', num2str(handles.MatProfileTEMP))
    return
end

[m n] = size(handles.MatList);
i = 0;
while i ~= m && isempty(cell2mat(handles.MatList(i+1,2))) ~= 1
    i = i + 1;
end

if min(MatProfile) < 0
    errordlg('Minimum Material Number should be 0. 0 will make use of the default Stress - Strain curve.',...
        'Invalid Entry for Material.')
    set(handles.MatProfile, 'String', num2str(handles.MatProfileTEMP))
    return
end

if sum(rem(MatProfile,1)) ~= 0
    errordlg('Material Number should be an integer.')
    set(handles.MatProfile, 'String', num2str(handles.MatProfileTEMP))
    return
end

if max(MatProfile) > i
    errordlg(['Material Number is not defined for ', num2str(max(MatProfile)), '.', 'Check if the material is defined for the Material desired, and that there is no empty entries inbetween.'],...
        'Invalid Entry for Material.')
    set(handles.MatProfile, 'String', num2str(handles.MatProfileTEMP))
    return
end

if length(MatProfile) ~= length(handles.MatProfileTEMP) || ...
        sum(MatProfile == handles.MatProfileTEMP)/length(MatProfile) ~= 1
    
    handles = OutputOff(handles);
    
    if handles.analysis == 1
        
        set(handles.su, 'enable', 'off')
        set(handles.su, 'string', num2str(ones(1, length(handles.VsTEMP))))
        handles.suTEMP = ones(1, length(handles.VsTEMP));
        set(handles.GammaMax, 'enable', 'off')
        set(handles.GammaMax, 'string', num2str(zeros(1, length(handles.VsTEMP))))
        handles.GammaMaxTEMP = zeros(1, length(handles.VsTEMP));
        set(handles.suMulti, 'enable', 'off')
        set(handles.suMulti, 'string', '15')
        handles.suMultiTEMP = 15;
        set(handles.GammaMaxMulti, 'enable', 'off')
        set(handles.GammaMaxMulti, 'string', '0')
        handles.GammaMaxMultiTEMP = 0;
        set(handles.NYS, 'enable', 'off')
        set(handles.NYS, 'string', num2str(1))
        handles.NYSTEMP = 1;
        
    end
    
end


handles.MatProfileTEMP = MatProfile;
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function MatProfile_CreateFcn(hObject, eventdata, handles)
% hObject    handle to MatProfile (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 LayersMulti_Callback(hObject, eventdata, handles)
% hObject    handle to LayersMulti (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 LayersMulti as text
%        str2double(get(hObject, 'String')) returns contents of LayersMulti as a double
LayersMulti = str2num(get(hObject, 'String'));

if isempty(LayersMulti) == 1
    set(handles.LayersMulti, 'string', '0')
    LayersMulti=0;
end

if LayersMulti ~= handles.LayersMultiTEMP
    
    handles = OutputOff(handles);
    
end

handles.LayersMultiTEMP = LayersMulti;
guidata(hObject, handles);


%---Executes during object creation, after setting all properties.
function LayersMulti_CreateFcn(hObject, eventdata, handles)
% hObject    handle to LayersMulti (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 VsMulti_Callback(hObject, eventdata, handles)
% hObject    handle to VsMulti (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 VsMulti as text
%        str2double(get(hObject, 'String')) returns contents of VsMulti as a double
VsMulti = str2num(get(hObject, 'String'));

if isempty(VsMulti) == 1
    set(handles.VsMulti, 'string', '0')
    VsMulti=0;
end

if VsMulti ~= handles.VsMultiTEMP
    
    handles = OutputOff(handles);
    
end

handles.VsMultiTEMP = VsMulti;
guidata(hObject, handles);


%---Executes during object creation, after setting all properties.
function VsMulti_CreateFcn(hObject, eventdata, handles)
% hObject    handle to VsMulti (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 suMulti_Callback(hObject, eventdata, handles)
% hObject    handle to suMulti (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 suMulti as text
%        str2double(get(hObject, 'String')) returns contents of suMulti as a double
suMulti = str2num(get(hObject, 'String'));

if isempty(suMulti) == 1
    set(handles.suMulti, 'string', '0')
    suMulti=0;
end

if suMulti ~= handles.suMultiTEMP
    
    handles = OutputOff(handles);
    
end

handles.suMultiTEMP = suMulti;
guidata(hObject, handles);

%---Executes during object creation, after setting all properties.
function suMulti_CreateFcn(hObject, eventdata, handles)
% hObject    handle to suMulti (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 GammaMaxMulti_Callback(hObject, eventdata, handles)
% hObject    handle to GammaMaxMulti (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 GammaMaxMulti as text
%        str2double(get(hObject, 'String')) returns contents of GammaMaxMulti as a double
GammaMaxMulti = str2num(get(hObject, 'String'));

if isempty(GammaMaxMulti) == 1
    set(handles.GammaMaxMulti, 'string', '0')
    GammaMaxMulti=0;
end

if GammaMaxMulti ~= handles.GammaMaxMultiTEMP
    
    handles = OutputOff(handles);
    
end

handles.GammaMaxMultiTEMP = GammaMaxMulti;
guidata(hObject, handles);


%---Executes during object creation, after setting all properties.
function GammaMaxMulti_CreateFcn(hObject, eventdata, handles)
% hObject handle to GammaMaxMulti (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 denMulti_Callback(hObject, eventdata, handles)
% hObject  handle to denMulti (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 denMulti as text
% str2double(get(hObject, 'String')) returns contents of denMulti as a double
denMulti = str2num(get(hObject, 'String'));

if isempty(denMulti) == 1
    set(handles.denMulti, 'string', '0')
    denMulti=0;
end

if denMulti ~= handles.denMultiTEMP
    
    handles = OutputOff(handles);
    
end

handles.denMultiTEMP = denMulti;
guidata(hObject, handles);


%---Executes during object creation, after setting all properties.
function denMulti_CreateFcn(hObject, eventdata, handles)
% hObject  handle to denMulti (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 NodeOut_Callback(hObject, eventdata, handles)
% hObject handle to NodeOut (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 NodeOut as text
% str2double(get(hObject, 'String')) returns contents of NodeOut as a double
NodeOut = str2num(get(hObject, 'String'));

if isempty(NodeOut)
    NodeOut = 0;
    set(handles.NodeOut, 'string', '0')
end

handles.NodeOutTEMP = NodeOut;
guidata(hObject, handles);


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

% Hint: get(hObject, 'Value') returns toggle state of SeparatedOut
SeparatedOut = get(hObject, 'Value');

handles.SeparatedOutTEMP = SeparatedOut;
guidata(hObject, handles);


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

% Hint: get(hObject, 'Value') returns toggle state of EarthquakeOut
EarthquakeOut = get(hObject, 'Value');

handles.EarthquakeOutTEMP = EarthquakeOut;
guidata(hObject, handles);


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

% Hint: get(hObject, 'Value') returns toggle state of DisplacementOut
DisplacementOut = get(hObject, 'Value');

handles.DisplacementOutTEMP = DisplacementOut;
guidata(hObject, handles);


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

% Hint: get(hObject, 'Value') returns toggle state of VelocityOut
VelocityOut = get(hObject, 'Value');

handles.VelocityOutTEMP = VelocityOut;
guidata(hObject, handles);


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

% Hint: get(hObject, 'Value') returns toggle state of StressOut
StressOut = get(hObject, 'Value');

handles.StressOutTEMP = StressOut;
guidata(hObject, handles);


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

% Hint: get(hObject, 'Value') returns toggle state of AccelerationOut
AccelerationOut = get(hObject, 'Value');

handles.AccelerationOutTEMP = AccelerationOut;
guidata(hObject, handles);


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

% Hint: get(hObject, 'Value') returns toggle state of StrainOut
StrainOut = get(hObject, 'Value');

handles.StrainOutTEMP = StrainOut;
guidata(hObject, handles);


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

% Hint: get(hObject, 'Value') returns toggle state of StressOut2
StressOut2 = get(hObject, 'Value');

handles.StressOut2TEMP = StressOut2;
guidata(hObject, handles);


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

if isempty(EleOut)
    EleOut = 0;
    set(handles.EleOut, 'string', '0')
end

handles.EleOutTEMP = EleOut;
guidata(hObject, handles);


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

% Hint: get(hObject,'Value') returns toggle state of SpectralOut
SpectralOut = get(hObject, 'Value');

handles.SpectralOutTEMP = SpectralOut;
guidata(hObject, handles);



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

% Hint: get(hObject,'Value') returns toggle state of FourierOut
FourierOut = get(hObject, 'Value');

handles.FourierOutTEMP = FourierOut;
guidata(hObject, handles);


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

% Hint: get(hObject,'Value') returns toggle state of TransferOut
TransferOut = get(hObject, 'Value');

handles.TransferOutTEMP = TransferOut;
guidata(hObject, handles);


%---Executes on button press in OutputFig.
function OutputFig_Callback(hObject, eventdata, handles)

handles.output = hObject;

N = handles.NTEMP;
h = handles.hTEMP;
Base = handles.BaseTEMP;
T = handles.T;
AC = handles.AC;
Disp = handles.Disp;
Vel = handles.Vel;
ACC = handles.ACC;
Strain = handles.Strain;
Stress = handles.Stress;

NodeOut = handles.NodeOutTEMP;
ElementOut = handles.EleOutTEMP;
DisplacementOut = handles.DisplacementOutTEMP;
VelocityOut = handles.VelocityOutTEMP;
AccelerationOut = handles.AccelerationOutTEMP;
StressOut = handles.StressOutTEMP;
StressOut2 = handles.StressOut2TEMP;
StrainOut = handles.StrainOutTEMP;
EarthquakeOut = handles.EarthquakeOutTEMP;
SeparatedOut = handles.SeparatedOutTEMP;

SpectralOut = handles.SpectralOutTEMP;
FourierOut = handles.FourierOutTEMP;
TransferOut = handles.TransferOutTEMP;
AccelerationXLabel = handles.AccelerationXLabelTEMP;

if Base == 1
    AC = ACC(end,:)';
end

if sum(ElementOut) == 0 && sum(StrainOut + StressOut + StressOut2) ~= 0
    errordlg('Need to enter which elements you want a figure of. Can be entered as a vector. (ex. 1:3 or 1 2 3)',...
        'Element # missing')
elseif sum(NodeOut + ElementOut) == 0 && EarthquakeOut == 0
    errordlg('Need to enter node and/or elements you want a figure of. Can be entered as a vector. (ex. 1:3 or 1 2 3)',...
        'No element or node selected.')
elseif NodeOut > N
    errordlg('Invalid Node Selection NodeOut > N.', 'Invalid Node')
elseif ElementOut > N
    errordlg('Invalid Element Selection. For Rigid Base there are only N number of elements, and for Flexible Base there are N - 1',...
        'Invalid element selection')
elseif sum(NodeOut) == 0 && sum(DisplacementOut + VelocityOut + AccelerationOut + FourierOut + TransferOut) ~= 0
    errordlg('Need to enter which nodes you want a figure of. Can be entered as a vector. (ex. 1:3 or 1 2 3)',...
        'Node # missing')
elseif sum(rem(NodeOut,1)) ~= 0
    errordlg('Node numbers should all be integers. (ex. 1:3 or 1 2 3)', 'Invalid Input for Node')
elseif sum(rem(ElementOut,1))~=0
    errordlg('Element numbers should all be integers. (ex. 1:3 or 1 2 3)', 'Invalid Input for Element')
else
    
    k1 = 0;
    k2 = 0;
    width = 0.6;
    height = 0.6;
    count = 0;
    
    if handles.unitTEMP == 1
        A = ' (m/s^2)';
        B = ' (m/s)';
        C = ' (m)';
        D = ' (kPa)';
        g = 9.81;
    else
        A = ' (in/s^2)';
        B = ' (in/s)';
        C = ' (in)';
        D = ' (ksi)';
        g = 32.2*12;
    end
    
    if handles.out == 0
        handles.out = 1;
        handles.TrueSeparated =  handles.SeparatedOutTEMP;
        num1 = DisplacementOut + VelocityOut + AccelerationOut;
        num2 = StressOut + StressOut2 + StrainOut;
        lineWidth = 1.8;
        Tn1 = 0.02;
        Tn2 = 50;
        
        co = 0;
        c1 = 0;
        if EarthquakeOut == 1
            %% Time History
            co = co + 1;
            c1 = c1 + 1;
            h1 = figure(co);
            plot(T, AC, 'r','LineWidth',lineWidth)
            title('Input Motion')
            xlabel('Time (s)'); ylabel(['Acceleration ', A]);
            grid on;
            
            set(h1,'name', 'Input Motion', 'numbertitle', 'off')
            
            % Positions Figures
            set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
            
            k2 = k2 + 0.05;
            count = count + 1;
            
            if count == 10
                k1 = k1 + 0.2;
                k2 = 0;
                count = 0;
                
                if k1 >= .8
                    k1 = 0;
                end
            end
            
            %% Response Spectrum
            co = co + 1;
            c1 = c1 + 1;
            
            ug1 = AC;
            dt = T(2)-T(1);
            ug2 = interp1(dt:dt:length(AC)*dt,AC,dt:0.001:length(AC)*dt)';
            m = 1;
            dr = .05;
            gamma = 0.5;
            beta = 0;
            Tn = 10.^(log10(Tn1):(log10(Tn2)-log10(Tn1))/1000:log10(Tn2));
            wn = 2*pi./Tn;
            
            Sa = zeros(length(wn),1);
            Sv = zeros(length(wn),1);
            Sd = zeros(length(wn),1);
            
            for j = 1:length(wn)
                
                if wn(j) > 10
                    ug = ug2;
                    dt = 0.001;
                else
                    ug = ug1;
                    dt = T(2)-T(1);
                end
                c = 2 * wn(j) * m * dr; % Damping Coefficent
                k = m * wn(j)^2; % Stiffness Coefficent
                d = zeros(1,length(ug)); % Relative Displacement Vector (m)
                v = zeros(1,length(ug)); % Relative Velocity Vector (m/s)
                a = zeros(1,length(ug)); % Relative Acceleration Vector (m/s/s)
                r1 = zeros(1,length(ug)); % Restoring Force
                
                d(1) = 0; % Initial Displacement
                v(1) = 0; % Initial Velocity
                a(1) = (-m*ug(1)-c*v(1)-k*d(1))/m; % Initial Acceleration
                
                for i=1:length(ug)-1
                    d(i+1) = d(i) + dt*v(i) + dt^2*((1/2-beta)*a(i) + beta*a(i));
                    
                    f = -m*ug(i+1); % Force from ground motion
                    r1(i+1) = k * d(i+1); % Restoring Force
                    
                    a(i+1) = (f - r1(i+1) - c*(v(i) + dt*(1-gamma)*a(i)))/(m + dt*gamma*c);
                    v(i+1) = v(i) + dt*((1-gamma)*a(i) + gamma*a(i+1));
                end
                
                at = a' + ug; % Total Acceleration
                Sa(j) = max(abs(at)); % Spectral Acceleration
                Sv(j) = max(abs(v)); % Spectral Velocity
                Sd(j) = max(abs(d)); % Spectral Displacement
            end
            
            h1 = figure(co);
            
            if AccelerationXLabel == 1
                semilogx(1./Tn,Sa(:)/g, 'r','LineWidth',lineWidth)
                xlabel('Frequency (Hz)')
            elseif AccelerationXLabel == 2
                semilogx(Tn,Sa(:)/g, 'r','LineWidth',lineWidth)
                xlabel('Period (Sec)')
                
            end
            
            xlim([0.02 50])
            ylabel('S_A (g)')
            title(['Acceleration Response Spectrum for Base ',...
                '; Depth: ', num2str(h*N), C])
            grid on;
            
            set(h1,'name', 'Input Motion', 'numbertitle', 'off')
            
            % Positions Figures
            set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
            
            k2 = k2 + 0.05;
            count = count + 1;
            
            if count == 10
                k1 = k1 + 0.2;
                k2 = 0;
                count = 0;
                
                if k1 >= .8
                    k1 = 0;
                end
            end
            
            %% Fourier Analysis
            co = co + 1;
            c1 = c1 + 1;
            a1 = AC;
            a1 = a1.*hanning(length(a1));
            
            Fs = 1/dt;
            L = length(a1);
            
            NFFT = 2^nextpow2(L); % Next power of 2 from length of y
            f = Fs/2*linspace(0,1,NFFT/2+1);
            Y11 = fft(a1,NFFT)/L;
            
            y = (2*abs(Y11(1:NFFT/2+1)));
            
            for i = 1:10
                y(2:end-1) = 0.25*y(1:end-2) + 0.5*y(2:end-1) + 0.25*y(3:end);
            end
            h1 = figure(co);
            
            if AccelerationXLabel == 1
                semilogx(f,y, 'r','LineWidth',lineWidth)
                xlabel('Frequency (Hz)')
            elseif AccelerationXLabel == 2
                semilogx(1./f,y, 'r','LineWidth',lineWidth)
                xlabel('Period (Sec)')
            end
            
            xlim([0.02 50])
            ylabel(['FAS ', B])
            title(['Fourier Amplitude Spectrum for Base: ',...
                '; Depth: ', num2str(h*N), C])
            grid on;
            
            set(h1,'name', 'Input Motion', 'numbertitle', 'off')
            
            % Positions Figures
            set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
            
            k2 = k2 + 0.05;
            count = count + 1;
            
            if count == 10
                k1 = k1 + 0.2;
                k2 = 0;
                count = 0;
                
                if k1 >= .8
                    k1 = 0;
                end
            end
            
        end
        
        if NodeOut ~= 0
            
            if SpectralOut == 1
                ug1 = AC;
                dt = T(2)-T(1);
                ug2 = interp1(dt:dt:length(AC)*dt,AC,dt:0.005:length(AC)*dt)';
                m = 1;
                dr = .05;
                gamma = 0.5;
                beta = 0,25;
                Tn = 10.^(log10(Tn1):(log10(Tn2)-log10(Tn1))/1000:log10(Tn2));
                wn = 2*pi./Tn;
                
                Sa_AC = zeros(length(wn),1);
                for j = 1:length(wn)
                    
                    if wn(j) > 10
                        ug = ug2;
                        dt = 0.005;
                    else
                        ug = ug1;
                        dt = T(2)-T(1);
                    end
                    c = 2 * wn(j) * m * dr; % Damping Coefficent
                    k = m * wn(j)^2; % Stiffness Coefficent
                    d = zeros(1,length(ug)); % Relative Displacement Vector (m)
                    v = zeros(1,length(ug)); % Relative Velocity Vector (m/s)
                    a = zeros(1,length(ug)); % Relative Acceleration Vector (m/s/s)
                    r1 = zeros(1,length(ug)); % Restoring Force
                    
                    d(1) = 0; % Initial Displacement
                    v(1) = 0; % Initial Velocity
                    a(1) = (-m*ug(1)-c*v(1)-k*d(1))/m; % Initial Acceleration
                    
                    for i=1:length(ug)-1
                        d(i+1) = d(i) + dt*v(i) + dt^2*((1/2-beta)*a(i) + beta*a(i));
                        
                        f = -m*ug(i+1); % Force from ground motion
                        r1(i+1) = k * d(i+1); % Restoring Force
                        
                        a(i+1) = (f - r1(i+1) - c*(v(i) + dt*(1-gamma)*a(i)))/(m + dt*gamma*c);
                        v(i+1) = v(i) + dt*((1-gamma)*a(i) + gamma*a(i+1));
                    end
                    
                    at = a' + ug; % Total Acceleration
                    Sa_AC(j) = max(abs(at)); % Spectral Acceleration
                end
                
            end
            
            if FourierOut == 1
                %% Fourier Analysis
                a1 = AC;
                a1 = a1.*hanning(length(a1));
                
                L = length(a1);
                NFFT = 2^nextpow2(L); % Next power of 2 from length of y
                Y11 = fft(a1,NFFT)/L;
                
                y_AC = (2*abs(Y11(1:NFFT/2+1)));
                
                for i = 1:10
                    y_AC(2:end-1) = 0.25*y_AC(1:end-2) + 0.5*y_AC(2:end-1) + 0.25*y_AC(3:end);
                end
                
            end
            
            for n = NodeOut
                j = 0;
                c1 = c1 + 1;
                if DisplacementOut == 1
                    co = co + 1;
                    if SeparatedOut == 0
                        j = j + 1;
                        h1 = figure(c1);
                        subplot(num1, 1, j); plot(T, Disp(n, :), 'k','LineWidth',lineWidth);
                        title(['Relative Displacement for Node: ',...
                            num2str(n), '; Depth: ', num2str((n - 1)*h), C])
                    else
                        h1 = figure(co);
                        plot(T, Disp(n, :), 'k','LineWidth',lineWidth);
                        title(['Relative Displacement for Node: ',...
                            num2str(n), '; Depth: ', num2str((n - 1)*h), C])
                        
                        set(h1,'name', ['Relative Displacement for Node: ',...
                            num2str(n), '; Depth: ', num2str((n - 1)*h), C]...
                            , 'numbertitle', 'off')
                        
                        % Positions Figures
                        set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
                        
                        k2 = k2 + 0.05;
                        count = count + 1;
                        
                        if count == 10
                            k1 = k1 + 0.2;
                            k2 = 0;
                            count = 0;
                            
                            if k1 >= .8
                                k1 = 0;
                            end
                        end
                        
                    end
                    
                    xlabel('Time (s)'); ylabel(['Displacement ', C]);
                    grid on;
                    
                end
                
                
                if VelocityOut == 1
                    co = co + 1;
                    if SeparatedOut == 0
                        j = j + 1;
                        h1 = figure(c1);
                        subplot(num1, 1, j); plot(T, Vel(n, :), 'k','LineWidth',lineWidth); hold on;
                        title(['Relative Velocity for Node: ', num2str(n), '; Depth: ', num2str((n - 1)*h), C])
                    else
                        h1 = figure(co);
                        plot(T, Vel(n, :), 'k','LineWidth',lineWidth); hold on;
                        title(['Relative Velocity for Node: ', num2str(n), '; Depth: ', num2str((n - 1)*h), C])
                        
                        set(h1,'name', ['Relative Velocity for Node: ',...
                            num2str(n), '; Depth: ', num2str((n - 1)*h), C],...
                            'numbertitle', 'off')
                        
                        % Positions Figures
                        set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
                        
                        k2 = k2 + 0.05;
                        count = count + 1;
                        
                        if count == 10
                            k1 = k1 + 0.2;
                            k2 = 0;
                            count = 0;
                            
                            if k1 >= .8
                                k1 = 0;
                            end
                        end
                        
                    end
                    
                    xlabel('Time (s)'); ylabel(['Velocity ', B]);
                    grid on;
                    
                end
                
                if AccelerationOut == 1
                    co = co + 1;
                    if SeparatedOut == 0
                        j = j + 1;
                        h1 = figure(c1);
                        subplot(num1, 1, j); plot(T, ACC(n, :), 'k','LineWidth',lineWidth);
                        title(['Acceleration for Node: ', num2str(n), '; Depth: ', num2str((n - 1)*h), C])
                    else
                        h1 = figure(co);
                        plot(T, ACC(n, :), 'k','LineWidth',lineWidth);
                        title(['Acceleration for Node: ', num2str(n), '; Depth: ', num2str((n - 1)*h), C])
                        
                        set(h1,'name', ['Acceleration for Node: ', num2str(n),...
                            '; Depth: ', num2str((n - 1)*h), C], 'numbertitle',...
                            'off')
                        
                        % Positions Figures
                        set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
                        
                        k2 = k2 + 0.05;
                        count = count + 1;
                        
                        if count == 10
                            k1 = k1 + 0.2;
                            k2 = 0;
                            count = 0;
                            
                            if k1 >= .8
                                k1 = 0;
                            end
                        end
                        
                    end
                    
                    xlabel('Time (s)'); ylabel(['Acceleration ', A]);
                    grid on;
                    
                end
                
                if SeparatedOut == 0 && sum(DisplacementOut + VelocityOut +...
                        AccelerationOut) ~= 0
                    
                    set(h1,'name', ['Node: ',...
                        num2str(n), '; Depth: ', num2str((n - 1)*h), C]...
                        , 'numbertitle', 'off')
                    
                    % Positions Figures
                    set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
                    
                    k2 = k2 + 0.05;
                    count = count + 1;
                    
                    if count == 10
                        k1 = k1 + 0.2;
                        k2 = 0;
                        count = 0;
                        
                        if k1 >= .8
                            k1 = 0;
                        end
                    end
                    
                end
                
                if SpectralOut == 1
                    
                    co = co + 1;
                    c1 = c1 + 1;
                    ug1 = ACC(n, :)';
                    dt = T(2)-T(1);
                    ug2 = interp1(dt:dt:length(ug1)*dt,ug1,dt:0.005:length(ug1)*dt)';
                    m = 1;
                    dr = .05;
                    gamma = 0.5;
                    beta = 0;
                    Tn = 10.^(log10(Tn1):(log10(Tn2)-log10(Tn1))/1000:log10(Tn2));
                    
                    wn = 2*pi./Tn;
                    Sa = zeros(length(wn),1);
                    for j = 1:length(wn)
                        
                        if wn(j) > 10
                            ug = ug2;
                            dt = 0.005;
                        else
                            ug = ug1;
                            dt = T(2)-T(1);
                        end
                        c = 2 * wn(j) * m * dr; % Damping Coefficent
                        k = m * wn(j)^2; % Stiffness Coefficent
                        d = zeros(1,length(ug)); % Relative Displacement Vector (m)
                        v = zeros(1,length(ug)); % Relative Velocity Vector (m/s)
                        a = zeros(1,length(ug)); % Relative Acceleration Vector (m/s/s)
                        r1 = zeros(1,length(ug)); % Restoring Force
                        
                        d(1) = 0; % Initial Displacement
                        v(1) = 0; % Initial Velocity
                        a(1) = (-m*ug(1)-c*v(1)-k*d(1))/m; % Initial Acceleration
                        
                        for i=1:length(ug)-1
                            d(i+1) = d(i) + dt*v(i) + dt^2*((1/2-beta)*a(i) + beta*a(i));
                            
                            f = -m*ug(i+1); % Force from ground motion
                            r1(i+1) = k * d(i+1); % Restoring Force
                            
                            a(i+1) = (f - r1(i+1) - c*(v(i) + dt*(1-gamma)*a(i)))/(m + dt*gamma*c);
                            v(i+1) = v(i) + dt*((1-gamma)*a(i) + gamma*a(i+1));
                        end
                        
                        at = a' + ug; % Total Acceleration
                        Sa(j) = max(abs(at)); % Spectral Acceleration
                    end
                    
                    h1 = figure(co);
                    if AccelerationXLabel == 1
                        semilogx(1./Tn,Sa_AC(:)/g, 'r','LineWidth',lineWidth)
                        hold all
                        semilogx(1./Tn,Sa(:)/g, 'k','LineWidth',lineWidth)
                        xlabel('Frequency (Hz)')
                    elseif AccelerationXLabel == 2
                        semilogx(Tn,Sa_AC(:)/g, 'r','LineWidth',lineWidth)
                        hold all
                        semilogx(Tn,Sa(:)/g, 'k','LineWidth',lineWidth)
                        xlabel('Period (Sec)')
                    end
                    
                    xlim([0.02 50])
                    ylabel('S_A (g)')
                    title(['Acceleration Response Spectrum for Node: ', num2str(n),...
                        '; Depth: ', num2str((n - 1)*h), C])
                    legend('Base',['Node: ', num2str(n)])
                    grid on;
                    
                    set(h1,'name', ['Node: ',...
                        num2str(n), '; Depth: ', num2str((n - 1)*h), C]...
                        , 'numbertitle', 'off')
                    
                    % Positions Figures
                    set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
                    
                    k2 = k2 + 0.05;
                    count = count + 1;
                    
                    if count == 10
                        k1 = k1 + 0.2;
                        k2 = 0;
                        count = 0;
                        
                        if k1 >= .8
                            k1 = 0;
                        end
                    end
                    
                    
                end
                
                if FourierOut == 1
                    
                    co = co + 1;
                    c1 = c1 + 1;
                    a1 = ACC(n, :)';
                    a1 = a1.*hanning(length(a1));
                    
                    dt = T(2)-T(1);
                    Fs = 1/dt;
                    L = length(a1);
                    
                    NFFT = 2^nextpow2(L); % Next power of 2 from length of y
                    f = Fs/2*linspace(0,1,NFFT/2+1);
                    Y11 = fft(a1,NFFT)/L;
                    
                    %%
                    y = (2*abs(Y11(1:NFFT/2+1)));
                    
                    for i = 1:10
                        y(2:end-1) = 0.25*y(1:end-2) + 0.5*y(2:end-1) + 0.25*y(3:end);
                    end
                    h1 = figure(co);
                    
                    if AccelerationXLabel == 1
                        semilogx(f,y_AC, 'r','LineWidth',lineWidth)
                        hold all
                        semilogx(f,y, 'k','LineWidth',lineWidth)
                        xlabel('Frequency (Hz)')
                    elseif AccelerationXLabel == 2
                        semilogx(1./f,y_AC, 'r','LineWidth',lineWidth)
                        hold all
                        semilogx(1./f,y, 'k','LineWidth',lineWidth)
                        xlabel('Period (Sec)')
                    end
                    
                    xlim([0.02 50])
                    ylabel(['FAS ', B])
                    title(['Fourier Amplitude Spectrum for Node: ', num2str(n),...
                        '; Depth: ', num2str((n - 1)*h), C])
                    legend('Base',['Node: ', num2str(n)])
                    grid on;
                    
                    set(h1,'name', ['Node: ',...
                        num2str(n), '; Depth: ', num2str((n - 1)*h), C]...
                        , 'numbertitle', 'off')
                    
                    % Positions Figures
                    set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
                    
                    k2 = k2 + 0.05;
                    count = count + 1;
                    
                    if count == 10
                        k1 = k1 + 0.2;
                        k2 = 0;
                        count = 0;
                        
                        if k1 >= .8
                            k1 = 0;
                        end
                    end
                    
                end
                
                if TransferOut == 1
                    
                    co = co + 1;
                    c1 = c1 + 1;
                    base = AC;
                    a1 = ACC(n, :)';
                    a1 = a1.*hanning(length(a1));
                    base = base.*hanning(length(base));
                    
                    dt = T(2)-T(1);
                    Fs = 1/dt;
                    L = length(a1);
                    
                    NFFT = 2^nextpow2(L); % Next power of 2 from length of y
                    
                    f = Fs/2*linspace(0,1,NFFT/2+1);
                    
                    Y11 = fft(a1,NFFT)/L;
                    Y12 = fft(base,NFFT)/L;
                    
                    %%
                    y = (2*abs(Y11(1:NFFT/2+1)))./(2*abs(Y12(1:NFFT/2+1)));
                    
                    for i = 1:10
                        y(2:end-1) = 0.25*y(1:end-2) + 0.5*y(2:end-1) + 0.25*y(3:end);
                    end
                    
                    h1 = figure(co);
                    if AccelerationXLabel == 1
                        semilogx(f,y, 'k','LineWidth',lineWidth)
                        xlabel('Frequency (Hz)')
                    elseif AccelerationXLabel == 2
                        semilogx(1./f,y, 'k','LineWidth',lineWidth)
                        xlabel('Period (Sec)')
                    end
                    
                    xlim([0.02 50])
                    ylabel(['Amplitude (Node ', num2str(n),'/Base)'])
                    title(['Transfer Function for Node: ', num2str(n),...
                        '; Depth: ', num2str((n - 1)*h), C])
                    grid on;
                    
                    set(h1,'name', ['Node: ',...
                        num2str(n), '; Depth: ', num2str((n - 1)*h), C]...
                        , 'numbertitle', 'off')
                    
                    % Positions Figures
                    set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
                    
                    k2 = k2 + 0.05;
                    count = count + 1;
                    
                    if count == 10
                        k1 = k1 + 0.2;
                        k2 = 0;
                        count = 0;
                        
                        if k1 >= .8
                            k1 = 0;
                        end
                    end
                    
                    
                end
                
            end
            
        end
        
        
        if ElementOut ~= 0
            for n = ElementOut
                j = 0;
                c1 = c1 + 1;
                if StrainOut == 1
                    co = co + 1;
                    if SeparatedOut == 0
                        j = j + 1;
                        h1 = figure(c1);
                        subplot(num2, 1, j); plot(T, Strain(n, :), 'k','LineWidth',lineWidth);
                        title(['Strain for Element: ', num2str(n), '; Depth: ', num2str((n - 0.5)*h), C])
                    else
                        h1 = figure(co);
                        plot(T, Strain(n, :), 'k','LineWidth',lineWidth);
                        title(['Strain for Element: ', num2str(n), '; Depth: ', num2str((n - 0.5)*h), C])
                        
                        set(h1,'name', ['Strain for Element: ', num2str(n),...
                            '; Depth: ', num2str((n - 0.5)*h), C], 'numbertitle',...
                            'off')
                        
                        % Positions Figures
                        set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
                        
                        k2 = k2 + 0.05;
                        count = count + 1;
                        
                        if count == 10
                            k1 = k1 + 0.2;
                            k2 = 0;
                            count = 0;
                            
                            if k1 >= .8
                                k1 = 0;
                            end
                        end
                        
                    end
                    xlabel('Time (s)'); ylabel('Strain (%)');
                    grid on;
                    
                end
                
                if StressOut == 1
                    co = co + 1;
                    if SeparatedOut == 0
                        j = j + 1;
                        h1 = figure(c1);
                        subplot(num2, 1, j); plot(T, Stress(n, :), 'k','LineWidth',lineWidth);
                        title(['Stress for Element: ', num2str(n),...
                            '; Depth: ', num2str((n - 0.5)*h), C])
                    else
                        h1 = figure(co);
                        plot(T, Stress(n, :), 'k','LineWidth',lineWidth);
                        title(['Stress for Element: ', num2str(n),...
                            '; Depth: ', num2str((n - 0.5)*h), C])
                        
                        set(h1,'name', ['Stress for Element: ', num2str(n),...
                            '; Depth: ', num2str((n - 0.5)*h), C],...
                            'numbertitle', 'off')
                        
                        % Positions Figures
                        set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
                        
                        k2 = k2 + 0.05;
                        count = count + 1;
                        
                        if count == 10
                            k1 = k1 + 0.2;
                            k2 = 0;
                            count = 0;
                            
                            if k1 >= .8
                                k1 = 0;
                            end
                        end
                        
                    end
                    xlabel('Time (s)'); ylabel(['Stress ', D]);
                    grid on;
                    
                end
                
                if StressOut2 == 1
                    co = co + 1;
                    if SeparatedOut == 0
                        j = j + 1;
                        h1 = figure(c1);
                        subplot(num2, 1, j); plot(Strain(n, :),...
                            Stress(n, :), 'k','LineWidth',lineWidth);
                        title(['Strain vs Stress for Element: ',...
                            num2str(n), '; Depth: ', num2str((n - 0.5)*h), C])
                    else
                        h1 = figure(co);
                        plot(Strain(n, :), Stress(n, :), 'k','LineWidth',lineWidth);
                        title(['Strain vs Stress for Element: ',...
                            num2str(n), '; Depth: ', num2str((n - 0.5)*h), C])
                        
                        set(h1,'name', ['Strain vs Stress for Element: ',...
                            num2str(n), '; Depth: ', num2str((n - 0.5)*h), C],...
                            'numbertitle', 'off')
                        
                        % Positions Figures
                        set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
                        
                        k2 = k2 + 0.05;
                        count = count + 1;
                        
                        if count == 10
                            k1 = k1 + 0.2;
                            k2 = 0;
                            count = 0;
                            
                            if k1 >= .8
                                k1 = 0;
                            end
                        end
                        
                    end
                    xlabel('Strain (%)'); ylabel(['Stress ', D]);
                    grid on;
                    
                end
                
                if SeparatedOut == 0
                    
                    set(h1,'name', ['Node: ',...
                        num2str(n), '; Depth: ', num2str((n - 1)*h), C]...
                        , 'numbertitle', 'off')
                    
                    % Positions Figures
                    set(h1, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, width, height])
                    
                    k2 = k2 + 0.05;
                    count = count + 1;
                    
                    if count == 10
                        k1 = k1 + 0.2;
                        k2 = 0;
                        count = 0;
                        
                        if k1 >= .8
                            k1 = 0;
                        end
                    end
                    
                end
                
            end
        end
        set(handles.OutputFig, 'string', 'Close Figures')
    else
        handles.out = 0;
        nfigs = get(0, 'Children');
        for i = 2:length(nfigs)
            close(nfigs(i))
        end
        set(handles.OutputFig, 'string', 'Output Figures')
    end
    
    handles.co = 0;
    handles.c1 = 0;
end
guidata(hObject, handles);


%---Executes during object deletion, before destroying properties.
function figure1_DeleteFcn(hObject, eventdata, handles)
% hObject handle to figure1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)


%---Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject handle to figure1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure
choice = questdlg('Would you like to save your inputs?', 'Save Inputs', 'Yes', 'No', 'No');

if isempty(choice) == 1
else
    switch choice
        case 'Yes'
            
            FileName = handles.FileNameTEMP;
            
            EarthquakeFile = handles.EarthquakeFile;
            LocationControl_val = handles.LocationControl_valTEMP;
            g = handles.gTEMP;
            unit = handles.unitTEMP;
            dt = handles.dtTEMP;
            Scale = handles.ScaleTEMP;
            
            analysis = handles.analysis;
            NYS = handles.NYSTEMP;
            Tol = handles.TolTEMP;
            
            Base = handles.BaseTEMP;
            Incident = handles.IncidentTEMP;
            
            h = handles.hTEMP;
            N = handles.NTEMP;
            pr = handles.prTEMP;
            Vsr = handles.VsrTEMP;
            
            Layers = handles.LayersTEMP;
            Vs = handles.VsTEMP;
            su = handles.suTEMP;
            GammaMax = handles.GammaMaxTEMP;
            den = handles.denTEMP;
            LayersMulti = handles.LayersMultiTEMP;
            VsMulti = handles.VsMultiTEMP;
            suMulti = handles.suMultiTEMP;
            GammaMaxMulti = handles.GammaMaxMultiTEMP;
            denMulti = handles.denMultiTEMP;
            mode_i_frequency = handles.mode_i_frequencyTEMP;
            mode_j_frequency = handles.mode_j_frequencyTEMP;
            mode_i_damping = handles.mode_i_dampingTEMP;
            mode_j_damping = handles.mode_j_dampingTEMP;
            
            choice = 0;
            
            a = cd; %#ok<*MCCD>
            [SaveFile FileLocation] = uiputfile('.mat', 'Save File',...
                [a, '\Saved Inputs']);
            save(strcat(FileLocation, SaveFile))
        case 'No'
    end
    close all force
end


% --------------------------------------------------------------------
function File_Callback(hObject, eventdata, handles)
% hObject    handle to File (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function New_Callback(hObject, eventdata, handles)
% hObject    handle to New (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

load([handles.GUILOCATION, '\Saved Inputs\Default.mat'])

handles.FileNameTEMP = 'Untitled';
set(handles.FileSaveName, 'string', handles.FileNameTEMP)

handles.K = 0;
handles.M = 0;
handles.T = 0;
handles.Disp = 0;
handles.Vel = 0;
handles.ACC = 0;
handles.Strain = 0;
handles.Stress = 0;

handles.LocationControl_valTEMP = LocationControl_val;

handles.EarthquakeFile = EarthquakeFile;
set(handles.Location, 'string', handles.EarthquakeFile)

handles.EQ = EQ;
handles.TTEMP = TTEMP;
handles.gTEMP = g;
if g == 1
    set(handles.g, 'value', 0)
else
    set(handles.g, 'value', 1)
end

handles.unitTEMP = unit;
if unit == 1
    set(handles.meter, 'value', 1)
    set(handles.hstr, 'string', 'Height of ele (m)')
    set(handles.RockDenstr, 'string', 'Density of Rock (kg/m^3)')
    set(handles.Vsrstr, 'string', 'Vs of Rock (m/s)')
    set(handles.Vsstr, 'string', 'Shear Wave Velocity (m/s)')
    set(handles.sustr, 'string', 'Shear Strength (kPa)')
    set(handles.GammaMaxstr, 'string', 'Max Shear Strain (%)')
    set(handles.denstr, 'string', 'Mass Density (kg/m^3)')
else
    set(handles.foot, 'value', 1)
    set(handles.hstr, 'string', 'Height of ele (in)')
    set(handles.RockDenstr, 'string', 'Density of Rock (slug/in^3)')
    set(handles.Vsrstr, 'string', 'Vs of Rock (in/s)')
    set(handles.Vsstr, 'string', 'Shear Wave Velocity (in/s)')
    set(handles.sustr, 'string', 'Shear Strength (ksi)')
    set(handles.GammaMaxstr, 'string', 'Max Shear Strain (%)')
    set(handles.denstr, 'string', 'Mass Density (slug/in^3)')
end

handles.dtTEMP = dt;
set(handles.dt, 'string', num2str(handles.dtTEMP))

handles.ScaleTEMP = Scale;
set(handles.Scale, 'string', num2str(handles.ScaleTEMP))

handles.analysis = analysis;
if analysis == 1
    
    set(handles.Linear, 'value', 1)
    set(handles.su, 'enable', 'off')
    set(handles.GammaMax, 'enable', 'off')
    set(handles.MatProfile, 'enable', 'off')
    set(handles.suMulti, 'enable', 'off')
    set(handles.GammaMaxMulti, 'enable', 'off')
    set(handles.NYS, 'enable', 'off')
    
else
    
    set(handles.Nonlinear, 'value', 1)
    set(handles.su, 'enable', 'on')
    set(handles.GammaMax, 'enable', 'on')
    set(handles.MatProfile, 'enable', 'on')
    set(handles.suMulti, 'enable', 'on')
    set(handles.GammaMaxMulti, 'enable', 'on')
    set(handles.NYS, 'enable', 'on')
    
end

handles.NYSTEMP = NYS;
set(handles.NYS, 'string', num2str(handles.NYSTEMP))

handles.TolTEMP = Tol;
set(handles.Tol, 'string', num2str(handles.TolTEMP))

handles.BaseTEMP = Base;
if Base == 0
    
    set(handles.RigidBase, 'value', 1)
    set(handles.pr, 'enable', 'off')
    set(handles.Vsr, 'enable', 'off')
    set(handles.pr, 'string', '')
    set(handles.Vsr, 'string', '')
    set(handles.Incident, 'enable', 'off')
    
    if DampType == 1
        set(handles.RayleighDamping, 'Value', 1)
    else
        set(handles.LucoDamping, 'Value', 0)
    end
    
else
    
    set(handles.FlexibleBase, 'value', 1)
    set(handles.pr, 'enable', 'on')
    set(handles.Vsr, 'enable', 'on')
    set(handles.pr, 'string', num2str(pr))
    set(handles.Vsr, 'string', num2str(Vsr))
    set(handles.Incident, 'enable', 'on')
    set(handles.LucoDamping, 'enable', 'off')
    
end

handles.IncidentTEMP =  Incident;
set(handles.Incident, 'value', handles.IncidentTEMP)

handles.hTEMP = h;
set(handles.h, 'string', num2str(handles.hTEMP))

handles.NTEMP = N;
set(handles.N, 'string', num2str(handles.NTEMP))

handles.prTEMP = pr;
handles.VsrTEMP = Vsr;

handles.mode_i_frequencyTEMP = mode_i_frequency;
handles.mode_j_frequencyTEMP = mode_j_frequency;
handles.mode_i_dampingTEMP = mode_i_damping;
handles.mode_j_dampingTEMP = mode_j_damping;
handles.damping = damping;
handles.DampType = DampType;

handles.MatNumTEMP = 1;
set(handles.MatNum, 'String', handles.MatNumTEMP)

handles.MatTEMP = MatList(1,2);
set(handles.Mat, 'String', handles.MatTEMP)

handles.MatList = MatList;
set(handles.MatTable,'Data',MatList)

handles.LayersTEMP = Layers;
set(handles.Layers, 'string', num2str(handles.LayersTEMP))

handles.VsTEMP = Vs;
set(handles.Vs, 'string', num2str(handles.VsTEMP))

handles.suTEMP = su;
set(handles.su, 'string', num2str(handles.suTEMP))

handles.GammaMaxTEMP = GammaMax;
set(handles.GammaMax, 'string', num2str(GammaMax))

handles.denTEMP = den;
set(handles.den, 'string', num2str(den))

handles.MatProfileTEMP = MatProfile;
set(handles.MatProfile, 'string', num2str(handles.MatProfileTEMP))

handles.LayersMultiTEMP = LayersMulti;
set(handles.LayersMulti, 'string', num2str(LayersMulti))

handles.VsMultiTEMP = VsMulti;
set(handles.VsMulti, 'string', num2str(VsMulti))

handles.suMultiTEMP = suMulti;
set(handles.suMulti, 'string', num2str(handles.suMultiTEMP))

handles.GammaMaxMultiTEMP = GammaMaxMulti;
set(handles.GammaMaxMulti, 'string', num2str(GammaMaxMulti))

handles.denMultiTEMP = denMulti;
set(handles.denMulti, 'string', num2str(denMulti))

handles = OutputOff(handles);

plot(zeros(N, 1), 0:h:(N - 1)*h, 'o')
hold on
plot(zeros(N + 1, 1), 0:h:N*h, ' - ')
axis([ - 0.05 0.05 0 N*h]);
set(gca, 'YDir', 'reverse')
set(gca, 'xtick', []);

if handles.unitTEMP == 1
    ylabel('depth (m)')
else
    ylabel('depth (in.)')
end

if Base == 1
    plot(0, (N*h), 'o')
end
hold off

guidata(hObject, handles);


% --------------------------------------------------------------------
function Open_Callback(hObject, eventdata, handles)
% hObject    handle to Open (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[FileNameTEMP b] = uigetfile('.mat','Load File', [handles.GUILOCATION,'\Saved Inputs']);

% If FileNameTEMP ==0  then user did not select anything and just closed the
% Standard open file dialog box.
if FileNameTEMP == 0
else
    
    handles.FileNameTEMP = FileNameTEMP;
    c = strrep(handles.FileNameTEMP, '.mat', '');
    set(handles.FileSaveName, 'string', c)
    load([b handles.FileNameTEMP])
    
    handles.K = 0;
    handles.M = 0;
    handles.T = 0;
    handles.Disp = 0;
    handles.Vel = 0;
    handles.ACC = 0;
    handles.Strain = 0;
    handles.Stress = 0;
    
    handles.LocationControl_valTEMP = LocationControl_val;
    
    handles.EarthquakeFile = EarthquakeFile;
    set(handles.Location, 'string', handles.EarthquakeFile)
    
    handles.EQ = EQ;
    handles.TTEMP = TTEMP;
    handles.gTEMP = g;
    if g == 1
        set(handles.g, 'value', 0)
    else
        set(handles.g, 'value', 1)
    end
    
    handles.unitTEMP = unit;
    if unit == 1
        set(handles.meter, 'value', 1)
        set(handles.hstr, 'string', 'Height of ele (m)')
        set(handles.RockDenstr, 'string', 'Density of Rock (kg/m^3)')
        set(handles.Vsrstr, 'string', 'Vs of Rock (m/s)')
        set(handles.Vsstr, 'string', 'Shear Wave Velocity (m/s)')
        set(handles.sustr, 'string', 'Shear Strength (kPa)')
        set(handles.GammaMaxstr, 'string', 'Max Shear Strain (%)')
        set(handles.denstr, 'string', 'Mass Density (kg/m^3)')
    else
        set(handles.foot, 'value', 1)
        set(handles.hstr, 'string', 'Height of ele (in)')
        set(handles.RockDenstr, 'string', 'Density of Rock (slug/in^3)')
        set(handles.Vsrstr, 'string', 'Vs of Rock (in/s)')
        set(handles.Vsstr, 'string', 'Shear Wave Velocity (in/s)')
        set(handles.sustr, 'string', 'Shear Strength (ksi)')
        set(handles.GammaMaxstr, 'string', 'Max Shear Strain (%)')
        set(handles.denstr, 'string', 'Mass Density (slug/in^3)')
    end
    
    handles.dtTEMP = dt;
    set(handles.dt, 'string', num2str(handles.dtTEMP))
    
    handles.ScaleTEMP = Scale;
    set(handles.Scale, 'string', num2str(handles.ScaleTEMP))
    
    handles.analysis = analysis;
    if analysis == 1
        
        set(handles.Linear, 'value', 1)
        set(handles.su, 'enable', 'off')
        set(handles.GammaMax, 'enable', 'off')
        set(handles.MatProfile, 'enable', 'off')
        set(handles.suMulti, 'enable', 'off')
        set(handles.GammaMaxMulti, 'enable', 'off')
        set(handles.NYS, 'enable', 'off')
        
    else
        
        set(handles.Nonlinear, 'value', 1)
        set(handles.su, 'enable', 'on')
        set(handles.GammaMax, 'enable', 'on')
        set(handles.MatProfile, 'enable', 'on')
        set(handles.suMulti, 'enable', 'on')
        set(handles.GammaMaxMulti, 'enable', 'on')
        set(handles.NYS, 'enable', 'on')
        
    end
    
    handles.NYSTEMP = NYS;
    set(handles.NYS, 'string', num2str(handles.NYSTEMP))
    
    handles.TolTEMP = Tol;
    set(handles.Tol, 'string', num2str(handles.TolTEMP))
    
    handles.BaseTEMP = Base;
    if Base == 0
        
        set(handles.RigidBase, 'value', 1)
        set(handles.pr, 'enable', 'off')
        set(handles.Vsr, 'enable', 'off')
        set(handles.pr, 'string', '')
        set(handles.Vsr, 'string', '')
        set(handles.Incident, 'enable', 'off')
        
        set(handles.LucoDamping, 'enable', 'on')
        handles.DampType = DampType;
        
        if DampType == 1
            set(handles.RayleighDamping, 'Value', 1)
        else
            set(handles.LucoDamping, 'Value', 1)
        end
        
    else
        
        set(handles.FlexibleBase, 'value', 1)
        set(handles.pr, 'enable', 'on')
        set(handles.Vsr, 'enable', 'on')
        set(handles.pr, 'string', num2str(pr))
        set(handles.Vsr, 'string', num2str(Vsr))
        set(handles.Incident, 'enable', 'on')
        
        set(handles.RayleighDamping, 'Value', 1)
        set(handles.LucoDamping, 'enable', 'off')
        
        handles.DampType = 1;
        
    end
    
    handles.IncidentTEMP =  Incident;
    set(handles.Incident, 'value', handles.IncidentTEMP)
    
    handles.hTEMP = h;
    set(handles.h, 'string', num2str(handles.hTEMP))
    
    handles.NTEMP = N;
    set(handles.N, 'string', num2str(handles.NTEMP))
    
    handles.prTEMP = pr;
    handles.VsrTEMP = Vsr;
    
    handles.mode_i_frequencyTEMP = mode_i_frequency;
    handles.mode_j_frequencyTEMP = mode_j_frequency;
    handles.mode_i_dampingTEMP = mode_i_damping;
    handles.mode_j_dampingTEMP = mode_j_damping;
    handles.damping = damping;
    
    handles.MatNumTEMP = 1;
    set(handles.MatNum, 'String', handles.MatNumTEMP)
    
    handles.MatTEMP = MatList(1,2);
    set(handles.Mat, 'String', handles.MatTEMP)
    
    handles.MatList = MatList;
    set(handles.MatTable,'Data',MatList)
    
    handles.LayersTEMP = Layers;
    set(handles.Layers, 'string', num2str(handles.LayersTEMP))
    
    handles.VsTEMP = Vs;
    set(handles.Vs, 'string', num2str(handles.VsTEMP))
    
    handles.suTEMP = su;
    set(handles.su, 'string', num2str(handles.suTEMP))
    
    handles.GammaMaxTEMP = GammaMax;
    set(handles.GammaMax, 'string', num2str(GammaMax))
    
    handles.denTEMP = den;
    set(handles.den, 'string', num2str(den))
    
    handles.MatProfileTEMP = MatProfile;
    set(handles.MatProfile, 'string', num2str(handles.MatProfileTEMP))
    
    handles.LayersMultiTEMP = LayersMulti;
    set(handles.LayersMulti, 'string', num2str(LayersMulti))
    
    handles.VsMultiTEMP = VsMulti;
    set(handles.VsMulti, 'string', num2str(VsMulti))
    
    handles.suMultiTEMP = suMulti;
    set(handles.suMulti, 'string', num2str(handles.suMultiTEMP))
    
    handles.GammaMaxMultiTEMP = GammaMaxMulti;
    set(handles.GammaMaxMulti, 'string', num2str(GammaMaxMulti))
    
    handles.denMultiTEMP = denMulti;
    set(handles.denMulti, 'string', num2str(denMulti))
    
    handles = OutputOff(handles);
    
    plot(zeros(N, 1), 0:h:(N - 1)*h, 'o')
    hold on
    plot(zeros(N + 1, 1), 0:h:N*h, ' - ')
    axis([ - 0.05 0.05 0 N*h]);
    set(gca, 'YDir', 'reverse')
    set(gca, 'xtick', []);
    
    if handles.unitTEMP == 1
        ylabel('depth (m)')
    else
        ylabel('depth (in.)')
    end
    
    if Base == 1
        plot(0, (N*h), 'o')
    end
    hold off
end
guidata(hObject, handles);


% --------------------------------------------------------------------
function Save_Callback(hObject, eventdata, handles)
% hObject    handle to Save (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

FileName = strrep(handles.FileNameTEMP, '.mat', '');

% Check if user is trying to save over Default
while exist([handles.GUILOCATION, '\Saved Inputs\', FileName, '.mat'])
    
    % Request user to change name if trying to save as Default
    if strcmpi(FileName, 'Default') == 1
        
        prompt = {'Enter a Different Name:'};
        dlg_title = 'Error cannot modify Default file';
        num_lines = [1, 50];
        def = {strrep(handles.FileNameTEMP, '.mat', '')};
        options.Resize='on';
        options.WindowStyle='normal';
        answer = inputdlg(prompt,dlg_title,num_lines,def,options);
        
        if isempty(answer)
            return
        end
        
        FileName = cell2mat(answer);
        
    else
        
        choice = questdlg('A file with that name already exist are you sure you want to save?', 'Save Inputs',...
            'Yes', 'No', 'Change Name', 'No');
        
        switch choice
            case 'Yes'
                break
            case 'No'
                return
            case 'Change Name'
                
                % Prompt user to change name
                prompt = {'Enter a Different File Name:'};
                dlg_title = 'Enter a Different File Name';
                num_lines = [1, 50];
                
                def = {FileName};
                
                options.Resize='on';
                options.WindowStyle='normal';
                
                answer = inputdlg(prompt,dlg_title,num_lines,def,options);
                
                if isempty(answer)
                    return
                end
                
                FileName = cell2mat(answer(1));
                
        end
        
    end
end

handles.FileNameTEMP = FileName;
guidata(hObject, handles);

EarthquakeFile = handles.EarthquakeFile;
TTEMP = handles.TTEMP;
EQ = handles.EQ;
LocationControl_val = handles.LocationControl_valTEMP;
g = handles.gTEMP;
unit = handles.unitTEMP;
dt = handles.dtTEMP;
Scale = handles.ScaleTEMP;

analysis = handles.analysis;
NYS = handles.NYSTEMP;
Tol = handles.TolTEMP;

Base = handles.BaseTEMP;
Incident = handles.IncidentTEMP;

h = handles.hTEMP;
N = handles.NTEMP;
pr = handles.prTEMP;
Vsr = handles.VsrTEMP;

MatList = handles.MatList;

Layers = handles.LayersTEMP;
Vs = handles.VsTEMP;
su = handles.suTEMP;
GammaMax = handles.GammaMaxTEMP;
den = handles.denTEMP;
MatProfile = handles.MatProfileTEMP;
LayersMulti = handles.LayersMultiTEMP;
VsMulti = handles.VsMultiTEMP;
suMulti = handles.suMultiTEMP;
GammaMaxMulti = handles.GammaMaxMultiTEMP;
denMulti = handles.denMultiTEMP;
mode_i_frequency = handles.mode_i_frequencyTEMP;
mode_j_frequency = handles.mode_j_frequencyTEMP;
mode_i_damping = handles.mode_i_dampingTEMP;
mode_j_damping = handles.mode_j_dampingTEMP;
damping = handles.damping;
DampType = handles.DampType;

set(handles.FileSaveName, 'string', FileName)
GUILOCATION = handles.GUILOCATION;
clear hObject handles eventdata
save([GUILOCATION, '\Saved Inputs\', FileName])

msgbox('Complete')


% --------------------------------------------------------------------
function Material_Callback(hObject, eventdata, handles)
% hObject    handle to Material (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function DNM_Callback(hObject, eventdata, handles)
% hObject    handle to DNM (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

FileName = '';
Strain = '';
Modulus = '';
Strain2 = '';
Damping = '';
Selection = 2;

Type = 1;
plot1(handles.GUILOCATION, FileName, Strain, Modulus, Strain2, Damping, Selection, Type)


% --------------------------------------------------------------------
function MEM_Callback(hObject, eventdata, handles)
% hObject    handle to MEM (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[FileName, Strain, Modulus, Strain2, Damping, Selection] = MatRead(handles.GUILOCATION);

Type = 2;
plot1(handles.GUILOCATION, FileName, Strain, Modulus, Strain2, Damping, Selection, Type)


% --------------------------------------------------------------------
function USSDM_Callback(hObject, eventdata, handles)
% hObject    handle to USSDM (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[FileName, Strain, Modulus, Strain2, Damping, Selection] = MatRead(handles.GUILOCATION);

Type = 3;
plot1(handles.GUILOCATION, FileName, Strain, Modulus, Strain2, Damping, Selection, Type)


% --------------------------------------------------------------------
function UNMM_Callback(hObject, eventdata, handles)
% hObject    handle to UNMM (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[FileName, Strain, Modulus, Strain2, Damping, Selection] = MatRead(handles.GUILOCATION);

Type = 4;
plot1(handles.GUILOCATION, FileName, Strain, Modulus, Strain2, Damping, Selection, Type)


% --------------------------------------------------------------------
function VEM_Callback(hObject, eventdata, handles)
% hObject    handle to VEM (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[FileName, Strain, Modulus, Strain2, Damping, Selection] = MatRead(handles.GUILOCATION);
Shear = Modulus.*Strain/100;

% Plots values of selected material
f = figure;
set(f, 'units', 'normalized', 'outerposition', [0, 0.4, 0.5, 0.6])
set(f,'name','Strain vs. Modulus','numbertitle','off')
semilogx([0 Strain], [1 Modulus])
xlabel('Strain (%)')
ylabel('Modulus (G/Gmax)')
title('Strain (%) vs. Modulus (G/Gmax)')

f = figure;
set(f, 'units', 'normalized', 'outerposition', [0.5, 0.4, 0.5, 0.6])
set(f,'name','Strain vs. Shear','numbertitle','off')
plot([0 Strain], [0 Shear])
xlabel('Strain (%)')
ylabel('Stress (/Gmax)')
title('Strain (%) vs. Stress (/Gmax)')

% Create figure listing values from selected material file
f = figure;
size = get(0,'ScreenSize');
set(f, 'outerposition',[(size(3)-700)/2, 0, 700, 360]);
set(f, 'menubar', 'none');
set(f, 'name', 'View Material')
set(f, 'NumberTitle', 'off')

uicontrol('Style', 'text',...
    'String', num2str(Shear),...
    'Position', [10 20 680 30])

uicontrol('Style', 'text',...
    'String', 'Shear Stress Values (/Gmax)',...
    'Position', [10 60 680 20])

uicontrol('Style', 'text',...
    'String', num2str(Modulus),...
    'Position', [10 100 680 30])

uicontrol('Style', 'text',...
    'String', 'Modulus Reduction Values (G/Gmax)',...
    'Position', [10 140 680 20])

uicontrol('Style','text',...
    'String', num2str(Strain),...
    'Position', [10 180 680 30])

uicontrol('Style', 'text',...
    'String', 'Strain Values (%)',...
    'Position', [10 220 680 20])

uicontrol('Style','text',...
    'String', strrep(FileName, '.txt', ''),...
    'Position', [10 260 680 20])

uicontrol('Style', 'text',...
    'String', 'File Name',...
    'Position', [10 290 680 20])


% --------------------------------------------------------------------
function ModeShape_Callback(hObject, eventdata, handles)
% hObject    handle to ModeShape (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

h = handles.hTEMP;
Base = handles.BaseTEMP; % (Rigid Base = 0; Flexible Base = 1)
N = handles.NTEMP; % Number of Elements
pr = handles.prTEMP; % Density of Rock
Vsr = handles.VsrTEMP; % Shear Wave Velocity of Rock
str = handles.LayersTEMP; % Assigns Layer (used to define different Vs,
% den, etc.)

A = strread(str, '%s', 'delimiter', ' '); %#ok<REMFF1> % Layers is initially stored as a string.

if Base == 1
    errordlg('Mode Shapes are currently only available for rigid base.')
    return
end

for i = 1:length(A)
    Layers{i} = str2num(cell2mat(A(i))); %#ok<AGROW> % Changes Layers to numbers.
end

tot = 0;
for i = 1:length(A)
    tot = tot + length(Layers{i}); % Used to check length of layer
end

check = zeros(length(A), 1);
for i = 1:length(A)
    check(i) = max(Layers{i}); % Used to check if there is N number of elements is defined
end

n = 0;
for i = 1:length(A)
    test = Layers{i};
    for j = 1:length(test)
        n = n + 1;
        check2(n) = test(j); %#ok<AGROW> % Used to check for repeated values in Layers
    end
end

if length(check2) ~= length(unique(check2))
    errordlg('Please check Element vector for duplicates', ...
        'Duplicate Entries in Element Vector')
    return
elseif max(check) ~= N
    errordlg('Please change # of Ele/Node to match Elements or vice versa.', ...
        '# of Ele/Nodes do not match Elements.');
    return
elseif tot ~= N
    errordlg('Please check Element vector some values are missing. (ex. 1:3 5:6 Number 4 is missing).', ...
        'Missing Values in Element Vector')
    return
elseif sum(mod(check2, 1)) ~= 0
    errordlg('Please check that all entries in element vector are integers.', ...
        'Noninteger Entries in Element Vector')
    return
elseif length(Layers) ~= length(handles.VsTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Velocity should have 3 inputs.',...
        'Please check Velocity')
    return
elseif length(handles.VsTEMP) ~= sum(handles.VsTEMP>0)
    errordlg('Enter Velocity as positive values.', 'Error in Velocity')
    return
elseif length(Layers) ~= length(handles.denTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Density should have 3 inputs.', 'Please check Density')
    return
elseif length(handles.denTEMP) ~= sum(handles.denTEMP>0)
    errordlg('Density should all be positive.', 'Error in density')
    return
end

% Obtain Mode Number
answer = N + 1;
while answer > N
    
    prompt = {['Enter desired mode number between 1 and ', num2str(N), ':']};
    dlg_title = 'Mode Shape';
    num_lines = [1, 50];
    def = {''};
    options.Resize='on';
    options.WindowStyle='normal';
    answer = inputdlg(prompt,dlg_title,num_lines,def,options);
    
    if isempty(answer)
        return
    elseif isempty(str2num(cell2mat(answer))) || str2num(cell2mat(answer)) < 1 || str2num(cell2mat(answer)) > N || rem(str2num(cell2mat(answer)),1) ~= 0
        f = errordlg(['Invalid input, please enter an integer between 1 and ', num2str(N)]);
        uiwait(f);
        answer = N + 1;
    else
        answer = str2num(cell2mat(answer));
    end
    
end

ModeNum = answer;

h = handles.hTEMP*ones(N, 1); % Creates Vector of height (for each element)
Vs = handles.VsTEMP*10^handles.VsMultiTEMP; % Vector of Shear Wave Velocity
den = handles.denTEMP*10^handles.denMultiTEMP; % Vector Density

% Reorganizes Vs, and den according to how layers were
% defined

Vstemp = zeros(N-Base, 1);
dentemp = zeros(N-Base, 1);
for i = 1:length(Layers)
    Vstemp(Layers{i}) = Vs(i);
    dentemp(Layers{i}) = den(i);
end

Vs = Vstemp;
den = dentemp;

G = den.*Vs.^2; %Soil Shear Modulus Vector (Pa)
Rho = den; %mass Density Vector

M = zeros(N,N);
M(sub2ind([N, N], (1:N-Base), (1:N-Base))) = Rho(1:N-Base).*h(1:N-Base);
S(1:N-Base) = G(1:N-Base)./h(1:N-Base);

M(1, 1) = M(1, 1)/2;	%this adjust mass of first degree of freedom
%at the shear beam top (at ground surface) to half value as it should be

%Define Initial Tangent Stiffness for each spring
%S = G/h

%Assemble Global Tangent Stiffness Matrix KT (Subroutine 2)
K = Sub2(S, N);

[V, D] = eig((inv(M)*K)^.5);
D(~D) = nan;

for i = 1:answer-1
    [a b] = find(D == min(min(D)));
    D(a, b) = inf;
end

[a b] = find(D == min(min(D)));
x = V(:, a)';

x(end + 1) = 0;

if x(1)<0
    x = -x;
end

clear str

str(1) = {['Mode ', num2str(ModeNum)]};
str(2) = {['f = ', num2str(D(a, b)/2/pi), ' Hz.']};

plot(x, 0:h:N*h)

dif = abs(max(x) - min(x))*0.05;

Llim = min(x) - dif;
Rlim = max(x) + dif;

axis([Llim Rlim 0 N*h(1)]);

text(Llim + dif, 0.05*(N)*h(1), str)

set(gca, 'YDir', 'reverse')

if handles.unitTEMP == 1
    ylabel('depth (m)')
else
    ylabel('depth (in.)')
end


% --------------------------------------------------------------------
function Analysis_Callback(hObject, eventdata, handles)
% hObject    handle to Analysis (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function Run_Callback(hObject, eventdata, handles)
% hObject    handle to Run (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%% Error checking
disp('Error Checking')

Base = handles.BaseTEMP; % (Rigid Base = 0; Flexible Base = 1)
N = handles.NTEMP; % Number of Elements
pr = handles.prTEMP; % Density of Rock
Vsr = handles.VsrTEMP; % Shear Wave Velocity of Rock

str = handles.LayersTEMP; % Assigns Layer (used to define different Vs,
% den, etc.)

A = strread(str, '%s', 'delimiter', ' '); %#ok<REMFF1> % Layers is initially stored as a string.

for i = 1:length(A)
    Layers{i} = str2num(cell2mat(A(i))); %#ok<AGROW> % Changes Layers to numbers.
end

tot = 0;
for i = 1:length(A)
    tot = tot + length(Layers{i}); % Used to check length of layer
end

check = zeros(length(A), 1);
for i = 1:length(A)
    check(i) = max(Layers{i}); % Used to check if there is N number of elements is defined
end

n = 0;
for i = 1:length(A)
    test = Layers{i};
    for j = 1:length(test)
        n = n + 1;
        check2(n) = test(j); %#ok<AGROW> % Used to check for repeated values in Layers
    end
end

if max(check) ~= N
    errordlg('Please change # of Ele/Node to match Elements or vice versa.', ...
        '# of Ele/Nodes do not match Elements.');
elseif tot ~= N
    errordlg('Please check Element vector some values are missing. (ex. 1:3 5:6 Number 4 is missing).', ...
        'Missing Values in Element Vector')
elseif length(check2) ~= length(unique(check2))
    errordlg('Please check Element vector for duplicates', ...
        'Duplicate Entries in Element Vector')
elseif sum(mod(check2, 1)) ~= 0
    errordlg('Please check that all entries in element vector are integers.', ...
        'Noninteger Entries in Element Vector')
elseif pr == 0&&Base == 1
    errordlg('Density of Rock should be a positive value', ...
        'Check Density of Rock')
elseif Vsr == 0&&Base == 1
    errordlg('Velocity of Rock should be a positive value', ...
        'Check Velocity of Rock')
elseif length(Layers) ~= length(handles.VsTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Velocity should have 3 inputs.',...
        'Please check Velocity')
elseif length(handles.VsTEMP) ~= sum(handles.VsTEMP>0)
    errordlg('Enter Velocity as positive values.', 'Error in Velocity')
elseif length(Layers) ~= length(handles.suTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Shear Strength should have 3 inputs.',...
        'Please check Shear Strength')
elseif length(Layers) ~= length(handles.GammaMaxTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Max Shear Strain should have 3 inputs.',...
        'Please check Max Shear Strain')
elseif length(Layers) ~= length(handles.denTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Density should have 3 inputs.', 'Please check Density')
elseif length(handles.denTEMP) ~= sum(handles.denTEMP>0)
    errordlg('Density should all be positive.', 'Error in density')
elseif length(Layers) ~= length(handles.MatProfileTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Material should have 3 inputs. Can enter 0 to use default Stress Strain logic.',...
        'Please check Material')
elseif handles.EQ == 0
    errordlg('Please select an earthquake.')
else
    
    %% Defining Variables
    disp('Defining Variables')
    
    EQ = handles.EQ; % EQ acceleration Data
    dt = handles.dtTEMP; % Time Step
    g = handles.gTEMP; % gravity
    Scale = handles.ScaleTEMP; % Used to scale the EQ acceleration Data
    Incident = handles.IncidentTEMP; % If EQ data was not incident will divide by 2
    
    Tol = handles.TolTEMP; % Tolerance
    NYS = handles.NYSTEMP; % Number of Yield Surfaces
    
    N = N + Base; % Increases N by 1 if flexible otherwise N stays the same
    h = handles.hTEMP*ones(N, 1); % Creates Vector of height (for each element)
    
    Vs = handles.VsTEMP*10^handles.VsMultiTEMP; % Vector of Shear Wave Velocity
    
    den = handles.denTEMP*10^handles.denMultiTEMP; % Vector Density
    su = handles.suTEMP*10^handles.suMultiTEMP*1000; % Vector of Shear Strength Converts kPa/ksi to Pa/psi to keep units consistent
    GammaMax = handles.GammaMaxTEMP*10^handles.GammaMaxMultiTEMP/100; % Vector of Max Shear Strain Converts (%) to decimal
    mode_i_frequency = handles.mode_i_frequencyTEMP; % Frequency 1 (for damping)
    mode_j_frequency = handles.mode_j_frequencyTEMP; % Frequency 2 (for damping)
    mode_i_damping = handles.mode_i_dampingTEMP; % Damping at Frequency 1
    mode_j_damping = handles.mode_j_dampingTEMP; % Damping at Frequency 2
    
    % Reorganizes Vs, den, su, and GammaMax according to how layers were
    % defined
    
    Vstemp = zeros(N-Base, 1);
    dentemp = zeros(N-Base, 1);
    sutemp = zeros(N-Base, 1);
    GammaMaxtemp = zeros(N-Base, 1);
    for i = 1:length(Layers)
        Vstemp(Layers{i}) = Vs(i);
        dentemp(Layers{i}) = den(i);
        sutemp(Layers{i}) = su(i);
        GammaMaxtemp(Layers{i}) = GammaMax(i);
    end
    
    Vs = Vstemp;
    den = dentemp;
    su = sutemp;
    GammaMax = GammaMaxtemp;
    
    % If EQ data was not Incident, then Incident = 1 and will divide EQ by 2
    % otherwise Incident = 0 and will divide EQ by 1
    AC = EQ*Scale*g/(1 + Incident);
    
    % For flexible base
    if Base == 1
        % Use the bottom most Fmax and Fmin for the extra node
        su(end + 1)=su(end);
        GammaMax(end + 1)=GammaMax(end);
        
        % Creates velocity vector using EQ acceleration
        utr=zeros(length(AC), 1);
        %utr(1)=(AC(1))/2*dt;
        
        for i=2:length(AC) - 1
            utr(i) = utr(i-1) + ((AC(i-1) + AC(i))/2*dt)/2;
        end
        
    end
    
    % Initial Constants
    d0 = 0;
    v0 = 0;
    
    %Implied Declarations
    NTS = size(AC, 1); %Number of Time Steps
    
    G = den.*Vs.^2; %Soil Shear Modulus Vector (Pa)
    
    Rho = den; %mass Density Vector
    
    %Vector/Matrix Initialization
    
    d = zeros(N, 1); %Relative displacement
    v = zeros(N, 1); %Relative velocity
    a = zeros(N, 1); %Relative acceleration
    
    CF = zeros(N, 1); %Damping Force
    KF = zeros(N, 1); %Stiffness Force
    
    DeltaD = zeros(N, 1); %Delta Displacement Vector
    DeltaV = zeros(N, 1); %Delta Velocity Vector
    
    M = zeros(N, N); %Mass Matrix
    mstar = zeros(N, N); %Effective Mass Matrix
    MS = zeros(N, N); %Inverted Mass Matrix
    
    CT = zeros(N, N); %Tangent Damping Matrix
    KT = zeros(N, N); %Tangent Stiffness Matrix
    
    S = zeros(N, 1); %Spring Stiffness Vector
    ST = zeros(N, 1); %Tangent Stiffness for each Spring
    
    Fspring = zeros(N, 1); %Spring Force Vector
    f = zeros(N, 1); %External Force Vector
    
    T = zeros(1, NTS); %Time Array Vector
    
    AA = zeros(N, NTS); %Relative Acceleration Plot
    Vel = zeros(N, NTS); %Relative Velocity Plot
    Disp = zeros(N, NTS); %Relative Displacement Plot
    ACC = zeros(N, NTS); %Absolute Acceleration Plot
    Fspr = zeros(N, NTS); %Spring Force Plot
    
    VP = zeros(N, 1); %Predictor Velocity
    DP = zeros(N, 1); %Predictor Displacement
    
    DA = zeros(N, 1); %Delta Acceleration
    DF = zeros(N, 1); %Delta Force
    DFone = zeros(N, 1); %Delta Force of First Iteration
    
    AN = zeros(N, 1); %Acceleration at Current Time Step
    VN = zeros(N, 1); %Velocity at Current Time Step
    DN = zeros(N, 1); %Displacement at Current Time Step
    
    FspringN = zeros(N, 1); %Spring force at current Time Step
    
    DFm = zeros(NTS, 1);
    
    m = zeros(N, 1); %Active yield surface number
    tau = zeros(N, 1); %Current shear stress
    deltaGamma = zeros(NYS, N);%New strain increment
    Gamma_r = zeros(N, 1); %Model Constant
    alpha = zeros(NYS, N); %Yield surface coordinates
    k = zeros(NYS, N); %Yield surface sizes
    H = zeros(NYS, N); %Tangent moduli
    GT = zeros(N, 1); %Current Tangent moduli
    
    layerDeltaD = zeros(N, 1); %Nodal Displacement
    
    % If the EQ data was single column handles.TTEMP = 0, will need to
    % construct time vector using dt. If EQ data was double column, meaning it
    % has both acceleration and time, then handles.TTEMP will equal that time
    % vector.
    if handles.TTEMP == 0; % handles.TTEMP=0 when EQ data does not provide time vector
        T(2:NTS) = dt:dt:dt*(NTS - 1);
    else
        T = handles.TTEMP'; % handles.TTEMP is the Time vector for EQ data if it was provided
    end
    
    %Initial Velocity and Displacement Vectors
    d = d0*ones(N, 1);
    v = v0*ones(N, 1);
    
    %Define initial Delta Displacement and Delta Velocity
    DeltaD = d;
    DeltaV = v;
    
    %Build Mass Matrix and Spring Stiffness Vector
    
    M(sub2ind([N, N], (1:N-Base), (1:N-Base))) = Rho(1:N-Base).*h(1:N-Base);
    S(1:N-Base) = G(1:N-Base)./h(1:N-Base);
    
    M(1, 1) = M(1, 1)/2;	%this adjust mass of first degree of freedom
    %at the shear beam top (at ground surface) to half value as it should be
    
    % Flexible Base
    if Base == 1
        M(N, N) = M(N - 1, N - 1)/2; %this adjust mass for the extra node
        S(N) = G(N - 1)/h(N - 1);
        G(N) = G(N-1);
    end
    
    %Define Initial Tangent Stiffness for each spring
    %S = G/h
    ST(find(m == 0)) = S(find(m == 0)); %#ok<F%#ok<*MSNU> NDSB>
    ST(find(m > 0)) = H(sub2ind(size(H), m(find(m > 0)), find(m > 0))); %#ok<FND%#ok<*MSNU> SB>
    
    %Assemble Global Tangent Stiffness Matrix KT (Subroutine 2)
    KT = Sub2(ST, N);
    
    if handles.DampType == 1
        
        % Assemble Globtal tangent Damping matrix CT
        w_i = 2*pi*mode_i_frequency;
        w_j = 2*pi*mode_j_frequency;
        
        W_ray = 1/2*[1/w_i, w_i;
            1/w_j, w_j];
        
        % Right now w_i and w_j are being divided by 100 so that they can be
        % entered as percentage (like it is in Cyclic 1D).
        damp_ratio_ray = [mode_i_damping/100; mode_j_damping/100];
        
        a_damping = W_ray\damp_ratio_ray;
        
        a0 = a_damping(1);
        a1 = a_damping(2);
        
        %End of that aside, continue assembling global tangent Damping matrix
        CT = Sub3(a0, a1, M, KT);
        
    else
        
        damping = handles.damping/100;
        CT = 2*damping*M^.5*(M^-.5*KT*M^-.5)^.5*M^.5;
        
        if isreal(CT) ~= 1
            errordlg('Error assembling Damping Matrix.')
            return
        end
        
    end
    
    % Store K and M to calculate natural frequency later
    % For the first time step KT is the unyielded stiffness
    if Base == 0;
        % Rigid Base
        handles.K = KT;
        handles.M = M;
    else
        % Flexible Base
        % For flexible base last stiffness element is defined differently since
        % there is only one spring and not two.
        KT(N, N)= -KT(N-1, N);
        handles.K = KT(1:N-1, 1:N-1);
        handles.M = M(1:N-1, 1:N-1);
        
        % The last Node is being damped by the rock and is represented as a
        % dashpot (cr).
        cr = pr*Vsr;
        CT(N, N) = CT(N, N) + cr;
    end
    
    % Defining Parameters for Stresspoint (Stress - Strain relationship)
    Gamma_r = su.*GammaMax./(G.*GammaMax - su);
    
    if sum(Gamma_r > 0)/length(Gamma_r) ~= 1&& handles.analysis == 2
        errordlg('Check Shear Velocity, Density, Max Shear Strain, and Shear Strength.      (G = den*Vs^2). Please satisfy: G*Max Shear Strain > Shear Strength.',...
            'Check Inputs')
        return
    end
    
    ind = ones(1, NYS);
    suMAT(1:NYS, 1:N) =repmat(su, 1, NYS)';
    Gamma_rMAT = repmat(Gamma_r',NYS,1);
    GMAT = repmat(G',NYS,1);
    
    % Stress for each yield surface
    
    k = meshgrid(1:NYS, ones(N, 1))'.*suMAT/NYS;
    
    % Strain for each yield surface
    deltaGamma(2:NYS, 1:N) = k(2:NYS, 1:N).*Gamma_rMAT(2:NYS, 1:N)./...
        (GMAT(2:NYS, 1:N).*Gamma_rMAT(2:NYS, 1:N) - k(2:NYS, 1:N)) - ...
        k(1:NYS - 1, 1:N).*Gamma_rMAT(2:NYS, 1:N)./(GMAT(2:NYS, 1:N).*...
        Gamma_rMAT(2:NYS, 1:N) - k(1:NYS - 1, 1:N));
    
    % Slope of Stress Strain Curve for each yield surface
    H(1:NYS - 1, :) = (k(2:NYS, :) - k(1:NYS - 1, :))./deltaGamma(2:NYS, :);
    
    NYS = NYS*ones(N,1);
    MatProfile = handles.MatProfileTEMP;
    MatList = handles.MatList;
    
    for count = 1:length(MatProfile)
        
        if MatProfile(count)== 0
        else
            
            a = textread([handles.GUILOCATION, '\Material\', cell2mat(MatList(MatProfile(count),2)), '.txt'], '%s'); %#ok<RE%#ok<*MSNU> MFF1>
            
            count3 = str2num(strcat(cell2mat(a(2))));
            
            if isempty(count3)
                errordlg('Invalid format for material.')
                return
            end
            
            Strain = zeros(count3,1);
            Modulus = zeros(count3,1);
            
            % Assign values to appropriate vectors
            try
                
                for i = 1:count3
                    Strain(i) = str2num(cell2mat(a(2*i+3)));
                    Modulus(i) = str2num(cell2mat(a(2*i+4)));
                end
                
                count4 = str2num(strcat(cell2mat(a(6+count3*2))));
                Strain2 = zeros(count4,1);
                Damping = zeros(count4,1);
                for i = 1:count4
                    Strain2(i) = str2num(cell2mat(a(2*i+7+count3*2)));
                    Damping(i) = str2num(cell2mat(a(2*i+8+count3*2)));
                end
                
                Strain = Strain';
                Modulus = Modulus';
                Strain2 = Strain2';
                Damping = Damping';
                
                StrainOri = Strain;
                ModulusOri = Modulus;
                Strain2Ori = Strain2;
                DampingOri = Damping;
                
                Shear = Strain.*Modulus/100;
                ShearOri = Strain.*Modulus/100;
                ShearStrength = Shear(end);
                Slope = zeros(1, length(Shear));
                
            catch
                
                errordlg('Invalid format for material.')
                return
                
            end
            
            Gmax = G(Layers{count});
            Gmax = Gmax(1,1);
            
            Strain = Strain/100;
            
            Shear = Shear*Gmax;
            
            ShearStrength = su(Layers{count});
            ShearStrength = ShearStrength(1,1);
            
            Slope = zeros(1, length(Shear));
            Slope(1) = Shear(1)/Strain(1);
            for i = 1:length(Shear)-1
                Slope(end-i+1) = (Shear(end-i+1) - Shear(end-i))/(Strain(end-i+1) - Strain(end-i));
            end
            
            Adj1 = 0;
            Adj2 = 0;
            ShearTrial = Shear(end);
            i = 0;
            if ShearTrial > ShearStrength
                
                while Shear(end-i) > ShearStrength
                    Shear(end-i) = ShearStrength;
                    i = i + 1;
                    Adj1 = 1;
                end
                
            elseif ShearStrength > ShearTrial
                
                while ShearStrength > ShearTrial
                    
                    i = i + 1;
                    
                    if i == length(Shear)
                        errordlg('Max Shear Strength too high.')
                        break
                    end
                    
                    ShearTrial = Shear(end-i) + (Strain(end) - Strain(end-i))*Slope(end-i);
                    Adj2 = 1;
                    
                end
                
                if ShearTrial > ShearStrength && Adj2 == 1
                    Shear(end) = ShearStrength;
                elseif ShearTrial <= ShearStrength && Adj2 == 1
                    Shear(end) = ShearTrial;
                end
                
                for j = i-1:-1:1;
                    Shear(end-j) = Shear(end-j-1) + (Strain(end-j) - Strain(end-j-1))*Slope(end-i);
                    if Shear(end-j) > ShearStrength
                        Shear(end-j) = ShearStrength;
                    end
                end
                
            end
            
            Slope = zeros(1, length(Shear));
            for j = 1:length(Shear)-1
                Slope(j) = (Shear(j+1) - Shear(j))/(Strain(j+1) - Strain(j));
            end
            
            if Adj1 == 1
                NYS(Layers{count}) = length(Shear) - i + 1;
            elseif Adj2 == 1
                NYS(Layers{count}) = length(Shear);
            elseif Adj1 == 0 && Adj2 == 0
                i = 0;
                while Slope(end-i-1) == 0;
                    i = i + 1;
                end
                NYS(Layers{count}) = length(Shear) - i;
            end
            
            [m1, ~] = size(k);
            if length(Shear) - i + 1 > m1
                k2 = k;
                H2 = H;
                
                k = zeros(length(Shear) - i + 1,N);
                H = zeros(length(Shear) - i + 1,N);
                alpha = zeros(length(Shear) - i + 1,N);
                
                k(1:m1,1:n1) = k2;
                H(1:m1,1:n1) = H2;
            end
            k(1:NYS(Layers{count}), Layers{count}) = repmat(Shear(1:NYS(Layers{count}))',1,length(Layers{count}));
            H(1:NYS(Layers{count}), Layers{count}) = repmat(Slope(1:NYS(Layers{count}))',1,length(Layers{count}));
            k(NYS(Layers{count})+1:end, Layers{count}) = 0;
            H(NYS(Layers{count})+1:end, Layers{count}) = 0;
        end
        
    end
    
    if Base == 1
        k(:,end) = k(:,end-1);
        H(:,end) = H(:,end-1);
        NYS(end) = NYS(end-1);
    end
    
    for j = 1:N
        [alpha(1:NYS(j),j), m(j), tau(j)] = Stresspoint(alpha(1:NYS(j),j), 0, G(j), H(1:NYS(j),j), k(1:NYS(j),j), m(j), NYS(j), tau(j));
    end
    
    
    %% Begin the solution portion.
    % This first bit is for the very first time step (so called i = 1). It
    % will be treated seperately before the main round.
    
    disp('Time Step 1')
    
    i = 1;
    
    if Base==0
        %Calculate external force vector for rigid base
        f = Sub4(M, AC(i), N);
    else
        %Calculate external force vector for flexible base
        f(N, 1) = 2*pr*Vsr*utr(i);
    end
    
    %Calculate the inverse mass matrix and store
    MS = inv(M);
    
    %Solve for a (Subroutine 5)
    a = Sub5(MS, f, CF, CT, DeltaV, KF, KT, DeltaD);
    
    %Update Damping force and stiffness force for each DOF
    CF = CF + CT*DeltaV;
    KF = KF + KT*DeltaD;
    
    CT_Init = CT;
    KT_Init = KT;
    
    
    %% First time step
    Stress = zeros(N, NTS); %Shear Stress Plot
    Strain = zeros(N, NTS); %Shear Strain Plot
    
    i = 1;
    % Store Initial Time Step (Subroutine 6)
    [AA(:, i), ACC(:, i), Vel(:, i), Disp(:, i), FSpr(:, i), Stress(:, i), ...
        Strain(:, i)] = Sub6(a, AC(i), v, d, tau(:), h, N, Base);
    
    %% Time for main for loop
    
    disp('Main Loop')
    
    % Keep Track of how far along the calculation you are
    mark = ceil(NTS/10);
    count = 1;
    
    msgbox('0%', 'replace')
    
    for i = 2:NTS
        
        if i == count*mark
            msgbox([num2str(count*10) '%'], 'replace')
            count = count + 1;
        end
        
        if Base == 0
            f = Sub4(M, AC(i), N); % Force for rigid base
        else
            f(N, 1) = 2*pr*Vsr*utr(i); % Force for flexible base
        end
        
        %Predictor Phase
        
        VP = v + (dt/2)*a;
        DP = d + v*dt + (dt^2/4)*a;
        
        VN = VP;
        DN = DP;
        AN = zeros(N, 1);
        
        nit = 0; %Set Iterations to 0 before ending while loop
        Bal = 100;
        
        while true
            % The change in nodal displacement
            DeltaD = DN - d;
            % The change in nodal velocity
            DeltaV = VN - v;
            
            % Reseting Stiffness
            
            ST = S;
            
            % Defining vector of nodal displacement
            layerDeltaD(N) = DeltaD(N);
            layerDeltaD(1:end - 1) = DeltaD(1:end - 1) - DeltaD(2:end);
            layerDeltaD = layerDeltaD./h;
            
            % Calculating Stress and location along the stress strain curve
            % The calculated values are trial values, when the loop converges
            % the trial values for when it converged will be stored and used
            % for the next time step
            
            for j = 1:N
                [alphaTEMP(1:NYS(j),j), mTEMP(j), tau(j)] = Stresspoint(alpha(1:NYS(j),j), layerDeltaD(j), G(j),...
                    H(1:NYS(j),j), k(1:NYS(j),j), m(j), NYS(j), Stress(j, i-1));
            end
            
            % Calculating Tangent Stiffness
            % (only for elements that do not have layerDeltaD=0)
            % if layerDeltaD is 0 then use the value of ST from the last time
            % step
            STTemp=ST;
            STTemp(find(layerDeltaD~=0))=(tau(find(layerDeltaD~=0)) - ...
                Stress(find(layerDeltaD~=0), i - 1))./...
                layerDeltaD(find(layerDeltaD~=0))./h(find(layerDeltaD~=0)); %#ok<*FNDSB>
            
            % Get the tangent stiffness matrix
            KT = Sub2(STTemp, N);
            
            % Too be updated later
            %CT = Sub3(a0, a1, M, KT);
            
            if Base==1
                KT(N, N)= - KT(N - 1, N);
                % Too be updated later
                % CT(N, N)=CT(N, N) + cr;
            end
            
            % Calculate the change in force
            DF = f - M*AN - CF - CT*DeltaV - KF - KT*DeltaD;
            
            if nit == 0
                DFone = DF;
                
                SumSQDFone = sqrt(sum(DFone.^2));
            end
            
            SumSQDF = sqrt(sum(DF.^2));
            
            Bal = SumSQDF/SumSQDFone;
            
            % SumSQDFone==0 is to solve the problem of dividing by zero
            if Bal < Tol||SumSQDFone==0
                break
            end
            
            mstar = M + CT_Init*(dt/2) + KT_Init*(dt^2/4);
            MS = inv(mstar);
            
            % The change in acceleration
            DA = MS*DF;
            
            % Update the acceleration for the current step
            AN = AN + DA;
            % Update the velocity for the current step
            VN = VP + (dt/2)*AN;
            % Update the displacment for the current step
            DN = DP + (dt^2/4)*AN;
            
            if nit == 5000
                msgbox(['Failed to Converge after 5000 iterations. Time: '...
                    num2str(T(i))], 'replace')
                return
            end
            
            % Increment the itterator
            nit = nit + 1;
            
        end
        
        % Loop converged. Storing values.
        % if nit > 1
        % fprintf('%s%i%s%i%s\n', 'Itterated on step ', i, ' ', nit, ' times')
        % end
        
        % Updating Damping force (CF), spring force (KF), velocity (v),
        % displacement (d), and acceleration (a). These values will be used as
        % parameters for the next time step.
        [CF, KF, v, d, a] = Sub7(CF, CT, DeltaV, KF, KT, DeltaD, ...
            VN, DN, AN);
        
        % Storing the trial values. (These values are the values for when the
        % loop converged. alpha (yield surface coordinates), m (active yield
        % surface), ST (tangent stiffness). ST is being stored only to deal
        % with the case the case of dividing by 0. If layerDeltaD is zero, then
        % you do not update that part of ST, and use the value from the last
        % time step.
        alpha=alphaTEMP;
        m=mTEMP;
        ST=STTemp;
        
        [AA(:, i), ACC(:, i), Vel(:, i), Disp(:, i), FSpr(:, i), Stress(:, i), ...
            Strain(:, i)] = Sub6(a, AC(i), v, d, tau, h, N, Base);
        
    end
    
    % If flexible base, change velocity and displacement from total to relative
    if Base == 1
        for i = 1:N - Base
            Disp(i, :) = Disp(i, :) - Disp(N, :);
            Vel(i, :) = Vel(i, :) - Vel(N, :);
        end
    end
    
    %% Analysis Completed. Saving variables, and enabling options for user.
    disp('Anaylsis Completed/Saving Outputs')
    
    handles.T = T;
    handles.AC = AC;
    handles.Disp = Disp;
    handles.Vel = Vel;
    handles.ACC = ACC;
    handles.Strain = Strain*100;
    handles.Stress = Stress/1000;
    handles.out = 0;
    guidata(hObject, handles);
    
    handles = OutputOn(handles);
    
    msgbox('Complete', 'replace')
    
    guidata(hObject, handles);
end


% --------------------------------------------------------------------
function SaveOutput_Callback(hObject, eventdata, handles)
% hObject    handle to SaveOutput (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

T = handles.T';
AC = handles.AC;
Disp = handles.Disp';
Vel = handles.Vel';
ACC = handles.ACC';
Strain = handles.Strain';
Stress = handles.Stress';
N = handles.NTEMP + handles.BaseTEMP;
FileName = strrep(handles.FileNameTEMP, '.mat', '');
Base = handles.BaseTEMP;

[SaveFile b] = uiputfile('.xlsx', 'Save File',...
    [handles.GUILOCATION, '\Saved Outputs']);
File = [b, SaveFile];

if b == 0
else
    
    A0 = [T AC];
    A = zeros(length(T), 2);
    A = {'Time' 'Input Motion'};
    
    for i = 1:length(T)
        for j = 1:2
            A{i + 1, j} = A0(i, j);
        end
    end
    
    if handles.unitTEMP == 1;
        unit{1} = ' (m)';
        unit{2} = ' (m/s)';
        unit{3} = ' (m/s^2)';
        unit{4} = ' (%)';
        unit{5} = ' (kPa)';
        unit{6} = ' (s)';
    else
        unit{1} = ' (in)';
        unit{2} = ' (in/s)';
        unit{3} = ' (in/s^2)';
        unit{4} = ' (%)';
        unit{5} = ' (ksi';
        unit{6} = ' (s)';
    end
    
    for j = 1:N
        if Base == 1&&j == N
            
            B(1, :, j) = {'', 'Node', ' ', ' ', ' ', ' '}; %#ok<*AGROW>
            B(2, :, j) = {'Time (s)' ['Displacement', unit{1}]...
                ['Velocity', unit{2}] ['Acceleration', unit{3}], ' ', ' '};
            n = 1:4;
            
        else
            
            B(1, :, j) = {'', 'Node', '', '', 'Element', ''};
            B(2, :, j) = {'Time (s)', ['Displacement', unit{1}], ['Velocity',...
                unit{2}], ['Acceleration', unit{3}], 'Strain (%)',...
                ['Stress', unit{5}]};
            n = 1:6;
            
        end
        
        B0 = [T, Disp(:, j), Vel(:, j), ACC(:, j), Strain(:, j), Stress(:, j)];
        
        for i = 1:length(T)
            for k = n
                B{i + 2, k, j} = B0(i, k); %#ok<AGROW>
            end
        end
    end
    
    Excel = actxserver('Excel.Application');
    
    set(Excel, 'Visible', 1);
    Workbooks = Excel.Workbooks;
    Workbook = invoke(Workbooks, 'Add');
    Worksheets = Excel.Worksheets;
    
    Sheets = Excel.ActiveWorkBook.Sheets;
    
    for i = 1:N + 1
        
        if Base == 1 && i == N + 1
            k = 1:3;
        else
            k = 1:6;
        end
        
        if i == 1
            sheet = get(Sheets, 'Item', i);
            invoke(sheet, 'Activate');
            Activesheet = Workbook.Activesheet;
            ran=strcat('A1:B', num2str(length(T) + 1));
            
            try
                Chart = Excel.ActiveSheet.Shapes.AddChart;
            catch
                Chart = Excel.ActiveSheet.Shapes.AddChart2;
            end
            Chart.Name = 'NewChart';
            NewChart = Excel.ActiveSheet.ChartObjects('NewChart');
            NewChart.Activate;
            try
                NewChart.Chart.Legend.Delete;
            end
            
            Chart.Width = 400;
            Chart.Height= 200;
            Chart.Left = 400;
            Chart.Top = 0;
            
            ActivesheetRange = get(Activesheet, 'Range', ran);
            set(ActivesheetRange, 'Value', A);
            Sheets.Item(1).Name = 'Input Motion';
            
            NewChart.Activate;
            NewSeries = invoke(Excel.ActiveChart.SeriesCollection, 'NewSeries');
            NewSeries.XValues = Excel.ActiveSheet.Range(strcat('A2:A',...
                num2str((length(T) + 1))));
            NewSeries.Values = Excel.ActiveSheet.Range(strcat('B2:B',...
                num2str((length(T) + 1))));
            
            Excel.ActiveChart.ChartType = 'xlXYScatterLinesNoMarkers';
            Axes = invoke(Excel.ActiveChart, 'Axes', 1);
            set(Axes, 'HasTitle', 1);
            set(Axes.AxisTitle, 'Caption', ['Time', unit{6}])
            Axes = invoke(Excel.ActiveChart, 'Axes', 2);
            set(Axes, 'HasTitle', 1);
            set(Axes.AxisTitle, 'Caption', ['Acceleration', unit{3}])
            Excel.ActiveChart.HasTitle = 1;
            Excel.ActiveChart.ChartTitle.Characters.Text = 'Input Motion';
            
        elseif i == 2||i == 3
            
            try
                sheet = get(Sheets, 'Item', i);
                invoke(sheet, 'Activate');
                Activesheet = Workbook.Activesheet;
                ran = strcat('A1:F', num2str(length(T) + 2));
                
                for j = 1:6
                    
                    try
                        Chart = Excel.ActiveSheet.Shapes.AddChart;
                    catch
                        Chart = Excel.ActiveSheet.Shapes.AddChart2;
                    end
                    Chart.Name = strcat('Chart', j);
                    NewChart = Excel.ActiveSheet.ChartObjects(strcat('Chart', j));
                    NewChart.Activate;
                    try
                        NewChart.Chart.Legend.Delete;
                    end
                    
                    if j < 4
                        left = 400;
                    else
                        left = 800;
                    end
                    
                    if j == 1||j == 4
                        top = 0;
                    elseif j == 2||j == 5
                        top = 200;
                    else
                        top = 400;
                    end
                    
                    Chart.Width = 400;
                    Chart.Height= 200;
                    Chart.Left = left;
                    Chart.Top = top;
                    
                end
                
                ActivesheetRange = get(Activesheet, 'Range', ran);
                set(ActivesheetRange, 'Value', B(:, :, i - 1));
                Sheets.Item(i).Name=['Node and Element ', num2str(i - 1)];
                
                for j = 1:6
                    
                    if j == 1
                        let = 'B';
                        let2 = 'A';
                        type = ' Node :';
                        type2 = 'Displacement';
                    elseif j == 2
                        let = 'C';
                        let2 = 'A';
                        type = ' Node :';
                        type2 = 'Velocity';
                    elseif j == 3
                        let = 'D';
                        let2 = 'A';
                        type = ' Node :';
                        type2 = 'Acceleration';
                    elseif j == 4
                        let = 'E';
                        let2 = 'A';
                        type = ' Element :';
                        type2 = 'Strain';
                    elseif j == 5
                        let = 'F';
                        let2 = 'A';
                        type = ' Element :';
                        type2 = 'Stress';
                    else
                        let = 'F';
                        let2 = 'E';
                        type = 'Element :';
                        type2 = 'Stress vs. Strain';
                    end
                    
                    Chart.Name = strcat('Chart', j);
                    NewChart = Excel.ActiveSheet.ChartObjects(strcat('Chart', j));
                    NewChart.Activate;
                    
                    NewSeries = invoke(Excel.ActiveChart.SeriesCollection, 'NewSeries');
                    NewSeries.XValues = Excel.ActiveSheet.Range(strcat(let2, '3:',...
                        let2, num2str((length(T) + 2))));
                    NewSeries.Values = Excel.ActiveSheet.Range(strcat(let, '3:',...
                        let, num2str((length(T) + 2))));
                    
                    Excel.ActiveChart.ChartType = 'xlXYScatterLinesNoMarkers';
                    Axes = invoke(Excel.ActiveChart, 'Axes', 1);
                    
                    if j ~= 6
                        set(Axes, 'HasTitle', 1);
                        set(Axes.AxisTitle, 'Caption', strcat('Time', unit{6}))
                        Axes = invoke(Excel.ActiveChart, 'Axes', 2);
                        set(Axes, 'HasTitle', 1);
                        set(Axes.AxisTitle, 'Caption', strcat(type2, unit{j}))
                    else
                        set(Axes, 'HasTitle', 1);
                        set(Axes.AxisTitle, 'Caption', 'Strain')
                        Axes = invoke(Excel.ActiveChart, 'Axes', 2);
                        set(Axes, 'HasTitle', 1);
                        set(Axes.AxisTitle, 'Caption', strcat('Stress', unit{5}))
                    end
                    
                    Excel.ActiveChart.HasTitle = 1;
                    
                    Excel.ActiveChart.ChartTitle.Characters.Text = ...
                        [type2, ' for', type, ' ', num2str(i - 1)];
                end
            catch
                Worksheets.Add([], Worksheets.Item(Worksheets.Count));
                Activesheet = Workbook.Activesheet;
                ran=strcat('A1:F', num2str(length(T) + 2));
                
                for j = k
                    
                    try
                        Chart = Excel.ActiveSheet.Shapes.AddChart;
                    catch
                        Chart = Excel.ActiveSheet.Shapes.AddChart2;
                    end
                    
                    Chart.Name = strcat('Chart', j);
                    NewChart = Excel.ActiveSheet.ChartObjects(strcat('Chart', j));
                    NewChart.Activate;
                    
                    try
                        NewChart.Chart.Legend.Delete;
                    end
                    
                    if j < 4
                        left = 400;
                    else
                        left = 800;
                    end
                    
                    if j == 1||j == 4
                        top = 0;
                    elseif j == 2||j == 5
                        top = 200;
                    else
                        top = 400;
                    end
                    
                    Chart.Width = 400;
                    Chart.Height = 200;
                    Chart.Left = left;
                    Chart.Top = top;
                    
                end
                
                ActivesheetRange = get(Activesheet, 'Range', ran);
                set(ActivesheetRange, 'Value', B(:, :, i - 1));
                Sheets.Item(i).Name = ['Node and Element', ' ', num2str(i - 1)];
                
                for j = k
                    
                    if j == 1
                        let = 'B';
                        let2 = 'A';
                        type = ' Node :';
                        type2 = 'Displacement';
                    elseif j == 2
                        let = 'C';
                        let2 = 'A';
                        type = ' Node :';
                        type2 = 'Velocity';
                    elseif j ==  3
                        let = 'D';
                        let2 = 'A';
                        type = ' Node :';
                        type2 = 'Acceleration';
                    elseif j == 4
                        let = 'E';
                        let2 = 'A';
                        type = ' Element :';
                        type2 = 'Strain';
                    elseif j == 5
                        let = 'F';
                        let2 = 'A';
                        type = ' Element :';
                        type2 = 'Stress';
                    else
                        let = 'F';
                        let2 = 'E';
                        type = ' Element :';
                        type2 = 'Stress vs. Strain';
                    end
                    
                    Chart.Name = strcat('Chart', j);
                    NewChart = Excel.ActiveSheet.ChartObjects(strcat('Chart', j));
                    NewChart.Activate;
                    
                    NewSeries = invoke(Excel.ActiveChart.SeriesCollection, 'NewSeries');
                    NewSeries.XValues = Excel.ActiveSheet.Range(strcat(let2, '3:', let2, num2str((length(T) + 2))));
                    NewSeries.Values = Excel.ActiveSheet.Range(strcat(let, '3:', let, num2str((length(T) + 2))));
                    
                    Excel.ActiveChart.ChartType = 'xlXYScatterLinesNoMarkers';
                    Axes = invoke(Excel.ActiveChart, 'Axes', 1);
                    
                    if j ~= 6
                        set(Axes, 'HasTitle', 1);
                        set(Axes.AxisTitle, 'Caption', strcat('Time', unit{6}))
                        Axes = invoke(Excel.ActiveChart, 'Axes', 2);
                        set(Axes, 'HasTitle', 1);
                        set(Axes.AxisTitle, 'Caption', strcat(type2, unit{j}))
                    else
                        set(Axes, 'HasTitle', 1);
                        set(Axes.AxisTitle, 'Caption', 'Strain')
                        Axes = invoke(Excel.ActiveChart, 'Axes', 2);
                        set(Axes, 'HasTitle', 1);
                        set(Axes.AxisTitle, 'Caption', strcat('Stress', unit{5}))
                    end
                    
                    Excel.ActiveChart.HasTitle = 1;
                    
                    Excel.ActiveChart.ChartTitle.Characters.Text = [type2, ' for', type, ' ', num2str(i - 1)];
                    
                end
            end
            
        else
            Worksheets.Add([], Worksheets.Item(Worksheets.Count));
            Activesheet = Workbook.Activesheet;
            ran=strcat('A1:F', num2str(length(T) + 2));
            
            for j = k
                
                try
                    Chart = Excel.ActiveSheet.Shapes.AddChart;
                catch
                    Chart = Excel.ActiveSheet.Shapes.AddChart2;
                end
                Chart.Name = strcat('Chart', j);
                NewChart = Excel.ActiveSheet.ChartObjects(strcat('Chart', j));
                NewChart.Activate;
                try
                    NewChart.Chart.Legend.Delete;
                end
                
                if j < 4
                    left = 400;
                else
                    left = 800;
                end
                
                if j == 1||j == 4
                    top = 0;
                elseif j == 2||j == 5
                    top = 200;
                else
                    top = 400;
                end
                
                Chart.Width = 400;
                Chart.Height = 200;
                Chart.Left = left;
                Chart.Top = top;
                
            end
            
            ActivesheetRange = get(Activesheet, 'Range', ran);
            set(ActivesheetRange, 'Value', B(:, :, i - 1));
            Sheets.Item(i).Name = ['Node and Element', ' ', num2str(i - 1)];
            
            for j = k
                
                if j == 1
                    let = 'B';
                    let2 = 'A';
                    type = ' Node :';
                    type2 = 'Displacement';
                elseif j == 2
                    let = 'C';
                    let2 = 'A';
                    type = ' Node :';
                    type2 = 'Velocity';
                elseif j ==  3
                    let = 'D';
                    let2 = 'A';
                    type = ' Node :';
                    type2 = 'Acceleration';
                elseif j == 4
                    let = 'E';
                    let2 = 'A';
                    type = ' Element :';
                    type2 = 'Strain';
                elseif j == 5
                    let = 'F';
                    let2 = 'A';
                    type = ' Element :';
                    type2 = 'Stress';
                else
                    let = 'F';
                    let2 = 'E';
                    type = ' Element :';
                    type2 = 'Stress vs. Strain';
                end
                
                Chart.Name = strcat('Chart', j);
                NewChart = Excel.ActiveSheet.ChartObjects(strcat('Chart', j));
                NewChart.Activate;
                
                NewSeries = invoke(Excel.ActiveChart.SeriesCollection, 'NewSeries');
                NewSeries.XValues = Excel.ActiveSheet.Range(strcat(let2, '3:', let2, num2str((length(T) + 2))));
                NewSeries.Values = Excel.ActiveSheet.Range(strcat(let, '3:', let, num2str((length(T) + 2))));
                
                Excel.ActiveChart.ChartType = 'xlXYScatterLinesNoMarkers';
                Axes = invoke(Excel.ActiveChart, 'Axes', 1);
                
                if j ~= 6
                    set(Axes, 'HasTitle', 1);
                    set(Axes.AxisTitle, 'Caption', strcat('Time', unit{6}))
                    Axes = invoke(Excel.ActiveChart, 'Axes', 2);
                    set(Axes, 'HasTitle', 1);
                    set(Axes.AxisTitle, 'Caption', strcat(type2, unit{j}))
                else
                    set(Axes, 'HasTitle', 1);
                    set(Axes.AxisTitle, 'Caption', 'Strain')
                    Axes = invoke(Excel.ActiveChart, 'Axes', 2);
                    set(Axes, 'HasTitle', 1);
                    set(Axes.AxisTitle, 'Caption', strcat('Stress', unit{5}))
                end
                
                Excel.ActiveChart.HasTitle = 1;
                
                Excel.ActiveChart.ChartTitle.Characters.Text = [type2, ' for', type, ' ', num2str(i - 1)];
                
            end
        end
    end
    
    Workbook.SaveAs(File);
    Workbook.Saved = 1;
    invoke(Workbook, 'Close');
    
    % Quit Excel
    invoke(Excel, 'Quit');
    
    % End process
    delete(Excel);
    
end


% --------------------------------------------------------------------
function Plot_Callback(hObject, eventdata, handles)
% hObject    handle to Plot (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function MAV_Callback(hObject, eventdata, handles)
% hObject    handle to MAV (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function MAVA_Callback(hObject, eventdata, handles)
% hObject    handle to MAVA (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Checks if user ran analysis
if handles.ACC == 0
    errordlg('Please run Analysis first.')
    return
end

% Obtain results from analysis
ACC = handles.ACC;
Base = handles.BaseTEMP;
h = handles.hTEMP;
N = handles.NTEMP;

x = zeros(N, 1);

for i = 1:N + Base
    x(i) = max(abs(ACC(i, :)));
end

% Add the acc of the last node (Last node acc is equal to input motion)
if Base == 0
    x(end + 1) = max(handles.AC);
end

% Plot result
plot(x, 0:h:N*h)

% Set limit on axis
dif = abs(max(x) - min(x))*0.05;
Llim = min(x) - dif;
Rlim = max(x) + dif;
axis([Llim Rlim 0 N*h]);

set(gca, 'YDir', 'reverse')

% Label axis
if handles.unitTEMP == 1
    ylabel('depth (m)')
    xlabel('(m/s/s)')
else
    ylabel('depth (in.)')
    xlabel('(in./s/s)')
end


% --------------------------------------------------------------------
function MAVV_Callback(hObject, eventdata, handles)
% hObject    handle to MAVV (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Checks if user ran analysis
if handles.Vel == 0
    errordlg('Please run Analysis first.')
    return
end

% Obtain results from analysis
Vel = handles.Vel;
Base = handles.BaseTEMP;
h = handles.hTEMP;
N = handles.NTEMP;

x = zeros(N, 1);
for i = 1:N
    x(i) = max(abs(Vel(i, :)));
end

x(end + 1) = 0;

% Plot output
plot(x, 0:h:N*h)

% Set limit on axis
dif = abs(max(x) - min(x))*0.05;
Llim = min(x) - dif;
Rlim = max(x) + dif;
axis([Llim Rlim 0 N*h]);

set(gca, 'YDir', 'reverse')

% Labels axis
if handles.unitTEMP == 1
    ylabel('depth (m)')
    xlabel('(m/s)')
else
    ylabel('depth (in.)')
    xlabel('(in./s)')
end


% --------------------------------------------------------------------
function MAVD_Callback(hObject, eventdata, handles)
% hObject    handle to MAVD (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Checks if user ran analysis
if handles.Disp == 0
    errordlg('Please run Analysis first.')
    return
end

% Obtains results from analysis
Disp = handles.Disp;
Base = handles.BaseTEMP;
h = handles.hTEMP;
N = handles.NTEMP;

x = zeros(N,1);
for i = 1:N
    x(i) = max(abs(Disp(i, :)));
end

x(end + 1) = 0;

% Plot output
plot(x, 0:h:N*h)

% Sets limit on axis
dif = abs(max(x) - min(x))*0.05;
Llim = min(x) - dif;
Rlim = max(x) + dif;
axis([Llim Rlim 0 N*h]);

set(gca, 'YDir', 'reverse')

% Labels axis
if handles.unitTEMP == 1
    ylabel('depth (m)')
    xlabel('(m)')
else
    ylabel('depth (in.)')
    xlabel('(in.')
end


% --------------------------------------------------------------------
function MaxMin_Callback(hObject, eventdata, handles)
% hObject    handle to MaxMin (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function MMA_Callback(hObject, eventdata, handles)
% hObject    handle to MMA (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Checks if user ran analysis
if handles.ACC == 0
    errordlg('Please run Analysis first.')
    return
end

% Obtain result of analysis
ACC=handles.ACC;
Base=handles.BaseTEMP;
h=handles.hTEMP;
N=handles.NTEMP;

x1 = zeros(N, 1);
x2 = zeros(N, 1);
for i=1:N + Base
    x1(i) = max(ACC(i, :));
    x2(i) = min(ACC(i, :));
end

% If rigid base add the base node (which is zero acc)
if Base == 0
    x1(end + 1) = max(handles.AC);
    x2(end + 1) = min(handles.AC);
end

% Plots output
plot(x1, 0:h:N*h)
hold on
plot(x2, 0:h:N*h)
hold off

% Sets limit on axis
dif = abs(max(x1) - min(x2))*0.05;
Llim = min(x2) - dif;
Rlim = max(x1) + dif;
axis([Llim Rlim 0 N*h]);

set(gca, 'YDir', 'reverse')

% Labels axis
if handles.unitTEMP == 1
    ylabel('depth (m)')
    xlabel('(m/s/s)')
else
    ylabel('depth (in.)')
    xlabel('(in./s/s)')
end


% --------------------------------------------------------------------
function MMV_Callback(hObject, eventdata, handles)
% hObject    handle to MMV (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Checks if user ran analysis
if handles.Vel == 0
    errordlg('Please run Analysis first.')
    return
end

% Obtain result of analysis
Vel = handles.Vel;
Base = handles.BaseTEMP;
h = handles.hTEMP;
N = handles.NTEMP;

x1 = zeros(N, 1);
x2 = zeros(N, 1);
for i = 1:N
    x1(i)=max(Vel(i, :));
    x2(i)=min(Vel(i, :));
end

x1(end + 1)=0;
x2(end + 1)=0;

% PLots output
plot(x1, 0:h:N*h)
hold on
plot(x2, 0:h:N*h)
hold off

% Sets limits on axis
dif = abs(max(x1) - min(x2))*0.05;
Llim = min(x2) - dif;
Rlim = max(x1) + dif;
axis([Llim Rlim 0 N*h]);

set(gca, 'YDir', 'reverse')

% Labels axis
if handles.unitTEMP == 1
    ylabel('depth (m)')
    xlabel('(m/s)')
else
    ylabel('depth (in.)')
    xlabel('(in./s)')
end


% --------------------------------------------------------------------
function MMD_Callback(hObject, eventdata, handles)
% hObject    handle to MMD (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Checks if user ran analysis
if handles.Disp == 0
    errordlg('Please run Analysis first.')
    return
end

% Obtains results from analysis
Disp = handles.Disp;
Base = handles.BaseTEMP;
h = handles.hTEMP;
N = handles.NTEMP;

x1 = zeros(N, 1);
x2 = zeros(N, 1);
for i = 1:N
    x1(i)=max(Disp(i, :));
    x2(i)=min(Disp(i, :));
end

x1(end + 1) = 0;
x2(end + 1) = 0;

% Plots output
plot(x1, 0:h:N*h)
hold on
plot(x2, 0:h:N*h)
hold off

% Sets limits on axis
dif = abs(max(x1) - min(x2))*0.05;
Llim = min(x2) - dif;
Rlim = max(x1) + dif;
axis([Llim Rlim 0 N*h]);

set(gca, 'YDir', 'reverse')

% Labels axis
if handles.unitTEMP == 1
    ylabel('depth (m)')
    xlabel('(m)')
else
    ylabel('depth (in.)')
    xlabel('(in.)')
end


% --------------------------------------------------------------------
function Animation_Callback(hObject, eventdata, handles)
% hObject    handle to Animation (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function Create_Callback(hObject, eventdata, handles)
% hObject    handle to Create (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Checks if user ran analysis
if handles.Disp == 0
    errordlg('Please run Analysis first.')
    return
end

% Obtain desired file name for animation

FileName = strrep(handles.FileNameTEMP, '.mat', '');
[SaveFile b] = uiputfile('.avi', 'Save File',...
    [handles.GUILOCATION '\Saved Animations']);

if b == 0
    return
end

% Loads values from interface and from analysis
Base = handles.BaseTEMP;
dt = handles.dtTEMP;
N = handles.NTEMP;
height = handles.hTEMP;
Disp = handles.Disp;
AC = handles.AC;
unit = handles.unitTEMP;

Disp1 = Disp';
Disp2 = Disp1;

if Base == 0
    Disp2(:, end + 1) = 0;
else
    Disp1(:, end) = 0;
    Disp2(:, end) = 0;
end

figure
h = gcf;
set(figure(h), 'visible', 'off')
set(figure(h), 'resize', 'off')

% Formats avi
fps = 1/dt;

n = 1;
while fps>30
    n = n + 1;
    fps = (1/dt)/n;
end

if exist([b SaveFile], 'file')
    delete([b SaveFile]);
end

aviobj = avifile([b SaveFile], 'fps', fps, 'compression', 'FFDS');
aviobj.Quality = 0;

lim = max(abs(1.05*min(min(Disp1))), abs(1.05*max(max(Disp1))));

% Select unit for labeling axis
if unit == 1
    unit1 = '(m)';
    unit2 = '(m/s/s)';
else
    unit1 = '(ft)';
    unit2 = '(ft/s/s)';
end

try
    % Create animation
    for i = 1:n:length(Disp)
        
        figure(h)
        figure(h)
        subplot(2, 1, 2, 'position', [.1 (1/3*.3) .8 1/3*.6]);
        plot(handles.T, AC)
        xlabel('Time (sec)')
        ylabel(['Acc ', unit2])
        
        hold on
        g = plot(handles.T(i), AC(i), 'or');
        set(g, 'MarkerEdgeColor', 'none', 'MarkerFaceColor', 'r')
        
        subplot(2, 1, 1, 'position', [.1 1/3*1.2 .8 2/3*.8]);
        
        plot(Disp2(i, :), 0:height:height*(N), ' - ')
        axis([-lim lim -.25*height height*(N)+height*0.25]);
        set(gca,'YDir','reverse');
        xlabel(['Disp ', unit1])
        ylabel(['Depth ', unit1])
        
        str1(1) = {strcat('Time : ', num2str(dt*i), ' sec')};
        title(str1)
        
        aviobj = addframe(aviobj, h); %#ok<NOD%#ok<*MSNU> EF>
        
    end
    
    close(h)
    aviobj = close(aviobj); %#ok<NODEF,RE%#ok<*MSNU,*MSNU> DEF>
    
    msgbox('Complete')
    
catch
    
    errordlg('Error creating .avi file. Potential source of error 1) closing figures   2) Codec is not installed. (Need ffdshow.)')
    try
        aviobj = close(aviobj);
    end
    
end


% --------------------------------------------------------------------
function Play_Callback(hObject, eventdata, handles)
% hObject    handle to Play (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

v1 = 0;

% Obtain list of file names of saved animations
d = dir([pwd, '\Saved Animations\']);
str = {d.name};
str = str(3:end);

if isempty(str)
    errordlg('Animations found. Please create animation first.')
    return
end

% Prompts user to select one of the saved animations
[s1,v1] = listdlg('PromptString','Select a file:',...
    'SelectionMode','single',...
    'ListString',str);

if v1 ~= 0
    
    try
        % Plays selected animation
        implay([pwd, '\Saved Animations\', cell2mat(str(s1))]);
    catch
        errordlg('Error playing video.')
    end
end


% --------------------------------------------------------------------
function Close_Callback(hObject, eventdata, handles)
% hObject    handle to Close (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function CloseFigures_Callback(hObject, eventdata, handles)
% hObject    handle to CloseFigures (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% Obtain number of figure
nfigs = get(0, 'Children');

% Close all figures
for i = 1:length(nfigs)
    if round(nfigs(i)) == nfigs(i)
        close(nfigs(i))
    end
end


% --------------------------------------------------------------------
function Help_Callback(hObject, eventdata, handles)
% hObject    handle to Help (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function OUM_Callback(hObject, eventdata, handles)
% hObject    handle to OUM (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

open([handles.GUILOCATION, '\User Manual\User Manual.pdf'])


% --------------------------------------------------------------------
function About_Callback(hObject, eventdata, handles)
% hObject    handle to About (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

msgbox('     SB2013 is a free software, that is developed to be capable of performing 1D earthquake analysis of soil.', 'About SB2013')


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

Base = handles.BaseTEMP; % (Rigid Base = 0; Flexible Base = 1)
N = handles.NTEMP; % Number of Elements

str = handles.LayersTEMP; % Assigns Layer (used to define different Vs,
% den, etc.)

A = strread(str, '%s', 'delimiter', ' '); %#ok<REMFF1> % Layers is initially stored as a string.

for i = 1:length(A)
    Layers{i} = str2num(cell2mat(A(i))); %#ok<AGROW> % Changes Layers to numbers.
end

tot = 0;
for i = 1:length(A)
    tot = tot + length(Layers{i}); % Used to check length of layer
end

check = zeros(length(A), 1);
for i = 1:length(A)
    check(i) = max(Layers{i}); % Used to check if there is N number of elements is defined
end

n = 0;
for i = 1:length(A)
    test = Layers{i};
    for j = 1:length(test)
        n = n + 1;
        check2(n) = test(j); %#ok<AGROW> % Used to check for repeated values in Layers
    end
end

if max(check) ~= N
    errordlg('Please change # of Ele/Node to match Elements or vice versa.', ...
        '# of Ele/Nodes do not match Elements.');
    return
elseif tot ~= N
    errordlg('Please check Element vector some values are missing. (ex. 1:3 5:6 Number 4 is missing).', ...
        'Missing Values in Element Vector')
    return
elseif length(check2) ~= length(unique(check2))
    errordlg('Please check Element vector for duplicates', ...
        'Duplicate Entries in Element Vector')
    return
elseif sum(mod(check2, 1)) ~= 0
    errordlg('Please check that all entries in element vector are integers.', ...
        'Noninteger Entries in Element Vector')
    return
elseif length(Layers) ~= length(handles.VsTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Velocity should have 3 inputs.',...
        'Please check Velocity')
    return
elseif length(handles.VsTEMP) ~= sum(handles.VsTEMP>0)
    errordlg('Enter Velocity as positive values.', 'Error in Velocity')
    return
elseif length(Layers) ~= length(handles.suTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Shear Strength should have 3 inputs.',...
        'Please check Shear Strength')
    return
elseif length(Layers) ~= length(handles.GammaMaxTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Max Shear Strain should have 3 inputs.',...
        'Please check Max Shear Strain')
    return
elseif length(Layers) ~= length(handles.denTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Density should have 3 inputs.', 'Please check Density')
    return
elseif length(handles.denTEMP) ~= sum(handles.denTEMP>0)
    errordlg('Density should all be positive.', 'Error in density')
    return
elseif length(Layers) ~= length(handles.MatProfileTEMP)
    errordlg('The inputs for the vector should match Elements defined. (ex. If Elements was 1:10 11:20 21 then Material should have 3 inputs. Can enter 0 to use default Stress Strain logic.',...
        'Please check Material')
    return
end

NYS = handles.NYSTEMP; % Number of Yield Surfaces

N = N + Base; % Increases N by 1 if flexible otherwise N stays the same
h = handles.hTEMP*ones(N, 1); % Creates Vector of height (for each element)

Vs = handles.VsTEMP*10^handles.VsMultiTEMP; % Vector of Shear Wave Velocity

den = handles.denTEMP*10^handles.denMultiTEMP; % Vector Density
su = handles.suTEMP*10^handles.suMultiTEMP*1000; % Vector of Shear Strength Converts kPa/ksi to Pa/psi to keep units consistent
GammaMax = handles.GammaMaxTEMP*10^handles.GammaMaxMultiTEMP/100; % Vector of Max Shear Strain Converts (%) to decimal
mode_i_frequency = handles.mode_i_frequencyTEMP; % Frequency 1 (for damping)
mode_j_frequency = handles.mode_j_frequencyTEMP; % Frequency 2 (for damping)
mode_i_damping = handles.mode_i_dampingTEMP; % Damping at Frequency 1
mode_j_damping = handles.mode_j_dampingTEMP; % Damping at Frequency 2

% Reorganizes Vs, den, su, and GammaMax according to how layers were
% defined

Vstemp = zeros(N-Base, 1);
dentemp = zeros(N-Base, 1);
sutemp = zeros(N-Base, 1);
GammaMaxtemp = zeros(N-Base, 1);
for i = 1:length(Layers)
    Vstemp(Layers{i}) = Vs(i);
    dentemp(Layers{i}) = den(i);
    sutemp(Layers{i}) = su(i);
    GammaMaxtemp(Layers{i}) = GammaMax(i);
end

Vs = Vstemp;
den = dentemp;
su = sutemp;
GammaMax = GammaMaxtemp;

if Base == 1
    % Use the bottom most Fmax and Fmin for the extra node
    su(end + 1)=su(end);
    GammaMax(end + 1)=GammaMax(end);
end

G = den.*Vs.^2; %Soil Shear Modulus Vector (Pa)

m = zeros(N, 1); %Active yield surface number
tau = zeros(N, 1); %Current shear stress
deltaGamma = zeros(NYS, N);%New strain increment
Gamma_r = zeros(N, 1); %Model Constant
alpha = zeros(NYS, N); %Yield surface coordinates
k = zeros(NYS, N); %Yield surface sizes
H = zeros(NYS, N); %Tangent moduli
GT = zeros(N, 1); %Current Tangent moduli
% Flexible Base
if Base == 1
    G(N) = G(N-1);
end

Gamma_r = su.*GammaMax./(G.*GammaMax - su);

if sum(Gamma_r > 0)/length(Gamma_r) ~= 1&& handles.analysis == 2
    errordlg('Check Shear Velocity, Density, Max Shear Strain, and Shear Strength. (G = den*Vs^2) G*Max Shear Strain should be greater than Shear Strength.',...
        'Check Inputs')
    return
end

ind = ones(1, NYS);
suMAT(1:NYS, 1:N) =repmat(su, 1, NYS)';
Gamma_rMAT = repmat(Gamma_r',NYS,1);
GMAT = repmat(G',NYS,1);

% Stress for each yield surface

k = meshgrid(1:NYS, ones(N, 1))'.*suMAT/NYS;

% Strain for each yield surface
deltaGamma(2:NYS, 1:N) = k(2:NYS, 1:N).*Gamma_rMAT(2:NYS, 1:N)./...
    (GMAT(2:NYS, 1:N).*Gamma_rMAT(2:NYS, 1:N) - k(2:NYS, 1:N)) - ...
    k(1:NYS - 1, 1:N).*Gamma_rMAT(2:NYS, 1:N)./(GMAT(2:NYS, 1:N).*...
    Gamma_rMAT(2:NYS, 1:N) - k(1:NYS - 1, 1:N));

% Slope of Stress Strain Curve for each yield surface
H(1:NYS - 1, :) = (k(2:NYS, :) - k(1:NYS - 1, :))./deltaGamma(2:NYS, :);

NYS = NYS*ones(N,1);
MatProfile = handles.MatProfileTEMP;
MatList = handles.MatList;
k1 = 0;
k2 = 0;

if handles.unitTEMP == 1
    unit1 = '(Pa)';
else
    unit1 = '(psi)';
end

for count = 1:length(MatProfile)
    
    if MatProfile(count)== 0
        
        index = Layers{count};
        index = index(1);
        
        Shear = k(1:NYS(index), index);
        
        Strain(1:NYS(index)) = k(1:NYS(index), index).*Gamma_rMAT(1:NYS(index), index)./...
            (GMAT(1:NYS(index), index).*Gamma_rMAT(1:NYS(index), index) - k(1:NYS(index), index));
        
        % Plots values of selected material
        h = figure;
        set(h, 'units', 'normalized', 'outerposition', [0.5, 0.4, 0.5, 0.6])
        set(h,'name','Strain vs. Shear','numbertitle','off')
        plot([0 Strain(1:NYS(index))], [0 Shear(1:NYS(index))'])
        xlabel('Strain (%)')
        ylabel(['Stress ' unit1])
        title(['Strain (%) vs. Stress ' unit1])
        
        % Changes name of figures
        set(h,'name',['Layers ', num2str(Layers{count})],'numbertitle','off')
        
        % Positions Figures
        set(h, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, 0.4, 0.6])
        
        i = i + 1;
        k2 = k2 + 0.05;
        count = count + 1;
        
        if count == 10
            k1 = k1 + 0.2;
            k2 = 0;
            count = 0;
            
            if k1 >= .6
                k1 = 0;
            end
        end
        
    else
        
        a = textread([handles.GUILOCATION, '\Material\', cell2mat(MatList(MatProfile(count),2)), '.txt'], '%s'); %#ok<RE%#ok<*MSNU> MFF1>
        
        count3 = str2num(strcat(cell2mat(a(2))));
        
        if isempty(count3)
            errordlg('Invalid format for material.')
            return
        end
        
        Strain = zeros(count3,1);
        Modulus = zeros(count3,1);
        
        % Assign values to appropriate vectors
        try
            for i = 1:count3
                Strain(i) = str2num(cell2mat(a(2*i+3)));
                Modulus(i) = str2num(cell2mat(a(2*i+4)));
            end
            
            count4 = str2num(strcat(cell2mat(a(6+count3*2))));
            Strain2 = zeros(count4,1);
            Damping = zeros(count4,1);
            for i = 1:count4
                Strain2(i) = str2num(cell2mat(a(2*i+7+count3*2)));
                Damping(i) = str2num(cell2mat(a(2*i+8+count3*2)));
            end
            
            Strain = Strain';
            Modulus = Modulus';
            Strain2 = Strain2';
            Damping = Damping';
            
            StrainOri = Strain;
            ModulusOri = Modulus;
            Strain2Ori = Strain2;
            DampingOri = Damping;
            
            Shear = Strain.*Modulus/100;
            ShearOri = Strain.*Modulus/100;
            ShearStrength = Shear(end);
            Slope = zeros(1, length(Shear));
            
        catch
            
            errordlg('Invalid format for material.')
            return
            
        end
        
        Gmax = G(Layers{count});
        Gmax = Gmax(1,1);
        
        Strain = Strain/100;
        
        Shear = Shear*Gmax;
        
        ShearStrength = su(Layers{count});
        ShearStrength = ShearStrength(1,1);
        
        Slope = zeros(1, length(Shear));
        Slope(1) = Shear(1)/Strain(1);
        for i = 1:length(Shear)-1
            Slope(end-i+1) = (Shear(end-i+1) - Shear(end-i))/(Strain(end-i+1) - Strain(end-i));
        end
        
        Adj1 = 0;
        Adj2 = 0;
        ShearTrial = Shear(end);
        i = 0;
        
        if ShearTrial > ShearStrength
            
            while Shear(end-i) > ShearStrength
                Shear(end-i) = ShearStrength;
                i = i + 1;
                Adj1 = 1;
            end
            
        elseif ShearStrength > ShearTrial
            
            while ShearStrength > ShearTrial
                
                i = i + 1;
                
                if i == length(Shear)
                    errordlg('Max Shear Strength too high.')
                    break
                end
                
                ShearTrial = Shear(end-i) + (Strain(end) - Strain(end-i))*Slope(end-i);
                Adj2 = 1;
                
            end
            
            if ShearTrial > ShearStrength && Adj2 == 1
                Shear(end) = ShearStrength;
            elseif ShearTrial <= ShearStrength && Adj2 == 1
                Shear(end) = ShearTrial;
            end
            
            for j = i-1:-1:1;
                Shear(end-j) = Shear(end-j-1) + (Strain(end-j) - Strain(end-j-1))*Slope(end-i);
                if Shear(end-j) > ShearStrength
                    Shear(end-j) = ShearStrength;
                end
            end
            
        end
        
        for j = 1:i
            Slope(end-j+1) = (Shear(end-j+1) - Shear(end-j))/(Strain(end-j+1) - Strain(end-j));
        end
        
        i = 1;
        while Slope(i) ~= 0 && i < count3
            i = i + 1;
        end
        
        if Slope(i) ~= 0
            i = i + 1;
            Slope(i) = 0;
        end
        
        NYS(Layers{count}) = i - 1;
        
        [m1 n1] = size(k);
        if i > m1
            k(m1+1:i,:) = zeros(i-m1,N);
            H(m1+1:i,:) = zeros(i-m1,N);
            deltaGamma(m1+1:i,:) = zeros(i-m1,N);
            alpha(m1+1:i,:) = zeros(i-m1,N);
        end
        
        k(1:i-1, Layers{count}) = repmat(Shear(1:i-1)',1,length(Layers{count}));
        
        H(1:i-1, Layers{count}) = repmat(Slope(2:i)',1,length(Layers{count}));
        
        deltaGamma(2:i-1, Layers{count}) = repmat(Strain(2:i-1)',1,length(Layers{count}))...
            -repmat(Strain(1:i-2)',1,length(Layers{count}));
        
        index = Layers{count};
        index = index(1);
        
        % Plots values of selected material
        h = figure;
        set(h, 'units', 'normalized', 'outerposition', [0.5, 0.4, 0.5, 0.6])
        set(h,'name','Strain vs. Shear','numbertitle','off')
        plot([0 Strain(1:NYS(index))], [0 Shear(1:NYS(index))])
        xlabel('Strain (%)')
        ylabel(['Stress ' unit1])
        title(['Strain (%) vs. Stress ' unit1])
        
        % Changes name of figures
        set(h,'name',['Layers ', num2str(Layers{count})],'numbertitle','off')
        
        % Positions Figures
        set(h, 'units', 'normalized', 'outerposition', [k1, 0.4 - k2, 0.4, 0.6])
        
        i = i + 1;
        k2 = k2 + 0.05;
        count = count + 1;
        
        if count == 10
            k1 = k1 + 0.2;
            k2 = 0;
            count = 0;
            
            if k1 >= .6
                k1 = 0;
            end
        end
        
    end
    
end


% --------------------------------------------------------------------
function ModeShapeTab_Callback(hObject, eventdata, handles)
% hObject    handle to ModeShapeTab (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)



% --- Executes when selected object is changed in Acceleration.
function Acceleration_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in Acceleration
% eventdata  structure with the following fields (see UIBUTTONGROUP)
%	EventName: string 'SelectionChanged' (read only)
%	OldValue: handle of the previously selected object or empty if none was selected
%	NewValue: handle of the currently selected object
% handles    structure with handles and user data (see GUIDATA)
handles.output = hObject;
switch get(eventdata.NewValue, 'Tag')
    case 'Frequency'
        AccelerationXLabel = 1;
    case 'Period'
        AccelerationXLabel = 2;
end
handles.AccelerationXLabelTEMP = AccelerationXLabel;
guidata(hObject, handles);


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

FileName = strrep(handles.FileNameTEMP,'.mat','');
handles.output = hObject;
N = handles.NTEMP;
h = handles.hTEMP;
Base = handles.BaseTEMP;
T = handles.T;
AC = handles.AC;
Disp = handles.Disp;
Vel = handles.Vel;
ACC = handles.ACC;
Strain = handles.Strain;
Stress = handles.Stress;

NodeOut = handles.NodeOutTEMP;
ElementOut = handles.EleOutTEMP;
DisplacementOut = handles.DisplacementOutTEMP;
VelocityOut = handles.VelocityOutTEMP;
AccelerationOut = handles.AccelerationOutTEMP;
StressOut = handles.StressOutTEMP;
StressOut2 = handles.StressOut2TEMP;
StrainOut = handles.StrainOutTEMP;
EarthquakeOut = handles.EarthquakeOutTEMP;

SpectralOut = handles.SpectralOutTEMP;
FourierOut = handles.FourierOutTEMP;
TransferOut = handles.TransferOutTEMP;
AccelerationXLabel = handles.AccelerationXLabelTEMP;

if Base == 1
    AC = ACC(end,:)';
end

mkdir([pwd,'\Saved Outputs\'],FileName)
[b]=uigetdir([pwd,'\Saved Outputs\',FileName]);

if b==0
else
    
    if sum(ElementOut) == 0 && sum(StrainOut + StressOut + StressOut2) ~= 0
        errordlg('Need to enter which elements you want to save the data for. Can be entered as a vector. (ex. 1:3 or 1 2 3)',...
            'Element # missing')
    elseif sum(NodeOut + ElementOut) == 0 && EarthquakeOut == 0
        errordlg('Need to enter node and/or elements you want to save the data for. Can be entered as a vector. (ex. 1:3 or 1 2 3)',...
            'No element or node selected.')
    elseif NodeOut > N
        errordlg('Invalid Node Selection NodeOut > N.', 'Invalid Node')
    elseif ElementOut > N
        errordlg('Invalid Element Selection. For Rigid Base there are only N number of elements, and for Flexible Base there are N - 1',...
            'Invalid element selection')
    elseif sum(NodeOut) == 0 && sum(DisplacementOut + VelocityOut + AccelerationOut + FourierOut + TransferOut) ~= 0
        errordlg('Need to enter which nodes you want a figure of. Can be entered as a vector. (ex. 1:3 or 1 2 3)',...
            'Node # missing')
    elseif sum(rem(NodeOut,1)) ~= 0
        errordlg('Node numbers should all be integers. (ex. 1:3 or 1 2 3)', 'Invalid Input for Node')
    elseif sum(rem(ElementOut,1))~=0
        errordlg('Element numbers should all be integers. (ex. 1:3 or 1 2 3)', 'Invalid Input for Element')
    else
        
        if handles.unitTEMP == 1
            g = 9.81;
        else
            g = 32.2*12;
        end
        
        if EarthquakeOut == 1
            msgbox('Saving Output for Earthquake', 'replace')
            %% Time History
            fid = fopen([b,'\Input_Time_History.txt'],'w');
            
            for i = 1:length(T)
                fprintf(fid,'%7e %7e\n',T(i),AC(i));
                [i T(i) AC(i)]
            end
            fclose(fid);
            
            
            %% Response Spectrum
            ug1 = AC;
            dt = T(2)-T(1);
            ug2 = interp1(dt:dt:length(AC)*dt,AC,dt:0.001:length(AC)*dt)';
            m = 1;
            dr = .05;
            gamma = 0.5;
            beta = 0;
            Tn1 = 0.02;
            Tn2 = 25;
            Tn = 10.^(log10(Tn1):(log10(Tn2)-log10(Tn1))/1000:log10(Tn2));
            wn = 2*pi./Tn;
            
            Sa = zeros(length(wn),1);
            for j = 1:length(wn)
                
                if wn(j) > 10
                    ug = ug2;
                    dt = 0.001;
                else
                    ug = ug1;
                    dt = T(2)-T(1);
                end
                c = 2 * wn(j) * m * dr; % Damping Coefficent
                k = m * wn(j)^2; % Stiffness Coefficent
                d = zeros(1,length(ug)); % Relative Displacement Vector (m)
                v = zeros(1,length(ug)); % Relative Velocity Vector (m/s)
                a = zeros(1,length(ug)); % Relative Acceleration Vector (m/s/s)
                r1 = zeros(1,length(ug)); % Restoring Force
                
                d(1) = 0; % Initial Displacement
                v(1) = 0; % Initial Velocity
                a(1) = (-m*ug(1)-c*v(1)-k*d(1))/m; % Initial Acceleration
                
                for i=1:length(ug)-1
                    d(i+1) = d(i) + dt*v(i) + dt^2*((1/2-beta)*a(i) + beta*a(i));
                    
                    f = -m*ug(i+1); % Force from ground motion
                    r1(i+1) = k * d(i+1); % Restoring Force
                    
                    a(i+1) = (f - r1(i+1) - c*(v(i) + dt*(1-gamma)*a(i)))/(m + dt*gamma*c);
                    v(i+1) = v(i) + dt*((1-gamma)*a(i) + gamma*a(i+1));
                end
                
                at = a' + ug; % Total Acceleration
                Sa(j) = max(abs(at)); % Spectral Acceleration
            end
            
            fid = fopen([b,'\Input_Spectral_Acceleration.txt'],'w');
            if AccelerationXLabel == 1
                for i = 1:length(Tn)
                    fprintf(fid,'%7e %7e\n',1/Tn(i),Sa(i)/g);
                end
            elseif AccelerationXLabel == 2
                for i = 1:length(Tn)
                    fprintf(fid,'%7e %7e\n',Tn(i),Sa(i)/g);
                end
            end
            fclose(fid);
            
            %% Fourier Analysis
            a1 = AC;
            a1 = a1.*hanning(length(a1));
            
            Fs = 1/dt;
            L = length(a1);
            
            NFFT = 2^nextpow2(L); % Next power of 2 from length of y
            f = Fs/2*linspace(0,1,NFFT/2+1);
            Y11 = fft(a1,NFFT)/L;
            
            y = (2*abs(Y11(1:NFFT/2+1)));
            
            for i = 1:10
                y(2:end-1) = 0.25*y(1:end-2) + 0.5*y(2:end-1) + 0.25*y(3:end);
            end
            
            fid = fopen([b,'\Input_Fourier_Spectrum.txt'],'w');
            if AccelerationXLabel == 1
                for i = 1:length(f)
                    fprintf(fid,'%7e %7e\n',f(i),y(i));
                end
            elseif AccelerationXLabel == 2
                for i = 1:length(f)
                    fprintf(fid,'%7e %7e\n',1/f(i),y(i));
                end
            end
            fclose(fid);
            
        end
        
        if NodeOut ~= 0
            
            if SpectralOut == 1
                ug1 = AC;
                dt = T(2)-T(1);
                ug2 = interp1(dt:dt:length(AC)*dt,AC,dt:0.005:length(AC)*dt)';
                m = 1;
                dr = .05;
                gamma = 0.5;
                beta = 0;
                Tn1 = 0.02;
                Tn2 = 25;
                Tn = 10.^(log10(Tn1):(log10(Tn2)-log10(Tn1))/1000:log10(Tn2));
                wn = 2*pi./Tn;
                
                Sa_AC = zeros(length(wn),1);
                for j = 1:length(wn)
                    
                    if wn(j) > 10
                        ug = ug2;
                        dt = 0.005;
                    else
                        ug = ug1;
                        dt = T(2)-T(1);
                    end
                    c = 2 * wn(j) * m * dr; % Damping Coefficent
                    k = m * wn(j)^2; % Stiffness Coefficent
                    d = zeros(1,length(ug)); % Relative Displacement Vector (m)
                    v = zeros(1,length(ug)); % Relative Velocity Vector (m/s)
                    a = zeros(1,length(ug)); % Relative Acceleration Vector (m/s/s)
                    r1 = zeros(1,length(ug)); % Restoring Force
                    
                    d(1) = 0; % Initial Displacement
                    v(1) = 0; % Initial Velocity
                    a(1) = (-m*ug(1)-c*v(1)-k*d(1))/m; % Initial Acceleration
                    
                    for i=1:length(ug)-1
                        d(i+1) = d(i) + dt*v(i) + dt^2*((1/2-beta)*a(i) + beta*a(i));
                        
                        f = -m*ug(i+1); % Force from ground motion
                        r1(i+1) = k * d(i+1); % Restoring Force
                        
                        a(i+1) = (f - r1(i+1) - c*(v(i) + dt*(1-gamma)*a(i)))/(m + dt*gamma*c);
                        v(i+1) = v(i) + dt*((1-gamma)*a(i) + gamma*a(i+1));
                    end
                    
                    at = a' + ug; % Total Acceleration
                    Sa_AC(j) = max(abs(at)); % Spectral Acceleration
                end
                
            end
            
            if FourierOut == 1
                %% Fourier Analysis
                a1 = AC;
                a1 = a1.*hanning(length(a1));
                
                L = length(a1);
                NFFT = 2^nextpow2(L); % Next power of 2 from length of y
                Y11 = fft(a1,NFFT)/L;
                
                y_AC = (2*abs(Y11(1:NFFT/2+1)));
                
                for i = 1:10
                    y_AC(2:end-1) = 0.25*y_AC(1:end-2) + 0.5*y_AC(2:end-1) + 0.25*y_AC(3:end);
                end
                
            end
            
            for n = NodeOut
                msgbox(['Saving Output for Node: ', num2str(n)], 'replace')
                if DisplacementOut == 1
                    fid = fopen([b,'\Node_',num2str(n),'_Displacement.txt'],'w');
                    for i = 1:length(T)
                        fprintf(fid,'%7e %7e\n',T(i),Disp(n, i));
                    end
                    fclose(fid);
                end
                
                if VelocityOut == 1
                    fid = fopen([b,'\Node_',num2str(n),'_Velocity.txt'],'w');
                    for i = 1:length(T)
                        fprintf(fid,'%7e %7e\n',T(i),Vel(n, i));
                    end
                    fclose(fid);
                end
                
                if AccelerationOut == 1
                    fid = fopen([b,'\Node_',num2str(n),'_Acceleration.txt'],'w');
                    for i = 1:length(T)
                        fprintf(fid,'%7e %7e\n',T(i),ACC(n, i));
                    end
                    fclose(fid);
                end
                
                if SpectralOut == 1
                    ug1 = ACC(n, :)';
                    dt = T(2)-T(1);
                    ug2 = interp1(dt:dt:length(ug1)*dt,ug1,dt:0.005:length(ug1)*dt)';
                    m = 1;
                    dr = .05;
                    gamma = 0.5;
                    beta = 0;
                    Tn1 = 0.02;
                    Tn2 = 25;
                    Tn = 10.^(log10(Tn1):(log10(Tn2)-log10(Tn1))/1000:log10(Tn2));
                    wn = 2*pi./Tn;
                    
                    Sa = zeros(length(wn),1);
                    for j = 1:length(wn)
                        
                        if wn(j) > 10
                            ug = ug2;
                            dt = 0.005;
                        else
                            ug = ug1;
                            dt = T(2)-T(1);
                        end
                        c = 2 * wn(j) * m * dr; % Damping Coefficent
                        k = m * wn(j)^2; % Stiffness Coefficent
                        d = zeros(1,length(ug)); % Relative Displacement Vector (m)
                        v = zeros(1,length(ug)); % Relative Velocity Vector (m/s)
                        a = zeros(1,length(ug)); % Relative Acceleration Vector (m/s/s)
                        r1 = zeros(1,length(ug)); % Restoring Force
                        
                        d(1) = 0; % Initial Displacement
                        v(1) = 0; % Initial Velocity
                        a(1) = (-m*ug(1)-c*v(1)-k*d(1))/m; % Initial Acceleration
                        
                        for i=1:length(ug)-1
                            d(i+1) = d(i) + dt*v(i) + dt^2*((1/2-beta)*a(i) + beta*a(i));
                            
                            f = -m*ug(i+1); % Force from ground motion
                            r1(i+1) = k * d(i+1); % Restoring Force
                            
                            a(i+1) = (f - r1(i+1) - c*(v(i) + dt*(1-gamma)*a(i)))/(m + dt*gamma*c);
                            v(i+1) = v(i) + dt*((1-gamma)*a(i) + gamma*a(i+1));
                        end
                        
                        at = a' + ug; % Total Acceleration
                        Sa(j,n) = max(abs(at)); % Spectral Acceleration
                    end
                    
                    fid = fopen([b,'\Node_',num2str(n),'_Spectral_Acceleration.txt'],'w');
                    if AccelerationXLabel == 1
                        for i = 1:length(Tn)
                            fprintf(fid,'%7e %7e\n',1/Tn(i),Sa(i)/g);
                        end
                    elseif AccelerationXLabel == 2
                        for i = 1:length(Tn)
                            fprintf(fid,'%7e %7e\n',Tn(i),Sa(i)/g);
                        end
                    end
                    fclose(fid);
                    
                end
                
                if FourierOut == 1
                    a1 = ACC(n, :)';
                    a1 = a1.*hanning(length(a1));
                    
                    Fs = 1/dt;
                    L = length(a1);
                    
                    NFFT = 2^nextpow2(L); % Next power of 2 from length of y
                    f = Fs/2*linspace(0,1,NFFT/2+1);
                    Y11 = fft(a1,NFFT)/L;
                    
                    y = (2*abs(Y11(1:NFFT/2+1)));
                    
                    for i = 1:10
                        y(2:end-1) = 0.25*y(1:end-2) + 0.5*y(2:end-1) + 0.25*y(3:end);
                    end
                    
                    fid = fopen([b,'\Node_',num2str(n),'_Fourier_Spectrum.txt'],'w');
                    if AccelerationXLabel == 1
                        for i = 1:length(f)
                            fprintf(fid,'%7e %7e\n',f(i),y(i));
                        end
                    elseif AccelerationXLabel == 2
                        for i = 1:length(f)
                            fprintf(fid,'%7e %7e\n',1/f(i),y(i));
                        end
                    end
                    fclose(fid);
                    
                end
                
                if TransferOut == 1
                    
                    base = AC;
                    a1 = ACC(n, :)';
                    a1 = a1.*hanning(length(a1));
                    base = base.*hanning(length(base));
                    
                    Fs = 1/dt;
                    L = length(a1);
                    
                    NFFT = 2^nextpow2(L); % Next power of 2 from length of y
                    
                    f = Fs/2*linspace(0,1,NFFT/2+1);
                    
                    Y11 = fft(a1,NFFT)/L;
                    Y12 = fft(base,NFFT)/L;
                    
                    y = (2*abs(Y11(1:NFFT/2+1)))./(2*abs(Y12(1:NFFT/2+1)));
                    
                    for i = 1:10
                        y(2:end-1) = 0.25*y(1:end-2) + 0.5*y(2:end-1) + 0.25*y(3:end);
                    end
                    
                    fid = fopen([b,'\Node_',num2str(n),'_Transfer_Function.txt'],'w');
                    if AccelerationXLabel == 1
                        for i = 1:length(f)
                            fprintf(fid,'%7e %7e\n',f(i),y(i));
                        end
                    elseif AccelerationXLabel == 2
                        for i = 1:length(f)
                            fprintf(fid,'%7e %7e\n',1/f(i),y(i));
                        end
                    end
                    fclose(fid);
                end
                
            end
            
        end
        
        if ElementOut ~= 0
            for n = ElementOut
                msgbox(['Saving Output for Element: ', num2str(n)], 'replace')
                if StrainOut == 1
                    fid = fopen([b,'\Node_',num2str(n),'_Strain.txt'],'w');
                    for i = 1:length(T)
                        fprintf(fid,'%7e %7e\n',T(i),Strain(n, i));
                    end
                    fclose(fid);
                end
                
                if StressOut == 1
                    fid = fopen([b,'\Node_',num2str(n),'_Stress.txt'],'w');
                    for i = 1:length(T)
                        fprintf(fid,'%7e %7e\n',T(i),Stress(n, i));
                    end
                    fclose(fid);
                end
                
                if StressOut2 == 1
                    fid = fopen([b,'\Node_',num2str(n),'_StressStrain.txt'],'w');
                    for i = 1:length(T)
                        fprintf(fid,'%7e %7e\n',Strain(n, i),Stress(n, i));
                    end
                    fclose(fid);
                end
            end
        end
        fclose('all');
        msgbox('Complete', 'replace')
    end
    
end
guidata(hObject, handles);


% --- Executes on button press in Reset.
function Reset_Callback(hObject, eventdata, handles)
% hObject    handle to Reset (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(handles.OutputFig, 'string', 'Output Figures')

set(handles.DisplacementOut, 'value', 0)
set(handles.VelocityOut, 'value', 0)
set(handles.AccelerationOut, 'value', 1)
set(handles.NodeOut, 'string', num2str(1))
set(handles.StrainOut, 'value', 0)
set(handles.StressOut, 'value', 0)
set(handles.StressOut2, 'value', 0)
set(handles.EleOut, 'string', num2str(0))
set(handles.SeparatedOut, 'value', 0)
set(handles.EarthquakeOut, 'value', 0)
set(handles.SpectralOut, 'value', 0)
set(handles.FourierOut, 'value', 0)
set(handles.TransferOut, 'value', 0)
set(handles.Frequency, 'value', 1)

handles.DisplacementOutTEMP = 0;
handles.VelocityOutTEMP = 0;
handles.AccelerationOutTEMP = 1;
handles.NodeOutTEMP = 1;
handles.StrainOutTEMP = 0;
handles.StressOutTEMP = 0;
handles.StressOut2TEMP = 0;
handles.EleOutTEMP = 0;
handles.SeparatedOutTEMP = 0;
handles.EarthquakeOutTEMP = 0;
handles.SpectralOutTEMP = 0;
handles.FourierOutTEMP = 0;
handles.TransferOutTEMP = 0;
handles.AccelerationXLabelTEMP = 1;
guidata(hObject, handles);
