### $Id: muestreador.R, v 1.01 2004/10/22 
### 
### Muestreador de poblaciones artificiales 
### 
### Copyright 2002-2004 José A. Palazón & José F. Calvo
### <palazon@um.es> <jfcalvo@um.es> 
### 
### This file is made available under the terms of the GNU General 
### Public License, version 2, or at your option, any later version, 
### incorporated herein by reference. 
### 
### This program is distributed in the hope that it will be 
### useful, but WITHOUT ANY WARRANTY; without even the implied 
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
### PURPOSE. See the GNU General Public License for more 
### details. 
### 



fe<-100/465


lam2<-list(
sp1=matrix(c(
x1=c(337,361,28,161,221,319,133,161,213,133,321,276,10,317,374,374,149,26,270,44,131,367,400,130,388,94,28,116,314,23,224,53,389,36,324,21,65,28,119,320,47,54,432,367,163,400,393,113,13,366,94,383,14,317,17,30,39,276,20,275,267,133,27,23,169,414,280,228,397,318,46,268,137,346,99,388,132,141,302,40,73,356,58,153,390,371,141,270,63,14,327,334,25,149,116,374,363,201),
x2=c(195,270,99,312,325,438,308,330,293,407,205,155,379,446,442,284,306,158,149,99,357,34,313,404,285,418,173,352,450,98,310,421,442,348,201,404,359,408,360,443,95,362,109,16,311,132,300,51,396,263,423,282,268,196,246,398,415,133,425,55,150,233,376,415,317,92,64,325,118,207,86,58,409,330,36,125,366,253,429,417,369,40,358,401,142,436,317,140,357,176,448,310,394,399,408,129,266,306),
x3=c(4,5,5,7,4,4,7,10,8,5,6,7,5,6,6,7,9,7,9,6,6,6,6,3,9,7,8,7,9,8,4,10,5,7,4,4,5,3,6,4,4,9,10,8,10,8,9,4,5,5,6,4,6,7,6,10,4,5,7,7,5,10,6,4,9,4,11,10,7,10,4,7,7,4,4,4,6,6,9,9,10,6,3,6,8,8,6,9,9,6,11,8,7,9,9,11,7,7)
),

98,3)*fe,

sp2=matrix(c(
x3=c(363,253,368,230,110,201,228,380,331,364,104,386,216,94,237,93,96,112,98,225,352,226,363,341,262,346,68,238,234,116,90,385,373,100,350,259,344,227,131,229,104,92,225,107,81,271,84,64,81,410,222,135,97,190,224,401,238,239,90,125,95,234,74,252,120,76,367,91,74,230,341,204,343,229,225,265,244),
y3=c(96,113,327,360,55,123,360,275,75,312,216,189,192,261,354,258,249,229,92,240,208,108,232,111,322,362,65,228,66,372,336,81,84,112,93,115,238,104,346,85,394,260,229,216,211,110,72,60,118,95,98,369,224,393,225,367,108,89,238,275,97,339,244,230,78,194,99,99,243,78,330,118,132,236,128,65,368),
z3=c(6,6,5,8,6,7,6,9,5,6,8,7,10,5,9,4,4,9,9,9,10,5,9,5,7,6,6,3,7,9,6,10,7,6,7,7,6,7,4,9,4,5,9,4,7,9,9,10,6,7,9,5,6,7,7,4,9,5,8,9,4,11,4,7,3,8,9,7,4,6,4,4,6,11,4,6,6)
),
77,3)*fe,

sp3=matrix(c(
x5=c(61,253,305,227,153,262,301,167,198,128,180,229,208,86,40,102,53,226,169,248),
y5=c(171,356,321,278,422,338,363,205,405,30,340,237,385,181,33,438,137,258,449,23),
z5=c(7,11,6,8,7,6,9,6,9,9,3,4,10,4,6,9,9,4,8,7)
),
20,3)*fe,


sp4=matrix(c(
x4=c(245,302,191,32,226,135,311,421,342,79,308,316,270,259),
y4=c(149,356,26,49,51,53,338,45,145,211,337,216,216,55),
z4=c(9,10,5,4,9,9,7,5,9,11,9,9,9,8)
),
14,3)*fe,


sp5=matrix(c(
x2=c(212,238,24,85,353),
y2=c(455,451,283,367,133),
z2=c(4,11,7,8,6)
),
5,3)*fe,


n.especies=5, n.especies<-5,
nom.especies=paste(rep("Esp",n.especies),1:n.especies),
col.especies=c("olivedrab3","olivedrab4","burlywood4","orange3","brown3")
)




