\documentclass[pnumromarab,normaltoc]{abnt} 

\usepackage[brazil]{babel}
\usepackage[utf8]{inputenc}
\usepackage{abnt-alf}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{alltt}
\usepackage{caption}
\usepackage{textcomp}
\usepackage[T1]{fontenc} 
\usepackage{aurical}    
\usepackage[dvips]{graphicx}
\usepackage{epic,curves}
\usepackage{exscale,relsize}
\usepackage{enumerate}

\makeatletter 
\begin{document}
\noindent \textbf{SCRIPTS PARA O ESTUDO DE SIMULAÇÃO:}\\

\begin{verbatim}
## Função para simular um conjunto de dados
simu.dados <- function(theta,n.simul){
  mu1 <- theta[1]
  mu2 <- theta[2]
  s1 <- theta[3]
  s2 <- theta[4]
  tau1 <- theta[5]
  tau2 <- theta[6]
  phi <- theta[7]
  rho <- theta[8]
  ## Gerando erros de uma normal bivariada com mu=0, 
  ## variâncias tau1, tau2 e correlação rho 
  v <- matrix(c(tau1^2,tau1*tau2*rho,tau1*tau2*rho,tau2^2),2,2)
  Sc <- kronecker(diag(n.simul/2),v)
  set.seed(3333)
  e <- mvrnorm(n.simul,mu=c(0,0),v)
  ## Gerando um campo aleatório gaussiano com média 0, variância 1 e 
  ## parâmetros de alcance phi 
  set.seed(3333)
  U <- grf(n.simul,grid='irreg',cov.model='exp',cov.pars=c(1,phi),
           mean=0,messages=FALSE)
  S1 <- s1*U$data
  S2 <- s2*U$data
  ## Calculando e organizando os dados
  y1 <-as.geodata(data.frame(U$coords[,1],U$coords[,2],mu1+S1+e[,1]))
  y2 <- as.geodata(data.frame(U$coords[,1],U$coords[,2],mu2+S2+e[,2]))
  seq1 <- seq(1,length(y1$data)*2,by=2)
  seq2 <- seq(2,length(y2$data)*2,by=2)
  y <- c()
  y[seq1] <- y1$data
  y[seq2] <- y2$data
  coords <- data.frame(U$coords[,1],U$coords[,2])
  ## Preparando a saída
  retorna <- list()
  retorna[1][[1]] <- y1$data
  retorna[2][[1]] <- y2$data
  retorna[3][[1]] <- y
  retorna[4][[1]] <- coords
  return(retorna)
}
\end{verbatim}

\noindent \underline{\textbf{Script referente a 1a configuração de diagrama ternário}}  

\begin{verbatim}
## Script rodado para a 1a configuração de diagrama ternário 
## Gera valores preditos por quadratura e simulação pra fazer mapas 
require(MASS)
require(geoR)
require(geoComp)
require(compositions)
require(statmod)
source("simula.R")
## Função que gera diagrama ternário
gerater <- function(theta,n.simul=100){
  dados <- simu.dados(theta,n.simul=n.simul)
  y1 <- dados[1][[1]]
  y2 <- dados[2][[1]]
  y12 <- data.frame(y1,y2)
  comp <- agl(y12)
  comp1 <- acomp(comp)
  names(comp1) <- c("Areia", "Silte", "Argila")
  plot(comp1)
  plot(mean(comp1),add=T,pch=20,col='red')
  ellipses(mean(comp1),var(comp1),col='red',r=2)
  ellipses(mean(comp1),var(comp1),col='red',r=4)
  retorna <- data.frame(comp,dados[[4]][,1],dados[[4]][,2])
  return(retorna)
}
###### 1a configuração ################
theta <- c(-0.2,-0.5,1,1.5,0.3,0.3,0.6,0.9)
###### 2a configuração : theta <- c(1,1,1.2,1.5,0.9,1,0.6,0.5)
###### 3a configuração : theta <- c(-0.2,-1,0.45,0.13,0.3,0.3,0.6,0.95)
#set.seed(1333)
dados1 <- gerater(theta,n.simul=100)
names(dados1) <- c("Areia","Silte","Argila","Coord.X", "Coord.Y")
write.table(dados1,"dad_conf1_by0.03k7ns1000.txt")
dados1 <- as.geoComp(dados1)
write.table(dados1[[3]],"points1_by0.03k7ns1000.txt")
estima <- mec(dados1)
write.table(estima[[1]],"estima1par_by0.03k7ns1000.txt")
write.table(estima[[2]],"estima1log_by0.03k7ns1000.txt")
borda <- cbind(c(0,0,1,1,0),c(0,1,1,0,0))
gr <- pred_grid(borda, by=0.03)
source("cokrigagem.R")
md.cov.ck <- cokrigagem(estima[[1]]$Estimativas,loc=gr,
                        dados.comp=dados1)
preditos.gh <- volta.quad(md.cov.ck,n.pontos=7,Variancia=FALSE)
write.table(preditos.gh,"dt1predgh_by0.03k7ns1000.txt")
preditos.simu <- volta.cokri(md.cov.ck,num.simu=1000,int.conf=0.95)
write.table(preditos.simu,"dt1predsimu_by0.03k7ns1000.txt")
\end{verbatim}
\end{document}

