function [I]=gauss_jacobi_integration_I(N,alpha,beta,a,b,f)
% 
% This function computes the integral of function f with gaussian formulas.
% 
% INPUTS:
% N - number of nodes
% alpha,beta: w(x)=((1-x)^alpha)*((1+x)^beta)
% a,b - if alpha=0, beta=0, they are the integration extremes, 
%    else a=-1, b=1.
% f - integrand with weight function w.
% 
% OUTPUTS:
% I - integral of f on [a,b]

ab=r_jacobi(N,alpha,beta); % recursive terms
%printf('ab=%f\n',ab)
xw=gauss(N,ab); % matrix of nodes/weights
%printf('xw=%f\n',xw)
x=xw(:,1);
w=xw(:,2);%Gauss-Jacobi nodes/weights on [-1,1]
%printf('x=%f \nw=%f\n',x,w)
if (alpha == 0) & (beta == 0) & ((a~=-1)|(b~=1))
 x=((a+b)/2)+((b-a)/2)*x;%Gauss-Jacobi nodes on [a,b]
 w=((b-a)/2)*w;%Gauss-Jacobi weights on [a,b]
end
fx=feval(f,x);%function evaluating
I=w'*fx;%integral