lam1<-list(
sp1=matrix(c(
c(109,110,112,112,113,115,119,124,125,131,134,136,138,140,141,145,150,151,159,160,171,187,190,198,200,202,204,208,212,216,217,219,220,223,238,241,245,25,261,27,270,270,271,272,272,280,282,287,287,287,290,291,292,294,300,304,312,313,314,320,327,328,333,34,342,344,349,359,364,368,372,372,389,389,401,401,405,407,413,413,418,422,428,433,442,445,445,445,446,447,52,64,69,71,77,80,82,83,86,89),
c(306,181,119,127,201,281,188,376,107,141,329,379,378,411,294,27,349,198,432,244,338,234,245,25,205,288,49,412,278,309,112,168,186,362,409,337,52,243,131,358,152,343,89,158,404,99,276,113,267,316,376,394,50,288,236,197,107,447,180,334,320,302,406,369,395,343,63,393,220,75,293,76,348,393,169,41,67,221,109,441,178,71,38,345,276,206,22,43,382,182,80,408,356,113,260,139,284,261,340,220),
c(5,8,14,11,6,5,8,12,9,12,11,10,6,12,8,11,6,10,13,12,7,10,9,12,7,8,8,14,9,12,10,13,11,12,11,10,14,6,7,6,7,11,7,13,10,14,8,6,9,9,8,7,9,5,10,14,9,8,14,11,13,12,10,5,14,5,12,5,13,9,8,10,11,6,12,10,12,5,13,11,10,9,13,6,6,7,11,10,6,7,6,14,14,9,14,9,13,8,12,8)
),
100,3)*fe,

sp2=matrix(c(
x2=c(111,114,139,141,143,15,15,151,151,151,154,167,168,184,185,188,189,209,214,228,231,241,242,247,253,262,265,276,277,295,295,306,306,310,315,320,329,33,335,340,345,381,387,45,58,75,81,83,92,94),
y2=c(342,172,359,42,151,251,425,345,355,371,131,271,26,257,268,429,272,366,414,427,360,33,213,270,404,119,132,400,198,104,32,213,263,136,330,343,312,176,350,314,338,232,93,248,397,359,142,439,234,331),
z2=c(14,5,6,5,8,5,7,8,8,10,7,5,10,5,8,9,8,14,11,9,8,8,13,12,13,6,9,11,8,11,9,11,7,11,13,9,10,13,9,5,13,12,5,10,8,12,10,10,13,10)
),
50,3)*fe,

sp3=matrix(c(
x3=c(114,127,140,18,180,185,200,234,237,247,264,286,290,306,31,322,344,346,35,360,369,38,396,414,433,441,45,70,83,99),
y3=c(403,155,195,275,399,235,178,269,203,21,308,260,284,364,89,277,112,67,281,72,433,71,255,149,215,211,145,170,107,174),
z3=c(7,8,7,5,9,14,10,5,8,8,10,6,8,5,11,14,6,14,11,9,7,9,12,11,11,7,6,7,12,9)
),
30,3)*fe,

sp4=matrix(c(
x4=c(111,131,177,179,182,190,234,310,315,322,331,344,417,430,45),
y4=c(331,258,195,231,118,164,79,118,177,110,250,173,95,400,194),
z4=c(10,14,5,5,13,7,6,5,12,13,8,7,13,8,9)
),
15,3)*fe,

sp5=matrix(c(
x5=c(175,333,364,52,90),
y5=c(81,407,360,428,27),
z5=c(11,14,6,12,11)
),
5,3)*fe,

n.especies=5, n.especies<-5,
nom.especies=paste(rep("Esp",n.especies),1:n.especies),
col.especies=c("olivedrab3","olivedrab4","burlywood4","orange3","brown3")
)


# -----------------------------------------------------------
v.lamina<-function(loc=lam1){
# -----------------------------------------------------------

par(pty="s")

plot(0,xlim=c(0,100),ylim=c(0,100),type="n",xlab=" ",ylab=" ")

for(i in 1:5){
   symbols(loc[[i]][,1], loc[[i]][,2], loc[[i]][,3],inches=F,add=T,bg=loc$col.especies[i],fg="pink")
   symbols(loc[[i]][,1], loc[[i]][,2], loc[[i]][,3],inches=F,add=T,fg=1)
}

for(i in 1:5)points( loc[[i]][,1], loc[[i]][,2],pch=".",col="white")
#for(i in 1:5)points(loc[[i]][,1:2],cex=loc[[i]][,3]/5,col=1+i,pch=1)
#rect(0,0,100,100)
}


