PRO PlotSpectra

; =============================================================================
; Code to plot horizontal spectra as a function of wavenumber and vertical mode
; =============================================================================

;.compile PlotSpectra.pro
;.compile BoxPlot.pro



;-----------------------------------------------------------
; Set-up 16 colours + 0=black, 1=white
;-----------------------------------------------------------
r = [0,255,0,0,0,0,66,184,245,255,255,255,230,191,148,105,161,184]
g = [0,255,84,199,255,255,255,255,255,209,135,28,0,0,0,0,0,0]
b = [0,255,255,255,255,178,0,0,0,0,0,0,0,0,0,20,135,186]
TVLCT, r,g,b



;-----------------------------------------------------------
; Set-up data and axis arrays
;-----------------------------------------------------------
NoOfWavenumbers = 143
NoOfVertModes   = 38

Data        = FLTARR(NoOfWavenumbers+1,NoOfVertModes,4,2)
PlotData    = FLTARR(NoOfWavenumbers+1,NoOfVertModes)
WNarray     = FINDGEN(NoOfWavenumbers+1)
VMarray     = FINDGEN(NoOfVertModes)
StartWn     = FLTARR(5)
NumModes    = FLTARR(5)

;Different parameters have different numbers of wavenumbers
StartWn(1)  = 1   ;for psi
StartWn(2)  = 1   ;for chi
StartWn(3)  = 0   ;for Ap
StartWn(4)  = 0   ;for mu

;Different parameters have different number of vertical modes
NumModes(1) = 38  ;for psi
NumModes(2) = 38  ;for chi
NumModes(3) = 38  ;for Ap
NumModes(4) = 28  ;for mu


FOR loop = 0, 2 DO BEGIN
  ;loop=0 is for standard data
  ;loop=1 is for PV data
  ;loop=2 is for difference PV-standard

  PRINT, '================================='
  PRINT, loop
  PRINT, '================================='


  CASE loop OF
    0 : BEGIN
          ;Deal with standard data
          Readfile = 1
          file     = 'Spectrum0_St'
          title    = 'Log spectrum (standard scheme)'
          PSfile   = 'LogSpectrum_St.ps'
        END
    1 : BEGIN
          ;Deal with PV data
          Readfile = 1
          file     = 'Spectrum0_PV'
          title    = 'Log spectrum (PV scheme)'
          PSfile   = 'LogSpectrum_PV.ps'
        END
    2 : BEGIN
          ;Deal with difference PV-standard data
          Readfile = 0
          file     = ''
          title    = 'Spectrum (PV-standard)'
          PSfile   = 'Spectrum_Diff.ps'
        END
  ENDCASE



  linestring = ''
  a          = 0
  b          = 0.0

  ;Set-up postscript file for output
  SET_PLOT, 'ps'
  DEVICE, FILE=PSfile, /PORTRAIT, /COLOR, XSIZE=21.0, YSIZE=29.5, $
          XOFFSET=2.0, YOFFSET=2.0, FONT_SIZE=8

  IF (Readfile EQ 1) THEN BEGIN
    OPENR, unit, file, /GET_LUN
  ENDIF

  FOR para = 1, 4 DO BEGIN
    ;Loop over parameter (psi, chi, Ap, mu)
    PRINT, 'Parameter No. ',para

    IF (Readfile EQ 1) THEN BEGIN
      ;Read-in data for this parameter into an array
      READF, unit, linestring
      FOR mode = 1, NumModes(para) DO BEGIN
        READF, unit, linestring
        READF, unit, linestring
        READF, unit, linestring
        FOR wn = StartWn(para), NoOfWavenumbers DO BEGIN
          READF, unit, a, b
          Data(wn,mode-1,para-1,loop) = b       ;Store for later use
          PlotData(wn,mode-1)         = ALOG(b) ;This is the data array that will be plotted
        ENDFOR
      ENDFOR
    ENDIF ELSE BEGIN
      ;Calculate the difference between the previous two loops
      PlotData(*,*) = Data(*,*,para-1,1) - Data(*,*,para-1,0)
    ENDELSE

    ;Plot these data
    BoxPlot, PlotData(StartWn(para):NoOfWavenumbers, 0:NumModes(para)-1), 0,        $
             WNarray(StartWn(para):NoOfWavenumbers), VMarray(0:NumModes(para)-1),   $
             0.1, 0.95-0.22*para, 0.6, 0.9-0.22*(para-1),                           $
             0.0, title, 0, 0, 1, 0, 0, 0

  ENDFOR

  IF (Readfile EQ 1) THEN BEGIN
    CLOSE, unit
    FREE_LUN, unit
  ENDIF

ENDFOR


DEVICE, /CLOSE

END