c
c
      subroutine FFTPLOT ( t, cmplxts, ndp, delt, nf, freq, cwrk1,
c
     1                     cwrk2, cwrk3, ftreal, 
c
     2                     ftimag, title, 
c
     3                     preal, pimag, power, pwrdb, ts )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      character*(*) title
c
      logical ts, skip1, skip2
c
      dimension  t(ndp) 
c
      complex  cmplxts(ndp), cwrk1(nf), cwrk2(nf), cwrk3(nf)
c
      dimension  freq(nf), ftreal(nf), ftimag(nf), f0(3), zero(3),
c
     1           preal(nf), pimag(nf), power(nf), pwrdb(nf)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c        t = array of time points for the input time series
c
c  cmplxts = complex array containing the times series to be analyzed
c
c      ndp = number of terms in the time series to be analyzed
c
c     delt = sample spacing for the time series ( delt = t(i+1) - t(i) )
c
c       nf = number of frequencies desired for the Fourier transform
c            (must be a power of two and .ge. ndp)
c
c     freq = array of frequencies for the Fourier transform (computed herein)
c
c    cwrk1 = complex vector of working storage
c
c    cwrk2 = complex vector of working storage
c
c    cwrk3 = complex vector of working storage
c
c   ftreal = real vector containing real components of DFT }  computed herein
c                                                          }        and
c   ftimag = real vector containing imaginary components   }      returned
c
c       ts = logical variable indicating that the data are a time series
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c
      skip1 = .false.
      skip2 = .false.
c
      delf = 1.0 / ( FLOAT( nf ) * delt )
c
      nfo2 = nf / 2
c
      do j = 1, nf
         freq(j) = FLOAT( j - 1 - nfo2 ) * delf
      end do
c
c    ------------------------------------
      call PREFFT ( nf, 0, nexp, cwrk1 )
c    ------------------------------------
c
      do i = 1, ndp
         cwrk2(i) = cmplxts(i)
      end do
c
      if ( ndp .lt. nf ) then
                              do i = ndp+1, nf
                                 cwrk2(i) = CMPLX( 0.0, 0.0 )
                              end do
      end if
c
c    ----------------------------------------------
      call FFT ( nf, 0, delt, nexp, cwrk1, cwrk2 )
c    ----------------------------------------------
c                             FFT replaces zero-padded complex time series in
c                                                    array cwrk2 with its DFT
c ---  Sort DFT into proper frequency order
c
      do j = 1, nfo2
         cwrk3(j) = cwrk2(nfo2+j)
      end do
c
      do j = nfo2+1, nf
         cwrk3(j) = cwrk2(j-nfo2)
      end do
c
      do j = 1, nf
         ftreal(j) = REAL ( cwrk3(j) )
         ftimag(j) = AIMAG( cwrk3(j) )
      end do
c
      f0(1) = freq(1)
      f0(2) = 0.0
      f0(3) = freq(nf)
c
      zero(1) = 0.0
      zero(2) = 0.0
      zero(3) = 0.0
c
c   -------------------------------------------------------
      call TWOPLT ( freq, ftreal, nf, f0, zero, 3, 0, 1,
c
     1              'Real Part{ DFT[ '// title //' ] }',
c
     2              'Frequency [kHz]', 'Re{ DFT }',
c
     3              freq, ftimag, nf, f0, zero, 3, 0, 1,
c
     4              'Imag. Part{ DFT[ ' // title //' ] }',
c
     5              'Frequency [kHz]', 'Im{ DFT }', ts )
c   -------------------------------------------------------
c
c  Compute and plot separate power spectra for Real and Imaginary parts of DFT
c
      do j = 1, nf
         preal(j) = ftreal(j)**2 / float( ndp )
         pimag(j) = ftimag(j)**2 / float( ndp )
         power(j) = preal(j) + pimag(j) 
      end do
c
c   -----------------------------------------------------
      call TWOPLT ( freq, preal, nf, f0, zero, 3, 0, 1,
c
     1              'Power[ Real Part{ DFT } ]',
c
     2              'Frequency [kHz]', 'Power',
c
     3              freq, pimag, nf, f0, zero, 3, 0, 1,
c
     4              'Power[ Imag. Part{ DFT } ]',
c
     5              'Frequency [kHz]', 'Power', ts )
c   -------------------------------------------------------
c
      prmax = 0.0
      pimax = 0.0
      pwmax = 0.0
c
      do j = 1, nf
c
         if ( (preal(j) .eq. 0.0) .or. (pimag(j) .eq. 0.0) ) then
                print *, ' preal(', j, ') =', preal(j)
                print *, ' pimag(', j, ') =', pimag(j)
                skip1 = .true.
         end if
c
         if ( preal(j) .gt. prmax )  prmax = preal(j)
         if ( pimag(j) .gt. pimax )  pimax = pimag(j)
c
         if ( power(j) .eq. 0.0 ) then
                print *, ' power(', j, ') =', power(j)
                skip2 = .true.
         end if
c
         if ( power(j) .gt. pwmax )  pwmax = power(j)
c
      end do
c
      do j = 1, nf
         if ( .not. skip1 ) preal(j) = 10.0 * LOG10( preal(j) / prmax )
         if ( .not. skip1 ) pimag(j) = 10.0 * LOG10( pimag(j) / pimax )
         if ( .not. skip2 ) pwrdb(j) = 10.0 * LOG10( power(j) / pwmax )
      end do
c
c
      if ( .not. skip1 ) then
c           ----------------------------------------------------------
             call TWOPLT ( freq, preal, nf, f0, zero, 3, 0, 1,
c
     1                     'Power[ Real Part{ DFT } ]',
c
     2                     'Frequency [kHz]', 'Rel. Power [db]',
c
     3                     freq, pimag, nf, f0, zero, 3, 0, 1,
c
     4                     'Power[ Imag. Part{ DFT } ]',
c
     5                    'Frequency [kHz]', 'Rel. Power [db]', ts )
c           -----------------------------------------------------------
      end if
c
c
      if ( .not. skip2 ) then
c           -----------------------------------------------------------
             call TWOPLT ( freq, power, nf, f0, zero, 3, 0, 1,
c
     1                     'Periodogram [ ' // title // ' ]',
c
     2                     'Frequency [kHz]', 'Power',
c
     3                     freq, pwrdb, nf, f0, zero, 3, 0, 1,
c   
     4                     'Periodogram [ ' // title // ' ]',
c
     5                     'Frequency [kHz]', 'Rel. Power [db]', ts )
c           -------------------------------------------------------------
      end if
c
c    ++++++++
      return
c    ++++++++
c
      end


