C* PS.F General routines for power spectra
C 
C+
      SUBROUTINE PSPEC(DATA,NDAT,NFFT,WINDOW,TSAMP,WORK,STATUS)
      IMPLICIT NONE
      INTEGER NDAT,NFFT,STATUS
      REAL*8 DATA(NFFT),TSAMP,WORK(NFFT)
      CHARACTER*(*) WINDOW
C
C Calculates the power spectrum of a function, multiplying first by
C a window function.
C Arguments:
C  Input:
C    NDAT   No of data samples
C    NFFT   Size of FFT. If NFFT>NDAT the data is padded with zeroes 
C           (there must be enough space in DATA to do this)
C           If NFFT<NDAT, it is truncated.
C    WINDOW Type of window to multiply data before Fourier transforming.
C           Can be one of "PARZEN", "TOPHAT" or "HANNING". In each case
C           the window goes to zero outside the range given by 
C           MIN(NDAT,NFFT) and has unit rms amplitude, the average being
C           taken over NFFT samples.
C    TSAMP  Sample time in seconds.
C  Altered:
C    DATA   On input, contains the data samples. On output it contains 
C           the power spectrum W(0:NFFT/2), normalised s.t. 
C              \sum_{I=-(NFFT/2-1)}^{NFFT/2} W(I)*DELTA_F = <DATA(I)**2>
C            where
C                DELTA_F=1/(NFFT*TSAMP)
C            i.e. it is normalised in Nyingi**2/Hz if
C            DATA is sampled in units of Nyingi at intervals of TSAMP seconds
C
C            Note that this is a TWO-SIDED power spectrum definition, even
C            though only one side is returned. Note also that this works
C            for NFFT not a multiple of 2, by changing NFFT/2 -> (NFFT-1)/2
C            and (NFFT/2-1) -> (NFFT-1)/2 in the above definitions. 
C    WORK    Array of size at least NFFT for use as workspace
C    STATUS  Status return. Zero on exit if ok. Returns without doing
C            anything if non-zero on entry.
C
C History:
C   ??-???-1990 DFB Created
C   27-Mar-1991 DFB Double precision version using Num. Recipes FFT
C                   Consequently, only powers of 2 are allowed for NFFT.
C-
      INTEGER NMAX
      INTEGER LIM(2),I,IFAIL
      REAL*8 SCALE,WSS

      IF(STATUS.NE.0)RETURN

      LIM(1)=1
      LIM(2)=MIN(NDAT,NFFT)
      CALL GENDWIND(WINDOW,WORK,LIM)
      WSS=0.0
      DO I=1,LIM(2)
         DATA(I)=DATA(I)*WORK(I)
         WSS=WSS+WORK(I)**2
      END DO
      DO I=LIM(2)+1,NFFT
         DATA(I)=0.
      END DO
      IFAIL=-1
      CALL DREALFT(DATA,NFFT,1)
      SCALE=(TSAMP/WSS)
      DO I=0,NFFT/2
        DATA(I+1)=(DATA(2*I+1)**2+DATA(2*I+2)**2)*SCALE
      END DO

 999  END


      SUBROUTINE GENDWIND(WINDOW,WF,XR)
