function [y,a,x,Z,b]=ieq(ctrl,x,y)
% Solves the integral equation
%   1/4 Int_0^1(sin^2(x-u)/(x-u)^2*y(u))du = y(x)-x
% by the Method of Moments.
% The GUI can be used to choose the basis functions, method, and number of
% terms to be used. The function outputs are sent to MATLAB workspace.
% y : The solution y at points x.
% a : The solved coefficients for the basis functions.
% x : The points at which y is evaluated.
% Z : The coefficient matrix for MOM.
% b : The excitation matrix for MOM.
% The output can be plotted by the command
%       plot(x,y);

% Dr. Sencer Koç, 12 March 2008

if nargin==0        % This part is to create the GUI
    figure(1)
    set(gcf,'toolbar','figure')
    htbf = uicontrol('Style','Frame');
    set(htbf,'units','normalized','position',[0.02 0.9 0.25 0.05])
    htbasis = uicontrol('Style','Text','String','Choose Basis Functions :');
    set(htbasis,'units','normalized','position',[0.03 0.905 0.23 0.04],...
        'HorizontalAlignment','left')
    hbasis = uicontrol('Style', 'popup','tag','btype','String',...
        'Full-domain, polynomial|Full-domain, exponential|Sub-domain, rectangular|Sub-domain, triangular');
    set(hbasis,'units','normalized','position',[0.02 0.8 0.25 0.1])

    htmf = uicontrol('Style','Frame');
    set(htmf,'units','normalized','position',[0.3 0.9 0.18 0.05])
    htmethod = uicontrol('Style','Text','String','Choose Method :');
    set(htmethod,'units','normalized','position',[0.31 0.905 0.16 0.04],...
        'HorizontalAlignment','left')
    hmethod = uicontrol('Style','popup','tag','mtype','String',...
        'Point matching|Galerkin|');
    set(hmethod,'units','normalized','position',[0.3 0.8 0.18 0.1])

    hnf = uicontrol('Style','Frame');
    set(hnf,'units','normalized','position',[0.5 0.9 0.18 0.05])
    htn = uicontrol('Style','Text','String','Number of terms :');
    set(htn,'units','normalized','position',[0.51 0.905 0.16 0.04],...
        'HorizontalAlignment','left')
    hn = uicontrol('Style', 'edit','String','10','tag','N',...
        'HorizontalAlignment','left');
    set(hn,'units','normalized','position',[0.5 0.85 0.18 0.05],'Callback','ieq(1);')

    hrun = uicontrol('style','pushbutton','string','RUN',...
        'tag','run','Callback','[y,a,x,Z,b]=ieq(4);');
    set(hrun,'units','normalized','position',[0.7 0.9 0.1 0.05])
    hexit = uicontrol('style','pushbutton','string','EXIT',...
        'tag','exit','Callback','close all,clear all,clc');
    set(hexit,'units','normalized','position',[0.7 0.85 0.1 0.05])
    
    hplt = uicontrol('style','pushbutton','string','PLOT',...
        'tag','plot','Callback','ieq(2,x,y);');
	set(hplt,'units','normalized','position',[0.83 0.9 0.15 0.05])
    hclrplt = uicontrol('style','pushbutton','string','CLEAR PLOT',...
        'tag','clrplt','Callback','ieq(3);');
	set(hclrplt,'units','normalized','position',[0.83 0.85 0.15 0.05])
    
    haxes = axes('position',[.1  .1  .8  .7]);
    hold on
