function [integr] = int_CompSimpsonRule(func,model)
% Intergration with Composite Simpson's Rule
% To get an integral quantity on the Gauss grid without losing of accuracy.
% func: parameter on Gauss grid to integrate.
% intgr: result of integration on gauss grid.

% Define number of intervals to integrate
nint = 6;
% Gauss grid
rhogauss = model.rgrid.rhogauss;
% Number of Gauss point
ngauss = model.rgrid.nrhogauss;
% Initialization
integr = zeros(ngauss,1);

% Grid with intermediate points
jj = 1;
for ii=1:ngauss-1
    rhofull(jj:jj+nint) =  linspace(rhogauss(ii),rhogauss(ii+1),nint+1);
    jj = jj+nint;
end
% Get function values for the new grid with intermediate points
ff = interp1(rhogauss,func,rhofull,'spline');
    
% Weights
weights = [41. 216. 27. 272. 27. 216. 21.];
for ii=2:ngauss    
    % Normalization
    norm = (rhogauss(ii) - rhogauss(ii-1))/sum(weights);
    % Integration
    int0 = 1 + (ii-2)*nint;
    inte = nint + 1 + (ii-2)*nint;
    x_sum = norm*sum(weights.*ff(int0:inte));
    integr(ii) = integr(ii-1) + x_sum;
end

return
% % Another way: need even numbe of integrals
% % Integral calculation
% for ii=2:ngauss
%     % Define function values on Gauss interval
%     fa = func(ii-1);
%     fb = func(ii);    
%     % Set spacing step
%     h = (rhogauss(ii) - rhogauss(ii-1))/nint;
%     % Initialization for even and odd parts of integration
%     x_even = 0.0;
%     x_odd = 0.0;
%     % Integration
%     for jj=1:nint-1
%         if mod(jj,2) % even jj            
%             x_odd = x_odd + ff(jj+1 + (ii-2)*nint);
%         else % odd jj
%             x_even = x_even + ff(jj+1 + (ii-2)*nint);
%         end        
%     end
%     x_int = h*(fa + 2.*x_even + 4.*x_odd + fb)./3.;
%     integr(ii) = integr(ii-1) + x_int;
% end 
    