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

% Last Modified by GUIDE v2.5 22-Aug-2012 11:19:50

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

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


getUserData(handles);
% Initialise plot
updatePlot(handles);

% Update handles structure
guidata(hObject, handles);

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


% --- Outputs from this function are returned to the command line.
function varargout = orth_poly_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 update.
function update_Callback(hObject, eventdata, handles)
% hObject    handle to update (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
updatePlot(handles); % Assume getUserData already called

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

contents=cellstr(get(hObject,'String'));
index=get(hObject,'Value');

switch(index)
    case 1
        % do nothing

    case 2
        % Legendre
        set(handles.edit5,'String',1);    
        set(handles.edit6,'String','-1');
        set(handles.edit7,'String','1');
        
    case 3
        % Chebyshev (1st kind)
        set(handles.edit5,'String','1/sqrt(1-x^2)');    
        set(handles.edit6,'String','-1');
        set(handles.edit7,'String','1');
        
    case 4
        % Laguerre
        set(handles.edit5,'String','exp(-x)');    
        set(handles.edit6,'String','0');
        set(handles.edit7,'String','Inf');
    case 5
        % Hermite
        set(handles.edit5,'String','exp(-x^2)');    
        set(handles.edit6,'String','-Inf');
        set(handles.edit7,'String','Inf');       
end

getUserData(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 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
getUserData(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 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
set(handles.popupmenu1,'Value',1);
getUserData(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 edit6_Callback(hObject, eventdata, handles)
% hObject    handle to edit6 (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 edit6 as text
%        str2double(get(hObject,'String')) returns contents of edit6 as a double
getUserData(handles);

% --- Executes during object creation, after setting all properties.
function edit6_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit6 (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 edit7_Callback(hObject, eventdata, handles)
% hObject    handle to edit7 (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 edit7 as text
%        str2double(get(hObject,'String')) returns contents of edit7 as a double
getUserData(handles);

% --- Executes during object creation, after setting all properties.
function edit7_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit7 (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 updatePlot(handles)
    cla(handles.axes1,'reset');
    set(handles.axes1,'Visible','on');

    % Assume we've already got valid UserData
    ud=get(handles.figure1,'UserData');

    a=ud.a;
    plot_a=a;
    b=ud.b;
    plot_b=b;
    n=ud.n;
    w=ud.w;

    if( isinf(-a) )
        plot_a=-1;
        % go sufficiently negative for plotting purposes
        while(abs(w(plot_a))>1e-6) 
            plot_a=plot_a-1;
        end
    end

    if( isinf(b) )
        plot_b=1;
        % go sufficiently negative for plotting purposes
        while(abs(w(plot_b))>1e-6) 
            plot_b=plot_b+1;
        end
    end

    w=@(x)(w(x));
    w=@(x)(arrayfun(w,x));


    p=orthpoly(w, a, b, n);
    [r,c]=size(p);
    entries=cell(1,r);

    co=get(handles.axes1,'ColorOrder');
    [mco,~]=size(co);

    styles={'-','--','-.'};

    hold on;
    for l=1:r
        plot_fun=@(x)(polyval(p(l,:),x));
        h=ezplot(plot_fun,[plot_a plot_b]);
        set(h,'Color',co(mod(l-1,mco)+1,:));    

        set(h,'LineStyle',styles{mod(l-1,3)+1});
        entries{l}=['Deg. ' sprintf('%d',l-1)];
    end
    title('Orthogonal Polynomials');
    xlabel('x');
    legend(entries,'Location','SouthEast');
    axis([plot_a plot_b -1.5*max(abs(plot_a),abs(plot_b)) 1.5*max(abs(plot_a),abs(plot_b))]);
    axis square
    hold off;

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

function getUserData(handles)

try
    ud=get(handles.figure1,'UserData');
    
    
    %***Check degree***
    n=str2double(get(handles.edit4,'String'));
    assert(~isnan(n) && ~isinf(n) && ceil(n) == n && isreal(n) &&...
        and((n >= 1),(n <= 5)),...
        'Degree must be an integer between 1 and 5.');
    ud.n=n;

    %***Check interval***
    a=str2double(get(handles.edit6,'String'));
    b=str2double(get(handles.edit7,'String'));
    
    assert(~or(isnan(a),isnan(b)) && and(isreal(a),isreal(b)) ...
        && a<b, 'Interval limits must be real and increasing.');
    ud.a=a;
    ud.b=b;
    
    
    %***Check weight w***
    w = inline((get(handles.edit5,'String')),'x');
    x = sym('x');
    feval(w,x);
    
    %Check w(x) is [approximately] [almost everywhere] positive on interval

    integr = @(x) (abs(w(x)) - w(x));
    integr = @(x)arrayfun(integr,x);
    integral = quadgk(integr,a,b);
    assert(integral < 10^(-12),'Quadrature differs for w(x) and abs(w(x)). Is w(x) [almost-everywhere] positive?');
    
    ud.w=w;
    
    set(handles.figure1,'UserData',ud);
    
    set(handles.update,'Enable','on');
    set(handles.latex,'Enable','on');
catch err
    cla(handles.axes1,'reset');
    set(handles.axes1,'Visible','off')
    set(handles.update,'Enable','off');
    set(handles.latex,'Enable','off');
    h=msgbox(err.message,'Error','error','modal');
    waitfor(h);
end

function updateLatex(handles)
cla(handles.axes1,'reset');

getUserData(handles);

ud=get(handles.figure1,'UserData');

a=ud.a;
b=ud.b;
n=ud.n;
w=ud.w;

w=@(x)(w(x));
w=@(x)(arrayfun(w,x));

p=orthpoly(w, a, b, n);

rationals=get(handles.checkbox1,'Value');

str=poly2latex(p,rationals);

text(0.5,0.5,str,...
    'Interpreter','latex',...
    'FontSize',16.0,...
    'HorizontalAlignment','Center',...
    'VerticalAlignment','Middle');
set(handles.axes1,'XTick',[]);
set(handles.axes1,'YTick',[]);
axis off


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

function [ string ] = poly2latex( p, rational )
%poly2latex Prints a set of polynomials stored as an array in LaTeX form.
%   Input arguments:
%       p - an nxn matrix storing n monic polynomials (in the rows) with
%           degree equal to the row index.
%       rational - true if rational approximations are to be used, false
%                  otherwise
%   Output arguments:
%       string - the LaTeX markup

string='';
rattol=1e-9;
tol=1e-3;

[rownum,colnum]=size(p);

for row=1:rownum
    for col=1:colnum
        coeff=p(row,col);
        exponent=colnum-col;
        
        if(rational)
            [a,b]=rat(coeff,rattol);
            if( a ~= 0 )
                if( sign(a)*sign(b) > 0 && col ~= (colnum-row+1) )
                    string=[string '+'];
                elseif( sign(a)*sign(b) < 0 )
                    string=[string '-'];
                end
                if( b~=1 )
                    string=[string '\frac' sprintf('{%d}{%d}',abs(a),abs(b)) ];
                elseif( abs(a) ~= 1)
                    string=[string sprintf('%d',abs(a))];
                elseif( exponent == 0)
                    string=[string '1'];
                end
                
                if( exponent > 0 )
                string = [string ' x'];
                if( exponent > 1 )
                    string = [string sprintf('^%d',exponent)];
                end
        end
            end
        else
            if( abs(coeff) > tol ) % nonzero up to tolerance
                if( row ~= 1 && col ~= colnum-row+1 )
                    if( abs(coeff-1) < tol)
                        string=[string '+'];
                        if( exponent==0)
                            string=[string '1'];
                        end
                    elseif(abs(coeff+1) < tol)
                        string=[string '-'];
                        if( exponent==0)
                            string=[string '1'];
                        end
                    else
                        string=[string sprintf('%+0.4g',coeff)];
                    end
                else
                    if( exponent == 0 )
                        string=[string '1'];
                    end
                end
                                   
                if( exponent > 0 )
                    string = [string ' x'];
                    if( exponent > 1 )
                        string = [string sprintf('^%d',exponent)];
                    end
                end
            end                    
        end
        if( col ~= colnum)
            string=[string '&'];
        end
    end
    string = [string '\\'];
end

columnspec=repmat(['c'],1,colnum);

string = ['$\begin{array}{' columnspec '}' string '\end{array}$'];

function [ p ] = orthpoly(w, a, b, n)
%orthpoly Returns the orthogonal polynomials of order no more than n
%
%   Orthogonal polynomials are considered on the interval [a,b] with 
%   non-negative weight function w(x) supported almost-everywhere in [a,b].
%
%   Input Arguments:
%       w - the weight function (as a function handle @(x)... )
%       a - the lower limit of the interval
%       b - the upper limit of the interval
%       n - the highest degree of a polynomial in the set
%   Output Arguments:
%       p - the set of polynomials, returned as an (n+1)x(n+1) matrix with
%           the degree of the polynomial in row N having degree N-1
%           (polynomials are stored in rows)

if( n < 1 || a >= b ) % incorrect usage...note that w is not checked!!!
    fprintf('ERROR: Incorrect a,b or n passed to orthpoly\n');
    p=NaN;
else
    p=zeros(n+1,n+1);
    
    % Recursion coefficients
   
    % n = 0

    % alpha0 = <1,x>/<1,1>
    f1=@(x)(w(x).*x); 
    pn_pn=quadgk(w,a,b);
    alpha_n=quadgk(f1,a,b)/pn_pn; % <1,x>/<1,1>
    
    p(1,n+1)=1; % p_0
    
    % p_1
    p(2,n)=1;    
    p(2,n+1)=-alpha_n;
    
    % <p_1,p_1>
    c=conv(p(2,:),p(2,:));
    f2=@(x)(w(x).*polyval(c,x));
    pn1_pn1=quadgk(f2,a,b);
    
    % advance to n = 1
    beta_n=pn1_pn1/pn_pn;
    pn_pn=pn1_pn1;
        
    % <p_1,x*p_1>
    c=conv(p(2,:),circshift(p(2,:),[0 -1]));
    f1=@(x)(w(x).*polyval(c,x));
    
    % alpha_1
    alpha_n=quadgk(f1,a,b)/pn_pn;
    

    for l=2:n % now do the rest
        c=conv(p(l,:),[1 -alpha_n]);
        c=c(2:end); % need to make the size n+1 by trimming zeros on left
        p(l+1,:)= c-beta_n*p(l-1,:); % three-term recurrence relation
        
        % integrand for <p_n,p_n>
        c=conv(p(l+1,:),p(l+1,:));
        f2=@(x)(w(x).*polyval(c,x));
        pn1_pn1=quadgk(f2,a,b); %  <p_(n+1),p_(n+1)> = <p_1,p_1>

        % advance n -> n+1
        beta_n=pn1_pn1/pn_pn;
        pn_pn=pn1_pn1;

        % integrand for <p_n,x*p_n>
        c=circshift(p(l+1,:),[0 -1]);   
        pn_xpn=conv(p(l+1,:),c);
        f1=@(x)(w(x).*polyval(pn_xpn,x));
        alpha_n=quadgk(f1,a,b)/pn_pn;
     end
end
