hohenPr<- read.table("C:/CLIM/HohenP.txt",header=T)        ###    
attach(hohenPr)
postscript(file="C:/CLIM/Hpgram.ps",height=6,width=15,horiz=F)  
quot<-"Hohenpeissenberg, Precip. Winter (Resid.) 1879-2010";quot

plotP<- function(k,z,nh,tylab,yli,bc,bct,lte,xte){
plot(k,z,type="l",lty=1,xlim=c(0,nh),ylim=yli,ylab=tylab,xlab="")
segments(1,bc,nh-4,bc,lty=2)        #Drawing simultaneous bounds
text(nh,bc,bct,cex=0.7)  
#T-values as text on the bottom margin (side=1,line=2)
mtext(lte,side=1,line=2,at=xte[1:7])
mtext(c("k","T"),side=1,line=c(1,2),at=xte[8])       
}         

#----------Preparation of the vector Y, to be analyzed ---------
Y1<- (dcly+jan+feb)/1000    #Amount of precipitation winter [dm]
Ja<- Year-1900; Ja2<- Ja*Ja; Ja3<- Ja2*Ja; Ja4<- Ja2*Ja2
#Removal of polynomial trend
Y<- Y1 - predict(lm(Y1~Ja+Ja2+Ja3+Ja4))
n<- length(Y); nh<- round(n/2); SD2<- var(Y)   
 
#Calculate periodogram via Fourier-coefficients a,b
seq<-1:nh; pg<-seq     #vectors of dim nh
for (k in 1:nh) 
{a<- 0; b<- 0; omk<- 2*pi*k/n
for (i in 1:n) 
{a<- a+ Y[i]*cos(omk*i)
 b<- b+ Y[i]*sin(omk*i)
 pg[k]<- (a*a+b*b)/(n*pi)}
}

#-------------Plotting the periodogram--------------------------
Pgr<- pg/SD2                                      #Standardizing  
tylab<- "Periodogram R(k)^2*n/(4 pi s^2)"       #R^2 = a^2 + b^2
#Simultaneous bounds b_l=-ln(alpha/l)/pi, l=1,4,12
bc<- -log(0.05/c(1,4,12))/pi              #3 bounds b_1,b_4,b_12
bct<- c("b_1","b_4","b_12")   
lte<- c("26.4","13.2","6.6","4.4","3.3","2.6","2.2")
xte<- c(5,10,20,30,40,50,60,70)  
yli<- c(0.0,1.8)  
      
plotP(seq,Pgr,nh,tylab,yli,bc,bct,lte,xte)              
title(main=quot)
 
dev.off()  

detach(hohenPr)                                             ###
rm(list=objects())                                          ###