# CALCULO da DIMENSÃO DE CORRELAÇÃO PARA UMA DADA DIMENSÃO DE IMERSÃO

# OBS: Os dados devem estar dispostos ao longo de uma coluna com header X


setwd("C:/Users/Sérgio/Documents/FisicaAmbiental/Dimensionalidade")



data <- read.csv('T_Sinop_Met_Jun2002.txt',header = TRUE, sep="\t")

ndat <- nrow(data) - 1

print(ndat)



#TEMPO DE DEFASAGEM E DIMENSÃO DE IMERSÃO:

TAO <- 10
Dim <- 10

LOGR_Min <- -2.0
LOGR_Max <- 2.0

deltalr <- (LOGR_Max - LOGR_Min)/400    # intervalo de log r dividido em 200 intervalos

RES <- matrix(0.0,nrow=401,ncol=3)


# CALCULA DA MATRIZ DE DISTÂNCIAS para N dimensões do Espaço de Imersão

DIST2N <- matrix(0.0,nrow=ndat,ncol=ndat)
RData <- matrix(0.0,nrow=401,ncol=1)

imin <- TAO*Dim
jmin <- TAO*Dim
imax <- ndat - TAO*Dim
jmax <- ndat - TAO*Dim


for(i in imin:imax){
print(i)
for(j in jmin:jmax){
	

		for(m in 1:Dim){

			DIST2N[i,j] <- DIST2N[i,j] + (data$X[i+(m-1)*TAO] - data$X[j+(m-1)*TAO])^2

			} #for m

	}
	}

#CALCULO DO COEFICIENTE DE CORRELAÇÃO DE HEAVISIDE:


for(i in imin:imax){
print(i)
for(j in jmin:jmax){
	
	if(DIST2N[i,j] != 0.0){	
			LOG_DIST <- log10(sqrt(DIST2N[i,j]))
			indice <- as.integer((LOG_DIST - LOGR_Min)/deltalr)
			if(indice > 0 && indice < 401){RData[indice,1] <- RData[indice,1] + 1}

		}

	} # for j
	} # for i

LOGR <- LOGR_Min

for(k in 1:401){

	RES[k,1] <- LOGR
	
	LOGR <- LOGR + deltalr

	for(m in 1:k){

		RES[k,2] <- RES[k,2] + RData[m,1]
		} # for m

	RES[k,2] <- log10((RES[k,2] + ndat)/(ndat^2))   # Acrescenta-se ndat ao RES(ultado) para levar em conta os próprios pontos no cômputo da dim corr


	} # for k

	

	plot(RES[,1],RES[,2],xlab="log r",ylab="log C", type="b")


# CALCULO DA RAZÃO logC / logr

ratio <- rep(0.0,times = 401)

ratioMAX <- 0.0

for(k in 2:400){

	ratio[k] <- ((RES[k,2] - RES[k-1,2])/(RES[k,1] - RES[k-1,1]) + (RES[k+1,2] - RES[k,2])/(RES[k+1,1] - RES[k,1]))/2
	RES[k,3] <- ratio[k]
	if(ratio[k] > ratioMAX) ratioMAX <- ratio[k]
	}

plot(RES[,1],ratio,xlab="log r",ylab="logC/logr", type="b")

write.table(RES,file = "Resultado", sep=" ", eol="\n")

print(ratioMAX)