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

% Last Modified by GUIDE v2.5 07-Sep-2012 09:54:36

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

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

% Update handles structure
guidata(hObject, handles);

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


% --- Outputs from this function are returned to the command line.
function varargout = advection_gui_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 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)
M = get_gridsize(handles);
f = get_function(handles);
method = get_method(handles);
mu = get_mu(handles);
N = get_steps(handles);

[u v] = advection(f, M, mu, method, N);

u(:,1);
p = plot(handles.axes1,1/(M+1)*(1:M), u(:,1),1/(M+1)*(1:M), v(:,1));
show_latex(M, mu, method);
axis ([0 1 -3 3])
t = 2;
while t < N && max(abs(u(:,t-1))) <= 3 && ishandle(p(1)) && ishandle(p(2))
    set(p(1), 'YData', u(:,t));
    set(p(2), 'YData', v(:,t));
    t = t+1;
    drawnow;
end

% --- 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


% --- 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


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

% --- 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 axes1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to axes1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes1


% --- 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);

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

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



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

function method = get_method(handles)
methodcontents = cellstr(get(handles.popupmenu1,'String'));
switch methodcontents{get(handles.popupmenu1,'Value')}
    case 'Upwind'
        method = 'upwind';
    case 'Downwind'
        method = 'downwind';
    case 'Lax-Wendroff'
        method = 'lax-wendroff';
    case 'Leapfrog'
        method = 'leapfrog';
    case 'Angled Derivative'
        method = 'angled-derivative';
    otherwise
        method = 'upwind';
end

function f = get_function(handles)
f = str2func(['@(x) double(' get(handles.edit1,'String') ')']);
arrayfun(f,0:0.1:1);

