function varargout = Sturm_Liouville(varargin)
% STURM_LIOUVILLE M-file for Sturm_Liouville.fig
%      STURM_LIOUVILLE, by itself, creates a new STURM_LIOUVILLE or raises the existing
%      singleton*.
%
%      H = STURM_LIOUVILLE returns the handle to a new STURM_LIOUVILLE or the handle to
%      the existing singleton*.
%
%
% How to use the GUI:
%   - Choose a family of eigenfunctions by clicking on one of the 
%   buttons in the panel entitled "Types of orthogonal functions".
%   Sines and cosines are chosen by default when the GUI initializes.
%   - Use the sliders in the third panel to select and plot two modes
%   in the chosen family. Twenty one modes are available in each family.
%   - The bottom panel shows the weighted product of the two selected
%   modes, as well as a numerical estimate of their weighted dot
%   product.
%   - Use the slider in the panel entitled "Domain of the Bessel
%   functions" to rescale the independent variable in the Bessel
%   functions J_n, n = 1, 2, or 3.
%   
% Developed by Benjamin Dyhr, Mathematics Department, University
% of Arizona, Spring 2007.
% 

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

global orth_functions

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

% Update handles structure
guidata(hObject, handles);

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

set(handles.figure1,'Name',' Sturm-Liouville Eigenfunctions');

orth_functions = 1;
[x,r,f1,f2,x_min,x_max,y_min,y_max]= modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);

