gene.div <- function(Data=Loci, Pop, Ind.Nr.Pop)
{
# Calculate allele frequency
# 1) number of alleles in each population at each locus
Nr.ind.Pop <- tapply(Pop, Pop, length)
allele.frequency  <- list()
test <- list()
for(i in 1:ncol(Loci))
	{
	b <- list(Loci[,i], Pop)
	allele.frequency [[i]] <- tapply(Loci[,i], b, length)
	}
allele.frequency 
# 2) calculate allele frequency
for(i in 1:ncol(Loci))
	{
	for(j in 1:length(Nr.ind.Pop))
		{
		allele.frequency[[i]][,j] <-allele.frequency [[i]][,j]/Nr.ind.Pop[j]
		}
	}
allele.frequency	# allele frequency per locus and population 	
# 3) replace the NAs with zero
b <- list()
for(i in 1:ncol(Loci))
	{
	b <- allele.frequency[[i]]
	b[is.na(b)]=0
	allele.frequency[[i]] <- b
	}
test <-allele.frequency 
test # allele.frequency at each locus, list.
# allele.frequency: allele frequency for each locus
# for haploid data, for diploid data only if the two alleles are in rows

### Calculate Jx per locus und Pop
C.nr <-(length(unique(Pop))^2 - length(unique(Pop)))/2;C.nr

Jx <- matrix(NA,ncol(Loci),length(unique(Pop)),dimnames=list(names(Loci),unique(Pop)) );Jx
for(x in 1:ncol(Loci))
{
test1 <-test[[x]]

# Calculate Jx ("average homozygosity within pop")
n <- 1
b.1 <- list()
for(j in 1:length(unique(Pop)))
	{
	b.1[[n]] <- (test1[,j]^2) # Quadrat of allele frequencies xij2 
	n <- n+1
	}
b.1
f <- list()
for(i in 1:length(unique(Pop)))
	{
	f[[i]] <- (sum(b.1[[i]])) # Sum products b into list a[[i]]
	}
f
# Write f[[i]] into matrix Jx

n <- 1
for(i in 1:length(unique(Pop)))
	{
	Jx[x,i] <- f[[n]]
	n <- n+1
	}
Jx
}
Jx

# Calculate heterozygosity from homozygosity Jx:
Jx <- 1-Jx

################First route
#unbiased estimate of "heterozygosity" of each locus, h:
#multiply Jx with n/(n-1),
#where n is the number of thalli investigated per Population.

#h <- matrix(NA,ncol(Loci),length(unique(Pop)),dimnames=list(names(Loci),unique(Pop)) )
#for(i in 1:length(unique(Pop)))
#	{
#	for(j in 1:ncol(Loci))
#		{
#		h[j,i] <- (Jx[j,i]*Ind.Nr.Pop[i])/(Ind.Nr.Pop[i]-1)
#		}
#	}
#h

##Unbiased gene diversity H
#Sum of h over r loci:
#sum.h <-rep(NA,length(unique(Pop)))
#n <- 1
#for(i in 1:length(unique(Pop)))
#	{
#	sum.h[n] <- sum(h[,i])
#	n <- n+1
#	}
#sum.h

##########################
#Alternative: calculate sum hk, multiply with n/(n-1)
h <- matrix(NA,ncol(Loci),length(unique(Pop)),dimnames=list(names(Loci),unique(Pop)) )
for(i in 1:length(unique(Pop)))
	{
	for(j in 1:ncol(Loci))
		{
		h[j,i] <- Jx[j,i]
		}
	}
h
#Sum of h over loci
sum.h <-rep(NA,length(unique(Pop)))
n <- 1
for(i in 1:length(unique(Pop)))
	{
	sum.h[n] <- sum(h[,i])
	n <- n+1
	}
sum.h

#Multiply sum.h with n/(n-1)
a.1 <- c()
for(i in 1:length(unique(Pop)))
	{
	a.1[i] <- (sum.h[i]*Ind.Nr.Pop[i])/(Ind.Nr.Pop[i]-1)
	}

#Divide by number of loci
GeneDiv <- matrix(NA,length(unique(Pop)),1, dimnames=list(unique(Pop),c("H")))
for(i in 1:length(unique(Pop)))
	{
	GeneDiv[i] <- a.1[i]/ncol(Loci)
	}
GeneDiv
GeneDiv[is.nan(GeneDiv)]<-0

#########################

#dimnames=list(c("H"),unique(Loci))
#Teile sum.h durch Anzahl Loci
#GeneDiv <- matrix(NA,length(unique(Pop)),1, dimnames=list(unique(Pop),c("H")))
#for(i in 1:length(unique(Pop)))
#	{
#	GeneDiv[i] <- sum.h[i]/ncol(Loci)
#	}
#GeneDiv
Result <- list()
Result$GeneDiversity.unbiased <- GeneDiv
Result
}