###
###    source("C:/File_Pathname/17_Periodogram_Script.txt", echo=FALSE)
###

#      Clear everything from memory.
rm(list=ls())
plot.new()

#      Set the working directory
setwd("C:/File_Pathname/")


####################################################################################################
###                                                                                              ###
###    Initialization Section.                                                                   ###
###    All possible variables used in this program are initialized in this section.              ###
###                                                                                              ###
####################################################################################################

#         Guidelines for setting interval_fct, depends on the number of samples in the time-series.
#         Make the integer "interval_fct" larger if more frequency-intervals are needed.
#                   If the number of samples in the time-series < 100, set interval_fct = 6
#                   100 < sample-size < 300, set interval_fct = 4
#                   300 < sample-size < 500, set interval_fct = 3
#                   500 < sample-size <1000, set interval_fct = 2
#                         sample-size >1000, set interval_fct = 1
interval_fct= 2
LS_s_title = "The Unified Cycle Theory,  http://www.uct-news.com                  Data Source:  Moberg et al., 2005          "
LS_size  = 1.00
LS_font  = 2
LS_color = "black"
LS_s_size  = .75
LS_s_font  = 2
LS_s_color = "black"

#       Determine if the spectrum-data is to be written to the hard-drive.
Write_Output_file = FALSE

draw_EWUS_primary = TRUE
draw_EWUS_secondary = FALSE
draw_solar_lines= FALSE

#      If numb_of_plots =  1 - Only the Smoothed Periodogram prints.
#                       =  2 - The Time-Series Plot is eliminated.
#                       =  3 - All 3 plots get displayed.
numb_of_plots = 2
ts_y_label  = "Delta Oxy-18"
ts_descr    = "2000-Yr Global Climate"

#******************************  Tests Begin Here (6 tests)  ****************************


test_no = 1



#    Test # 1: Global Climate
if(test_no == 1) { 
   file_name = "Data/Climate_a_19-yr.txt"
   period_EUWS = 19.0953697945932   #  The wavelength being investigated.
   period_solar= 700.441 
   search_adj_at_peaks = .012
   p_min  = 5                       #  Analyze wavelengths beginning with this minimum period.
   p_max  = 50                      #  Analyze wavelengths ending with this maximum period.
   start_time  = 0                  #  The time of the first observation.
   interval    = 2.1217077549548    #  Delta-t -- Interval between observations (in units = to timescale).
   Daniell_smooth = c(7,9)
   Daniell_null   = c(75,125)
   filter_meth = "Raw Data"         #  Describe the filtering method, or "Raw Data" if filtering was not used.
   #        Placement of text on Lomb-Scargle Periodogram.
   x_LS_wave = .055
   y_LS_wave = 8.25
   x_LS_confid = x_LS_wave
   y_LS_confid = 0.92 * y_LS_wave
   LS_text_size = .8
   x_LS_alt = .105
   y_LS_alt = 6.5
   LS_alt_message = "8.8814-yr"
   prd_multiplier = 1             #  Multiplier for reporting maximum.
   prd_label = "yr"               #  Unit of measure for reporting maximum. 
   #        Placement of text on Smoothed Periodogram.
   ylim_min_smooth = 0
   ylim_max_smooth = .095
   x_smooth_wave = 1.07 * x_LS_wave
   y_smooth_wave = 0.982 * ylim_max_smooth
   x_smooth_confid = x_smooth_wave
   y_smooth_confid = 0.905 * ylim_max_smooth
   confid_indc = "Confidence Level exceeds 95%"
   smooth_text_size = .8
   x_smooth_alt = 0.101
   y_smooth_alt = .056
   smooth_alt_message = "9.1188-yr"
   #     x-axis and y-axis labels for periodograms.
   LS_x_title = "Frequency (Cycles per Year)"
   LS_y_title = "Spectral Density"
   ts_x_label  = "Time (Years)"     #  x-axis label for time-series plot.
   time_scale  = "Yr"               #  Unit-of-Time for the analysis
   time_plural = "Yr"
}