else
	figure(1)
    ch = get(gcf,'children');
    for ii=1:length(ch)
        if strcmp(get(ch(ii),'tag'),'btype')
            hbtype = ii;
        end
        if strcmp(get(ch(ii),'tag'),'mtype')
            hmtype = ii;
        end
        if strcmp(get(ch(ii),'tag'),'N')
            hN = ii;
        end
        if strcmp(get(ch(ii),'tag'),'run')
            hrun = ii;
        end
        if strcmp(get(ch(ii),'tag'),'exit')
            hexit = ii;
        end
        if strcmp(get(ch(ii),'tag'),'plot')
            hplot = ii;
        end
        if strcmp(get(ch(ii),'tag'),'clrplt')
            hclrplt = ii;
        end
    end
    switch ctrl
        case 1      % Checking the input for N
            N = str2num(get(ch(hN),'string'));
            if numel(N)~=1 || ~isnumeric(N) || N<1
                set(ch(hN),'string','10');
            end
        case 2      % Plot function
            k = length(get(gca,'children'));
            hh = plot(x,y);
            cm = get(gca,'colororder');
            kc = mod(k+1,7);
            if kc==0, kc=7; end
            set(hh,'color',cm(kc,:))
            btype = get(ch(hbtype),'value');
            btype_string = get(ch(hbtype),'string');
            btype_string = deblank(btype_string(btype,:));
            mtype = get(ch(hmtype),'value');
            mtype_string = get(ch(hmtype),'string');
            mtype_string = deblank(mtype_string(mtype,:));
            N = str2num(get(ch(hN),'string'));
            if k==0
                legend(hh,[btype_string '; ' mtype_string '; N = ' num2str(N)],...
                    'location','best');
            else
                [legend_h,object_h,plot_h,text_strings] = legend;
                legend(text_strings,[btype_string '; ' mtype_string '; N = ' num2str(N)])
                legend('location','best');
            end
        case 3      % Clear plot function
            delete(get(gca,'children'));
            legend('off')
        case 4      % Execution
            ch = get(gcf,'children');
            set(ch(hrun),'enable','off')
            set(ch(hexit),'enable','off')
            set(ch(hplot),'enable','off')
            set(ch(hclrplt),'enable','off')
            btype = get(ch(hbtype),'value');
            mtype =  get(ch(hmtype),'value');
            N = str2num(get(ch(hN),'string'));
            if btype==1 && mtype==1     % Full-domain polynomial basis, point matching
                [y,a,x,Z,b]=ieqpfp(N);
            end
            if btype==2 && mtype==1     % Full-domain exponential basis, point matching
                [y,a,x,Z,b]=ieqefp(N);
            end
            if btype==3 && mtype==1     % Sub-domain rectangular basis, point matching
                [y,a,x,Z,b]=ieqrsp(N);
            end
            if btype==4 && mtype==1     % Sub-domain triangular basis, point matching
                [y,a,x,Z,b]=ieqtsp(N);
            end
            if btype==1 && mtype==2     % Full-domain polynomial basis, Galerkin
                [y,a,x,Z,b]=ieqpfg(N);
            end
            if btype==2 && mtype==2     % Full-domain exponential basis, Galerkin
                [y,a,x,Z,b]=ieqefg(N);
            end
            if btype==3 && mtype==2     % Sub-domain rectangular basis, Galerkin
                [y,a,x,Z,b]=ieqrsg(N);
            end
            if btype==4 && mtype==2     % Sub-domain triangular basis, Galerkin
                [y,a,x,Z,b]=ieqtsg(N);
            end
            set(ch(hrun),'enable','on')
            set(ch(hexit),'enable','on')
            set(ch(hplot),'enable','on')
            set(ch(hclrplt),'enable','on')
    end
end

function [y,a,x,A,b]=ieqpfg(N)
% Polynomial basis, Full domain, Galerkin
b = zeros(N+1,1);
A = zeros(N+1);
% Define the kernel as an anonymous function
K = @(x,ksi) sinc((x-ksi)/pi).^2/4;
hwait = waitbar(0,'Calculating! Please wait...');
for m=0:N
    % right hand side is calculated analytically
    b(m+1) = -1/(m+2);
    for n=0:N
        % Define the integrand as an anonymous function
        afun = @(x,ksi) K(x,ksi).*ksi.^n.*x.^m;
        % Double integration for testing
        A(m+1,n+1) = dblquad(afun,0,1,0,1);
        A(m+1,n+1) = A(m+1,n+1)-1/(n+m+1);
    end
    waitbar(m/(N+1),hwait)
