*
* @MESA( options )   series   start  end
* Computes and optionally graphs a spectrum using Maximum Entropy Method
*
* Options:
*   ORDINATES=Number of ordinates [depends upon length,seasonal]
*   POLES=Number of poles [depends upon number of observations]
*   [GRAPH]/NOGRAPH
*   HEADER=Header for graph
*   FOOTER=Footer for graph
*   [LOGSCALE]/NOLOGSCALE Graph spectrum in log scale or not.
*
* Revision Schedule
*   05/1992 Written by T. Doan
*   10/2000 Correction to RTOC instruction to deal with passing reduced
*           range of series
*   07/2009 Add FOOTER and LOGSCALE options. Switch to METHOD=BURG.
*           Use SCATTER rather than GRAPH.
*
procedure mesa series start end
type series     series
type integer    start end
*
option integer  ordinate
option integer  poles
option switch   graph       1
option string   header
option string   footer
option switch   logscale    1

local integer  nobs startl endl span freqord polesl half
local series   icorrs  spect  harmonic frequencies

if .not.%defined(series) {
   display "Syntax: @MESA series  start end"
   return
}

*
*  Determine the range of SERIES to use
*

inquire(series=series) startl<<start endl<<end
compute nobs=endl-startl+1

*
* Determine the ordinates. If the user provides a number, use it. If
* not, use the number of observations, if the user has requested the
* periodogram, or the recommended number, if a more general spectral
* density estimate is requested.
*

if %defined(ordinate).and.ordinate>=nobs
   compute freqord=ordinate
else
   compute freqord=%freqsize(nobs)

compute half = freqord/2+1

if .not.%defined(poles).or.poles<=0
  compute polesl=fix(sqrt(nobs))
else
  compute polesl=poles

*
* Determines where to show harmonics in the graph
*

if %clock(freqord,12)==12
   compute span=12
else
   compute span=4

*
* Compute inverse autocovariances
*

corr(inverse=icorrs,covar,noprint,number=polesl,method=burg) series startl endl

frequency 1 freqord
rtoc 1 poles+1
# icorrs
#   1

fft 1
cmult 1 1
cset 1 = icorrs(1)/%z(t,1)
*
if .not.graph
   return

ctor 1 half
#    1
#  spect
set harmonic    1 half = %clock(t,freqord/span)==1
set frequencies 1 half = (t-1.0)/half
if logscale {
   scatter(style=lines,header=header,footer=footer,hlabel="Fractions of pi",vlog=10.0)
   # frequencies spect 1 half
}
else {
   scatter(style=lines,header=header,footer=footer,hlabel="Fractions of pi")
   # frequencies spect 1 half
}

end mesa