#    Test # 2: Global Climate
if(test_no == 2) { 
   file_name = "Data/Climate_b_19-yr.txt"
   period_EUWS = 19.0953697945932   #  The wavelength being investigated.
   period_solar= 700.441 
   search_adj_at_peaks = .012
   p_min  = 5                       #  Analyze wavelengths beginning with this minimum period.
   p_max  = 50                      #  Analyze wavelengths ending with this maximum period.
   start_time  = 0                  #  The time of the first observation.
   interval    = 2.1217077549548    #  Delta-t -- Interval between observations (in units = to timescale).
   Daniell_smooth = c(7,9)
   Daniell_null   = c(75,125)
   filter_meth = "Raw Data"         #  Describe the filtering method, or "Raw Data" if filtering was not used.
   #        Placement of text on Lomb-Scargle Periodogram.
   x_LS_wave = .058
   y_LS_wave = 13.85
   x_LS_confid = x_LS_wave
   y_LS_confid = 0.92 * y_LS_wave
   LS_text_size = .8
   x_LS_alt = .101
   y_LS_alt = 7
   LS_alt_message = "9.2721-yr"
   prd_multiplier = 1             #  Multiplier for reporting maximum.
   prd_label = "yr"               #  Unit of measure for reporting maximum. 
   #        Placement of text on Smoothed Periodogram.
   ylim_min_smooth = 0
   ylim_max_smooth = .125
   x_smooth_wave = 1.00 * x_LS_wave
   y_smooth_wave = 0.982 * ylim_max_smooth
   x_smooth_confid = x_smooth_wave
   y_smooth_confid = 0.905 * ylim_max_smooth
   confid_indc = "Confidence Level exceeds 99.9%"
   smooth_text_size = .8
   x_smooth_alt = 0.097
   y_smooth_alt = .066
   smooth_alt_message = "9.7927-yr"
   #     x-axis and y-axis labels for periodograms.
   LS_x_title = "Frequency (Cycles per Year)"
   LS_y_title = "Spectral Density"
   ts_x_label  = "Time (Years)"     #  x-axis label for time-series plot.
   time_scale  = "Yr"               #  Unit-of-Time for the analysis
   time_plural = "Yr"
}

#    Test # 3: Global Climate
if(test_no == 3) { 
   file_name = "Data/Climate_c_19-yr.txt"
   period_EUWS = 19.0953697945932   #  The wavelength being investigated.
   period_solar= 700.441 
   search_adj_at_peaks = .012
   p_min  = 5                       #  Analyze wavelengths beginning with this minimum period.
   p_max  = 50                      #  Analyze wavelengths ending with this maximum period.
   start_time  = 0                  #  The time of the first observation.
   interval    = 2.1217077549548    #  Delta-t -- Interval between observations (in units = to timescale).
   Daniell_smooth = c(7,9)
   Daniell_null   = c(75,125)
   filter_meth = "Raw Data"         #  Describe the filtering method, or "Raw Data" if filtering was not used.
   #        Placement of text on Lomb-Scargle Periodogram.
   x_LS_wave = .066
   y_LS_wave = 24.7
   x_LS_confid = x_LS_wave
   y_LS_confid = 0.92 * y_LS_wave
   LS_text_size = .8
   x_LS_alt = .029
   y_LS_alt = 7.5
   LS_alt_message = "23.7852-yr                       10.9005-yr                                       6.0544-yr"
   prd_multiplier = 1             #  Multiplier for reporting maximum.
   prd_label = "yr"               #  Unit of measure for reporting maximum. 
   #        Placement of text on Smoothed Periodogram.
   ylim_min_smooth = 0
   ylim_max_smooth = .16
   x_smooth_wave = 1.03 * x_LS_wave
   y_smooth_wave = 0.982 * ylim_max_smooth
   x_smooth_confid = x_smooth_wave
   y_smooth_confid = 0.905 * ylim_max_smooth
   confid_indc = "Confidence Level exceeds 99.9%"
   smooth_text_size = .8
   x_smooth_alt = 0.019
   y_smooth_alt = .07
   smooth_alt_message = "22.6954-yr                                         10.8719-yr                            6.0522-yr"
   #     x-axis and y-axis labels for periodograms.
   LS_x_title = "Frequency (Cycles per Year)"
   LS_y_title = "Spectral Density"
   ts_x_label  = "Time (Years)"     #  x-axis label for time-series plot.
   time_scale  = "Yr"               #  Unit-of-Time for the analysis
   time_plural = "Yr"
}

#    Test # 4: Global Climate
if(test_no == 4) { 
   file_name = "Data/Climate_d_57-yr.txt"
   period_EUWS = 57.2861093837796   #  The wavelength being investigated.
   period_solar= 700.441 
   search_adj_at_peaks = .012
   p_min  = 20                  #  Analyze wavelengths beginning with this minimum period.
   p_max  = 200                 #  Analyze wavelengths ending with this maximum period.
   start_time  = 0                 #  The time of the first observation.
   interval    = 6.3651232648644    #  Delta-t -- Interval between observations (in units = to timescale).
   Daniell_smooth = c(5,7)
   Daniell_null   = c(65,105)
   filter_meth = "Raw Data"       #  Describe the filtering method, or "Raw Data" if filtering was not used.
   #        Placement of text on Lomb-Scargle Periodogram.
   x_LS_wave = .0215
   y_LS_wave = 9.9
   x_LS_confid = x_LS_wave
   y_LS_confid = 0.92 * y_LS_wave
   LS_text_size = .8
   x_LS_alt = .012
   y_LS_alt = 6.9
   LS_alt_message = "54.2777-yr                            33.293-yr"
   prd_multiplier = 1             #  Multiplier for reporting maximum.
   prd_label = "yr"               #  Unit of measure for reporting maximum. 
   #        Placement of text on Smoothed Periodogram.
   ylim_min_smooth = 0
   ylim_max_smooth = .24
   x_smooth_wave = 1.00 * x_LS_wave
   y_smooth_wave = 0.982 * ylim_max_smooth
   x_smooth_confid = x_smooth_wave
   y_smooth_confid = 0.905 * ylim_max_smooth
   confid_indc = "Confidence Level exceeds 99%"
   smooth_text_size = .8
   x_smooth_alt = 0.028
   y_smooth_alt = .17
   smooth_alt_message = "33.574-yr"
   #     x-axis and y-axis labels for periodograms.
   LS_x_title = "Frequency (Cycles per Year)"
   LS_y_title = "Spectral Density"
   ts_x_label  = "Time (Years)"     #  x-axis label for time-series plot.
   time_scale  = "Yr"               #  Unit-of-Time for the analysis
   time_plural = "Yr"
}