% --- Outputs from this function are returned to the command line.
function varargout = Sturm_Liouville_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 during object creation, after setting all properties.
function Phi_m_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Phi_m (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 Phi_m


% --- Executes on slider movement.
function slider1_Callback(hObject, eventdata, handles)
% hObject    handle to slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global orth_functions
% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

sl1 = floor(get(handles.slider1,'Value'));
set(handles.slider1,'Value',sl1);
set(handles.text1,'String',['Mode # ' int2str(sl1)]);

[x,r,f1,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);

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

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end


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

% Hint: place code in OpeningFcn to populate Phi_n


% --- Executes on slider movement.
function slider2_Callback(hObject, eventdata, handles)
% hObject    handle to slider2 (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,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
sl2 = floor(get(handles.slider2,'Value'));
set(handles.slider2,'Value',sl2);
set(handles.text2,'String',['Mode # ' int2str(sl2)]);


[x,r,f1,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);

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

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end


% --- Executes on slider movement.
function slider3_Callback(hObject, eventdata, handles)
% hObject    handle to slider3 (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,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

sl3 = ceil(get(handles.slider3,'Value'));
if sl3 == 0
    sl3 = 1;
end
set(handles.slider3,'Value',sl3);
set(handles.text3,'String',['[0, ' int2str(sl3) ']']);

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

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
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)
global orth_functions

orth_functions = 1;
[x,r,f1,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);
set(handles.uipanel2,'Title','Sines and cosines')

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

orth_functions = 2;
[x,r,f1,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);
set(handles.uipanel2,'Title','Legendre polynomials')

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

orth_functions = 3;
[x,r,f1,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);
set(handles.uipanel2,'Title','Bessel functions (n=0)')

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

orth_functions = 4;
[x,r,f1,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);
set(handles.uipanel2,'Title','Bessel functions (n=1)')

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

orth_functions = 5;
[x,r,f1,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);
set(handles.uipanel2,'Title','Bessel functions (n=3)')

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

orth_functions = 6;
[x,r,f1,f2,x_min,x_max,y_min,y_max] = modes(hObject, eventdata, handles);
plotmodes(hObject, eventdata, handles);
set(handles.uipanel2,'Title','Chebyshev polynomials')

function [x,r,f1,f2,x_min,x_max,y_min,y_max]=modes(hObject, eventdata, handles);

global orth_functions rst
global F

sl1 = floor(get(handles.slider1,'Value'));
sl2 = floor(get(handles.slider2,'Value'));
sl3 = floor(get(handles.slider3,'Value'));

if orth_functions == 1
    
    rst = 'Weight function: r(x)= 1';
    
    x_min = -pi;
    x_max = pi;
    y_min = -1.1;
    y_max = 1.2;
    
    x=[x_min:(x_max-x_min)/max(sl1*sl2*3,100):x_max];
    
    siz=size(x);
    r=ones(1,siz(2));
    
    if sl1-2*floor(sl1/2) == 0
        f1 = cos(sl1*x/2);
        fun1 = @(z)cos(sl1*z/2);
    else
        f1 = sin((sl1+1)*x/2);
        fun1 = @(z)sin((sl1+1)*z/2);
    end
    if sl2-2*floor(sl2/2) == 0
        f2 = cos(sl2*x/2);
        fun2 = @(z)cos(sl2*z/2);
    else
        f2 = sin((sl2+1)*x/2);
        fun2 = @(z)sin((sl2+1)*z/2);
    end
    F = @(z)fun1(z).*fun2(z);
    
elseif orth_functions == 2
    
    rst = 'Weight function: r(x)= 1';
    
    x_min = -1;
    x_max = 1;
    y_min = -1;
    y_max = 1;
    
    x=[x_min:(x_max-x_min)/max(sl1*sl2*3,100):x_max];
    
    siz=size(x);
    r=ones(1,siz(2));
  
    f1 = lgndr(sl1,x);
    f2 = lgndr(sl2,x);
    
    F = @(z) lgndr(sl1,z).*lgndr(sl2,z);
    
elseif orth_functions == 3
    
    rst = 'Weight function: r(x)= x';
    
    Bessel0_zeros = [2.40482 5.52007 8.65372 11.79153 14.93091 18.07106 21.21163 24.35247 27.49347 30.63460 33.77582 36.91709 40.05842 43.19979 46.34118 49.48260 52.62405 55.76551 58.90698 62.04846 65.18996];
    
   
    x_min = 0;
    x_max = sl3;
    y_min = -1;
    y_max = 1;
    
    x=[x_min:(x_max-x_min)/max(sl1*sl2*3,100):x_max];
    
    r=x;
    f1=besselj(0,Bessel0_zeros(1,sl1+1)/x_max.*x);
    f2=besselj(0,Bessel0_zeros(1,sl2+1)/x_max.*x);
    
    F = @(z)besselj(0,Bessel0_zeros(1,sl1+1)/x_max.*z)...
        .*besselj(0,Bessel0_zeros(1,sl2+1)/x_max.*z).*z;

elseif orth_functions == 4
    
    rst = 'Weight function: r(x)= x';
    
    Bessel1_zeros = [3.83170 7.015586 10.17346 13.32369 16.47063 19.61585 22.76008438059277 25.90367 29.04682 32.18967 35.33230 38.47476 41.61709 44.75931 47.90146 51.04353 54.18555 57.32752 60.46945 63.61135 66.75323];
    
   
    x_min = 0;
    x_max = sl3;
    y_min = -1;
    y_max = 1;
    
    x=[x_min:(x_max-x_min)/max(sl1*sl2*3,100):x_max];
    
    r=x;
    f1=besselj(1,Bessel1_zeros(1,sl1+1)/x_max.*x);
    f2=besselj(1,Bessel1_zeros(1,sl2+1)/x_max.*x);

    F = @(z)besselj(1,Bessel1_zeros(1,sl1+1)/x_max.*z)...
        .*besselj(1,Bessel1_zeros(1,sl2+1)/x_max.*z).*z;
    
elseif orth_functions == 5

    rst = 'Weight function: r(x)= x';
    
    Bessel2_zeros = [5.135622301840683 8.417244140399855 11.61984117214906 14.79595178235126 17.95981949498783 21.11699705302185 24.27011231357310 27.42057354998456 30.56920449551640 33.71651950922270 36.86285651128381 40.00844673347819 43.15345377837146 46.29799667723692 49.44216411041687 52.58602350681596 55.72962705320114 58.87301577261216 62.01622235921765 65.15927319075780 68.30219];
    
   
    x_min = 0;
    x_max = sl3;
    y_min = -1;
    y_max = 1;
    
    x=[x_min:(x_max-x_min)/max(sl1*sl2*3,100):x_max];
    
    r=x;
    f1=besselj(2,Bessel2_zeros(1,sl1+1)/x_max.*x);
    f2=besselj(2,Bessel2_zeros(1,sl2+1)/x_max.*x);
    
    F = @(z)besselj(2,Bessel2_zeros(1,sl1+1)/x_max.*z)...
        .*besselj(2,Bessel2_zeros(1,sl2+1)/x_max.*z).*z;
    
elseif orth_functions == 6
    
%     set(handles.text17,'String',['r(x)= (1-x^2)^{-1/2}' ]);
    rst = 'Weight function: r(x)= (1-x^2)^{-1/2}';
    
    x_min = -1;
    x_max = 1;
    y_min = -2;
    y_max = 2;
    
    x=[x_min:(x_max-x_min)/max(sl1*sl2*6,1000):x_max];
    
    r=sqrt(1-x.^2).^(-1);
    f1=cos(sl1.*acos(x));
    f2=cos(sl2.*acos(x));
    
    F = @(z)cos(sl1.*acos(z)).*cos(sl2.*acos(z))...
        .*sqrt(1-z.^2).^(-1);
    
end



function plotmodes(hObject, eventdata, handles)

global orth_functions rst
global F

[x,r,f1,f2,x_min,x_max,y_min,y_max]=modes(hObject, eventdata, handles);

axes(handles.Phi_m)
plot(x,f1,'LineWidth',2)
axis([x_min x_max y_min y_max])

axes(handles.Phi_n)
plot(x,f2,'LineWidth',2)
axis([x_min x_max y_min y_max])

axes(handles.PHI_mn)
plot(x,r.*f1.*f2,x,0,'-m','LineWidth',2)
xlim([x_min x_max])
title(rst)

if orth_functions == 6
    x_mi=-0.99999999999;
    x_ma=0.99999999999;
else
    x_mi=x_min;
    x_ma=x_max;
end
qq=quadl(F,x_mi,x_ma,0.0000000001);
if abs(qq) <= 0.00001
    qq = 0;
end
set(handles.text9,'String',num2str(qq,'%12.4f'));


function V = lgndr(n,x)

U=legendre(n,x);
V=U(1,:);