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

% Last Modified by GUIDE v2.5 07-Sep-2012 16:23:09

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @stability_OpeningFcn, ...
                   'gui_OutputFcn',  @stability_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 stability is made visible.
function stability_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to stability (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = stability_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on selection change in popupmenu1.
function popupmenu1_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (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 popupmenu1 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu1
update_grid(handles);



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

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



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

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (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 pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
N = 1000;
M = get_gridsize(handles);
mu = get_courant(handles);
f = get_function(handles);
method = get_method(handles);
g = get_g(handles);
c = get_coef(handles);
u = finite_differences(f, M, mu, c, g, method, N);

p = plot(handles.axes1,1/(M+1)*(0:M+1), u(:,1));

show_steps(mu, M, method);

axis ([0 1 -3 3])
t = 2;
while t <= N && max(abs(u(:,t-1))) <= 3 && ishandle(p)
    set(p, 'YData', u(:,t));
    t = t+1;
    drawnow;
end
    

% --- Executes on selection change in popupmenu2.
function popupmenu2_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu2 (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 popupmenu2 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu2
update_grid(handles);

        

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

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


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

% --- Executes during object creation, after setting all properties.
function edit2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit2 (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 selection change in popupmenu3.
function popupmenu3_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu3 (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 popupmenu3 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu3
update_grid(handles);

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

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



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

% --- Executes during object creation, after setting all properties.
function edit3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit3 (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 edit4_Callback(hObject, eventdata, handles)
% hObject    handle to edit4 (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 edit4 as text
%        str2double(get(hObject,'String')) returns contents of edit4 as a double
update_grid(handles);

% --- Executes during object creation, after setting all properties.
function edit4_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit4 (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 text6_CreateFcn(hObject, eventdata, handles)
% hObject    handle to text6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called


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

% --- Executes during object creation, after setting all properties.
function edit5_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit5 (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 update_grid(handles)
try
    f = get_function(handles);
    M = get_gridsize(handles);
    mu = get_courant(handles);
    method = get_method(handles);
    g = get_g(handles);
    c = get_coef(handles);
    u = finite_differences(f, M, mu, c, g, method, 1);

    p = plot(handles.axes1,1/(M+1)*(0:M+1), u(:,1));
    axis([0 1 -3 3])
    show_steps(mu, M, method);
    set(handles.pushbutton1,'Enable','on');
catch exception
     if (strcmp(exception.identifier,'MATLAB:inlineeval:InlineExprError')) || (strcmp(exception.identifier,'MATLAB:arrayfun:NotAScalarOutput'))
         waitfor(msgbox('The function you have entered cannot be understood by MATLAB(R)', 'Error', 'error'));
     elseif strcmp(exception.identifier, 'MATLAB:m_missing_operator')
         waitfor(msgbox('The function you have typed in appears to me missing an operator', 'Error', 'error'));
     else
         waitfor(msgbox(exception.message,'Error','error'));
     end
     set(handles.pushbutton1,'Enable','off');
end

function mu = get_courant(handles)
mu = str2double(get(handles.edit1, 'String'));   
assert(~isnan(mu) && isreal(mu) && ~isinf(mu) && mu > 0,...
    'The Courant number needs to be positive!');

function M = get_gridsize(handles)
M = str2double(get(handles.edit2,'String'));
assert(~isnan(M) && isreal(M) && ~isinf(M) && floor(M)==M...
    && and(M >= 1,M<=100),...
    'Mesh size must be an integer between 1 and 100.');

function method = get_method(handles)
methodcontents = cellstr(get(handles.popupmenu1,'String'));
switch methodcontents{get(handles.popupmenu1,'Value')};
    case 'Wave Equation'
        method = 'wave';
    case 'Diffusion Equation (Euler Method)'
        method = 'diffusion';
    case 'Diffusion Equation (Crank Nicolson)'
        method = 'crank-nicolson';
end

function f = get_function(handles)
f = inline(get(handles.edit5,'String'), 'x');
f = @(x)(f(x));
f = @(x)(arrayfun(f,x));
f(0:0.1:1);

function g = get_g(handles)
g = str2num(get(handles.edit4,'String'));
assert(ndims(g)==2 && size(g,1)==1 && size(g,2)==2 && ...
    all(isreal(g)) && all(~isinf(g)),...
    '[g(0) g(1)] must be a real row vector of length 2');

function c = get_coef(handles)
c = eval(get(handles.edit3,'String'));
assert(~any(isnan(c))&&  ndims(c)==2 && size(c,1)==1 && size(c,2)==2 && ...
    isreal(c) && ~any(isinf(c))&& any(c~=0),...
    '[a b] must be a real, non-zero row vector of length 2')

function show_steps(mu, M, method)
text('Units','Normalized','Position', [0.19, -0.18], 'VerticalAlignment','Top','HorizontalAlignment','Right','String', '$$au+b \frac{\partial u}{\partial n} =  g~;\quad[a~b]= $$', 'Interpreter', 'latex');
text('Units','Normalized','Position', [0.5, -0.207],'HorizontalAlignment','Right','String', '$$[g(0), g(1)] = $$', 'Interpreter', 'latex');
switch method
    case 'wave'
        dt = sqrt(mu)/(M+1);
    otherwise
        dt = mu/(M+1)^2;
end
[a1 a2] = sci_note(1/(M+1));
[a3 a4] = sci_note(dt);
t1 = sprintf('$$ \\Delta x = %0.4f \\times 10^{ %i }$$', a1, a2);
t2 = sprintf('$$ \\Delta t = %0.4f \\times 10^{ %i }$$', a3, a4);
text('Units','Normalized','Position', [0.4 -0.3], 'VerticalAlignment','Top','HorizontalAlignment','Left','String', {t1,t2}, 'Interpreter', 'latex');

function u = finite_differences ( f, M, mu, coef, g, method, N )
h = 1/(M+1);

u = zeros (M+2, N+2);
u(1,1) = 0;
u(2:(M+1),1) = f(h*(1:M));
u(M+2,1) = 0;

e = ones(M+1,1);
   

switch method

    case 'wave'
        Ain = spdiags( [mu*e (2-2*mu)*e mu*e], 0:2, M, M+2);
        b = zeros(M+2,1);
        S = speye(M+2);
        if coef(2) == 0
            l = zeros(1,M+2);
            r = zeros(1,M+2);
            l(1) = 1;
            r(M+2) = 1;
            u(1,1) = g(1)/coef(1);
            u(M+2,1) = g(2)/coef(1);
            S(1,1) = 0;
            S(M+2,M+2) = 0;
        else
           l = zeros(1,M+2);
           r = zeros(1,M+2);
           l(1) = 2-2*mu-2*mu*h*coef(1)/coef(2);
           l(2) = 2*mu;
           r(M+2) = 2-2*mu+2*mu*h*coef(1)/coef(2);
           r(M+1) = 2*mu;
           b(1) = 2*mu*h*g(1)/coef(2);
           b(M+2) = -2*mu*h*g(2)/coef(2); 
        end
        A = [l; Ain; r];
        u(:,2) = [l; 0.5*Ain; r]*u(:,1);
        for n = 3:N
            u(:,n) = A*u(:,n-1)-S*u(:,n-2)+b;
        end
    case 'diffusion'
        Ain = spdiags( [mu*e (1-2*mu)*e mu*e], 0:2, M, M+2);
        l = zeros(1,M+2);
        r = zeros(1,M+2);
        b = zeros(M+2,1);
        if coef(2) == 0
            l(1) = 0;
            r(M+2) = 0;
            b(1) = g(1)/coef(1);
            b(M+2) = g(2)/coef(1);
        else   
            l(1) = 1-2*mu-mu*2*h*coef(1)/coef(2);
            l(2) = 2*mu;
            r(M+2) = 1-2*mu+mu*2*h*coef(1)/coef(2);
            r(M+1) = 2*mu;
            b(1) = 2*mu*h*g(1)/coef(2);
            b(M+2) = 2*mu*h*g(2)/coef(2);
        end
        A = [l; Ain; r];
        for n = 2:N
            u(:,n) = A*u(:,n-1)+b;
        end
    case 'crank-nicolson'
        Bin = spdiags( [-0.5*mu*e (1+mu)*e -0.5*mu*e], 0:2, M, M+2);
        Cin = spdiags( [0.5*mu*e (1-mu)*e 0.5*mu*e], 0:2, M, M+2);
        l1 = zeros(1,M+2);
        l2 = l1;
        r1 = zeros(1,M+2);
        r2 = r1;
        b = zeros(M+2,1);
        if coef(2) == 0
            l1(1) = 1;
            l2 = l1;
            r1(M+2) = 1;
            r2 = r1;
            u(1,1) = g(1)/coef(1);
            u(M+2,1) = g(2)/coef(1);
        else
            l1(1) = 1-0.5*mu*(-2+2*h*coef(1)/coef(2));
            l1(2) = -mu;
            l2(1) = 1+0.5*mu*(-2+2*h*coef(1)/coef(2));
            l2(2) = mu;
            
            r1(M+2) = 1-0.5*mu*(-2-2*h*coef(1)/coef(2));
            r1(M+1) = -mu;
            r2(M+2) = 1+0.5*mu*(-2-2*h*coef(1)/coef(2));
            r2(M+1) = mu;
            
            b(1) = -2*mu*h*g(1)/coef(2);
            b(M+2) = 2*mu*h*g(2)/coef(2);
        end
        B = [l1; Bin; r1];
        C = [l2; Cin; r2];
        for n = 2:N
            u(:,n) = B\(C*u(:,n-1)+b);
        end
end

function [n, e] = sci_note (x)
e = floor((log(x)/log(10)));
n = x/10^e;