end
close(hwait)
% Solve the system of equations
a = A\b;
% Define plot points
x = linspace(0,1,100);
% Calculate the function at plot points
y = polyval(flipud(a),x);

function [y,a,x,A,b]=ieqefg(N)
% Exponential basis, Full domain, Galerkin
jey = sqrt(-1);
b = zeros(2*N+1,1);
A = zeros(2*N+1);
% Define the kernel as an anonymous function
K = @(x,ksi) sinc((x-ksi)/pi).^2/4;
hwait = waitbar(0,'Calculating! Please wait...');
for m=-N:N
    % Define the integrand for right hand side as anonymous function
    bfun = @(x) -x.*exp(jey*2*pi*m*x);
    % Evaluate the integral for the right hand side
    b(N+m+1) = quad(bfun,0,1);
    for n=-N:N
        % Define the integrand as an anonymous function
        zfun = @(x,ksi) K(x,ksi).*exp(jey*2*pi*n*ksi).*exp(jey*2*pi*m*x);
        % Double integration for testing
        A(N+m+1,N+n+1) = dblquad(zfun,0,1,0,1);
        % Evaluate the integral of cross-products of test functions
        afun = @(x) exp(jey*2*pi*n*x).*exp(jey*2*pi*m*x);
        A(N+m+1,N+n+1) = A(N+m+1,N+n+1)-quad(afun,0,1);
        waitbar(((N+m)*(2*N+1)+n)/(2*N+1)^2,hwait)
    end
end
close(hwait)
% Solve the system of equations
a = A\b;
% Define plot points
x = linspace(0,1,100);
% Calculate the function at plot points
y = zeros(size(x));
for m=-N:N
    y = y+real(a(N+m+1)*exp(jey*2*pi*m*x));
end

function [y,a,x,A,b]=ieqrsg(N)
% rectangular pulse basis, Subdomain, Galerkin
b = zeros(N,1);
A = zeros(N);
% Define the kernel as an anonymous function
K = @(x,ksi) sinc((x-ksi)/pi).^2/4;
% Define intervals
xi = [0:N]/N;
hwait = waitbar(0,'Calculating! Please wait...');
for m=1:N
    % right hand side is calculated analytically
    b(m) = -(xi(m+1)^2-xi(m)^2)/2;
    for n=1:N
        % Double integration for testing
        % With rectangular sub-domain basis only integration limits are
        % affected
        A(m,n) = dblquad(K,xi(m),xi(m+1),xi(n),xi(n+1));
        if m==n
            A(m,n) = A(m,n)-(xi(n+1)-xi(n));
        end
        waitbar(((m-1)*N+n)/N^2,hwait)
    end
end
close(hwait)
% Solve the system of equations
a = A\b;
% Define plot points
x = linspace(0,1,100);
% Calculate the function at plot points
y = zeros(size(x));
for m=1:N
    y = y+a(m)*(x>=xi(m)).*(x<=xi(m+1));
end