*     ==============================
* Double precision version of GENWIND.
* Generates the window function WF - this version accepts 'PARZEN' 
* (pyramid profile) 'HANNING' (vaguely Gaussian) `GAUSS' (Gaussian 
* with sigma=1/4 width of dataset) or 'TOPHAT'
      IMPLICIT UNDEFINED (A-Z)
      INTEGER XR(2),IX,N,J
      CHARACTER*(*)WINDOW
      REAL*8 WF(*)
      REAL*8 CENTRE,CONST

  1   IF(WINDOW.EQ.'PARZEN')THEN
        N=XR(2)-XR(1)+1
        DO IX=XR(1),XR(2)
           J=IX-XR(1)
           WF(IX)=1.0-ABS(J-0.5*(N-1))/(0.5*(N+1))
        END DO
      ELSE IF(WINDOW.EQ.'HANNING')THEN
        N=XR(2)-XR(1)+1
        DO IX=XR(1),XR(2)
           J=IX-XR(1)
           WF(IX)=0.5*(1.0-COS(2.0*3.141592654*FLOAT(J)/(N-1)))
        END DO
      ELSE IF(WINDOW.EQ.'GAUSS')THEN
        CENTRE=(XR(1)+XR(2))/2D0
        CONST=0.5D0/((XR(2)-XR(1))/4D0)**2
        DO IX=XR(1),XR(2)
           WF(IX)=EXP(-CONST*(DBLE(IX)-CENTRE)**2)
        END DO
      ELSE IF(WINDOW.EQ.'TOPHAT')THEN
        DO IX=XR(1),XR(2)
           WF(IX)=1.0
        END DO
      ELSE
        WRITE(*,1000)WINDOW
 1000   FORMAT(' ** Unknown window function: ',A,' - using PARZEN')
        WINDOW='PARZEN'
        GOTO 1
      ENDIF
      END

C+
      SUBROUTINE LINBIN(X,Y,N,BINFACT,NBIN)

      IMPLICIT NONE
      INTEGER N,BINFACT,NBIN
      REAL*8 X(N),Y(N)
C
C Bins data in-place by BINFACT. If last bin has less than a full
C complement of samples, it is not computed. No of computed bins 
C returned in NBIN
C-
      INTEGER I,J
      REAL SUMX,SUMY

      NBIN=N/BINFACT
      DO I=1,NBIN
        SUMX=0.
        SUMY=0.
        DO J=1,BINFACT
          SUMX=SUMX+X((I-1)*BINFACT+J)
          SUMY=SUMY+Y((I-1)*BINFACT+J)
        END DO
        X(I)=SUMX/BINFACT
        Y(I)=SUMY/BINFACT
      END DO
      END


C+
      SUBROUTINE LOGBIN(X,Y,N,UNITBIN,NBIN)

      IMPLICIT NONE
      INTEGER N,NBIN
      REAL*8 X(N),Y(N),UNITBIN
C
C Bins data in logarithmically-spaced bins. UNITBIN is the x-coord at
C which the bin size is equal to the initial bin size. For samples
C up to this ordinate, no binning is done, and after this the bin size 
C is the nearest integer to X/UNITBIN where X is the x-coord of the 
C leftmost initial bin. If last bin has less than a full
C complement of samples, it IS computed. The number of computed bins 
C returned in NBIN.
C-
      INTEGER INBIN,OUTBIN,J,BINFACT
      REAL*8 SUMX,SUMY

      OUTBIN=0
      INBIN=1
  10  CONTINUE
        IF(INBIN.GT.N)GOTO 90
        BINFACT=MAX(1,NINT(X(INBIN)/UNITBIN))
        BINFACT=MIN(BINFACT,N-INBIN+1)
        SUMX=0.
        SUMY=0.
        DO J=1,BINFACT
          SUMX=SUMX+X(INBIN+J-1)
          SUMY=SUMY+Y(INBIN+J-1)
        END DO
        INBIN=INBIN+BINFACT
        OUTBIN=OUTBIN+1
        X(OUTBIN)=SUMX/BINFACT
        Y(OUTBIN)=SUMY/BINFACT
      GOTO 10
 90   CONTINUE
      NBIN=OUTBIN
      END


      SUBROUTINE FIRSTDIFF(DATA,NDAT)
*     ===============================
* Takes first-differences of a dataset. The last value is set to
* zero so that the number of output items is the same as the number input.
* Thus one piece of information is lost, i.e. the D.C. level of the data.
* The advantage is that the dynamic range of the data is reduced for data
* with steep spectra.
      IMPLICIT NONE
      INTEGER NDAT
      REAL*8 DATA(NDAT)

      INTEGER I

      DO I=1,NDAT-1
         DATA(I)=DATA(I+1)-DATA(I)
      END DO
      DATA(NDAT)=0D0
      END


      SUBROUTINE POSTCOMP(PSPEC,FREQ,NVAL,TSAMP)
*     =======================================
* Post-compensates a power spectrum for the effects of taking 
* first-differences of the sample data. The D.C. value cannot be 
* compensated and so is left as-is.
      IMPLICIT NONE
      INTEGER NVAL
      REAL*8 TSAMP,PSPEC(NVAL),FREQ(NVAL)

      REAL*8 PI
      PARAMETER (PI=6.28318530717958647692D0/2D0)
      INTEGER I

      DO I=1,NVAL
         IF(FREQ(I).NE.0D0)THEN
            PSPEC(I)=PSPEC(I)/(2D0*SIN(PI*TSAMP*FREQ(I)))**2
         ENDIF
      END DO
      END


      SUBROUTINE DECONVOLVE(DATA,MASK,NFFT)
*     =====================================
* Input contains power spectrum of the data and of the data sampling.
* The theory is that DATA=(TRUE DATA)*(MASK) where * denotes convolution.
* The idea is to return some reasonable approximation of the true data.
* Try one is simple deconvolution-by-division, with some suitably-chosen
* apodisation, in this case PARZEN.

      IMPLICIT NONE
      INTEGER NFFT
      REAL*8 DATA(0:NFFT+1),MASK(0:NFFT+1)

      REAL*8 SMALL
      PARAMETER (SMALL=1E-6)
      INTEGER I,IZERO

C - Expand arrays to complex conjugate-symmetric, in place.
      DO I=NFFT/2,0,-1
         DATA(2*I)=DATA(I)
         DATA(2*I+1)=0
         MASK(2*I)=MASK(I)
         MASK(2*I+1)=0
      END DO
C - FT to get ACFs
      CALL DREALFT(DATA,NFFT,-1)
      CALL DREALFT(MASK,NFFT,-1)
C - Find first zero of mask
      DO I=1,NFFT/2
         IF(MASK(I).LT.MASK(0)*SMALL)GOTO 100
      END DO
      I=NFFT/2
 100  IZERO=I
      PRINT*,'First ACF zero at pixel # ',IZERO

C - Divide mask into data, apodise and make conj-symmetric all in one go!
      DO I=IZERO-1,0,-1
         DATA(2*I)=DATA(I)/MASK(I)*(IZERO-I)
         DATA(2*I+1)=0D0
      END DO
C - Pad
      DO I=IZERO,NFFT/2
         DATA(2*I)=0D0
         DATA(2*I+1)=0D0
      END DO
C - Transform back
      CALL DREALFT(DATA,NFFT,-1)
      END