#    Test # 5: Global Climate
if(test_no == 5) { 
   file_name = "Data/Climate_e_172-yr.txt"
   period_EUWS = 171.858328151339    #  The wavelength being investigated.
   period_solar= 700.441 
   search_adj_at_peaks = .012
   p_min  = 50                        #  Analyze wavelengths beginning with this minimum period.
   p_max  = 500                      #  Analyze wavelengths ending with this maximum period.
   start_time  = 0                   #  The time of the first observation.
   interval    = 19.0953697945932    #  Delta-t -- Interval between observations (in units = to timescale).
   Daniell_smooth = c(5,7)
   Daniell_null   = c(35,55)
   filter_meth = "Raw Data"          #  Describe the filtering method, or "Raw Data" if filtering was not used.
   #        Placement of text on Lomb-Scargle Periodogram.
   x_LS_wave = .0099
   y_LS_wave = 6.8
   x_LS_confid = x_LS_wave
   y_LS_confid = 0.92 * y_LS_wave
   LS_text_size = .8
   x_LS_alt = .0177
   y_LS_alt = 3.7
   LS_alt_message = "53.7939-yr"
   prd_multiplier = 1               #  Multiplier for reporting maximum.
   prd_label = "yr"                 #  Unit of measure for reporting maximum. 
   #        Placement of text on Smoothed Periodogram.
   ylim_min_smooth = 0
   ylim_max_smooth = .42
   x_smooth_wave = .97 * x_LS_wave
   y_smooth_wave = 0.982 * ylim_max_smooth
   x_smooth_confid = x_smooth_wave
   y_smooth_confid = 0.905 * ylim_max_smooth
   confid_indc = "Confidence Level exceeds 99.9%"
   smooth_text_size = .8
   x_smooth_alt = 0.0177
   y_smooth_alt = .24
   smooth_alt_message = "52.2937-yr"
   #     x-axis and y-axis labels for periodograms.
   LS_x_title = "Frequency (Cycles per Year)"
   LS_y_title = "Spectral Density"
   ts_x_label  = "Time (Years)"     #  x-axis label for time-series plot.
   time_scale  = "Yr"               #  Unit-of-Time for the analysis
   time_plural = "Yr"
}

#    Test # 6: Global Climate
if(test_no == 6) { 
   file_name = "Data/Climate_f_516-yr.txt"
   period_EUWS = 515.574984454017   #  The wavelength being investigated.
   period_solar= 700.441 
   search_adj_at_peaks = .012
   p_min  = 200                      #  Analyze wavelengths beginning with this minimum period.
   p_max  = 1000                     #  Analyze wavelengths ending with this maximum period.
   start_time  = 0                  #  The time of the first observation.
   interval    = 57.2861093837796   #  Delta-t -- Interval between observations (in units = to timescale).
   Daniell_smooth = c(5)
   Daniell_null   = c(9,19)
   filter_meth = "Raw Data"         #  Describe the filtering method, or "Raw Data" if filtering was not used.
   #        Placement of text on Lomb-Scargle Periodogram.
   x_LS_wave = .0028
   y_LS_wave = 6.4
   x_LS_confid = x_LS_wave
   y_LS_confid = 0.92 * y_LS_wave
   LS_text_size = .8
   x_LS_alt = .006
   y_LS_alt = -1000000
   LS_alt_message = "Alt: x-yr"
   prd_multiplier = 1               #  Multiplier for reporting maximum.
   prd_label = "yr"                 #  Unit of measure for reporting maximum. 
   #        Placement of text on Smoothed Periodogram.
   ylim_min_smooth = 0
   ylim_max_smooth = 1.0
   x_smooth_wave = .94 * x_LS_wave
   y_smooth_wave = 0.982 * ylim_max_smooth
   x_smooth_confid = x_smooth_wave
   y_smooth_confid = 0.905 * ylim_max_smooth
   confid_indc = "Confidence Level exceeds 75%"
   smooth_text_size = .8
   x_smooth_alt = 0.0022
   y_smooth_alt = -1000000
   smooth_alt_message = "Alt: x-yr (x% CL)"
   #     x-axis and y-axis labels for periodograms.
   LS_x_title = "Frequency (Cycles per Year)"
   LS_y_title = "Spectral Density"
   ts_x_label  = "Time (Years)"     #  x-axis label for time-series plot.
   time_scale  = "Yr"               #  Unit-of-Time for the analysis
   time_plural = "Yr"
}






