*
*  @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 coFREQORDs]
*    [GRAPH]/NOGRAPH
*    HEADER=Header for high-resolution graph
*
*  Revision Schedule
*    05/92 Written by T. Doan
*    10/00 Correction to RTOC instruction to deal with passing reduced range of series
*
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

local integer  nobs startl endl span freqord polesl
local series   icorrs  spectrum  harmonic

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 FREQORDs
*

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

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

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

corr(inverse=icorrs,covar,noprint,number=polesl) 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 graph
   {
   cln 1
   ctor 1 (freqord+1)/2
   #    1
   #  spectrum
   set harmonic 1 (freqord+1)/2 = %clock(t,freqord/span)==1
   graph(noticks,grid=harmonic,header=header) 1
   # spectrum 1 (freqord+1)/2
   }
end mesa