function [y,a,x,A,b]=ieqtsg(N)
% Triangular pulse basis, Subdomain, Galerkin
b = zeros(N+1,1);
A = zeros(N+1);
G = A;
% Define the kernel as an anonymous function
K = @(x,ksi) sinc((x-ksi)/pi).^2/4;
% Define a generic triangular function over [a,b]
T = @(x,a,b) (((x<(a+b)/2).*(x-a)+(x>=(a+b)/2).*(b-x))*2/(b-a)).*(x>=a).*(x<=b);
% Define intervals
xi = [0:N]/N;
hwait = waitbar(0,'Calculating! Please wait...');
for m=0:N
    % Define m'th basis function
    % First and last intervals should be treated separately
    if m==0
        Tm = @(x) T(x,-xi(2),xi(2)).*(x>=0);
        startm = 0;
        stopm = xi(2);
    elseif m==N
        Tm = @(x) T(x,xi(N),2-xi(N)).*(x<=1);
        startm = xi(N);
        stopm = 1;
    else
        Tm = @(x) T(x,xi(m),xi(m+2));
        startm = xi(m);
        stopm = xi(m+2);
    end
    Tmx = @(x) x.*Tm(x);
    % Evaluate the integral for the right hand side
    b(m+1) = -quad(Tmx,startm,stopm);
    for n=0:N
        % Define n'th basis function
        % First and last intervals should be treated separately
        if n==0
            Tn = @(x) T(x,-xi(2),xi(2)).*(x>=0);
            startn = 0;
            stopn = xi(2);
        elseif n==N
            Tn = @(x) T(x,xi(N),2-xi(N)).*(x<=1);
            startn = xi(N);
            stopn = 1;
        else
            Tn = @(x) T(x,xi(n),xi(n+2));
            startn = xi(n);
            stopn = xi(n+2);
        end
        % Define the integrand as an anonymous function
        zfun = @(x,ksi) K(x,ksi).*Tn(ksi).*Tm(x);
        % Double integration for testing
        A(m+1,n+1) = dblquad(zfun,startm,stopm,startn,stopn);
        afun = @(x) Tm(x).*Tn(x);
        A(m+1,n+1) = A(m+1,n+1)-quad(afun,min(startm,startn),max(stopm,stopn));
        waitbar((m*N+n)/(N+1)^2,hwait)
    end
end
close(hwait)
% Solve the system of equations
a = A\b;
% Define plot points
x = linspace(0,1,100);
% Calculate the function at plot points
y = zeros(size(x));
for m=0:N
    if m==0
        y = y+a(m+1)*T(x,-xi(2),xi(2));
    elseif m==N
        y = y+a(m+1)*T(x,xi(N),2-xi(N));
    else
        y = y+a(m+1)*T(x,xi(m),xi(m+2));
    end
end

function [y,a,x,A,b]=ieqpfp(N)
% Polynomial basis, Full domain, Point matching
b = zeros(N+1,1);
A = zeros(N+1);
% Define the kernel as an anonymous function
K = @(x,ksi) sinc((x-ksi)/pi).^2/4;
% Define test points as equally spaced over [0,1]
xtest = linspace(0,1,N+1);
hwait = waitbar(0,'Calculating! Please wait...');
for m=0:N
    b(m+1) = -xtest(m+1);
    for n=0:N
        % Define the integrand at the test point
        afun = @(ksi) K(xtest(m+1),ksi).*ksi.^n;
        % Evaluate the integrand at the test point
        A(m+1,n+1) = quad(afun,0,1);
        A(m+1,n+1) = A(m+1,n+1)-xtest(m+1)^n;
    end
    waitbar(m/(N+1),hwait)
end
close(hwait)
% Solve the system of equations
a = A\b;
% Define plot points
x = linspace(0,1,100);
% Calculate the function at plot points
y = polyval(flipud(a),x);

function [y,a,x,A,b]=ieqefp(N)
% Exponential basis, Full domain, point matching
jey = sqrt(-1);
b = zeros(2*N+1,1);
A = zeros(2*N+1);
% Define the kernel as an anonymous function
K = @(x,ksi) sinc((x-ksi)/pi).^2/4;
% Define test points as equally spaced over [0,1]
xtest = linspace(0.1,0.9,2*N+1);
% xtest = rand(2*N+1,1);
hwait = waitbar(0,'Calculating! Please wait...');
for m=-N:N
    b(N+m+1) = -xtest(N+m+1);
    for n=-N:N
        % Define the integrand at the test point
        zfun = @(ksi) K(xtest(N+m+1),ksi).*exp(jey*2*pi*n*ksi);
        % Evaluate the integrand at the test point
        A(N+m+1,N+n+1) = quad(zfun,0,1);
        A(N+m+1,N+n+1) = A(N+m+1,N+n+1)-exp(jey*2*pi*n*xtest(N+m+1));
    end
    waitbar((N+m+1)/(2*N+1),hwait)
