function[res]=spec2d(x,y,sigma,V,k)

%this function correspond to section 3.5
% it calls function spec2d-inner that gives the integrand
%It compute the density  of specular point in two dimensions 
% it refers to the paper "Some applications of Rice formulas to waves"
% by Azaïs Leon and Wschebor
% Jean-Marc Azaïs  is the author of this function
%MVD 2 2009
%(x,y)are the space variables 
% sigma is the 3x3 variance matrix of the paper
% V  is the 2x2 variance matrix of the gradient which is not assumed to be
% diagonal. 
%k is the parameter of specular points

%k=0.1; x=1; y=1;% mise au point 

global  delt1 delt2 delt3 s1 s2 s3 k2 d int;
k2 =k*k;
a= 0.001;  % limit of integration, stands for zero
b=50000; %limit of integration, stands for infinity 

%density
u= [k*x,k*y];
Vi =inv(V);
dd = (1/(2*pi))* (det(V) ) ^(-1/2)* exp(-0.5*((u*Vi*u')));


% preliminary calculations for the determinant
A=[ 0 , 0.5, 0
    0.5,0,0
    0, 0,-1];
sigmad = sqrtm(sigma);
adiag = sigmad*A*sigmad;
[v,d]=eig(adiag);

delt1= d(1,1);
delt2= d(2,2);
delt3= d(3,3);

SSS = sigmad * v';
s1= SSS(1,1) + SSS(2,1);
s2= SSS(1,2) + SSS(2,2);
s3= SSS(1,3) + SSS(2,3);




 int = quad('spec2d_inner',a,b);
 
 res = 2/pi *int*dd;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function[z]=spec2d_inner(t);
     
 % the function called by quad 
 % uses delt1 delt2 delt3 s1 s2 s3  kas global variables 
  %it refers to the paper "Some applications of Rice formulas to waves"
% by Azaïs Leon and Wschebor
% Jean-Marc Azaïs  is the author of this function

% 
global  delt1 delt2 delt3 s1 s2 s3 k2 ;
 
 t2=t.*t;
z= 1./(1+t.*t);
 d1 = exp(-k2*t2*0.5 .*  s1*s1/( 1 + 4*delt1^2*t2));
 d2 = exp(-k2*t2*0.5 .*  s2*s2/( 1 + 4*delt2^2*t2));
  d3 = exp(-k2*t2*0.5 .*  s3*s3/( 1 + 4*delt3^2*t2));
  
  
  phi1 = 0.5 * atan( 2*delt1*t);
   phi2 = 0.5 * atan( 2*delt2*t);
    phi3 = 0.5 * atan( 2*delt2*t);
    
    
    
    psi1 = 1/3 - t2 *s1*s1*delt1 ./ ( 1+4* delt1* delt1*t2);
    psi2 = 1/3 - t2 *s2*s2*delt2 ./ ( 1+4* delt2* delt2*t2);
    psi3 = 1/3 - t2 *s3*s3*delt3 ./ ( 1+4* delt3* delt3*t2);
    
    
    h =  d1./sqrt( 1+4* delt1* delt1*t2);
    h = h.*d2./sqrt( 1+4* delt2* delt2*t2);
    h = h.* d3./sqrt( 1+4* delt3* delt3*t2);
    h=h.* cos( (phi1+ k2*t.*psi1)+(phi2+ k2*t.*psi2)+(phi3+ k2*t.*psi3));
    
    
    z= (1-h) ./t2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Exemple of utilization

%sigma= [ 9 ,  3 ,  0
%    3  , 11,  0
 %   0,  0, 3];
%sigma= 0.0001*sigma; 

% V= 0.0001*[114 0
 %    0 81];
 %xx= -40:40;
 %yy=-40:40 ;
 %%
 %for i=1:81
 %  for j=1:81
 %      x= xx(i); y=yy(j);
%z= spec2d(x,y,sigma,V,0.01);
%xxx(i,j) = x;
%yyy(i,j) = y;
%zz(i,j) = z; 

  % end;
