#########first row from wave/power is smallest period

#units="temp(C)"
#units="precip(mm)"
units="acre-ft"
#units="index"


source("fcontour.R")

coi=fin$coi



divf=mean(fin$power)/5

DIVF=signif(divf, digits=2)

pow=fin$power/(DIVF)



my.levels <- quantile(pow, probs = seq(0, 1, 0.2))



#my.levels=c(0, .068, .19, .4, 1.4, 2)
    my.cols <- c("white", "blue", "green","yellow", "red")


p=t(pow)
	for( i in 1:length(p[,1])){
	p[i,]=rev(p[i,])
	}	

xper=subset(fin$period, fin$period<(.5*length(data[,1])))

y=1:length(xper)

del=length(fin$period)-length(xper)

P=p[,(1+del):(length(fin$period))]


###this takes every 3rd period to simplify y axis labeling

ylbl=1:floor(length(fin$period)/3)

pos=1

for(i in 1:length(ylbl)){

ylbl[i]=fin$period[pos]

pos=pos+3

}

source("local_noise.txt") #significance function

prob=.95 # local confidence level

xt=local_noise(prob ,XL, "w",  fin, length(data[,2]))



test=c(2,4,8,16,32,64,128,256,512)


l=length(fin$power[1,])

val=1:length(fin$period)

 for( i in 1:length(fin$period)){

 val[i]=mean(pow[i,])


 }

ts=1:l


VAL=val


#VAL[1]=VAL[2]


device = png
ext = "png"
device(paste("power_spectrum",ext, sep="."), width=1200, height=600)
par(mfrow=c(1,3), cex=1.3)
layout(matrix(c(1,2,3),ncol=3, nrow=1),widths = c(.575,.1,.325))
par(mar=c(4,5,2,1))

fcontour(data[,1], log(fin$period), t(pow), levels= my.levels,col=my.cols, ylim=c(log(2), max(log(fin$period))), nlevels=5,ylab="", xlab="",

plot.axes = {
           
            axis(1,cex.axis=1.5)
           # axis(2 , at = log(ylbl), labels = round(ylbl, digits=1), cex.axis=1.5)
	     axis(2 , at = log(test), labels = test, cex.axis=1.5)
	      	
	}
)


contour(data[,1], log(fin$period), t(pow), nlevels=1,levels=(xt$sig[1]/DIVF), add=T, drawlabels=F)
lines(data[,1], log(fin$coi), lty=2)
mtext("Period (Years)", side = 2, line = 3.5, cex=1.1)

mtext("Year", side = 1, line = 2.75, cex=1.1)


cbar(nlevels=5,col=my.cols, levels=c(0,1,2,3,4,5), key.title = "",key.axes = {axis(2,at =c(0,1,2,3,4,5) , labels = round(my.levels, digits=1), cex.axis=1.5)


})

mtext(parse(text=paste("(",units,"^2)")), side = 2, adj=.65,line = 4) #.53
mtext("Power/", side = 2, adj=.23,line = 4) #.27
mtext(paste(DIVF), side = 2, adj=.4,line = 4)

par(usr=c(min(VAL), max(VAL), min(log(fin$period)), max(log(fin$period))))
plot(VAL,log(fin$period), ylab="", xlab="", type="l", ylim=range(log(fin$scale)), xlim=range(VAL,xx$sig/DIVF),lwd=2,main=, axes=F, yaxs="i")


lines(xx$sig/DIVF, log(fin$period), col="grey")
lines(zz$sig/DIVF, log(fin$period), col="grey", lty=5)



axis(1, cex.axis=1.5)
#axis(2, at = log(ylbl), labels=round(ylbl, digits=1),las=2, cex.axis=1.5)
axis(2 , at = log(test), labels = test,las=2, cex.axis=1.5)

box()

#mtext(paste("sig=",lvl), side = 3, line = .7, cex=1.1)
mtext("Period (Years)", side = 2, line = 4, cex=1.1)

mtext("Variance/", side = 1, adj=.01,line = 2.75) #.1
mtext(parse(text=paste("(",units,"^2)")), side = 1, adj=.99,line =2.75) #.8
mtext(paste(DIVF), side = 1, adj=.46,line = 2.75) #.47

dev.off()



device = postscript
ext = "eps"
device(paste("power_spectrum",ext, sep="."),width=6, height=3,pointsize=8, onefile=F, horizontal=F)
par(mfrow=c(1,3), cex=1.3)
layout(matrix(c(1,2,3),ncol=3, nrow=1),widths = c(.575,.1,.325))
par(mar=c(4,5,2,1))

fcontour(data[,1], log(fin$period), t(pow), levels=my.levels,col=my.cols,ylim=c(log(2), max(log(fin$period))), nlevels=5,ylab="", xlab="",

plot.axes = {
           
            axis(1,cex.axis=1.5)
           # axis(2 , at = log(ylbl), labels = round(ylbl, digits=1), cex.axis=1.5)
	     axis(2 , at = log(test), labels = test, cex.axis=1.5)
	      	
	}
)

mtext("Period (Years)", side = 2, line = 3.5, cex=1.1)

mtext("Year", side = 1, line = 2.75, cex=1.1)

contour(data[,1], log(fin$period), t(pow), nlevels=1,levels=(xt$sig[1]/DIVF), add=T, drawlabels=F)

lines(data[,1], log(fin$coi), lty=2)


cbar(nlevels=5,col=my.cols, levels=c(0,1,2,3,4,5), key.title = "",key.axes = {axis(2,at =c(0,1,2,3,4,5) , labels = round(my.levels, digits=1), cex.axis=1.5)


})


mtext(parse(text=paste("(",units,"^2)")), side = 2, adj=.65,line = 4) #.53
mtext("Power/", side = 2, adj=.23,line = 4)  #.27
mtext(paste(DIVF), side = 2, adj=.4,line = 4)

par(usr=c(min(VAL), max(VAL), min(log(fin$period)), max(log(fin$period))))
plot(VAL,log(fin$period), ylab="", xlab="", type="l", ylim=range(log(fin$scale)), xlim=range(VAL,xx$sig/DIVF),lwd=2,main=, axes=F, yaxs="i")

lines(xx$sig/DIVF, log(fin$period), col="grey")
lines(zz$sig/DIVF, log(fin$period), col="grey", lty=5)


axis(1, cex.axis=1.5)
#axis(2, at = log(ylbl), labels=round(ylbl, digits=1),las=2, cex.axis=1.5)
axis(2 , at = log(test), labels = test,las=2, cex.axis=1.5)

box()

#mtext(paste("sig=",lvl), side = 3, line = .75, cex=1.1)
mtext("Period (Years)", side = 2, line = 4, cex=1.1)

mtext("Variance/", side = 1, adj=.01,line = 2.75) #.1
mtext(parse(text=paste("(",units,"^2)")), side = 1, adj=.99,line =2.75)  #.8
mtext(paste(DIVF), side = 1, adj=.46,line = 2.75)  
dev.off()