      SUBROUTINE AUTOCP(Q,TQ,NQ,KD,XNU,GG,FF,ST,MODE,FREMAX)
      COMMON /PLOTFL/PLNAME
      CHARACTER PLNAME*12
      DIMENSION Q(NQ),XNU(1)   
      DIMENSION FF(1),GG(1),ST(1),DUM(4)
      DOUBLE PRECISION TQ(NQ),ARG,ARGINC,TWOPI,A,B,DF,DFDF
C
C     THIS ROUTINE ASSUMES THE INPUT DATA IS A FOURIER TRANSFORM
C     WITH THE ORDINATES IN Q AND THE FREQUENCIES, MEASURED IN
C     MICROHERTZ, IN TQ. IT THEN CALCULATES THE AUTOCORRELATION
C     OF THIS FT WITH ITSELF FOR VARIOUS 'LAGS' WHICH ARE
C     EQUALLY SPACE IN PERIOD
C
C     READ THE HIGHERST PERIOD LAG AND THE STEP (FROM ZERO)
C
      READ *,PLAGHI,DPLAG
C
C     CALCULATE THE ZERO LAG VALUE
C
C
C     FIRST CALCULATE THE MEAN
C
      QMEAN=0.0
      DO J=1,NQ
        QMEAN=QMEAN+Q(J)
      END DO
      QMEAN=QMEAN/NQ
C
C     NOW CALCULATE THE MEAN SQUARE DEVIATION
C
      C0=0.0
      DO J=1,NQ
        C0=C0+(Q(J)-QMEAN)*(Q(J)-QMEAN)
      END DO
      C0=C0/NQ

C
C     NOW LOOP OVER THE VARIOUS LAGS AND CALCULATE THE COVARIANCE
C     ACCORDING TO THE RECIPE GIVEN IN 'THE ANALYSIS OF TIME SERIES:
C     AN INTRODUCTION' BY C. CHATFIELD IN THE ASTRONOMY LIBRARY ON
C     SHELF HH.
C
C
C     DNU IS THE SPACING BETWEEN FREQUENCIES IN THE INPUT SPECTRUM IN HZ
C
      DNU = (TQ(NQ)-TQ(1)) / (NQ-1) * 1.0E-6
      NLAG=PLAGHI/DPLAG

      PRINT 5,NLAG,PLAGHI,DPLAG,QMEAN,C0,DNU
    5 FORMAT(12X,'PERIOD AUTOCORRELATION CALCULATION:'/
     *       12X,I5,' LAGS FROM 0 LAG TO ',F7.2,' IN STEPS OF ',F7.2/
     *       12X,'DATA MEAN IS ',1PE12.4,'   C0 IS ',E12.4/
     *       12X,'DNU IS ',1PE12.4)

      DO I=1,NLAG
        PLAG=I*DPLAG
        NTERM=0
        CK=0.0
        DO J=1,NQ
C
C         CALCULATE THE INDEX OF THE TERM IN THE FOURIER TRANSFORM WHICH
C         HAS A FREQUENCY WHOSE PERIOD IS PLAG SECONDS LARGER THAN
C         THE FREQUENCY IN TQ(J)
C
          JJ=NINT(J/(1.0-DNU*J*PLAG))
          IF(J.EQ.1)LAG1=JJ-J
          IF(JJ.LT.1)THEN
            PRINT 10,JJ,J,DNU,PLAG
   10       FORMAT(12X,'*** INDEX ERROR IN AUTOCP ***'/
     *             12X,'JJ= ',I6,'  J= ',I6,'   DNU= ',1PE12.4,
     *             'PLAG= ',0PF7.3)
            STOP
          END IF
          IF(JJ.GT.NQ)GO TO 20
          NTERM=NTERM+1
          CK=CK+(Q(J)-QMEAN)*(Q(JJ)-QMEAN)
          LAGLAS=JJ-J
        END DO
   20   CK=CK/NTERM
        XNU(I)=PLAG
        FF(I)=CK/C0
        KD=I
        IF(MOD(I-1,20).EQ.0)PRINT 25,I,PLAG,LAG1,LAGLAS,NTERM
   25   FORMAT(12X,I5,'  LAG ',F7.1,'  JJ-J: 1ST ',I6,
     *         '  LAST ',I6,'  NTERM ',I6)
      END DO
   50 CONTINUE

C
C     INITIALISE THE PLOTTING
C
c     CALL SETCC
c     CALL PLFILE(PLNAME)
c     CALL PLOTS(IERR,IDUMMY,28)
c     IF(IERR.NE.0) STOP 'ERROR ON OPENING SPOOLTEK CALL'
      XLEN=6.0 
      YLEN=5.0 
      XHINU=0.0   
      XLONU=0.0
      DO J=1,KD
        IF(XHINU.LT.XNU(J))XHINU=XNU(J)   
        IF(XLONU.GT.XNU(J))XLONU=XNU(J)   
      END DO
      DELTAV=(XHINU-XLONU)/XLEN     
      XNU(KD+1)=XLONU   
      XNU(KD+2)=DELTAV  

c     CALL FACTOR(2.54) 
c     CALL FRAME(XLEN+3.5,YLEN+2.5)
c     CALL PLOT(2.,1.0,-3)   
      AMAXFF=0.0
      AMINFF=0.0
      DO J=1,KD    
        IF(FF(J).GT.AMAXFF)AMAXFF=FF(J) 
        IF(FF(J).LT.AMINFF)AMINFF=FF(J) 
      END DO
      DUM(1)=AMINFF
      DUM(2)=AMAXFF   
c     CALL SCALE(DUM,YLEN,2,1)    
      FF(KD+1)=DUM(3) 
      FF(KD+2)=DUM(4) 
      SCALEF=FF(KD+2)   
      PRINT 75,AMAXFF,AMINFF,XLONU,XHINU   
   75 FORMAT(12X,'PLOT SCALED TO A MAXIMUM OF ',F13.5/
     *       12X,'             AND MINIMUM OF ',F13.5/
     *       12X,'PLOTTING : LOW PERIOD LAG  ',1PE15.7/
     *       12X,'           HI  PERIOD LAG  ',E15.7)  
c     CALL FAXIS(0.,0.,'PERIOD LAG',-10,XLEN,XLONU,DELTAV)    
c     CALL AAXIS(0.,0.,'AMPLITUDE',9,YLEN,FF(KD+1),SCALEF)  
c     CALL PLOT(0.0,YLEN,3)
c     CALL PLOT(XLEN,YLEN,2)
c     CALL AAXIS(XLEN,0.,' ',-1,YLEN,FF(KD+1),SCALEF) 
      YYYMAX=0.0
      DO J=1,KD    
        XXX=(XNU(J)-XLONU)/DELTAV   
        YYY=(FF(J)-FF(KD+1))/SCALEF
        YYY=MIN(YYY,YLEN)   
        YYYMAX=MAX(YYY,YYYMAX)
c       IF(J.EQ.1)CALL PLOT(XXX,YYY,3)    
c       CALL PLOT(XXX,YYY,2)  
      END DO    
c     CALL PLOT(0.0,0.0,999) 

      RETURN
      END   