# -----------------------------------------------------------
ajustaventana<-function(){
# -----------------------------------------------------------
plot(500.5,500.5,xlim=c(0,100),ylim=c(0,100),cex=50)
rect(0,0,100,100)
}

# Estima visual de la cobertura
# -----------------------------------------------------------
mc.visual   <- function (sep=10,loc=lam1){
# -----------------------------------------------------------
v.lamina(loc)
abline(h=seq(0,100,sep),v=seq(0,100,sep))
}


# medida de la densidad unidad de muestreo cuadrada

# -----------------------------------------------------------
md.cuadrado <- function (
# -----------------------------------------------------------
n=10,
tam=10,
loc=lam1,
matdat=F
){

v.lamina(loc)
t<-tam/2
if(length(n)>1){
   enx<-n[1]                                                                     
   eny<-n[2]                                                                     
   n<-n[1]*n[2]                                                                  
   intx<-(100-tam)/(enx)                                                               
   inty<-(100-tam)/(eny)                                                               
   x0<-runif(1,0,intx*0.95)                                                      
   y0<-runif(1,0,inty*0.95)                                                      
   mx<-rep(t+x0+(0:(enx-1)*intx),eny)                                              
   my<-sort(t+rep(y0+(0:(eny-1)*inty),enx))
}else{
mx<-runif(n,0+t,100-t)
my<-runif(n,0+t,100-t)
}

rect(mx-t,my-t,mx+t,my+t)
text(mx,my,1:n)

matrix(rep(0,(n*5)),n,5)->hc
for(m in 1:n) for(i in 1:5){
   hc[m,i]<-
     sum((loc[[i]][,1]>=mx[m]-t & loc[[i]][,1]<=mx[m]+t) & 
         (loc[[i]][,2]>=my[m]-t & loc[[i]][,2]<=my[m]+t) )
   }

#apply(md.cuadrado()$matriz,2,var)

colnames(hc)<-loc$nom.especies                                                  
rownames(hc)<-1:n 

if(matdat) print(hc)

md.cuadrado<-hc

}

# medida de la densidad unidad de muestreo circular

# -----------------------------------------------------------
md.circular <- function (
# -----------------------------------------------------------
n=10,
tam=100,
loc=lam1,
matdat=F
){

#if(length(n)>1) print ("regular")

v.lamina(loc)

t<-sqrt(tam/pi)
tam<-t*2

if(length(n)>1){
   enx<-n[1]
   eny<-n[2] 
   n<-n[1]*n[2]       
   intx<-(100-tam)/(enx)
   inty<-(100-tam)/(eny) 
   x0<-runif(1,0,intx*0.95)         
   y0<-runif(1,0,inty*0.95)           
   mx<-rep(t+x0+(0:(enx-1)*intx),eny)           
   my<-sort(t+rep(y0+(0:(eny-1)*inty),enx))
}else{
   mx<-runif(n,0+t,100-t)
   my<-runif(n,0+t,100-t)
}


symbols(mx,my,rep(t,n),inches=F,add=T,fg=1)

text(mx,my,1:n)

matrix(rep(0,n*5),n,5)->hc
for(m in 1:n) for(i in 1:5){
   hc[m,i]<-
     sum( ( (loc[[i]][,1] - mx[m])^2 +
            (loc[[i]][,2] - my[m])^2   ) < t^2) 
	 
   }

#apply(md.cuadrado()$matriz,2,var)

colnames(hc)<-loc$nom.especies                                                  
rownames(hc)<-1:n 

if(matdat) print(hc)

md.cuadrado<-hc

}





