function intgauss = int_gauss(valgauss,model,varargin)
% integrate quantity on gauss points and return value on gauss points
% directly.
% intgauss = int_gauss(valgauss,model)
% intgauss = int_gauss(valgauss,model,irhogauss)
% 
%  valgauss: value on gauss points
%  model: model structx
%  irhogauss (optional): indices on which to return values

if nargin==3
  irhogauss = varargin{1};
else
  irhogauss = true(model.rgrid.nrhogauss,1);
end

  rho = model.rgrid.rho;
  rhogauss = model.rgrid.rhogauss;
  ngauss = model.rgrid.ngauss;
  wgauss = model.rgrid.wgauss;
  
  ndat = size(valgauss,2);
  intgaussall = zeros(ngauss,ndat);
  for idat=1:ndat % for several columns of data
    % standard Gaussian integration for rho grid
    intall = cumsum(wgauss.*valgauss(:,idat));
    intrho = [0;intall(ngauss:ngauss:end)];
    
    % fill bits in between for rhogauss points by cubic interpolation
    for ii=1:numel(rho)-1
      iig = (ii-1)*ngauss+(1:ngauss); % indices into full rhogauss grid
      
      % if any overlap with output grid point request
      if any(irhogauss(iig))
        rhog = rhogauss(iig); % rho points for this interval
        valg = valgauss(iig,idat); % values
        [intrhog] = interp_cubic(rhog,valg,[rho(ii);rhog],-1); % integrals
        intgaussall(iig,idat) = intrho(ii)+intrhog(2:end);
      end
    end
  end
  intgauss = intgaussall(irhogauss,:);
end