;------------------------------------------------------------------------------
;            -- Utility routine to output PS and GIF files --
;------------------------------------------------------------------------------

PRO make_ps_gif, file_name, file_tag, rotate

    DEVICE, /CLOSE
    psfile = file_name + file_tag + '.ps'
    giffile = file_name + file_tag + '.gif'
    SPAWN, 'mv idl.ps ' + psfile
    IF ( rotate EQ 1 ) THEN BEGIN
      SPAWN, 'convert -rotate 270 ' + psfile + ' ' + giffile
    ENDIF ELSE BEGIN
      SPAWN, 'convert ' + psfile + ' ' + giffile
    ENDELSE
    PRINT, FORMAT = '( 5x, "Files ", a, 2x, a, " written" )', psfile, giffile

END


;------------------------------------------------------------------------------
;                          -- Main routine --
;------------------------------------------------------------------------------

PRO plot_fft_test, orig_input, idlfft, realft, four1, $
                   ps = ps, scale = scale


; --------------------------
; Establish plotting colours
; --------------------------

  TVLCT, red_save, green_save, blue_save, /GET
  red   = [ 255B, 0, 0 ]
  green = [ 0B, 255, 0 ]
  blue  = [ 0B, 0, 255 ]
  TVLCT, red, green, blue, 1


; ---------------------------
; Setup for PostScript output
; ---------------------------

  IF ( KEYWORD_SET( ps ) ) THEN BEGIN
    IF ( NOT KEYWORD_SET( scale ) ) THEN BEGIN
      scale = 1.0
      file_tag = ''
      portrait = 0
      landscape = 1
    ENDIF ELSE BEGIN
      file_tag = '_tn'
      portrait = 1
      landscape = 0
    ENDELSE
    rotate = landscape
    SET_PLOT, 'PS'
    DEVICE, /COLOR, PORTRAIT = portrait, landscape = LANDSCAPE, $
            BITS_PER_PIXEL = 8, SCALE_FACTOR = scale
    !P.FONT = 0
  ENDIF


; ------------------------
; Plot original input data
; ------------------------

  xtickformat = '( i5 )'
  xtitle = 'Point number'
  ytitle = 'Counts'
  title  = 'AERI-TWP1 HBB forward scan average uncalibrated spectrum'
  xrange = [ 0, 16384 ]
  xticks = 8

  PLOT, orig_input, $
        XTICKFORMAT = xtickformat, $
        TITLE = title, $
        XTITLE = xtitle, $
        YTITLE = ytitle, $
        XRANGE = xrange, $
        /XSTYLE, $
        XTICKS = xticks, $
        /NODATA
  OPLOT, orig_input, COLOR = 1
  OPLOT, IMAGINARY( orig_input ), COLOR = 2

  XYOUTS, 0.5, 0.6, 'Real part', /NORM, COLOR = 1
  XYOUTS, 0.5, 0.55, 'Imaginary part', /NORM, COLOR = 2
  
  IF ( KEYWORD_SET( ps ) ) THEN BEGIN
    make_ps_gif, 'orig_input', file_tag, rotate
  ENDIF ELSE BEGIN
    PRINT, FORMAT = '( 5x, "Press <RETURN> for next plot...", $ )'
    result = GET_KBRD( 1 )
  ENDELSE


