function varargout = nma_advection_pde_1D(varargin)
% implement HW3, Math 228B, advection ODE solver

% by Nasser M. Abbasi
% UC Davis, winter 2011
%
% NMA_ADVECTION_PDE_1D M-file for nma_advection_pde_1D.fig
%      NMA_ADVECTION_PDE_1D, by itself, creates a new NMA_ADVECTION_PDE_1D or raises the existing
%      singleton*.
%
%      H = NMA_ADVECTION_PDE_1D returns the handle to a new NMA_ADVECTION_PDE_1D or the handle to
%      the existing singleton*.
%
%      NMA_ADVECTION_PDE_1D('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in NMA_ADVECTION_PDE_1D.M with the given input arguments.
%
%      NMA_ADVECTION_PDE_1D('Property','Value',...) creates a new NMA_ADVECTION_PDE_1D or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before nma_advection_pde_1D_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to nma_advection_pde_1D_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 nma_advection_pde_1D

% Last Modified by GUIDE v2.5 11-Feb-2012 03:12:50

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

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

set(handles.axes_2_1,'Visible','off');
set(handles.axes_2_2,'Visible','off');

set(handles.axes_3_1,'Visible','off');
set(handles.axes_3_2,'Visible','off');
set(handles.axes_3_3,'Visible','off');

set(handles.axes_4_1,'Visible','off');
set(handles.axes_4_2,'Visible','off');
set(handles.axes_4_3,'Visible','off');
set(handles.axes_4_4,'Visible','off');

set(handles.axes_5_1,'Visible','off');
set(handles.axes_5_2,'Visible','off');
set(handles.axes_5_3,'Visible','off');
set(handles.axes_5_4,'Visible','off');
set(handles.axes_5_5,'Visible','off');
set(handles.axes_1_1,'Visible','on');

data.state = 1;
set(handles.figure1, 'UserData',data);
set(handles.figure1,'Name','1D advection solver, By Nasser M. Abbasi, Version: feb 12,2012');

nma_set_figure_position(handles.figure1,0.15,0.1,0.7,0.78);

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = nma_advection_pde_1D_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 lax_wendroff_tag.
function lax_wendroff_tag_Callback(hObject, eventdata, handles)
% hObject    handle to lax_wendroff_tag (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 lax_wendroff_tag

% --- Executes on button press in upwind_tag.
function upwind_tag_Callback(hObject, eventdata, handles)
% hObject    handle to upwind_tag (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 upwind_tag

% --- Executes on button press in lax_fried_tag.
function lax_fried_tag_Callback(hObject, eventdata, handles)
% hObject    handle to lax_fried_tag (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 lax_fried_tag

% --- Executes on button press in leap_frog_tag.
function leap_frog_tag_Callback(hObject, eventdata, handles)
% hObject    handle to leap_frog_tag (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 leap_frog_tag

% --- Executes on button press in beam_warming_tag.
function beam_warming_tag_Callback(hObject, eventdata, handles)
% hObject    handle to beam_warming_tag (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 beam_warming_tag

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

% Hints: get(hObject,'String') returns contents of L_tag as text
%        str2double(get(hObject,'String')) returns contents of L_tag as a
%        double

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

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

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

% Hints: get(hObject,'String') returns contents of h_tag as text
%        str2double(get(hObject,'String')) returns contents of h_tag as a double

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

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

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

% Hints: get(hObject,'String') returns contents of a_tag as text
%        str2double(get(hObject,'String')) returns contents of a_tag as a double

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

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

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

% Hints: get(hObject,'String') returns contents of CFL_tag as text
%        str2double(get(hObject,'String')) returns contents of CFL_tag as a double

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

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

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

% Hints: get(hObject,'String') returns contents of user_function_btn_tag as text
%        str2double(get(hObject,'String')) returns contents of user_function_btn_tag as a double

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

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

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

% Hints: get(hObject,'String') returns contents of rect_width_tag as text
%        str2double(get(hObject,'String')) returns contents of rect_width_tag as a double

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

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

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

% Hints: get(hObject,'String') returns contents of rect_start_at_tag as text
%        str2double(get(hObject,'String')) returns contents of rect_start_at_tag as a double

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

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

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

% check refinement
data.doing_refinement=false;
data.num_refinement_study_runs=0;

%find method to run
LAX_WENDROFF=1; UPWIND=2; LAX_F=3; LEAPFROG=4; BEAM_WARMING=5; CN=6;
FTCS=7; LAX=8;

data.methods_mask = zeros(1,8);  %only 8 methods are supported now

if get(handles.lax_wendroff_tag,'Value')==1
   data.methods_mask(LAX_WENDROFF)=1;
end
if get(handles.upwind_tag,'Value')==1
   data.methods_mask(UPWIND)=1;
end
if get(handles.lax_fried_tag,'Value')==1
   data.methods_mask(LAX_F)=1;
end
if get(handles.leap_frog_tag,'Value')==1
   data.methods_mask(LEAPFROG)=1;
end
if get(handles.beam_warming_tag,'Value')==1
   data.methods_mask(BEAM_WARMING)=1;
end

if get(handles.CN_tag,'Value')==1
   data.methods_mask(CN)=1;
end

if get(handles.FTCS_tag,'Value')==1
   data.methods_mask(FTCS)=1;
end

if get(handles.LAX_tag,'Value')==1
   data.methods_mask(LAX)=1;
end

data.number_of_methods = sum(data.methods_mask);
if data.number_of_methods == 0
   uiwait(errordlg('Must select at least one method','Bad Input', 'modal'));
   uicontrol(handles.lax_wendroff_tag);
   status = false;
   return;
end

if data.number_of_methods > 1 && get(handles.ref_tag,'Value')==1
   uiwait(errordlg('Only one method can be used when refinement study is selected','Bad Input', 'modal'));
   uicontrol(handles.ref_tag);
   status = false;
   return;
end

if data.number_of_methods >5
   uiwait(errordlg('Limit exceeded. Maximum of 5 methods can be compared at the same time.','Bad Input', 'modal'));
   uicontrol(handles.lax_wendroff_tag);
   status = false;
   return;
end

if get(handles.ref_tag,'Value')==1
   str = get(handles.ref_number_runs_tag,'String');
   [data.num_refinement_study_runs, status] = verify_valid_positive_integer(str,handles.run_time_tag, 'number of refinement study runs');
   if not(status)
      return;
   end
   
   if data.num_refinement_study_runs>11
      uiwait(errordlg('number of runs for refinement study too large. Limit 11','Bad Input', 'modal'));
      uicontrol(handles.ref_number_runs_tag);
      status = false;
      return;
   end
   
   data.doing_refinement = true;
else
   data.doing_refinement = false;
end
% end check for refienemtn

str = get(handles.L_tag,'String');
[data.L, status] = nma_verify_valid_positive_numeric(str,handles.L_tag, 'length');
if not(status)
   return;
end

str = get(handles.h_tag,'String');
[data.h, status] = nma_verify_valid_positive_numeric(str,handles.h_tag, 'grid spacing');
if not(status)
   return;
end

if data.h>data.L
   uiwait(errordlg('grid spacing too large','Bad Input', 'modal'));
   uicontrol(handles.h_tag);
   status = false;
   return;
end

str = get(handles.a_tag,'String');
[data.a, status] = nma_verify_valid_positive_numeric(str,handles.h_tag, 'speed');
if not(status)
   return;
end

str = get(handles.CFL_tag,'String');
[data.CFL, status] = nma_verify_valid_positive_numeric(str,handles.CFL_tag, 'Courant number');
if not(status)
   return;
end

% if data.CFL>1
%    uiwait(errordlg('Courant number must be <= 1','Bad Input', 'modal'));
%    uicontrol(handles.CFL_tag);
%    status = false;
%    return;
% end

if get(handles.max_time_btn_tag,'Value')==1
   str = get(handles.run_time_tag,'String');
   [data.max_t, status] = nma_verify_valid_positive_numeric(str,handles.run_time_tag, 'run time');
   if not(status)
      return;
   end
else
   str = get(handles.number_steps_value_tag,'String');
   [data.nSteps, status] = verify_valid_positive_integer(str,handles.run_time_tag, 'number of time steps');
   if not(status)
      return;
   end
   delt = data.CFL*data.h/data.a;
   data.max_t = data.nSteps * delt;
end

data.X=0:data.h:data.L;


%find what initial data to use
data.x_start = 0;

if get(handles.user_function_btn_tag,'Value')==1
   ic= strtrim(get(handles.user_function_tag,'String'));
   if strcmp(ic, '0') || strcmp(ic, '')
      uiwait(errordlg('Missing input data function','Bad Input', 'modal'));
      uicontrol(handles.user_function_tag);
      status = false;
      return;
   end
   data.ic = ['@(X)' ic];
   
   try
      data.ic = str2func(data.ic);
      data.ic(.5);
   catch ME
      uiwait(errordlg('Invalid initial function, check for valid Matlab synatx...',...
         'Bad Input', 'modal'));
      uicontrol(handles.user_function_btn_tag);
      status = false;
      return
   end
   
   data.u_init = data.ic(data.X);
   
   %for plotting
   if abs(max(data.u_init)-0)<= eps
      data.max_y = 0.1;
   else
      data.max_y = max(data.u_init)+0.3*abs(max(data.u_init));
   end
   
   if abs(min(data.u_init)-0)<= eps
      data.min_y = -0.2;
   else
      data.min_y = min(data.u_init)-0.4*abs(max(data.u_init));
   end
elseif get(handles.rect_tag,'Value')==1
   str = get(handles.rect_width_tag,'String');
   [data.pulse_width, status] = nma_verify_valid_positive_numeric(str,handles.rect_width_tag, 'rectangular pulse width');
   if not(status)
      return;
   end
   
   if data.pulse_width>data.L
      uiwait(errordlg('rectangular pulse can not be wider than length...',...
         'Bad Input', 'modal'));
      uicontrol(handles.rect_width_tag);
      status = false;
      return
   end
   
   str = get(handles.rect_start_at_tag,'String');
   [data.x_start , status] = nma_verify_valid_positive_numeric(str,handles.rect_start_at_tag, 'rectangular start location');
   if not(status)
      return;
   end
   
   if data.x_start>data.L || data.x_start<0
      uiwait(errordlg('invalid staring position for rectangular pulse',...
         'Bad Input', 'modal'));
      uicontrol(handles.rect_start_at_tag);
      status = false;
      return
   end
   
   
   str = get(handles.rect_height_tag,'String');
   [data.pulse_height , status] = nma_verify_valid_positive_numeric(str,handles.rect_height_tag, 'rectangular height');
   if not(status)
      return;
   end
   
   data.pulse  = nma_rect_pulse_on_periodic_1D(data.pulse_width ,data.pulse_height,data.L,data.h);
   data.u_init = data.pulse.get(data.x_start);
   %for plotting
   data.max_y = 1.4*max(data.u_init);
   data.min_y = -0.5*max(data.u_init);
end
status = true;

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

disable_all(handles);
[status,data] = verify_all_input(handles);
if not(status)
   enable_all(handles);
   return;
end

if data.doing_refinement==true
   %
   % Do refinement study. Use refinement_study object
   %
   R = nma_refinement_study_manager(data.num_refinement_study_runs,1);
   current_refinement_run = 0;
   cla(handles.axes_3_3,'reset');
   cla(handles.axes_3_2,'reset');
   
   h=arrayfun(@(x) (1/2)^x, 1:10);  % Divide h by 2 each time
   
   for j = 1:data.num_refinement_study_runs
      set(handles.h_tag,'String',num2str(h(j),16));
      
      [status,data] = verify_all_input(handles); %do this again, to reload
      %data affected by changing h
      if status
         current_refinement_run         = current_refinement_run+1;
         data.current_refinement_run    =  current_refinement_run;
         [numerical_solution,k,u_exact] = run0(handles,data,h(j)*data.L);
         
         R.add_exact(numerical_solution(data.methods_mask==1,:), ...
            u_exact,h(j)*data.L,k);
         
         [hd,bdy] = R.format_table();
         cla(handles.axes_3_2,'reset');
         set(gcf,'CurrentAxes',handles.axes_3_2);
         
         text(-0.03,.40,bdy,'FontSize',7);
         set(handles.axes_3_2,'YTick',[]);
         set(handles.axes_3_2,'XTick',[]);
         text(-.03,.9,hd,'FontSize',8);
         title('result of refinement study');
         set(gca,'FontSize',7);
         drawnow();
         
         if current_refinement_run>= 4
            tbl  = R.get_table();
            set(0,'CurrentFigure',handles.figure1);
            cla(handles.axes_3_3,'reset');
            set(handles.figure1, 'CurrentAxes',handles.axes_3_3);
            set(0,'defaultaxesfontsize',8) ;
            loglog(tbl(4:end,3),tbl(4:end,4),'-d');
            xlabel('log(h)','FontSize',7);
            title('refinement study result, ','FontSize',7);
            ylabel('log(error norm)','FontSize',7);
            grid on;
            set(gca,'FontSize',7);
            drawnow();
         end
      end
   end
else
   disable_all(handles);
   run0(handles,data,data.h);
   enable_all(handles);
end

%---------------------------------
function [u,unew,u_n_minus_1,u_n]=...
   make_step(data,n,handles,u,unew,...
   A,B,CN_A,CN_B,u_n_minus_1,u_n,k,axes_handles,current_time)

LAX_WENDROFF=1; UPWIND=2; LAX_F=3; LEAPFROG=4; BEAM_WARMING=5; CN=6;
FTCS=7; LAX=8;



set(0,'CurrentFigure',handles.figure1);

if data.methods_mask(LAX_WENDROFF)==1
   i = 2:length(data.u_init)-1;
   U = u(LAX_WENDROFF,:);
   unew(LAX_WENDROFF,1)   = U(1)-(A/2)*(U(2)-U(end-1))+ ...
      B*(U(end-1)-2*U(1)+U(2));
   
   unew(LAX_WENDROFF,end) = U(end)-(A/2)*(U(2)-U(end-1))+...
      B*(U(end-1)-2*U(end)+U(2));
   
   unew(LAX_WENDROFF,i)   = U(i)-(A/2)*(U(i+1)-U(i-1))+...
      B*(U(i-1)-2*U(i)+U(i+1));
end

if data.methods_mask(UPWIND)==1
   i = 2:length(data.u_init)-1;
   if data.a>= 0
      U = u(UPWIND,:);
      unew(UPWIND,1)   = U(1)-A*(U(1)-U(end));
      unew(UPWIND,i)   = U(i)-A*(U(i)-U(i-1));
      unew(UPWIND,end) = U(end)-A*(U(end)-U(end-1));
   else %not used now, only a>0 supported in this version
      unew(UPWIND,1)   = u(1)-A*(u(2)-u(1));
      unew(UPWIND,i)   = u(i)-A*(u(i+1)-u(i));
      unew(UPWIND,end) = u(end)-A*(u(1)-u(end));
   end
end

if data.methods_mask(LAX_F) == 1
   i = 2:length(data.u_init)-1;
   U = u(LAX_F,:);
   unew(LAX_F,1)   = (1/2)*( U(end-1)+U(2))-A/2*(U(2)-U(end-1));
   unew(LAX_F,i)   = (1/2)*( U(i-1)+U(i+1))-A/2*(U(i+1)-U(i-1));
   unew(LAX_F,end) = (1/2)*( U(end-1)+U(1))-A/2*(U(2)-U(end-1));
end

if data.methods_mask(LEAPFROG) == 1
   i = 2:length(data.u_init)-1;
   
   if n == 1
      u_n_minus_1 = u(LEAPFROG,:);
      u_n = u_n_minus_1;
      
      %do half step using forward time centered space
      u_n(1)   = u_n_minus_1(1)+ data.a*k/(2*data.h)*...
         (u_n_minus_1(1)-u_n_minus_1(2));
      
      u_n(i)   = u_n_minus_1(i)+ data.a*k/(2*data.h)*...
         (u_n_minus_1(i+1)-u_n_minus_1(i));
      
      u_n(end) = u_n_minus_1(end)+ data.a*k/(2*data.h)*...
         (u_n_minus_1(2)-u_n_minus_1(end));
      
      u_n_plus_1 = u_n;
   else
      u_n_plus_1(1)   = u_n_minus_1(1) - A*(u_n(2)-u_n(end-1));
      u_n_plus_1(i)   = u_n_minus_1(i) - A*(u_n(i+1)-u_n(i-1));
      u_n_plus_1(end) = u_n_minus_1(end) - A*(u_n(2)-u_n(end-1));
   end
   
   unew(LEAPFROG,1:length(u_n_plus_1)) = u_n_plus_1;
   if n>=1
      u_n_minus_1 = u_n;
      u_n = u_n_plus_1;
   end
end

if data.methods_mask(BEAM_WARMING)==1
   U = u(BEAM_WARMING,:);
   unew(BEAM_WARMING,1) = U(1)-A/2*(3*U(1)-4*U(end-1)+U(end-2))+...
      B*(U(1)-2*U(end-1)+U(end-2));
   
   unew(BEAM_WARMING,2) =  U(2)-A/2*(3* U(2)-4* U(1)+ U(end-1))+...
      B*( U(2)-2* U(1)+ U(end-1));
   
   i = 3:size(u,2);
   unew(BEAM_WARMING,i)=  U(i)-A/2*(3* U(i)-4* U(i-1)+ ...
      U(i-2)) + B*( U(i)-2* U(i-1)+ U(i-2));
end

if data.methods_mask(CN)==1
   RHS        = CN_B*(u(CN,1:end-1)');
   unew(CN,1:end-1) = CN_A\RHS;
   unew(CN,end)     = unew(CN,1);
end

if data.methods_mask(FTCS)==1
   i = 2:length(data.u_init)-1;
   U = u(FTCS,:);
   unew(FTCS,1)   =  U(1)-A/2*(U(2)-U(end-1));
   unew(FTCS,i)   =  U(i)-A/2*(U(i+1)-U(i-1));
   unew(FTCS,end) =  U(end)-A/2*(U(2)-U(end-1));
end

u = unew;

if get(handles.user_function_btn_tag,'Value') == 1
   %u_exact = data.ic(mod(data.X-data.a*current_time,data.L));
   u_exact = data.ic(data.X-(data.a*current_time));
else
   u_exact = data.pulse.get(data.x_start+data.a*current_time);
end

update_plots(handles,data,u,u_exact,k,n,axes_handles,current_time);


%-----------------------------------------
function [u,k,u_exact] = run0(handles,data,h)
LAX_WENDROFF=1; UPWIND=2; LAX_F=3; LEAPFROG=4; BEAM_WARMING=5; CN=6;
FTCS=7; LAX=8;

turn_off_axes(data,handles);

if data.doing_refinement
   set(handles.axes_3_1,'Visible','on');
   set(handles.axes_3_2,'Visible','on');
   set(handles.axes_3_3,'Visible','on');
   
   axes_handles(1) = handles.axes_3_1;
else
   switch data.number_of_methods
      case 1
         set(handles.axes_1_1,'Visible','on');
         axes_handles = handles.axes_1_1;
      case 2
         set(handles.axes_2_1,'Visible','on');
         set(handles.axes_2_2,'Visible','on');
         
         axes_handles(1) = handles.axes_2_1;
         axes_handles(2) = handles.axes_2_2;
         
      case 3
         set(handles.axes_3_1,'Visible','on');
         set(handles.axes_3_2,'Visible','on');
         set(handles.axes_3_3,'Visible','on');
         
         axes_handles(1) = handles.axes_3_1;
         axes_handles(2) = handles.axes_3_2;
         axes_handles(3) = handles.axes_3_3;
      case 4
         set(handles.axes_4_1,'Visible','on');
         set(handles.axes_4_2,'Visible','on');
         set(handles.axes_4_3,'Visible','on');
         set(handles.axes_4_4,'Visible','on');
         
         axes_handles(1) = handles.axes_4_1;
         axes_handles(2) = handles.axes_4_2;
         axes_handles(3) = handles.axes_4_3;
         axes_handles(4) = handles.axes_4_4;
      case 5
         set(handles.axes_5_1,'Visible','on');
         set(handles.axes_5_2,'Visible','on');
         set(handles.axes_5_3,'Visible','on');
         set(handles.axes_5_4,'Visible','on');
         set(handles.axes_5_5,'Visible','on');
         
         axes_handles(1) = handles.axes_5_1;
         axes_handles(2) = handles.axes_5_2;
         axes_handles(3) = handles.axes_5_3;
         axes_handles(4) = handles.axes_5_4;
         axes_handles(5) = handles.axes_5_5;
   end
end

userData = get(handles.figure1, 'UserData');
userData.state = 1;
set(handles.figure1,'UserData',userData);
set(0,'CurrentFigure',handles.figure1);

k    = data.CFL*data.h/data.a;        % time step
unew = zeros(8,length(data.u_init));  % n+1 solution
u    = zeros(8,length(data.u_init));  % n   solution

for i = 1:8
   u(i,:) = data.u_init(:);    % initialize to initial data
end

% Make initial plot at t=0
update_plots(handles,data,u,u,k,0,axes_handles,0);  % initial plot, t=0

A            = data.a*k/h;   % parameters used below
B            = A^2/2;

u_exact = data.u_init;  %allocate storage for exact solution, updated below

CN_A=0;
CN_B=0;
u_n_minus_1=0;
u_n=0;

if data.methods_mask(CN)==1
   [CN_A,CN_B] = make_matrix_for_crank_niclson_solver(...
      length(data.u_init) ,data.a, k,h);
end

nSteps = data.max_t/k;
completeNumberOfSteps = floor(nSteps);

if (abs(nSteps-completeNumberOfSteps)>2*eps)
   %uiwait(errordlg('WARNING: time step calculated is too large to fit in the last interval.','Bad Input', 'modal'));
   %uicontrol(handles.CFL_tag);
   partialStep = (nSteps - completeNumberOfSteps)*k;
else
   partialStep = 0;
end

current_time = 0;
for n = 1:completeNumberOfSteps
    [u,unew,u_n_minus_1,u_n]=...
        make_step(data,n,handles,u,unew,A,B,CN_A,CN_B,u_n_minus_1,u_n,k,axes_handles,current_time);
    current_time = n*k;
    userData = get(handles.figure1, 'UserData');
    if userData.state == -1
        userData.state = 1;
        set(handles.figure1,'UserData',userData);
        return;
    end
end

if partialStep ~=0
   n = n+1;
   current_time = (n-1)*k + partialStep;
   [u,unew,u_n_minus_1,u_n]=...
   make_step(data,n,handles,u,unew,A,B,CN_A,CN_B,u_n_minus_1,u_n,partialStep,axes_handles,current_time);   
end

%----------------------------
% make plots
function make_plot(handles,data,u,u_exact,k,n,axes_handles,current_time,...
   plot_number,ID,name)

if data.doing_refinement
   ax = handles.axes_3_1;
else
   ax = axes_handles(plot_number);
end

set(gcf,'CurrentAxes',ax);
if get(handles.move_object_tag,'Value')==1
   plot(data.X,u(ID,1:length(data.X)),'r.-');
   hold on;
   plot(data.X,u_exact,'-');
   plot(data.X,data.u_init,'k-.');
   xlim([-data.h data.L+data.h]);   
else
   X = mod(data.X-data.a*current_time,data.L);
   
   set(gca,'FontSize',8);
   plot(X,u(ID,1:length(X)),'r.');
   hold on;
   plot(data.X,data.u_init,'-');
   xlim([-data.h data.L+data.h]);
end

ylim([data.min_y data.max_y]);
if data.doing_refinement
   title(ax,sprintf('%s, run# %d, time=%3.6f, time step=%3.6f, step number=%d',...
      name,data.current_refinement_run,current_time,k,n));
else
   title(sprintf('%s, time=%3.6f, time step=%3.6f, step number=%d, ||u||=%3.6f',...
      name,current_time,k,n,sqrt(data.h)*norm(u(ID,1:length(data.X)),2)));
end
grid on;
hold off;
drawnow();

%----------------------
function update_plots(handles,data,u,u_exact,k,n,axes_handles,current_time)

LAX_WENDROFF =1; UPWIND=2; LAX_F=3; LEAPFROG=4; BEAM_WARMING=5; CN=6;
FTCS=7; LAX=8;
plot_number  = 0;

if data.methods_mask(LAX_WENDROFF)==1
   plot_number = plot_number + 1;
   make_plot(handles,data,u,u_exact,k,n,axes_handles,current_time,plot_number,...
      LAX_WENDROFF,'Lax-Wendroff');
end

if data.methods_mask(UPWIND)==1
   plot_number = plot_number + 1;
   make_plot(handles,data,u,u_exact,k,n,axes_handles,current_time,plot_number,UPWIND,...
      'upwind');
end

if data.methods_mask(LAX_F)==1
   plot_number = plot_number + 1;
   make_plot(handles,data,u,u_exact,k,n,axes_handles,current_time,plot_number,LAX_F,...
      'Lax-Friedrich');
end

if data.methods_mask(LEAPFROG)==1
   plot_number = plot_number + 1;
   make_plot(handles,data,u,u_exact,k,n,axes_handles,current_time,plot_number,LEAPFROG,...
      'Leapfrog');
end

if data.methods_mask(BEAM_WARMING)==1
   plot_number = plot_number + 1;
   make_plot(handles,data,u,u_exact,k,n,axes_handles,current_time,plot_number,BEAM_WARMING,...
      'Beam-Warming');
end

if data.methods_mask(CN)==1
   plot_number = plot_number + 1;
   make_plot(handles,data,u,u_exact,k,n,axes_handles,current_time,plot_number,CN,...
      'Crank-Nicolson');
end

if data.methods_mask(FTCS)==1
   plot_number = plot_number + 1;
   make_plot(handles,data,u,u_exact,k,n,axes_handles,current_time,plot_number,FTCS,...
      'FTCS');
end

if data.methods_mask(LAX)==1
   plot_number = plot_number + 1;
   make_plot(handles,data,u,u_exact,k,n,axes_handles,current_time,plot_number,LAX,...
      'LAX');
end

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

data = get(handles.figure1, 'UserData');
data.state = -1;
enable_all(handles);
set(handles.figure1, 'UserData',data);


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

% Hints: get(hObject,'String') returns contents of run_time_tag as text
%        str2double(get(hObject,'String')) returns contents of run_time_tag as a double

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

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

% --- Executes when selected object is changed in initial_data_tag.
function initial_data_tag_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in initial_data_tag
% 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)

if get(handles.user_function_btn_tag,'Value')==1
   set(handles.user_function_tag,'Enable','on');
   set(handles.rect_width_tag,'Enable','off');
   set(handles.rect_start_at_tag,'Enable','off');
   set(handles.rect_height_tag,'Enable','off');
elseif get(handles.rect_tag,'Value')==1
   set(handles.user_function_tag,'Enable','off');
   set(handles.rect_width_tag,'Enable','on');
   set(handles.rect_start_at_tag,'Enable','on');
   set(handles.rect_height_tag,'Enable','on');
end

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

% Hints: get(hObject,'String') returns contents of rect_height_tag as text
%        str2double(get(hObject,'String')) returns contents of rect_height_tag as a double

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

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


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

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

% Hints: get(hObject,'String') returns contents of number_steps_value_tag as text
%        str2double(get(hObject,'String')) returns contents of number_steps_value_tag as a double

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

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

% --- Executes when selected object is changed in uipanel10.
function uipanel10_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in uipanel10
% 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)

if get(handles.max_time_btn_tag,'Value')==1
   set(handles.run_time_tag,'Enable','on');
   set(handles.number_steps_value_tag,'Enable','off');
elseif get(handles.number_steps_btn_tag,'Value')==1
   set(handles.run_time_tag,'Enable','off');
   set(handles.number_steps_value_tag,'Enable','on');
end

% --- Executes on button press in ref_tag.
function ref_tag_Callback(hObject, eventdata, handles)
% hObject    handle to ref_tag (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 ref_tag

if get(hObject,'Value')==1
   set(handles.ref_number_runs_tag,'Enable','on');
else
   set(handles.ref_number_runs_tag,'Enable','off');
end

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

% Hints: get(hObject,'String') returns contents of ref_number_runs_tag as text
%        str2double(get(hObject,'String')) returns contents of ref_number_runs_tag as a double

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

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

%-----------------------------------
function enable_all(handles)
set(handles.CN_tag,'Enable','on');
%set(handles.LAX_tag,'Enable','on');
set(handles.FTCS_tag,'Enable','on');

set(handles.lax_wendroff_tag,'Enable','on');
set(handles.upwind_tag,'Enable','on');
set(handles.lax_fried_tag,'Enable','on');
set(handles.leap_frog_tag,'Enable','on');
set(handles.beam_warming_tag,'Enable','on');

set(handles.user_function_btn_tag,'Enable','on');
if get(handles.user_function_btn_tag,'Value')==1
   set(handles.user_function_tag,'Enable','on');
end

set(handles.rect_tag,'Enable','on');
if get(handles.rect_tag,'Value')==1
   set(handles.rect_width_tag,'Enable','on');
   set(handles.rect_start_at_tag,'Enable','on');
   set(handles.rect_height_tag,'Enable','on');
end

set(handles.L_tag,'Enable','on');
set(handles.h_tag,'Enable','on');
set(handles.a_tag,'Enable','on');
set(handles.CFL_tag,'Enable','on');

set(handles.max_time_btn_tag,'Enable','on');
if get(handles.max_time_btn_tag,'Value')==1
   set(handles.run_time_tag,'Enable','on');
end

set(handles.number_steps_btn_tag,'Enable','on');
if get(handles.number_steps_btn_tag,'Value')==1
   set(handles.number_steps_value_tag,'Enable','on');
end

set(handles.ref_tag,'Enable','on');
if get(handles.ref_tag,'Value')==1
   set(handles.ref_number_runs_tag,'Enable','on');
end

set(handles.run_tag,'Enable','on');
%---------------------------------------
function disable_all(handles)
set(handles.CN_tag,'Enable','off');
%set(handles.LAX_tag,'Enable','off');
set(handles.FTCS_tag,'Enable','off');

set(handles.lax_wendroff_tag,'Enable','off');
set(handles.upwind_tag,'Enable','off');
set(handles.lax_fried_tag,'Enable','off');
set(handles.leap_frog_tag,'Enable','off');
set(handles.beam_warming_tag,'Enable','off');

set(handles.user_function_btn_tag,'Enable','off');
set(handles.user_function_tag,'Enable','off');

set(handles.rect_tag,'Enable','off');

set(handles.rect_width_tag,'Enable','off');
set(handles.rect_start_at_tag,'Enable','off');
set(handles.rect_height_tag,'Enable','off');

set(handles.L_tag,'Enable','off');
set(handles.h_tag,'Enable','off');
set(handles.a_tag,'Enable','off');
set(handles.CFL_tag,'Enable','off');

set(handles.max_time_btn_tag,'Enable','off');
set(handles.run_time_tag,'Enable','off');
set(handles.number_steps_btn_tag,'Enable','off');
set(handles.number_steps_value_tag,'Enable','off');
set(handles.ref_tag,'Enable','off');
set(handles.ref_number_runs_tag,'Enable','off');
set(handles.run_tag,'Enable','off');

%------------------------------
function turn_off_axes(data,handles)

cla(handles.axes_1_1,'reset');

cla(handles.axes_2_1,'reset');
cla(handles.axes_2_2,'reset');

cla(handles.axes_3_1,'reset');

if not(data.doing_refinement)
   cla(handles.axes_3_2,'reset');
   cla(handles.axes_3_3,'reset');
end

cla(handles.axes_4_1,'reset');
cla(handles.axes_4_2,'reset');
cla(handles.axes_4_3,'reset');
cla(handles.axes_4_4,'reset');

cla(handles.axes_5_1,'reset');
cla(handles.axes_5_2,'reset');
cla(handles.axes_5_3,'reset');
cla(handles.axes_5_4,'reset');
cla(handles.axes_5_5,'reset');

set(handles.axes_1_1,'Visible','off');

set(handles.axes_2_1,'Visible','off');
set(handles.axes_2_2,'Visible','off');

set(handles.axes_3_1,'Visible','off');
set(handles.axes_3_2,'Visible','off');
set(handles.axes_3_3,'Visible','off');

set(handles.axes_4_1,'Visible','off');
set(handles.axes_4_2,'Visible','off');
set(handles.axes_4_3,'Visible','off');
set(handles.axes_4_4,'Visible','off');

set(handles.axes_5_1,'Visible','off');
set(handles.axes_5_2,'Visible','off');
set(handles.axes_5_3,'Visible','off');
set(handles.axes_5_4,'Visible','off');
set(handles.axes_5_5,'Visible','off');

% --- Executes on button press in CN_tag.
function CN_tag_Callback(hObject, eventdata, handles)
% hObject    handle to CN_tag (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 CN_tag

%-----------------------
% This function makes the A and B matrix for use by the C-N solver
% since it is an implicit solver.
function [A,B] = make_matrix_for_crank_niclson_solver(n,a,k,h)

N  = n-1;
mu = a*k/h;

e  = ones(N,1);
D  = [-(mu/4)*e e (mu/4)*e];
A  = spdiags(D,[-1 0 1],N,N);
A(1,end) = -mu/4;
A(end,1) = mu/4;

B = A';

% --- Executes on button press in FTCS_tag.
function FTCS_tag_Callback(hObject, eventdata, handles)
% hObject    handle to FTCS_tag (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 FTCS_tag

% --- Executes on button press in LAX_tag.
function LAX_tag_Callback(hObject, eventdata, handles)
% hObject    handle to LAX_tag (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 LAX_tag