# ---------------------------------------------------------------------
# Program: MInfoPlot1.S  
#  Author: Steven M. Boker
#    Date: Thu Apr 22 09:38:41 EST 1999
#
# This program creates mutual information plots for the Lorenz System
# 
#
# ---------------------------------------------------------------------
# Revision History
#   Steven M. Boker -- Thu Apr 22 09:38:37 EST 1999
#      Create MInfoPlot1.S
#
# ---------------------------------------------------------------------


# -----------------------------
# Read in the The Lorenz attractor

lorenz <- scan("lor63.dat", list(x=0), skip=0)
maxSamples <- length(lorenz$x)
maxSamples
summary(data.frame(x=lorenz$x))

# -----------------------------
# Plot the first 1000 samples of the Lorenz time series

thePlotData <- lorenz$x[1:1000]
tepsfile <- paste("LorenzTimeSeries.eps",sep="")
tmain <- paste("Lorenz Time Series", sep="")
postscript(tepsfile,height=6.4,horizontal=F)
plot(c(0, 1000), c(-30, 30), xlab="Samples", ylab="Displacement", 
    main=tmain, type="n")
lines(seq(1, 1000),thePlotData, type="l")
dev.off()

# -----------------------------
# Plot the full Frequency Power Spectrum for Lorenz

tPgram <- pgram(lorenz$x)
maxSamples <- length(lorenz$x)

thePlotData <- (tPgram / sum(tPgram)) * 100
theFreq <- rep(0,maxSamples/2)
for(j in 1:(maxSamples/2)) {
	theFreq[j] <- (j/maxSamples) * (2*3.14159)
}
tepsfile <- paste("LorenzPgram.eps",sep="")
tmain <- paste("Frequency Power Spectrum of Lorenz", sep="")
postscript(tepsfile,height=6.4,horizontal=F)
plot(c(0, theFreq[maxSamples/2]), c(0,1), xlab="Frequency", 
      ylab="Percent of Total Power", main=tmain, type="n")
lines(theFreq,thePlotData[1:(maxSamples/2)], type="l")
dev.off()

# -----------------------------
# Calculate the Mutual Information (Nonlinear Dependency) for the Lorenz data

maxTau <- 100
write(t(format(round(lorenz$x,6))), file="temp.dat", ncol=1, append=F)
theNLDCommand <- paste("nldcalc -o temp.nld -taumax ", maxTau, " temp.dat", sep="")
unix(theNLDCommand) 
nldLorenz <- scan("temp.nld", list(nld=0), skip=0)
summary(nldLorenz$nld)

# -----------------------------
# Plot the Mutual Information (Nonlinear Dependency) for the Lorenz data

thePlotData <- nldLorenz$nld
tepsfile <- paste("NLDLorenz1.eps",sep="")
tmain <- paste("Mutual Information for Lorenz System", sep="")
postscript(tepsfile,height=6.4,horizontal=F)
plot(c(-maxTau, maxTau), c(0, 0.4), xlab="Tau", ylab="Mutual Information", 
    main=tmain, type="n")
lines(seq(-maxTau, maxTau), thePlotData, type="l")
dev.off()

# -----------------------------
# Plot the Mutual Information (Nonlinear Dependency) for the Lorenz data

thePlotData <- nldLorenz$nld
tepsfile <- paste("NLDLorenz2.eps",sep="")
tmain <- paste("Mutual Information for Lorenz System", sep="")
postscript(tepsfile,height=6.4,horizontal=F)
plot(c(-maxTau, maxTau), c(0, 0.05), xlab="Tau", ylab="Mutual Information", 
    main=tmain, type="n")
lines(seq(-maxTau, maxTau), thePlotData, type="l")
dev.off()


# ---------------------------------------------------------------------
# Quit the program here  

q()