function [ Q ] = integ( n,x,y,un,un1,xi1,xi2,a, mu,we,wn,ve,it,ht)
% computes the integral to use
%   on the wright hand side of the equations system
% using a gaussian guadrature
% un is the array with the solution values, ti is the present time
%-dimension n
% x,y are  given points
%we=weights(n); wn -gaussian nodes
% ve is the propagation speed
% if ve=0 we consider the non-delay case
% un1 is the long vector containing the last layer of un 
h=2/n;
Q=0;
%%% it is the index of the solution at time ti
for i=1:n
    for j=1:n
       ind=(i-1)*n+j;
       u=wn(i);
       v=wn(j);
%%%%%%%%%     nondelay case
       if(ve==0)
           umod=un1(ind);
           Q=Q+we(i)*we(j)*ker(x,y,u,v,xi1,xi2,a)*S(mu,umod);
       else
 %%%%%%%%%%%%%%%%%%%%%%  delay case (ve different from 0 )          
       %%% dist is the distance between (x,y) and (u,v)
       dist=sqrt((x-u)^2 + (y-v)^2);
       del= dist/ve;
       %%%%%%%%%%%%%% compute umod
       %%%% umod is the modified value of un, taking the delay into account
       %%%%%%  jd is the shift between ti and the delayed time
       jd= floor(del/ht);
       if (jd==0)  %%%%%%%%%%%%%%%%   delay <=ht
           umod1=un1(ind);
           ii= max(1,it-1);
           umod2=un(ii,ind); %%%% initial condition
       else  %%%% delay >ht
           ii=max(1,it-jd);
           umod1=un(ii,ind);
           ii=max(1,it-jd-1);
           umod2= un(ii,ind);
       end
       %%%%%% s is the coefficient for linear interpolation
       s=del/ht-jd;
       %%%%% linear interpolation
       umod=s*umod1+(1-s)*umod2;
       Q=Q+we(i)*we(j)*ker(x,y,u,v,xi1,xi2,a)*S(mu,umod);
       end
    end
end
end