####################################################################################################
###                                                                                              ###
###    Read in the time-series file.                                                             ###
###    Store the time-series to be analyzed in ExpressionProfile.                                ###
###                                                                                              ###
####################################################################################################

#      Read the input time series, and store it in x.

x = ts(scan(file_name), start=start_time, frequency=1/interval)

#      Calculate the 1st-order autocorrelation coefficient.
#      From Meko:  chapter 3 on Autocorrelation, formulas (2 & 3)
#      The 1st-order acc is used in the effective sample size calculation.

N  = length(x)
x1_mean = mean(x[1:N-1])
x2_mean = mean(x[2:N])
ds_mean = mean(x[1:N])
xnum = x
xden = x
xnum[N] = 0
xden[N] = 0
for (t in 1:N-1) {
  xnum[t] = ((x[t]-x1_mean) * (x[t+1]-x2_mean))
  xden[t] = ((x[t]-ds_mean) * (x[t]-ds_mean))
}
acc = mean(xnum) / mean(xden)

#      Calculate the effective sample size
#      From Meko:  chapter 3 on Autocorrelation, formulas (15)
#      Note:  For a highly autocorrelated time-series, as "acc" approaches 1,
#      the effective sample size approaches zero.

Eff_Samp_Size = N * (1-acc) / (1+acc)





####################################################################################################
###                                                                                              ###
###    Lomb-Scargle Periodogram function                                                         ###
###                                                                                              ###
####################################################################################################

LombScargle <- function(t, h, TestFrequencies, Nindependent)
{
  stopifnot( length(t) == length(h) )

  hResidual  = h - mean(h)

  SpectPowerDens = rep(0, length(TestFrequencies))

  for (i in 1:length(TestFrequencies))
  {
    Omega       = 2*pi*TestFrequencies[i]
    TwoOmegaT   = 2*Omega*t
    Tau         = atan2( sum(sin(TwoOmegaT)) , sum(cos(TwoOmegaT))) / (2*Omega)
    OmegaTMinusTau = Omega * (t - Tau)
    SpectPowerDens[i] = (sum (hResidual * cos(OmegaTMinusTau))^2) / sum( cos(OmegaTMinusTau)^2 ) +
                        (sum (hResidual * sin(OmegaTMinusTau))^2) / sum( sin(OmegaTMinusTau)^2 )
  }

  # The "normalized" spectral density refers to the variance term in the
  # denominator.  With this term the SpectPowerDens has an
  # exponential probability distribution with unit mean.

  SpectPowerDens = SpectPowerDens / ( 2 * var(h) )
  Probability = 1 - (1-exp(-SpectPowerDens))^Nindependent

  ##  print(SpectPowerDens)
  ##  print(1/TestFrequencies)

  PeakIndex    = match(max(SpectPowerDens), SpectPowerDens)
  LikelyPvalue = Probability[PeakIndex]

  #      Use Interpolation to improve on the estimated spectral peak.
  intvl        = (TestFrequencies[11]-TestFrequencies[1])/10
  pleft        = SpectPowerDens[PeakIndex-1]
  pcenter      = SpectPowerDens[PeakIndex]
  pright       = SpectPowerDens[PeakIndex+1]

  ##       Old 3-point interpolation method
  freq_est       = TestFrequencies[PeakIndex] + (.5*intvl*(pright-pleft)/(pcenter-(min(pleft,pright))))
  wavelength_peak= 1/TestFrequencies[PeakIndex]
  wavelength_est = 1/freq_est
  adjusted_est   = wavelength_est

  print("***********************************************************************************")
  print("Calculations from Lomb-Scargle function.")
  print("Spectral Peak")
  print(TestFrequencies[PeakIndex])
  print("Wavelength at Spectral Peak")
  print(prd_multiplier*wavelength_peak)
  print("Wavelength from 3-Point Linear Interpolation")
  print(prd_multiplier*wavelength_est)
  print(" ")

  ##       New Parabolic Interpolation method
  freq_est       =  TestFrequencies[PeakIndex] + intvl*(pright-pleft)/(2*(pcenter-pleft+pcenter-pright))
  wavelength_est =  1/freq_est
  number_cycles  =  (N-1) / (wavelength_est/interval)
  adj_factor     =  1 - ((1/(100*number_cycles^2))/2)
  adjusted_est   = adj_factor * wavelength_est

  print("*** Parabolic Interpolation ***")
  print("Frequency from Parabolic Interpolation")
  print(freq_est)
  print("Wavelength from Parabolic Interpolation")
  print(prd_multiplier*wavelength_est)
  print("Number of cycles in this time-series -- based on the above Wavelength estimate")
  print(number_cycles)
  print("Adjustment for reducing the estimate from Parabolic Interpolation")
  print(adj_factor)
  print("Final Adjusted Wavelength")
  print(prd_multiplier*adjusted_est)
  print("***********************************************************************************")

  LikelyPeriod = adjusted_est

  return( list( Probability=Probability,
                SpectPowerDens=SpectPowerDens,
                PeakIndex=PeakIndex,
                LikelyPeriod=LikelyPeriod,
                LikelyPvalue=LikelyPvalue
              ) )
}