; -------------------
; Plot FFT input data
; -------------------

  r = FLOAT( orig_input )
  i = IMAGINARY( orig_input )
  n = N_ELEMENTS( r )

  r = [ r, REVERSE( r[ 1 : n - 2 ] ) ]
  i = [ i, -1.0 * REVERSE( i[ 1 : n - 2 ] ) ]

  fft_input = COMPLEX( TEMPORARY( r ), TEMPORARY( i ) )

  xtickformat = '( i5 )'
  xtitle = 'Point number'
  ytitle = 'Counts'
  title  = 'Reflected data for input to complex FFT'
  xrange = [ 0, 32768 ]
  xticks = 8

  PLOT, fft_input, $
        XTICKFORMAT = xtickformat, $
        TITLE = title, $
        XTITLE = xtitle, $
        YTITLE = ytitle, $
        XRANGE = xrange, $
        /XSTYLE, $
        XTICKS = xticks, $
        /NODATA
  OPLOT, fft_input, COLOR = 1
  OPLOT, IMAGINARY( fft_input ), COLOR = 2

  XYOUTS, 0.5, 0.6, 'Real part', /NORM, COLOR = 1
  XYOUTS, 0.5, 0.55, 'Imaginary part', /NORM, COLOR = 2
  
  IF ( KEYWORD_SET( ps ) ) THEN BEGIN
    make_ps_gif, 'fft_input', file_tag, rotate
  ENDIF ELSE BEGIN
    PRINT, FORMAT = '( /5x, "Press <RETURN> for next plot...", $ )'
    result = GET_KBRD( 1 )
  ENDELSE


; --------------------
; Plot raw FFT results
; --------------------
 
  n = N_ELEMENTS( idlfft )
  n21 = n / 2 + 1L
  f = LINDGEN( N )
  f[ n21 ] = n21 - n + LINDGEN( n21 - 2 )

  xtickformat = '( i4 )'
  xtitle = 'Point number'
  ytitle = 'IFG magnitude'
  title  = 'Output interferograms (near ZPD) from FFT algorithms'
  xrange = [ -32, 32 ]
  xticks = 8

  PLOT, SHIFT( f, -n21 ), SHIFT( idlfft, -n21 ), $
        XTICKFORMAT = xtickformat, $
        TITLE = title, $
        XTITLE = xtitle, $
        YTITLE = ytitle, $
        XRANGE = xrange, $
        /XSTYLE, $
        XTICKS = xticks, $
        /NODATA
  OPLOT, SHIFT( f, -n21 ), SHIFT( idlfft, -n21 ), COLOR = 1
  OPLOT, SHIFT( f, -n21 ), SHIFT( realft, -n21 ), COLOR = 2, PSYM = 4, SYMSIZE = 1.5
  OPLOT, SHIFT( f, -n21 ), SHIFT( four1, -n21 ), COLOR = 3, PSYM = 1
        
  xsave = 0.15
  ysave = 0.9
  x = xsave
  y = ysave
  XYOUTS, x, y, 'IDL FFT', /NORM, COLOR = 1
  y = y - 0.05
  PLOTS, x, y + 0.005, /NORM, PSYM = 4, COLOR = 2
    XYOUTS, x + 0.02, y, 'Fortran REALFT', /NORM, COLOR = 2
  y = y - 0.05
  PLOTS, x, y + 0.005, /NORM, PSYM = 1, COLOR = 3
    XYOUTS, x + 0.02, y, 'Fortran FOUR1', /NORM, COLOR = 3

  IF ( KEYWORD_SET( ps ) ) THEN BEGIN
    make_ps_gif, 'fft_output', file_tag, rotate
  ENDIF ELSE BEGIN
    PRINT, FORMAT = '( /5x, "Press <RETURN> for next plot...", $ )'
    result = GET_KBRD( 1 )
  ENDELSE


