{

par(mfrow=c(3,1))

source("conflines.txt")
lf = matrix(scan("LeesFerry-monflows-1906-2006.txt"), ncol=13, byrow=T)
lf = lf[,2:13]
lf = lf*0.0004690502	# convert to cms


y=lf
y=array(t(y))
N=length(y)
cc=acf(y,type="covariance",lag.max=N-1,plot=FALSE)$acf
NN=length(cc)
NN1=NN-1
crev=rev(cc)
cfull=c(crev[2:NN],cc)
Nfull=length(cfull)

## number of frequencies..
freqs=2*pi*(1:round((N/2)))/N
plotfreqs=freqs/(2*pi*(1/12))
tordf=1:round((N/2))
tord=1:N
nfreq=length(freqs)


#pgrams1=Mod(fft((cfull/Nfull)))^2
#pgrams1=pgrams1/(pi)
#plot(plotfreqs,pgrams1[2:(nfreq+1)],xlab="freq. cy/yr",ylab="Pgram", type="h")

### Spectrum from sample ACF - 12.2.1b Wei
NN=length(cc)
for(i in 1:nfreq){
pgrams[i]=(cc[1]+(2*sum(cc[2:NN]*cos(freqs[i]*(1:(NN-1))))))/(2*pi)
}

plot(plotfreqs,pgrams,xlab="freq. cy/yr",ylab="Pgram", type="h")

## FFT of the ACF
cc=c(rev(cc[2:NN]),cc)
pgrams=1:nfreq
for(i in 1:nfreq){
  pgrams[i]=(Re(sum(cc*exp(-1i*freqs[i]*((-N+1):(N-1))))))/(2*pi)
}

plot(plotfreqs,pgrams,xlab="freq. cy/yr",ylab="Pgram", type="h")

}