####################################################################################################
###                                                                                              ###
###    Function to Draw Verical Frequency Lines on a Periodogram                                 ###
###                                                                                              ###
####################################################################################################

Draw_Vertical_Lines <- function(Solar_period, Solar_Indc, EUWS_Period, EUWS_P_Indc, EUWS_S_Indc, min_period, max_period, search_adj)
{

#        Calculate frequencies of EUWS harmonics
   e_period = 0

   e_period[1]  = (1/243)*EUWS_Period
   e_period[2]  = (2/243)*EUWS_Period
   e_period[3]  = (3/243)*EUWS_Period
   e_period[4]  = (6/243)*EUWS_Period
   e_period[5]  = (9/243)*EUWS_Period
   e_period[6]  = (18/243)*EUWS_Period

   e_period[7]  = (1/9)*EUWS_Period
   e_period[8]  = (2/9)*EUWS_Period
   e_period[9]  = (3/9)*EUWS_Period
   e_period[10] = (6/9)*EUWS_Period
   e_period[11] = 1*EUWS_Period
   e_period[12] = 2*EUWS_Period

   e_period[13] = 3*EUWS_Period
   e_period[14] = 6*EUWS_Period
   e_period[15] = 9*EUWS_Period
   e_period[16] = 18*EUWS_Period
   e_period[17] = 27*EUWS_Period
   e_period[18] = 54*EUWS_Period
   e_period[19] = 81*EUWS_Period
   e_period[20] = 162*EUWS_Period
   e_period[21] = 243*EUWS_Period

#        Calculate frequencies of Solar harmonics
#        Base solar period ~ 7.8075 days
   base_solar_period = .02137575
   s_period = 0
   s_period[1]  = .073922
   s_period[2]  = 1/8^4 * Solar_period
   s_period[3]  = 1/8^3 * Solar_period
   s_period[4]  = 1/8^2 * Solar_period
   s_period[5]  = 1/8^1 * Solar_period
   s_period[6]  = 8^0 * Solar_period
   s_period[7]  = 8^1 * Solar_period
   s_period[8]  = 8^2 * Solar_period
   s_period[9]  = 8^3 * Solar_period
   s_period[10] = 8^4 * Solar_period
   s_period[11] = 8^5 * Solar_period

#      Draw Primary EUWS harmonic lines.
   if(EUWS_P_Indc == TRUE) {
      for (i in 1:11) {
         k = (2*i) - 1
         x_line = c(1/e_period[k], 1/e_period[k])
         lines(x_line, c(-1000000000, 1000000000), lty=1, col="Black", lwd=2)
      }
   }

#      Draw Secondary EUWS harmonic lines.
   if(EUWS_S_Indc == TRUE) {
      for (i in 1:10) {
         j = (2*i)
         x_line = c(1/e_period[j], 1/e_period[j])
         lines(x_line, c(-1000000000, 1000000000), lty=1, col="Gray", lwd=2)
      }
   }

#      Draw Solar harmonic lines.
   if(Solar_Indc == TRUE) {
      for (i in 1:11) {
         x_line = c(1/s_period[i], 1/s_period[i])
         lines(x_line, c(-1000000000, 1000000000), lty=1, col="Yellow", lwd=2)
      }
   }

   output_s = NULL
   outs_max = 0
   for (i in 1:11) {
      if(s_period[i] > min_period) {
         if(s_period[i] < max_period) {
            FSP  = Find_Spectral_Peak(s_period[i],search_adj)
            output_s = c(output_s, FSP$out_line)
            outs_max = outs_max + 1
         }
      }
   }
   output_e = NULL
   oute_max = 0
   for (i in 1:21) {
      if(e_period[i] > min_period) {
         if(e_period[i] < max_period) {
            FSP  = Find_Spectral_Peak(e_period[i],search_adj)
            output_e = c(output_e, FSP$out_line)
            oute_max = oute_max + 1
         }
      }
   }

   return( list(output_s=output_s,
                outs_max=outs_max,
                output_e=output_e,
                oute_max=oute_max) )
}



