############################
# Fonction qui renvoie la distance terrestre en kilometres entre les n1 points 
# de coordonees long1[i] (longitude) and lat1[i] (latitude) et les n2 points 
# de coordonees long2[j] and lat2[j]. AUne matrice n1xn2.
lldist<- function(long1,lat1,long2,lat2){
pi180=pi/180;
earth_radius<-6378.137;
m=length(long1);
plat1=lat1*pi180;
distl = matrix(0,length(long1),length(long2))
for (i in 1:length(long2)){
 dlon=(matrix(long2[i],m,1)-long1)*pi180;
 plat2=matrix(lat2[i],m,1)*pi180;
 dlat=plat2-plat1;
 a=(sin(dlat/2))^2 + cos(plat1)*cos(plat2)*(sin(dlon/2))^2;
 angles=2 * atan2(sqrt(a),sqrt(1-a));
 distl[,i]=earth_radius*angles;
}
return(distl)
}
############################
# Fonction qui fait un trace couleur.
# Valeur z[i] au point (x[i],y[i]).
plotpoints = function(x,y,z){
L = 64 # Nombre de couleurs
cr=rainbow(L)
cr = cr[length(cr):1] # retourne le vecteur c
znor=(z-min(z))/(max(z)-min(z)) # normalisation sur [0,1]
ind= round( znor*(L-20))+20 # numero de couleur. Les 20 premieres ne sont pas utilisees
plot(x,y,col=cr[ind],pch=19)
}
############ Estimation non paramétrique de la fonction ############   
############ covariance=f(distance) ############
#
# En entrée: 
#  Dist: Le tableau carré des distances, 
#  y: Les valeurs
#  N: Le nombre de bin (defaut=12)
# En sortie
#  r: covariance
#  d: distance.
#
# r[i] est la covariance empirique de y, pour des paires de points 
# dont la distance est comprise entre q[i-1] et q[i], où q[i] est i-eme 
# quantile de Dist. d[i]= (q[i-1]+q[i])/2.
# (Le cas particulier i=1 correspond à d[1]=0, et donc r[1] est la variance.) 
# On fait donc des moyennes sur un nombre a peu près fixe d'échantillons

covariog=function(Dist,y,N=12){ 
  q=quantile(Dist, probs = seq(0, 1, 1./(N-1)))
  lq=length(q);
  r1=(1:lq); d1=r1; s1=r1;
  R1=outer(y,y);
  for (i in 2:lq){ 
    x=R1[(Dist>=q[i-1])&(Dist<q[i])];
    r1[i]=mean(x); 
    d1[i]=(q[i-1]+q[i])/2;
    s1[i]=sqrt(var(x)/length(x));
  }
  x=R1[Dist==0]; 
  r1[1]=mean(x);
  d1[1]=0;
  return(list(d=d1,r=r1))
}