end
close(hwait)
% Solve the system of equations
a = A\b;
% Define plot points
x = linspace(0,1,100);
% Calculate the function at plot points
y = zeros(size(x));
for m=-N:N
    y = y+real(a(N+m+1)*exp(jey*2*pi*m*x));
end

function [y,a,x,A,b]=ieqrsp(N)
% rectangular pulse basis, Subdomain, Point matching
b = zeros(N,1);
A = zeros(N);
% Define the kernel as an anonymous function
K = @(x,ksi) sinc((x-ksi)/pi).^2/4;
% Define intervals
xi = [0:N]/N;
% Define test points as centers of intervals
xtest = (xi(1:end-1)+xi(2:end))/2;
hwait = waitbar(0,'Calculating! Please wait...');
for m=1:N
    b(m) = -xtest(m);
    % Define the integrand at the test point
    afun = @(ksi) K(xtest(m),ksi);
    for n=1:N
        A(m,n) = quad(afun,xi(n),xi(n+1));
        if m==n
            A(m,n) = A(m,n)-1;
        end
    end
    waitbar(m/N,hwait)
end
close(hwait)
% Solve the system of equations
a = A\b;
% Define plot points
x = linspace(0,1,100);
% Calculate the function at plot points
y = zeros(size(x));
for m=1:N
    y = y+a(m)*(x>=xi(m)).*(x<=xi(m+1));
end

function [y,a,x,A,b]=ieqtsp(N)
% Triangular pulse basis, Subdomain, Point matching
b = zeros(N+1,1);
A = zeros(N+1);
G = A;
% Define the kernel as an anonymous function
K = @(x,ksi) sinc((x-ksi)/pi).^2/4;
% Define a generic triangular function over [a,b]
T = @(x,a,b) (((x<(a+b)/2).*(x-a)+(x>=(a+b)/2).*(b-x))*2/(b-a)).*(x>=a).*(x<=b);
% Define intervals
xi = [0:N]/N;
% Define test points as end points of intervals
xtest = xi;
hwait = waitbar(0,'Calculating! Please wait...');
for m=0:N
    % Define m'th basis function
    % First and last intervals should be treated separately
    if m==0
        Tm = @(x) T(x,-xi(2),xi(2)).*(x>=0);
        startm = 0;
        stopm = xi(2);
    elseif m==N
        Tm = @(x) T(x,xi(N),2-xi(N)).*(x<=1);
        startm = xi(N);
        stopm = 1;
    else
        Tm = @(x) T(x,xi(m),xi(m+2));
        startm = xi(m);
        stopm = xi(m+2);
    end
    Tmx = @(x) x.*Tm(x);
    b(m+1) = -xtest(m+1);
    for n=0:N
        % Define n'th basis function
        % First and last intervals should be treated separately
        if n==0
            Tn = @(x) T(x,-xi(2),xi(2)).*(x>=0);
            startn = 0;
            stopn = xi(2);
        elseif n==N
            Tn = @(x) T(x,xi(N),2-xi(N)).*(x<=1);
            startn = xi(N);
            stopn = 1;
        else
            Tn = @(x) T(x,xi(n),xi(n+2));
            startn = xi(n);
            stopn = xi(n+2);
        end
        % Define the integrand as an anonymous function
        afun = @(ksi) K(xtest(m+1),ksi).*Tn(ksi).*Tm(xtest(m+1));
        A(m+1,n+1) = quad(afun,startn,stopn);
        A(m+1,n+1) = A(m+1,n+1)-Tn(xtest(m+1));
    end
    waitbar((m+1)/(N+1),hwait)
end
close(hwait)
% Solve the system of equations
a = A\b;
% Define plot points
x = linspace(0,1,100);
% Calculate the function at plot points
y = zeros(size(x));
for m=0:N
    if m==0
        y = y+a(m+1)*T(x,-xi(2),xi(2));
    elseif m==N
        y = y+a(m+1)*T(x,xi(N),2-xi(N));
    else
        y = y+a(m+1)*T(x,xi(m),xi(m+2));
    end
end