####################################################################################################
###                                                                                              ###
###    Find the spectral peak and associated values within a specified range                     ###
###                                                                                              ###
####################################################################################################

Find_Spectral_Peak <- function(test_period, search_pct_adj)
{

   speak_start= (1/test_period) - (search_pct_adj*(1/test_period))
   speak_end  = (1/test_period) + (search_pct_adj*(1/test_period))

   Index1 = findInterval(speak_start,TestFrequencies)
   Index2 = findInterval(speak_end,TestFrequencies)

   p_Index = match(max(LS$SpectPowerDens[Index1:Index2]), LS$SpectPowerDens)

   speak_prob  = LS$Probability[p_Index]
   speak_power = LS$SpectPowerDens[p_Index]
   speak_freq  = TestFrequencies[p_Index]
   speak_theo  = prd_multiplier*test_period
   speak_period= prd_multiplier*(1/TestFrequencies[p_Index])
   speak_diff  = 100*((speak_period/speak_theo)-1)

   x1 = strtrim(speak_start, width=6)
   x2 = strtrim(speak_end,   width=6)
   x3 = strtrim(speak_freq,  width=6)
   x4 = strtrim(speak_theo,  width=7)
   x5 = strtrim(speak_period,width=7)
   x6 = strtrim(speak_diff,  width=6)
   x7 = strtrim(100*speak_power/LS$SpectPowerDens[LS$PeakIndex], width=5)
   x8 = strtrim(speak_prob,  width=6)
   out_line = paste(c(x1,x2,x3,x4,x5,x6,x7,x8))

   return( list(out_line = out_line) )
}




####################################################################################################
###                                                                                              ###
###    Set input parameters, then execute the Lomb-Scargle function.                             ###
###                                                                                              ###
####################################################################################################

ExpressionProfile = x

N        = length(ExpressionProfile) #  Number of observations in the time-series
tDelta   = interval                  #  Interval between sample points
t        = 0:(N-1) * tDelta          #  The time-series covers a period
M        = interval_fct * N          #  The number of frequency interval to create (ususlly either 2*N or 4*N)
frequency= 1/period_EUWS             #  The frequency
Nindependent = Eff_Samp_Size         #  Approx # of independent observations for autocorrelated series.

#     Set up the test frequencies (for the periods between p_min and p_max years).
TestFrequencies   = (1/p_max) + (1/p_min - 1/p_max)* 0:(M-1) / (M-1)

#     Call the Lomb-Scargle Periodogram function.
LS = LombScargle(t, ExpressionProfile, TestFrequencies, Nindependent)




####################################################################################################
###                                                                                              ###
###    Plot 2 graphs associated with the Lomb-Scargle Periodogram                                ###
###                                                                                              ###
###        1 -- The original time-series                                                         ###
###        2 -- Lomb-Scargle Periodogram                                                         ###
###                                                                                              ###
####################################################################################################

if(numb_of_plots == 1) {savepar = par(mfrow=c(1,1))}
if(numb_of_plots == 2) {savepar = par(mfrow=c(2,1))}
if(numb_of_plots == 3) {savepar = par(mfrow=c(3,1))}

#     Plot 1 -- the Original time-series
if(numb_of_plots > 2) {
   plot(t+start_time, ExpressionProfile,
   type = "o",
   col  = "Black",
   lwd  = 2,
   ylab = ts_y_label,
   xlab = ts_x_label,
   main = paste(c(ts_descr, "--", filter_meth, "( N =", N, ")"), collapse=" "))
}

#     Plot 2 -- Plot the Lomb-Scargle Periodogram
if(numb_of_plots < 4) {
   plot(TestFrequencies, LS$SpectPowerDens,
        type = "l",
        col  = "Blue",
        lwd  = 4,
        xlab = " ",
        ylab = " ",
        main = paste(c("Lomb-Scargle P'gram:", ts_descr, p_min, "to", p_max, time_scale), collapse=" "))

#  Draw a vertical line at the frequency being tested.
DVL = Draw_Vertical_Lines(period_solar, draw_solar_lines, period_EUWS, draw_EWUS_primary, draw_EWUS_secondary, p_min, p_max, search_adj_at_peaks)

#  Print text indicating the wavelength at the spectral peak.
   text(x_LS_wave, y_LS_wave,
        labels = paste("period =", round(LS$LikelyPeriod*prd_multiplier,5), prd_label),
        cex    = LS_text_size,
        font   =2,
        col    ="black",
        pos    = 4,
        offset = 0)

#  Print the p-value below the wavelength.
   text(x_LS_confid, y_LS_confid,
        labels = paste("p-value =", round(LS$LikelyPvalue,5)),
        cex    = LS_text_size,
        font   =2,
        col    ="black",
        pos    = 4,
        offset = 0)

#  Print the Alternative Message on the Lomb-Scargle Periodogram.
   text(x_LS_alt, y_LS_alt,
        labels = LS_alt_message,
        cex    = LS_text_size,
        font   =2,
        col    ="black",
        pos    = 4,
        offset = 0)

#  Use a red dot to mark the spectral peak.
   points(TestFrequencies[LS$PeakIndex], LS$SpectPowerDens[LS$PeakIndex],
        pch = 19,
        col = "Red")
}