; ----------------------------
; Plot massaged #1 FFT results
; ----------------------------
 
  xtickformat = '( i4 )'
  xtitle = 'Point number'
  ytitle = 'IFG magnitude'
  title  = 'Scaled output interferograms (near ZPD) from FFT algorithms'
  xrange = [ -32, 32 ]
  xticks = 8

  PLOT, SHIFT( f, -n21 ), SHIFT( idlfft, -n21 ) / 2.0, $
        XTICKFORMAT = xtickformat, $
        TITLE = title, $
        XTITLE = xtitle, $
        YTITLE = ytitle, $
        XRANGE = xrange, $
        /XSTYLE, $
        XTICKS = xticks, $
        /NODATA
  OPLOT, SHIFT( f, -n21 ), SHIFT( idlfft, -n21 ) / 2.0, COLOR = 1
  OPLOT, SHIFT( f, -n21 ), SHIFT( realft, -n21 ), COLOR = 2, PSYM = 4, SYMSIZE = 1.5
  OPLOT, SHIFT( f, -n21 ), SHIFT( four1, -n21 ) / 2.0, COLOR = 3, PSYM = 1
        
  x = xsave
  y = ysave
  XYOUTS, x, y, 'IDL FFT / 2', /NORM, COLOR = 1
  y = y - 0.05
  PLOTS, x, y + 0.005, /NORM, PSYM = 4, COLOR = 2
    XYOUTS, x + 0.02, y, 'Fortran REALFT', /NORM, COLOR = 2
  y = y - 0.05
  PLOTS, x, y + 0.005, /NORM, PSYM = 1, COLOR = 3
    XYOUTS, x + 0.02, y, 'Fortran FOUR1 / 2', /NORM, COLOR = 3

  IF ( KEYWORD_SET( ps ) ) THEN BEGIN
    make_ps_gif, 's_fft_output', file_tag, rotate
  ENDIF ELSE BEGIN
    PRINT, FORMAT = '( /5x, "Press <RETURN> for next plot..." )'
    result = GET_KBRD( 1 )
  ENDELSE


; ----------------------------
; Plot massaged #2 FFT results
; ----------------------------
 
  xtickformat = '( i4 )'
  xtitle = 'Point number'
  ytitle = 'IFG magnitude'
  title  = 'Reversed, scaled output interferograms (near ZPD) from FFT algorithms'
  xrange = [ -32, 32 ]
  xticks = 8

  PLOT, SHIFT( f, -n21 ), SHIFT( idlfft, -n21 ) / 2.0, $
        XTICKFORMAT = xtickformat, $
        TITLE = title, $
        XTITLE = xtitle, $
        YTITLE = ytitle, $
        XRANGE = xrange, $
        /XSTYLE, $
        XTICKS = xticks, $
        /NODATA
  OPLOT, SHIFT( f, -n21 ), SHIFT( idlfft, -n21 ) / 2.0, COLOR = 1
    temp_realft = [ realft[ 0 ], REVERSE( realft[ 1 : n - 1 ] ) ]
  OPLOT, SHIFT( f, -n21 ), SHIFT( temp_realft, -n21 ), COLOR = 2, PSYM = 4, SYMSIZE = 1.5
    temp_four1 = [ four1[ 0 ], REVERSE( four1[ 1 : n - 1 ] ) ]
  OPLOT, SHIFT( f, -n21 ), SHIFT( temp_four1, -n21 ) / 2.0, COLOR = 3, PSYM = 1
        
  x = xsave
  y = ysave
  XYOUTS, x, y, 'IDL FFT / 2', /NORM, COLOR = 1
  y = y - 0.05
  PLOTS, x, y + 0.005, /NORM, PSYM = 4, COLOR = 2
    XYOUTS, x + 0.02, y, 'Reversed Fortran REALFT', /NORM, COLOR = 2
  y = y - 0.05
  PLOTS, x, y + 0.005, /NORM, PSYM = 1, COLOR = 3
    XYOUTS, x + 0.02, y, 'Reversed Fortran FOUR1 / 2', /NORM, COLOR = 3

  IF ( KEYWORD_SET( ps ) ) THEN BEGIN
    make_ps_gif, 'rs_fft_output', file_tag, rotate
    SET_PLOT, 'X'
    !P.FONT = -1
  ENDIF ELSE BEGIN
    PRINT, FORMAT = '( /5x, "Press <RETURN> to exit..." )'
    result = GET_KBRD( 1 )
    WDELETE, !D.WINDOW
  ENDELSE


; ------------------------
; Restore old colour table
; ------------------------

  TVLCT, red_save, green_save, blue_save

END