mendel=function(loci, modelo, N=F)#Cria a função "mendel" com os argumentos loci, modelo e N
{
  install.packages("gplots")
  library(gplots)#Carregar pacote gplots para balloon plot
  #VERIFICANDO OS PARAMETROS
  if (modelo!= "t" && modelo != "s") #Se o argumento "modelo" for diferente de "s" ou "t"
  { 
    stop (" o argumento modelo deve ser 't' ou 's'") # Interrompe a funcao e exibe  a mensagem ao usuario
  }
  if (N!=FALSE)#Se N for diferente de FALSE
  {
    if (N != round (N) | N <= 0) #Se N nao for um numero iteiro e menor que 0
    {
      stop ("N deve ser um numero inteiro maior que 0.") #Interrompe a funcao e exibe a mensagem ao usuario
    }
  }
  if (loci !=1 && loci!= 2) #Se o argumento loci for diferente de 1 e diferente de 2
  {
    stop (" O argumento loci deve ser 1 ou 2") #Interrompe a funcao e exibe mensagem ao usuario.
  }
  #CRIANDO OBJETOS COMUNS PARA OS DOIS MODELOS
  n.gametas=2^loci#Para calcular a quantidade de gametas possiveis
  ger1a.1=c(1,1)#Genótipo do progenitor 1 da primeira geração de 1 locus (números se referindo à "AA") 
  ger1b.1=c(3,3)#Genótipo do progenitor 2 da primeira geração de 1 locus (números se referindo à "aa") 
  ger1a.2=c(10,10)#Genótipo do progenitor 1 da primeira geração do segundo locus (números se referindo à "BB") 
  ger1b.2=c(20,20)#Genótipo do progenitor 1 da primeira geração do segundo locus (números se referindo à "bb")
  AB=11#Cria objeto com o valor do gameta AB
  Ab=21#Cria objeto com o valor do gameta Ab
  aB=13#Cria objeto com o valor do gameta aB
  ab=23#Cria objeto com o valor do gameta ab
  ger2a.1= c(ger1a.1[1], ger1b.1[1])#genótipo do progenitor 1 oriundo da F1 (Aa), ou seja, um alelo de cada pai da geração parental
  ger2b.1= c(ger1a.1[2], ger1b.1[2])##genótipo do progenitor 2 oriundo da F1 (Aa)
  ger2a.2=c(AB, Ab, aB, ab)##Cria vetor com os gametas possiveis para o progenitor 1 da f2(2 loci)
  ger2b.2=c(AB, Ab, aB, ab)#Cria vetor com os gametas possiveis para o progenitor 2 da f2(2 loci)
  AA=0# objeto com valor zero para servir de contador de frequencia de 'AA'
  Aa=0# objeto com valor zero para servir de contador de frequencia de 'Aa'
  aa=0# objeto com valor zero para servir de contador de frequencia de 'aa'
  AABB=0#Cria objeto para contar a frequancia de AABB
  AaBb=0#Cria objeto para contar a frequancia de AaBb
  aabb=0#Cria objeto para contar a frequancia de aabb
  AABb=0#Cria objeto para contar a frequancia de AABb
  aaBB=0#Cria objeto para contar a frequancia de aaBB
  AaBB=0#Cria objeto para contar a frequancia de AaBB
  aaBb=0#Cria objeto para contar a frequancia de aaBb
  AAbb=0#Cria objeto para contar a frequancia de AAbb
  Aabb=0#Cria objeto para contar a frequancia de Aabb
  #MODELO TEORICO
  if (modelo=="t")#Se o argumento "modelo" for igual a "t", roda o modelo teorico
  {
  if (loci==1)#Se loci for igual 1, fazer cáculo de monohibridismo
  {
    punnet=matrix (NA, n.gametas,n.gametas, dimnames = list(c("a", "a"), c("A", "A")))#MOntar matriz com NAs para preencher com os valores, com número de linhas e colunas de acordo com número de gametas
    for (j in 1:nrow(punnet))#fluxo "for" com contatdor j para linhas da matriz "punnet"
    {
      for (i in  1:ncol(punnet))#fluxo "for" com contador i para colunas da matriz "punnet"
      {
        punnet[j,i]= ger1a.1[i]+ ger1b.1[j]#Somar valores correspondentes ao "A" do progenitor 1 com "a" do progenitor 2 e guardar em cada posição da matriz
      }
    }
    # #CONVERTER VALORES NÚMERICOS EM LETRAS PARA O USUARIO
    for (j in 1:nrow(punnet)) #Fluxo 'for' com contador j para cada linha de punnet
    {
      for (i in 1:nrow(punnet) )
      {
        if (punnet[j,i]==4)#Se a soma de alguma posição for igual a 4 
        {
          punnet[j,i]="Aa"#substitui por "Aa"
        }
        if (punnet[j,i]==6)#Se a soma de alguma posição for igual a 6
        {
          punnet[j,i]="aa"#substitui por "aa"
        }
        if (punnet[j,i]==2)#Se a soma de alguma posição for igual a 2
        {
          punnet[j,i]="AA"#substitui por "AA"
        }
      }
    }
#SEGUNDA GERAÇÃO
    punnet2=matrix (NA, n.gametas,n.gametas, dimnames = list(c("A", "a"), c("A", "a")))#Montar matriz com NAs para preencher com oos valores, com número de linhas e colunas de acordo com número de gametas
    for (j in 1:nrow(punnet2))#Fluxo "for" com contador j para linhas da matriz "punnet2"
    {
      for (i in 1:ncol(punnet2))#fluxo "for" com contador i para colunas da matriz "punnet2"
      {
        punnet2[j,i]=ger2a.1[i]+ger2b.1[j]#soma um alelo do progenitor 1 a um alelo do progeitor 2 e guarda em uma posição da matriz punnet2
      }
    }
#CONVERTER VALORES NÚMERICOS EM LETRAS PARA O USUARIO
    for (j in 1:nrow(punnet2)) #Fluxo 'for' com contador j para cada linha de punnet2
    {
      for (i in 1:nrow(punnet2) )
      {
        if (punnet2[j,i]==4)#Se a soma de alguma posição for igual a 4 
        {
          Aa = Aa + 1#soma 1 no objeto Aa
          punnet2[j,i]="Aa"#substitui por "Aa"
        }
        if (punnet2[j,i]==6)#Se a soma de alguma posição for igual a 6
        {
          aa= aa+1#soma 1 no objeto aa
          punnet2[j,i]="aa"#substitui por "aa"
        }
        if (punnet2[j,i]==2)#Se a soma de alguma posição for igual a 2
        {
          AA=AA+1#soma 1 no objeto AA
          punnet2[j,i]="AA" #substitui por "AA"
        }
      }
    }
    teorico.a=c(AA, Aa , aa)#Cria vetor com os valores contados nos objetos AA, Aa e aa durante o fluxo acima
    dominante=AA+Aa# Cria objeto para guardar os cotadores correspondentes aos fenotipos dominantes
    recessivo=aa#Cria objeto para guardar os cotadores correspondentes ao fenotipo recessivo
    fenotipos=c(dominante, recessivo)#Cria objeto que guarda ambos os fenótipos
    fenotipos.m=matrix(fenotipos, 1,2, dimnames = list(c("Freq"), c("Dominante(A_)", "Recessivo(aa)") ))#objeto  que cria uma matriz com os dados
    fenotipos.t=as.table(fenotipos.m)# objeto que guarda uma tabela com os fenotipos
    par(mfrow= c(1,2), tcl=0.8, family="serif")# Divide o dispositivo de tela em duas colunas.
    barplot(teorico.a, names.arg=c("AA", "Aa", "aa"), col="coral3", cex.lab=1.3, pch=19, main = "Proporção genotípica")#Barplot com as frequencias genotipicas da  f2
    balloonplot(fenotipos.t, label =T, show.margins = F, dotcolor = "coral3", main=" Proporção fenotípica")#balloonplot com as fraquencias fenotipicas
  }
#PARA DOIS LOCI
  if (loci==2)#Se loci for igual 2, fazer cáculo de diibridismo
  {
    punnet=matrix (NA,n.gametas,n.gametas, dimnames = list(c("AB", "Ab", "aB", "ab"), c("AB", "Ab", "aB", "ab")))#Montar matriz com NAs para preencher com os valores, com número de linhas e colunas de acordo com número de gametas
    for (j in 1:nrow(punnet))
    {
      for (i in  1:ncol(punnet))
      {
        punnet[j,i]=ger1a.1[2]+ger1b.1[1]+ger1a.2[1]+ger1b.2[2]#Soma de um alelo  de cada gene do genotipo do progenitor 1(ger1a.1 e ger1a.2) e do genotipo do progenitor2 (ger1b.1 e ger1b.2)
      }
    }
    #CONVERTER VALORES NÚMERICOS EM LETRAS PARA O USUARIO
    for (j in 1:nrow(punnet)) #Entra no fluxo for com contador j de 1 até numero de linhas da matriz punnet
    {
      for (i in 1:nrow(punnet))
      {
        if (punnet[j,i]==22) #Se a soma de alguma posição for igual a 22
        {
          punnet[j,i]="AABB"#substitui por "AABB"
        }
        if (punnet[j,i]==24)#Se a soma de alguma posição for igual a 24
        {
          punnet[j,i]="AaBB"#substitui por "AaBB"
        }
        if (punnet[j,i]==36)#Se a soma de alguma posição for igual a 36
        {
          punnet[j,i]="aaBb"#substitui por "aaBb"
        }
        if(punnet[j,i]==46)#Se a soma de alguma posição for igual a 46
        {
          punnet[j,i]="aabb"#substitui por "aabb"
        }
        if (punnet[j, i]==32)#Se a soma de alguma posição for igual a 32
        {
          punnet[j,i]="AABb"#substitui por "AABb"
        }
        if (punnet[j,i]==34)#Se a soma de alguma posição for igual a 34
        {
          punnet[j,i]="AaBb"#substitui por "AaBb"
        }
        if (punnet[j,i]==26)#Se a soma de alguma posição for igual a 26
        {
          punnet[j,i]="aaBB"#substitui por "aaBB"
        }
      }
    }
    #SEGUNDA GERAÇÃO
    punnet2=matrix (NA,n.gametas,n.gametas, dimnames = list(c("AB", "Ab", "aB", "ab"), c("AB", "Ab", "aB", "ab")))#Montar matriz com NAs para preencher com os valores, com número de linhas e colunas de acordo com número de gametas
    for (j in 1:nrow(punnet2)) #Entra no fluxo for com contador j de 1 até numero de linhas da matriz punnet2
    {
      for (i in  1:ncol(punnet2)) #Entra no fluxo for com contador j de 1 até numero de colunas da matriz punnet
      {
        punnet2[j,i]=ger2a.2[j]+ger2b.2[i]#Soma um gameta de um progenitor com do outro e guarda em y=uma posição da matriz punnet2
      }
    }
    #CONVERTER VALORES NÚMERICOS EM LETRAS PARA O USUARIO
    for (j in 1:nrow(punnet2)) #Entra no fluxo for com contador j de 1 até numero de linhas da matriz punnet2
    {
      for (i in 1:ncol(punnet2) )#Entra no fluxo for com contador i de 1 até numero de colunas da matriz punnet2
      {
        if (punnet2[j,i]==22)#Se a soma de alguma posição for igual a 22
        {
          AABB = AABB+1
          punnet2[j,i]="AABB"#substitui por "AABB"
        }
        if (punnet2[j,i]==24)#Se a soma de alguma posição for igual a 24
        {
          AaBB= AaBB+1
          punnet2[j,i]="AaBB"#substitui por "AaBB"
        }
        if (punnet2[j,i]==36)#Se a soma de alguma posição for igual a 36
        {
          aaBb=aaBb+1
          punnet2[j,i]="aaBb" #substitui por "aaBb"
        }
        if(punnet2[j,i]==46)#Se a soma de alguma posição for igual a 46
        {
          aabb=aabb+1
          punnet2[j,i]="aabb"#substitui por "aabb"
        }
        if (punnet2[j, i]==32)#Se a soma de alguma posição for igual a 32
        {
          AABb=AABb+1
          punnet2[j,i]="AABb"#substitui por "AABb"
        }
        if (punnet2[j,i]==34)#Se a soma de alguma posição for igual a 34
        {
          AaBb=AaBb+1
          punnet2[j,i]="AaBb"#substitui por "AaBb"
        }
        if (punnet2[j,i]==26)#Se a soma de alguma posição for igual a 26
        {
          aaBB=aaBB+1
          punnet2[j,i]="aaBB"#substitui por "aaBB"
        }
        if (punnet2[j,i]==42)#Se a soma de alguma posição for igual a 42
        {
          AAbb=AAbb+1
          punnet2[j,i]="AAbb"#substitui por "AAbb"
        }
        if (punnet2[j,i]==44)#Se a soma de alguma posição for igual a 42
        {
          Aabb=Aabb+1
          punnet2[j,i]="Aabb"#substitui por "Aabb"
        }
      }
    }
    contador.ab=c(AABB, AAbb, AABb, Aabb,AaBb, AaBB, aaBb, aaBB,  aabb)#Cria vetor que guarda todos os contadores
    A.B=AABB+AABb+AaBb+AaBB#Cria objeto que guarda a soma da frequencia de fenotipos duplamente dominantes
    A.bb=Aabb+AAbb#Cria objeto que guarda a soma da frequencia de fenotipos dominantes para A e recessivos para b
    aa.B=aaBb+aaBB#Cria objeto que guarda a soma da frequencia de fenotipos recessivos para a e dominantes para B
    aa.bb=aabb# Cria objeto que guarda a soma da frequencia dos fenotipos duplamente recessivos
    fenotipos=c(A.B, A.bb, aa.B, aa.bb)#Cria vetor que guarda as proporcoes fenotipicas
    fenotipos.m=matrix(fenotipos, 1,4, dimnames = list(c("freq"), c("A_B_", "A_bb", "aaB_", "aabb") ))#Cria objeto para gurdar fenotipos como matriz
    fenotipos.t=as.table(fenotipos.m)#cria objeto para guardar fenotipos como tabela
    par(mfrow= c(1,2), tcl=0.8, family="serif")#Divide o dispositivo de tela em duas colunas.
    barplot(contador.ab, names.arg=c("AABB", "AAbb", "AABb" , "Aabb","AaBb", "AaBB", "aaBb","aaBB",  "aabb"), col="coral2", cex.lab=1.3, pch=19, main = "Proporção genotípica")#Barplot com as frequencias genotipicas da F2
    balloonplot(fenotipos.t, label =T, show.margins = F, dotcolor = "coral3", main=" Proporção fenotípica")#balloonplot com as fraquencias fenotipicas
  }
    quadros=list(punnet, punnet2) #Cria lista com as matrizes punnet e punnet2
    names(quadros)=c("Geracao F1", "Geracao F2")#Nome para cada uma  das matrizes
    return(quadros)#Retornar lista com as matrizes
  }
  if (modelo=="s")#Se o argumento "modelo" dor igual a "s", roda a simulacao
  {
    if (N==F)#Se o argumento N nao for informado
    {
      stop("Para o argumento modelo='s' deve ser informado N")#interrompe e exibe a mensagem ao usuario
    }
    filho=c(NA, NA)#Cria objeto "filho", com duas posicoes com NA, para guardar os valores sorteados
    if (loci==1)#Se loci for igual a 1
    {
    for (i in 1:N)#Entra em um fluxo "for" com contador i de 1 ate N
    {
      filho[1]=sample(ger2a.1,1)#Sorteia um gameta do genotipo ger2a.1(A ou a)
      filho[2]=sample(ger2b.1,1)#Sorteia um gameta do genotipo ger2b.1(A ou a)
      if (filho[1]== 1 && filho[2]==1)#Se a posicao 1 de filho for igual a 1 e a posicao 2 for igual a 1
      {
        AA=AA+1#Soma 1 no contador AA
      }
      if (filho[1] != filho[2])#Se a posicao 1 de filho for diferente da posicao 2 
      {
        Aa=Aa+1#Soma 1 no contador Aa
      }
      
      if (filho[1]==3 && filho[2]==3)#Se a posicao 1 de filho for igual a 3 e a posicao 2 for igual a 3
      {
        aa=aa+1#Soma 1 no contador aa
      }
    }
    contador.a=c(AA,Aa,aa)#Cria vetor para guardar os tres contadores
    dominante=AA+Aa#Cria objeto para guardar os contadores de genotipos com fenotipo dominante
    recessivo=aa#Cria objeto para guardar contador de genotipos  com fenotipo recessivo
    fenotipos=c(dominante, recessivo)#Cria vetor para guardar de acordo com os fenotipos
    fenotipos.m=matrix(fenotipos, 1,2, dimnames = list(c("freq"), c("dominante(A_)", "recessivo(aa)") ))#Cria objeto para guardar fenotipo como matriz
    fenotipos.t=as.table(fenotipos.m)#Cria objeto para guardar fenotipo como tabela
    par(mfrow= c(1,2), tcl=0.8, family="serif")#Divide o dispositivo de tela em duas colunas.
    barplot(contador.a, names.arg=c("AA", "Aa", "aa"), col="coral3", cex.lab=1.3, pch=19, main = "Proporção genotípica")#Barplot com as frequencias genotipicas da F2 simulada
    balloonplot(fenotipos.t, label =T, show.margins = F, dotcolor = "coral3", main=" Proporção fenotípica")#balloonplot com as fraquencias fenotipicas simuladas
    esperado.d=0.75*N#Cria objeto para guardar a proporcao esperada de fenotipos dominantes para o N dado pelo usuario
    esperado.r=0.25*N#Cria objeto para guardar a proporcao esperada de fenotipos recessivos para o N dado pelo usuario
    chi.valor=sum(((dominante-esperado.d)^2/esperado.d)+((recessivo-esperado.r)^2/esperado.r))#calculando o chi-quadrado
    alfa=0.05#Cria objeto para o valor de alfa a ser consultado
    df=1#cria objeto com os graus de liberdade para essa simulacao (n-1)
    chi.critico=3.841#Objeto para guardar o valor critico de qui quadrado encontrado na tabela para o alfa e os df acima
    saida=list(chi.valor,alfa, df, chi.critico)#Cria uma lista para guardar todos os valores
    valores=as.data.frame(saida, col.names=c("Qhi obs","Alfa", "DF", "Qui crítico"), row.names = "valor")#tranforma a lista em data.frame
    return(valores)#retorna data.frame para o usuario
    }
    # DOIS LOCI
    if (loci==2)#se o argumento "loci" for igual a 2
    {
    for (i in 1:N)#Entra em um fluxo for com contador i de 1 ate N
    {
      filho[1]=sample(ger2a.2,1)#Sorteia um gameta do genotipo ger2a.2(AB, Ab,aB ou ab) e guarda na posicai 1 do objeto "filho"
      filho[2]=sample(ger2b.2,1)#Sorteia um gameta do genotipo ger2b.2(AB, Ab,aB ou ab)e guarda na posicai 2 do objeto "filho"
      filho3=filho[1]+filho[2]#Soma a posicao 1 e 2 e guarda no objeto "filho3"
      if (filho3==22)#Se filho3 for igual a 22
      {
        AABB=AABB+1#soma 1 no contador AABB
      }
      if(filho3==42)#Se filho3 for igual a 42
      {
        AAbb=AAbb+1#soma 1 no contador AAbb
      }
      if (filho3==46)#Se filho3 for igual a 46
      {
        aabb=aabb+1#soma 1 no contador aabb
      }
      if(filho3==26)#Se filho3 for igual a 26
      {
        aaBB=aaBB+1#soma 1 no contador aaBB
      }
      if (filho3==34)#Se filho3 for igual a 34
      {
        AaBb=AaBb+1#soma 1 no contador AaBb
      }
      if (filho3==36)#Se filho3 for igual a 36
      {
        aaBb=aaBb+1#soma 1 no contador aaBb
      }
      if(filho3==32)#Se filho3 for igual a 32
      {
        AABb=AABb+1#soma 1 no contador AABb
      }
      if(filho3==44)#Se filho3 for igual a 44
      {
        Aabb=Aabb+1#soma 1 no contador Aabb
      }
      if(filho3==24)#Se filho3 for igual a 24
      {
        AaBB=AaBB+1#soma 1 no contador AaBB
      }
    }
    simulado=c(AABB, AAbb, AABb, Aabb,AaBb, AaBB, aaBb, aaBB,  aabb)#Cria vetor que guarda todos os contadores
    A.B=AABB+AABb+AaBb+AaBB#Cria objeto que guarda a soma da frequencia de fenotipos duplamente dominantes
    A.bb=Aabb+AAbb#Cria objeto que guarda a soma da frequencia de fenotipos dominantes para A e recessivos para b
    aa.B=aaBb+aaBB#Cria objeto que guarda a soma da frequencia de fenotipos recessivos para a e dominantes para B
    aa.bb=aabb# Cria objeto que guarda a soma da frequencia dos fenotipos duplamente recessivos
    fenotipos=c(A.B, A.bb, aa.B, aa.bb)#Cria vetor para guardar todas as proporcoes fenotipicas
    fenotipos.m=matrix(fenotipos, 1,4, dimnames = list(c("freq"), c("A_B_", "A_bb", "aaB_", "aabb") ))#Cria objeto para guardar fenotipos em uma matriz
    fenotipos.t=as.table(fenotipos.m)##Cria objeto para guardar fenotipos em uma tabela
    par(mfrow= c(1,2), tcl=0.8, family="serif")#Divide o dispositivo de tela em duas colunas.
    barplot(simulado, names.arg=c("AABB", "AAbb", "AABb" , "Aabb","AaBb", "AaBB", "aaBb","aaBB",  "aabb"), col="coral2", cex.lab=1.3, pch=19, main = "Proporção genotípica simulada")#Barplot com frequencias genotipicas simuladas
    balloonplot(fenotipos.t, label =T, show.margins = F, dotcolor = "coral3", main=" Proporção fenotípica simulada")#Balloonplot com frequencias fenotipicas simuladas
    esperado.A.bb=0.1825*N#Cria objeto para guardar a proporcao esperada de fenotipo A_bb para o N dado pelo usuario
    esperado.A.B=0.5625*N#Cria objeto para guardar a proporcao esperada de fenotipo A_B_ para o N dado pelo usuario
    esperado.aa.B=0.1825*N#Cria objeto para guardar a proporcao esperada de fenotipo aaB_ para o N dado pelo usuario
    esperado.aa.bb=0.0625*N#Cria objeto para guardar a proporcao esperada de fenotipo aabb para o N dado pelo usuario
    chi.valor=sum((((A.B-esperado.A.B)^2)/esperado.A.B)+(((A.bb-esperado.A.bb)^2)/esperado.A.bb)+(((aa.bb-esperado.aa.bb)^2)/esperado.aa.bb)+(((aa.B-esperado.aa.B)^2)/esperado.aa.B))#calculando o chi-quadrado
    alfa=0.05#Cria objeto para o valor de alfa a ser consultado
    df=3#cria objeto com os graus de liberdade para essa simulacao (n-1)
    chi.critico=7.815#Objeto para guardar o valor critico de qui quadrado encontrado na tabela para o alfa e os df acima
    saida=list(chi.valor,alfa, df, chi.critico)#Cria uma lista para guardar todos os valores
    valores=as.data.frame(saida, col.names=c("Qhi obs","Alfa", "DF", "Qui crítico"), row.names = "valor")#tranforma a lista em data.frame
    return(valores)#Retorna o data.frame "valores"
    }
  }
}