####################################################################################################
###                                                                                              ###
###      Plot 3 -- The smoothed periodogram, null continuum, and confidence bands.               ###
###                                                                                              ###
####################################################################################################

if(numb_of_plots > 1) {

#      Calculate a smoothed periodogram.
smoothp = spec.pgram(x,
                  spans  = Daniell_smooth,
                  taper  = 0,
                  fast   = TRUE,
                  demean = TRUE,
                  detrend= TRUE,
                  plot   = FALSE,
                  log    = "no")

#      Calculate a red noise null-continuum.
nullcont = spec.pgram(x,
                  spans  = Daniell_null,
                  taper  = 0,
                  lwd    = 2,
                  fast   = TRUE,
                  demean = TRUE,
                  detrend= TRUE,
                  plot   = FALSE,
                  log    = "no")
nullMax = max(nullcont$spec)

#      Find the frequency-range to be plotted.
i_start = 1
i_end   = 20
testlen = length(TestFrequencies)
origlen = length(smoothp$spec)
for (i in 1:origlen) {
   if(smoothp$freq[i] < TestFrequencies[1]) {i_start = i}
   if(smoothp$freq[i] < TestFrequencies[testlen]) {i_end = i}
}
i_start = i_start + 1
if(i_start == 2)  {i_start = 1}

#      Use Interpolation to improve on the estimated spectral peak.
sIndex = 1
smooth_max = -1000000000000
for (i in i_start:i_end) {
   if(smoothp$spec[i] > smooth_max) {
      smooth_max = smoothp$spec[i]
      sIndex = i
   }
}
intvl    = (smoothp$freq[11]-smoothp$freq[1])/10
pleft    = smoothp$spec[sIndex-1]
pcenter  = smoothp$spec[sIndex]
pright   = smoothp$spec[sIndex+1]
##  freq_est = smoothp$freq[sIndex] + (.5*intvl*(pright-pleft)/(pcenter-(min(pleft,pright))))
freq_est = smoothp$freq[sIndex] + intvl*(pright-pleft)/(2*(pcenter-pleft+pcenter-pright))
wavelength_est =  1/freq_est
number_cycles  =  (N-1) / (wavelength_est/interval)
adj_factor     =  1 - ((1/(100*number_cycles^2))/2)
sWave          = adj_factor * wavelength_est
sFreq          = 1/sWave

#      Customize the Labels
LS_x_label = title(xlab=LS_x_title, cex.lab=LS_size, font.lab=LS_font, col.lab=LS_color)
LS_y_label = title(ylab=LS_y_title, cex.lab=LS_size, font.lab=LS_font, col.lab=LS_color)

#      Plot the smoothed periodogram.
plot(smoothp$freq[i_start:i_end], smoothp$spec[i_start:i_end],
           type = "l",
           col  = "Blue",
           lwd  = 4,
           ylim = c(ylim_min_smooth, ylim_max_smooth),
           xlab = " ",
           ylab = " ",
           sub  = LS_s_title,
           cex.sub=LS_s_size,
           font.sub = LS_s_font,
           main = paste(c("Smoothed P'gram with Confidence-Bands:", ts_descr), collapse=" "))

#      Draw a vertical line at the frequency being tested.
DVL = Draw_Vertical_Lines(period_solar, draw_solar_lines, period_EUWS, draw_EWUS_primary, draw_EWUS_secondary, p_min, p_max, search_adj_at_peaks)

#      Print the estimated wavelength at the spectral peak.
text(x_smooth_wave, y_smooth_wave,
           labels = paste("period =", round(sWave*prd_multiplier,5), prd_label),
           cex  = smooth_text_size,
           font =2,
           col  ="black",
           pos    = 4,
           offset = 0)

#      Use a red dot to mark the maximum spectral density.
points(smoothp$freq[sIndex], smoothp$spec[sIndex], pch = 19, col = "Red")

#      Plot the null-continuum.
lines(nullcont$freq,nullcont$spec, lty=1, col="Red", lwd=3)

#      Draw Confidence Bands using Meko's method.
adj_75 = smoothp$df/qchisq(.75, smoothp$df)
Meko_75 = ts((smoothp$spec)*adj_75)
lines(smoothp$freq, Meko_75, lty=1,col=11, lwd=2)
CI75_Max = max(Meko_75[i_start:i_end])

adj_95 = smoothp$df/qchisq(.95, smoothp$df)
Meko_95 = ts((smoothp$spec)*adj_95)
lines(smoothp$freq, Meko_95, lty=1,col="Green", lwd=2)
CI95_Max = max(Meko_95[i_start:i_end])

adj_99 = smoothp$df/qchisq(.99, smoothp$df)
Meko_99 = ts((smoothp$spec)*adj_99)
lines(smoothp$freq, Meko_99, lty=1,col="Yellow", lwd=2)
CI99_Max = max(Meko_99[i_start:i_end])

adj_999 = smoothp$df/qchisq(.999, smoothp$df)
Meko_999 = ts((smoothp$spec)*adj_999)
lines(smoothp$freq, Meko_999, lty=1,col="Pink", lwd=2)
CI999_Max = max(Meko_999[i_start:i_end])

adj_9999 = smoothp$df/qchisq(.9999, smoothp$df)
Meko_9999 = ts((smoothp$spec)*adj_9999)
lines(smoothp$freq, Meko_9999, lty=1,col="Purple", lwd=2)
CI9999_Max = max(Meko_9999[i_start:i_end])

#      Print the Confidence Level below the spectral peak.
text(x_smooth_confid, y_smooth_confid,
           labels = confid_indc,
           cex  = smooth_text_size,
           font =2,
           col  ="black",
           pos    = 4,
           offset = 0)

#      Print the Alternative Message on the Smoothed Periodogram.
text(x_smooth_alt, y_smooth_alt,
           labels = smooth_alt_message,
           cex  = smooth_text_size,
           font =2,
           col  ="black",
           pos    = 4,
           offset = 0)

#     Customize the Labels
LS_x_label = title(xlab=LS_x_title, cex.lab=LS_size, font.lab=LS_font, col.lab=LS_color)
LS_y_label = title(ylab=LS_y_title, cex.lab=LS_size, font.lab=LS_font, col.lab=LS_color)

#      Confidence Bands using Arunajadai method.
###  adj_75 = nullcont$df/qchisq(1-.75, nullcont$df)
###  Arun_75 = nullcont$spec * adj_75
###  lines(nullcont$freq, Arun_75, lty=2,col=11, lwd=2)

###  adj_95 = nullcont$df/qchisq(1-.95, nullcont$df)
###  Arun_95 = nullcont$spec * adj_95
###  lines(nullcont$freq, Arun_95, lty=2,col="Green", lwd=2)

###  adj_99 = nullcont$df/qchisq(1-.99, nullcont$df)
###  Arun_99 = nullcont$spec * adj_99
###  lines(nullcont$freq, Arun_99, lty=2,col="Yellow", lwd=2)

###  adj_999 = nullcont$df/qchisq(1-.999, nullcont$df)
###  Arun_999 = nullcont$spec * adj_999
###  lines(nullcont$freq, Arun_999, lty=2,col="Pink", lwd=2)

###  adj_9999 = nullcont$df/qchisq(1-.9999, nullcont$df)
###  Arun_9999 = nullcont$spec * adj_9999
###  lines(nullcont$freq, Arun_9999, lty=2,col="Purple", lwd=2)

}