# Estima de la cobertura: lineal
# -----------------------------------------------------------
mc.lineal <- function (
n=5,
loc=lam1,
azar=T,
matdat=F
){
# -----------------------------------------------------------

v.lamina(loc)

if(azar){
  my<-runif(n,0,100)
}else{
  inty<-100/n
  y0<-runif(1,0,inty*0.95)
  my<-my<-y0+(0:(n-1))*inty
}

my<-sort(my)

segments(0,my,100,my,col=1)
text(0,my,1:n,1)

matrix(rep(0,n*5),n,5)->hc

for(m in 1:n) for(i in 1:5){
   hc[m,i]<-sum( ((loc[[i]][,3] > (dist<-abs(loc[[i]][,2]-my[m])))+ 0) *
        sqrt(abs(loc[[i]][,3]^2 - dist^2)) * 2 )
}

#hc<-t(as.matrix(apply(hc,2,sum)))
colnames(hc)<-loc$nom.especies
rownames(hc)<-1:n

if(matdat) print(hc)

mc.lineal<-hc

}


# Estima de la cobertura: puntos
# -----------------------------------------------------------
mc.puntos  <- function (
n=10,
loc=lam1,
matdat=F,
tvp=3
){
# -----------------------------------------------------------

v.lamina(loc)

if(length(n)>1) {
  enx<-n[1]
  eny<-n[2]
  n<-n[1]*n[2]
  intx<-100/(enx)
  inty<-100/(eny)
  x0<-runif(1,0,intx*0.95)
  y0<-runif(1,0,inty*0.95)
  mx<-rep(x0+(0:(enx-1)*intx),eny) 
  my<-sort(rep(y0+(0:(eny-1)*inty),enx))
} else {
  mx<-runif(n,0,100)
  my<-runif(n,0,100)
}

#points(mx,my,col=1,pch=10,cex=tvp)
 points(mx,my,col="purple",pch=10,cex=tvp)
#text(mx,my,labels=1:n,col="purple",pos=3)

matrix(rep(0,n*5),n,5)->hc
for(m in 1:n) for(i in 1:5){
   hc[m,i]<-sum(
        (loc[[i]][,1] - mx[m])^2 + 
	(loc[[i]][,2] - my[m])^2 <= 
	 loc[[i]][,3]^2)
}

hc<-t(as.matrix(apply(hc,2,sum)))
colnames(hc)<-loc$nom.especies
rownames(hc)<-""
if(matdat) print(hc)
md.cuadrado<-hc
}

# Cobertura "oficial"
# 0.1505 0.0691 0.0386 0.0205 0.0092

# Cobertura real
# 0.145409 0.06682014 0.03731118 0.1983246 0.008979092 


# -----------------------------------------------------------
veoplot <-function(){
# -----------------------------------------------------------

plotF<-T
if(!par()$ask){par(ask<-T);plotF<-F}

as.vector(colors())[1:657]->colores
matrix(c(rep(1:30,22),sort(rep(1:22,30))),660)->x
rownames(x)<-c(colores,rep("white",3))
plot(x,col=rownames(x),cex=4,pch=15,xlab="",ylab="")
text(x[,1],x[,2],labels=1:660,cex=1,col=1)



if(!plotF)par(ask=F)
}

# -----------------------------------------------------------
valores.reales<-function(loc=lam1){
# -----------------------------------------------------------
Nom<-Tam<-Cob<-Col<-0
des<-data.frame(Nom,Tam,Cob,Col)
for (i in 1:5){
des[i,1]<-loc$nom.especies[i]
des[i,2]<-nrow(loc[[i]])
des[i,3]<-round(sum(loc[[i]][,3]^2*3)/10000,4)
des[i,4]<-loc$col.especies[i]
}
colnames(des)<-c("Nombre","Tamaño","Cobertura","Color")
(valores.reales<-des)
}

# -----------------------------------------------------------
cobertura.real<-function(loc=lam1){
# -----------------------------------------------------------
des<-NULL
for (i in 1:5) des<- c(des, round(sum(loc[[i]][,3]^2*3)/10000,4),)
des<-matrix(des,1)
rownames(des)<-" "
colnames(des)<-loc$nom.especies
(densidad.real<-des)
}

# -----------------------------------------------------------
densidad.real<-function(loc=lam1){
# -----------------------------------------------------------
des<-NULL
for (i in 1:5) des<- c(des,nrow(loc[[i]]))
des<-matrix(des,1)
rownames(des)<-" "
colnames(des)<-loc$nom.especies
(densidad.real<-des)
}


