
function V=eval_mod_op(deg,ab,a,b,t,c,type)

%--------------------------------------------------------------------------
% OBJECT.
%----------
% 
% 
%--------------------------------------------------------------------------
% FUNCTIONS USED.
%----------------
%
%--------------------------------------------------------------------------
%%
%% Copyright (C) 2010. Alvise Sommariva, Marco Vianello.
%%
%% This program is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation; either version 2 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
%%
%% Author:  Alvise Sommariva <alvise@euler.math.unipd.it>
%%          Marco Vianello   <marcov@euler.math.unipd.it>
%%
%% Date: October 18, 2010
%%
%--------------------------------------------------------------------------

if nargin < 3
    a=0; b=0;
end

L=size(t,1);
aa=ab(:,1);
bb=ab(:,2);

switch type
case 1 % MONIC JACOBI POLYNOMIALS.
    V=[zeros(L,1) ones(L,1)];
    cn=ones(size(ab,1),1);
    co=ones(size(ab,1),1);
case 2 % CLASSICAL JACOBI (SEE WOLFRAM HOMESITE)
    V=[zeros(L,1) ones(L,1)];
    [cn,co]=c_standard(deg,a,b);
    
case 3 % ORTHONORMAL JACOBI POLYNOMIALS.
    V=[zeros(L,1) (1/sqrt(ab(1,2)))*ones(L,1)]; % "ab(1,2)" IS "int_-1^1 w(x) dx".
    [cn,co]=c_orthnrm(deg,a,b);
    
case 4 % ORTHONORMAL JACOBI POLYNOMIALS USED BY DUBINER ORTHONORMAL BASIS OVER THE SYMPLEX.
    V=[zeros(L,1) (sqrt(8)/sqrt(ab(1,2)))*ones(L,1)]; % "ab(1,2)" IS "int_-1^1 w(x) dx".
    [cn,co]=c_dub(deg,a,b);
end

for ii=1:deg
    Vloc=cn(ii)*(t-c.*aa(ii)).*V(:,ii+1)-(c.^2)*co(ii)*bb(ii).*V(:,ii);
    V=[V Vloc];
end

V=V(:,2:end);


function [cn,co]=c_dub(ade,a,b)

[cn_orth,co_orth]=c_orthnrm(ade,a,b);
cn=2*cn_orth;
co=[1; cn(2:end,:).*cn(1:end-1,:)];

function [cn,co]=c_orthnrm(ade,a,b)

[cn_st,co_st]=c_standard(ade,a,b);

ab=r_jacobi(1,a,b);
an0=ab(1,2);

n=(1:ade)';

if length(n) > 0
    if abs(b) > 0 
        anloc=(2^(a+b+1)).*gamma(n+a+1).*gamma(n+b+1)./((2*n+a+b+1).*gamma(n+a+b+1).*gamma(n+1));
    else
        anloc=2^(a+b+1)./(2*n+a+b+1);
    end
else
    anloc=[];
end

an=[an0; anloc];

an=an.^(-1/2);
an_rat=an(2:end,:)./an(1:end-1,:);
cn=cn_st.*an_rat;
co=[1; cn(2:end,:).*cn(1:end-1,:)];





function [cn,co]=c_standard(deg,a,b)

s=a+b;

% INITIALIZATION TO AVOID SOME POSSIBLE OVERFLOW PROBLEMS FOR NEGATIVE s.
if (deg > 0)
    cn1=0.5*gamma(2+s+1)/gamma(s+2);
else
    cn1=[];
end

n=2:deg; n=n';

if (abs(s-floor(s)) > 0) 
    cnloc=(1./(2*n)).*(gamma(2*n+s+1)./gamma(2*n+s-1)).*(gamma(n+s)./gamma(n+s+1));
else
    cnloc=(1./(2*n)).*(2*n+s).*(2*n+s-1)./(n+s);
end
 
cn=[cn1; cnloc];

co=[1; cn(2:end,:).*cn(1:end-1,:)];