par(savepar)





####################################################################################################
###                                                                                              ###
###      Write the results to the file Spectrum_Output.txt                                       ###
###      To be imported into a spreadsheet at a later time.                                      ###
###                                                                                              ###
####################################################################################################

if(Write_Output_file == TRUE) {
   if(numb_of_plots > 1) {

      outlen = length(smoothp)
      cat(" ", file="Spectrum_Output.txt", sep="", fill=FALSE, append=FALSE, "\n")
      for(i in 1:outlen) cat(smoothp$freq[i],
                   1/smoothp$freq[i],
                   Meko_95[i], Meko_99[i], Meko_999[i], Meko_9999[i],
                   smoothp$spec[i], 
                   nullcont$spec[i],
                   file="Spectrum_Output.txt",
                   sep=", ",
                   fill=FALSE, 
                   append=TRUE, "\n")
      cat("  ", file="Spectrum_Output.txt", sep="", fill=FALSE, append=TRUE, "\n")

      outlen = length(TestFrequencies)
      for(i in 1:outlen) cat(TestFrequencies[i],
                   1/TestFrequencies[i],
                   LS$SpectPowerDens[i], 
                   file="Spectrum_Output.txt",
                   sep=", ",
                   fill=FALSE, 
                   append=TRUE, "\n")
      cat("  ", file="Spectrum_Output.txt", sep="", fill=FALSE, append=TRUE, "\n")
   }
}



####################################################################################################
###                                                                                              ###
###      Print a summary of the statistics calculated in this program.                           ###
###      Then exit.                                                                              ###
###                                                                                              ###
####################################################################################################

print("1st-Order Autocorrelation Coefficient")
print(acc)

print("Number of sample-points in the time-series")
print(N)

print("Effective Sample Size")
print(Eff_Samp_Size)