function mu = get_mu(handles)
mu = str2double(get(handles.edit2, '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.edit3,'String'));
assert(~isnan(M) && isreal(M) && ~isinf(M) && floor(M)==M...
    && and(M >= 1,M<=500),...
    'Time mesh must be an integer between 1 and 500.');

function N = get_steps(handles)
N = str2double(get(handles.edit4, 'String'));
assert(~isnan(N) && isreal(N) && ~isinf(N) && floor(N)==N...
    && and(N >= 1,N<=500),...
    'Space mesh must be an integer between 1 and 500.');

function update_grid(handles)
try
    M = get_gridsize(handles);
    N = get_steps(handles);
    mu = get_mu(handles);
    f = get_function(handles);
    method = get_method(handles);
    
    
    set(handles.pushbutton1,'Enable','on');

    [u v] = advection(f, M, mu, 'upwind', 3);

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

    show_latex(M, mu, method);
    % set(handles.axes1, 'YData', u(:,1));
    axis ([0 1 -3 3]);
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 show_latex(M, mu, method)
delete(findall(allchild(gca),'Type','Text'))
text('Units','Normalized','Position',[0.5 -0.1], 'String', '$$\mu=$$', 'Interpreter', 'latex');
text('Units','Normalized','Position',[0.5 -0.15], 'String', '$$M=$$', 'Interpreter', 'latex');
text('Units','Normalized','Position',[0.5 -0.2], 'String', '$$N=$$', 'Interpreter', 'latex');
space = '\qquad \; \;';
switch method
    case 'upwind'
        l1 = sprintf('$$ %0.2f \\ %0.2f $$',1-mu,mu);
        l2 = sprintf('$$ 0.00 \\ %2$0.2f \\ %3$0.2f $$',space,1-mu,mu);
        l3 = sprintf('$$ %1$s 0.00 \\ %2$0.2f \\ %3$0.2f $$',space, 1-mu,mu);
        l4 = sprintf('$$ %1$s %1$s \\ddots \\ \\ddots \\ \\ddots $$', space);
        
        spectral_radius = abs(1-mu);
        show_matrix({l1,l2,l3,l4},[1.01,0.8],spectral_radius)

    case 'downwind'
        l1 = sprintf('$$ \\ \\; \\; %0.2f \\ 0.00 $$', 1+mu);
        l2 = sprintf('$$ %1$0.2f \\ %2$0.2f \\ 0.00 $$',-mu,1+mu);
        l3 = sprintf('$$ %1$s %2$0.2f \\ %3$0.2f \\ 0.00 $$',space, -mu,1+mu);
        l4 = sprintf('$$ \\ \\; %1$s %1$s \\ddots \\ \\ \\ddots \\ \\ \\ddots $$', space);

        spectral_radius = abs(1+mu);
        show_matrix({l1,l2,l3,l4},[1.01,0.8],spectral_radius)
        
    case 'lax-wendroff'
        l1 = sprintf('$$ %0.2f \\ %0.2f $$', 1-mu^2, 0.5*mu*(1+mu));
        l2 = sprintf('$$ %0.2f \\ %0.2f \\ %0.2f $$',-0.5*mu*(1-mu),1-mu^2, 0.5*mu*(1+mu));
        l3 = sprintf('$$ %s %0.2f \\ %0.2f \\ %0.2f $$',space, -0.5*mu*(1-mu),1-mu^2, 0.5*mu*(1+mu));
        l4 = sprintf('$$ \\ \\; %1$s %1$s \\ddots \\ \\ \\ddots \\ \\ \\ddots $$', space);
               
        spectral_radius = abs(1-mu^2+1i*sqrt(mu^2*(1-mu^2)) * cos(pi/M+1));
        show_matrix({l1,l2,l3,l4},[1.01,0.8],spectral_radius);
end

function show_matrix(L, pos, spectral_radius)
% For some reason positioning doesn't work as expected. 
% This is a hack:
% pos0 = desired position
% pos1 = actual position
% corr = correction (actual minus desired) to correct future positions)
pos0=pos;
h=text('String', 'A = ', ...
'Units','Normalized','Position', pos0, 'HorizontalAlignment', 'Left', ...
    'Interpreter', 'latex', 'fontsize', 15);
pos1=get(h,'Extent');
corr=pos1(1:2)-pos0;
pos0=[pos1(1)+pos1(3),pos1(2)-1.5*pos1(4)]-corr;
h=text('String', '[', ...
    'Units','Normalized','Position', pos0, 'VerticalAlignment','Bottom', 'HorizontalAlignment', 'Left', ...
    'interpreter', 'latex', 'fontsize', 60);
pos1=get(h,'Extent');
corr=pos1(1:2)-pos0;
pos0=[pos1(1)+0.5*pos1(3),pos1(2)]-corr;
h=text('String', L, ...
    'Units','Normalized','Position', pos0, 'HorizontalAlignment','Left','VerticalAlignment','Bottom', ...
    'Interpreter', 'latex');
pos1=get(h,'Extent');
corr=pos1(1:2)-pos0;
pos0=[pos1(1)+pos1(3),pos1(2)]-corr;
h=text('String', ']', ...
    'Units','Normalized','Position', pos0, 'VerticalAlignment','Bottom', 'HorizontalAlignment', 'Left', ...
    'interpreter', 'latex', 'fontsize', 60);
pos1=get(h,'Extent');
corr=pos1(1:2)-pos0;
pos0=[1.05,pos1(2)-0.5*pos1(4)]-corr;
spectral_radius_text = sprintf('$$ \\rho(A) = %0.2f $$', spectral_radius);
text('String', spectral_radius_text,...
    'Units','Normalized','Position', pos0, 'HorizontalAlignment', 'Left', ...
     'interpreter', 'latex', 'fontsize', 15); 



function A = create_matrix(M, mu, method)
e = ones(M,1);
switch method
    case 'upwind'
        A = spdiags( [(1-mu)*e mu*e], 0:1, M, M);
    case 'downwind'
        A = spdiags( [-mu*e (1+mu)*e ], -1:0, M, M);
    case 'lax-wendroff'
        A = spdiags( [ -0.5*mu*(1-mu)*e (1-mu^2)*e 0.5*mu*(1+mu)*e], -1:1, M, M);
    case 'leapfrog'
        A = spdiags( [-mu*e 0*e mu*e], -1:1, M, M);
    case 'angled-derivative'
        A = spdiags( [(1-2*mu)*e -(1-2*mu)*e], 0:1, M, M);
end

function [u v] = advection (f, M, mu, method, N)
h = 1/(M+1);
dt = mu*h;

u = zeros (M, N+2);
v = zeros (M, N+2);

u(1:M,1) = arrayfun(f,h*(1:M));

T = floor(M/((M+1)*dt)+1);

for t = 1:T
    i = min(M,(1-dt*(t-1))*(M+1));
    i = floor(i);
    v(1:i,t) = arrayfun(f,h*(1:i)+dt*(t-1));
end


A = create_matrix(M, mu, method);

switch method
    case 'leapfrog'
        u(1:M,2) = arrayfun(f,h*(1:M)+dt);
    case 'angled-derivative'
        u(1:M,2) = arrayfun(f,h*(1:M)+dt);
    otherwise
        u(:,2) = A*u(:,1);
end

for t = 3:N
    switch method
        case 'leapfrog'
            u(:,t) = A*u(:,t-1)+u(:,t-2);
            u(1,t) = u(1,t) - mu*u(1,t-2);
        case 'angled-derivative'
            z = zeros(M,1);
            z(1:M-1) = u(2:M,t-2);
            u(:,t) = A*u(:,t-1)+z;
%             u(1,t) = u(1,t) - mu*u(1,t-2);
        otherwise
            u(:,t) = A*u(:,t-1);
    end
end


% --- Executes when figure1 is resized.
function figure1_ResizeFcn(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)
if(isstruct(handles))
    show_latex(get_gridsize(handles),get_mu(handles),get_method(handles));
end
