####################################################################
##Execute following on dos command prompt (omit the "##"):
##/usr/bin/R -q --no-restore --no-save < Spectra.r >& R.out
##"S:\\R\\rw2001\\bin\\Rterm.exe" -q --no-restore --no-save < Spectra_Test.r > R.out
####################################################################

options(echo = FALSE)

#--------------------------------------------------------------------
# All plots will be named plot001.png, plot002.png... from each
# successive call to a plotting function below.

png(file="plot%03d.png", width = 680, height = 450, pointsize = 12, bg="transparent", res = 200)
#bitmap(file="plot%03d.png", width=5, height=5, res=400)

#--------------------------------------------------------------------
# Simple harmonic series: sin[(pi/2)*t], sin[2*pi*t]; 

x1 <- sin((pi/2)*seq(1:136)/8);
x2 <- sin(2*pi*seq(1:136)/8);
x3 <- rnorm(136); # gaussian noise.

startyear <- 1988;
startperiod <- 1;
ts1 <- ts( x1, start=c(startyear, startperiod), frequency=4 );
ts2 <- ts( x2, start=c(startyear, startperiod), frequency=4 );

# plot timeseries:
ts.plot(ts1, ts2, gpars=list(ylab="Value", col=c(1:2)),
        main="sin[(pi/2)*t] (BLACK); sin[2*pi*t] (RED)");

#--------------------------------------------------------------------
# Spectra.

# frequency scale is fraction of a cycle per Delta_t sampling interval
# in input series. Max=0.5 cycle per Delta_t => 1 cycle=2*Delta_t=2 Qtrs:
# Period (time represented by 1 complete cycle) = "1/f" qtrs. (= 2pi/omega where
# omega = 2*pi*f)

spectrum_xd <- spectrum(cbind(x1,x2, x3), method="pgram", plot=TRUE);
spectrum_xd$freq
spectrum_xd$spec

dev.off()
q()