# -----------------------------------------------------------
l.especies <-function(loc=lam1){
# -----------------------------------------------------------
vactual<-dev.cur()
x11(width=5.25,height=6, pointsize = 12)
n<-(1:5)*100
x<-n/n
plot(x,-n,col=loc$col.especies,type="n",xlim=c(0.5,1.2),ylim=c(1,-max(n)*1.1),axes=F,ann=F,xlab=paste("Lámina ",loc))
text(x,-n,loc$nom.especies,1,cex=3)
symbols(x+0.07,-n,x/20,inches=F,add=T,bg=loc$col.especies,fg=1)
dev.set(vactual)
#par(ask=T)
#plot(1)
#par(ask=F)
#dev.off()
}


### MEDIDAS DE DISTANCIA PARA ESTIMAS DE DENSIDAD
## Método del individuo más próximo

# -----------------------------------------------------------
mdist.ind <- function(
n=10,
tam=10,
loc=lam1,
matdat=F,
lines=T,
tvp=2
){
# -----------------------------------------------------------

v.lamina(loc)

mx<-runif(n,10,90)
my<-runif(n,10,90)

points(mx,my,col="purple",pch=10,cex=tvp)
#text(mx,my,labels=1:n,col="purple",pos=3)

matrix(1:(n*5),n,5)->dist

for(i in 1:5) {
  if (lines) {
   matrix(999,n,length(loc[[i]][,1]))->vec
      for(m in 1:n) {
         for(j in 1:length(loc[[i]][,1])) {vec [m,j]<-  (sqrt ( (loc[[i]][j,1] - mx[m]) ^2 + (loc[[i]][j,2] - my[m])^2))}
      }
      for(m in 1:n) {
         for(j in 1:length(loc[[i]][,1])){
            if ( vec[m,j] == min (vec[m,] ) ) lines( c ( mx[m] , loc[[i]][j,1] ) , c ( my[m] , loc[[i]][j,2] ) )
         }
      }
   }
   for(m in 1:n) {dist[m,i]<-  min(sqrt ( (loc[[i]][,1] - mx[m]) ^2 + (loc[[i]][,2] - my[m])^2))}
}

rownames(dist)<-1:n
colnames(dist)<-loc$nom.especies

if(matdat) list (distancias=dist,densidad= n / apply((pi*dist^2),2,sum))
else mdist.ind= n / apply((pi*dist^2),2,sum)

}



## Método del vecino más próximo

# -----------------------------------------------------------
mdist.vec <- function(
n=10,
tam=10,
loc=lam1,
matdat=F,
lines=T,
tvp=2
){
# -----------------------------------------------------------

v.lamina(loc)

mx<-runif(n,10,90)
my<-runif(n,10,90)
ix<-runif(n,0,100)
iy<-runif(n,0,100)

points(mx,my,col="purple",pch=10,cex=tvp)
#text(mx,my,labels=1:n,col="purple",pos=3)

matrix(1:(n*5),n,5)->dist

for(i in 1:5) {
   matrix(999,n)->ix
   matrix(999,n)->iy
   matrix(999,n,length(loc[[i]][,1]))->vec
    for(m in 1:n) {
        for(j in 1:length(loc[[i]][,1])) {vec[m,j]<-(sqrt((loc[[i]][j,1]-mx[m])^2+(loc[[i]][j,2]-my[m])^2))}
   }
   for(m in 1:n) {
        for(j in 1:length(loc[[i]][,1])){
         if ( vec[m,j] == min (vec[m,] ) ) {
             if (lines) {lines( c ( mx[m] , loc[[i]][j,1] ) , c ( my[m] , loc[[i]][j,2] ) )}
            ix[m] <- loc[[i]][j,1]
            iy[m] <- loc[[i]][j,2] 
            if (lines) {points(ix,iy,col=1,pch=3,cex=tvp)}
         }
        }
   }

   matrix(999,n,length(loc[[i]][,1]))->vec1
    for(m in 1:n) {
        for(j in 1:length(loc[[i]][,1])) {vec1[m,j]<-(sqrt((loc[[i]][j,1]-ix[m])^2+(loc[[i]][j,2]-iy[m])^2))}
   }
   if (lines) {
      for(m in 1:n) {
           for(j in 1:length(loc[[i]][,1])){
            if ( vec1[m,j] == min ( (vec1[m,]==0)*99999+vec1[m,] ) ) {
                lines( c ( ix[m] , loc[[i]][j,1] ) , c ( iy[m] , loc[[i]][j,2] ) , col="red")
            }
           }
      }
   }

   for(m in 1:n) {
      dist[m,i]<-  min(((sqrt ( (loc[[i]][,1] - ix[m]) ^2 + (loc[[i]][,2] - iy[m])^2)) == 0)* 99999 + 
                  (sqrt ( (loc[[i]][,1] - ix[m]) ^2 + (loc[[i]][,2] - iy[m])^2)) )
   }
}

rownames(dist)<-1:n
colnames(dist)<-loc$nom.especies

if(matdat) list (distancias=dist,densidad= (2*n)/apply((pi*dist^2),2,sum))
else mdist.vec=(2*n)/apply((pi*dist^2),2,sum)

}


