% Specular points in two dimensions 
%This script 
%enters the parameters
%computes the parameter m_2 of formula (48)
%calls the function m2_inner
%makes the picture 
% 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 script 
%Tlse 3 2009

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

V= 0.0001*[114 0
     0 81];
 
 A=[ 0 , 0.5, 0
    0.5,0,0
    0, 0,-1];
sigmad = sqrtm(sigma);
adiag = sigmad*A*sigmad;
[v,d]=eig(adiag);
global  delt1 delt2 delt3; % to pass them to m2_inner
delt1= d(1,1);
delt2= d(2,2);
delt3= d(3,3);

%computation of the integral using quad
 
  a = 0.00001; %stands  for zero
  b= 100000;% stands for infinity
 int = quad('m2_inner',a,b);
 1000*int
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function[res]=m2_inner(t)

%this function computes the intergrand for m2

t2= t.*t;
global  delt1 delt2 delt3 
A1= (1+ 4* delt1*delt1*t2);
A1=A1.^(-0.5);
 B1 = sqrt((1-A1)/(1+A1));
 
 A2= (1+ 4* delt2*delt2.*t.*t).^(-0.5);
 B2 = sqrt((1-A2)/(1+A2));
 
 A3= (1+ 4* delt3*delt3.*t.*t).^(-0.5);
 B3 = sqrt((1-A3)./(1+A3));
  res = 2 ^(-1.5) * A1.* sqrt(1+A1).*A2.* sqrt(1+A2).*A3.* sqrt(1+A3);
   res = 1 - res.* (1-B1.*B2-B2.*B3-B1.*B3);
   res = res./(t.*t);
 