## Método de los cuadrantes

# -----------------------------------------------------------
mdist.4 <- function(
n=10,
tam=10,
loc=lam1,
matdat=F,
lines=T,
tvp=16
){
# -----------------------------------------------------------

v.lamina(loc)

mx<-runif(n,10,90)
my<-runif(n,10,90)

points(mx,my,col="purple",pch=3,cex=tvp)
#text(mx,my,labels=1:n,col="purple",pos=3.5,offset=0.4,cex=2)

matrix(1:(n*5),n,5)->dist

for(i in 1:5) {
    matrix(999,n,length(loc[[i]][,1]))->vec1
    matrix(999,n,length(loc[[i]][,1]))->vec2
    matrix(999,n,length(loc[[i]][,1]))->vec3
    matrix(999,n,length(loc[[i]][,1]))->vec4
    for(m in 1:n) {
      for(j in 1:length(loc[[i]][,1])) {
        if ( (loc[[i]][j,1] >= mx[m]) & (loc[[i]][j,2] >= my[m]) ) {
          vec1[m,j]<-  (sqrt ( (loc[[i]][j,1] - mx[m]) ^2 + (loc[[i]][j,2] - my[m])^2))
        }
        if ( (loc[[i]][j,1] >= mx[m]) & (loc[[i]][j,2] < my[m]) ) {
          vec2[m,j]<-  (sqrt ( (loc[[i]][j,1] - mx[m]) ^2 + (loc[[i]][j,2] - my[m])^2))
        }
        if ( (loc[[i]][j,1] < mx[m]) & (loc[[i]][j,2] < my[m]) ) {
          vec3[m,j]<-  (sqrt ( (loc[[i]][j,1] - mx[m]) ^2 + (loc[[i]][j,2] - my[m])^2))
        }
        if ( (loc[[i]][j,1] < mx[m]) & (loc[[i]][j,2] >= my[m]) ) {
          vec4[m,j]<-  (sqrt ( (loc[[i]][j,1] - mx[m]) ^2 + (loc[[i]][j,2] - my[m])^2))
        }
      }
    }
  if (lines) {
    for(m in 1:n) {
      for(j in 1:length(loc[[i]][,1])){
        if ( (vec1[m,j] == min (vec1[m,])) & (min(vec1[m,]) <999) ) {
        lines( c ( mx[m] , loc[[i]][j,1] ) , c ( my[m] , loc[[i]][j,2] ) )}
        if ( (vec2[m,j] == min (vec2[m,])) & (min(vec2[m,]) <999) ) {
        lines( c ( mx[m] , loc[[i]][j,1] ) , c ( my[m] , loc[[i]][j,2] ) )}
        if ( (vec3[m,j] == min (vec3[m,])) & (min(vec3[m,]) <999) ) {
        lines( c ( mx[m] , loc[[i]][j,1] ) , c ( my[m] , loc[[i]][j,2] ) )}
        if ( (vec4[m,j] == min (vec4[m,])) & (min(vec4[m,]) <999) ) {
        lines( c ( mx[m] , loc[[i]][j,1] ) , c ( my[m] , loc[[i]][j,2] ) )}
      }
    }
  }

    for(m in 1:n) {
    c(min(vec1[m,]),min(vec2[m,]),min(vec3[m,]),min(vec4[m,]))->temp
    NA->temp[temp==999]
    temp[!is.na(temp)]->temp
    dist[m,i]<-  sum(temp^2)
  }
}

rownames(dist)<-1:n
colnames(dist)<-loc$nom.especies

if(matdat) list (sumdistcuad=dist,densidad= (4*(4*n-1))/(pi*apply((dist),2,sum)))
else mdist.4=(4*(4*n-1))/(pi*apply((dist),2,sum))
}