      SUBROUTINE EQUALP(Q,TQ,NQ,KD,XNU,FF,FREMAX,TITLE,HELJD)
      COMMON /PLOTFL/PLNAME
      CHARACTER PLNAME*12
      CHARACTER TITLE*80,YTITLE*20,XTITLE*20,FNAME*30
      DIMENSION Q(NQ),XNU(1)   
      DIMENSION FF(1),GG(1),ST(1),DUM(4),OMC(2000),CYC(2000)
      DIMENSION OMCBES(2000),CYCBES(2000),PERIOD(2000)
      DOUBLE PRECISION TQ(NQ),ARG,ARGINC,TWOPI,A,B,DF,DFDF,HELJD(*)
C
C     THIS ROUTINE ASSUMES THE INPUT DATA ARE A SET OF NQ PERIODS  
C     (IN TQ IN SECONDS) WITH CORRESPONDING POWER (IN MICRO-OOMPS) 
C     IN ARRAY Q. IT THEN CALCULATES THE INVERSE VARIANCE
C     OF THIS SET OF PERIODS WITH RESPECT TO A DIRAC COMB OF
C     CONSTANT PERIOD SPACING AND STORES THE RESULTS FOR THE KD
C     DIFFERENT PERIODS SPACINGS IN ARRAYS XNU AND FF. XNU CONTAINS
C     THE INVERSE OF THE PERIOD SPACING AND FF CONTAINS THE INVERSE
C     VARIANCE FOR THAT PERIOD SPACING.
C
C     READ THE SHORTEST AND LONGEST PERIOD SPACING SEARCHED
C
      READ *,PLOW,PHIGH
 
      MODE=0
      IF(PLOW.LT.0.0)THEN
        PLOW=-PLOW
        MODE=1
      END IF

C
C     CALCULATE THE CORRESPONDING FREQUENCIES AND STEP LENGTH
C
      FLO=1./PHIGH
      FHI=1.0/PLOW
      DF=1.0/20.0/(TQ(NQ)-TQ(1))
      N=(FHI-FLO)/DF
C
C     FIRST CALCULATE THE HIGHEST PEAK
C
      QMAX=0.0
      DO J=1,NQ
        IF(QMAX.LT.Q(J))THEN
          QMAX=Q(J)
C         THIS REFERENCE PERIOD IS TO BE USED TO CALCULATE THE
C         'CYCLE' NUMBERS IN EACH TRIAL EPHEMERIS
          PREF=TQ(J)
        END IF
      END DO
C
C     NOW LOOP OVER THE VARIOUS PERIOD SPACINGS AND CALCULATE THE
C     INVERSE VARIANCE OF THE FIT OF OF THE GIVEN PERIOD SPACING TO 
C     THE OBSERVED PERIODS.
C
      PRINT 5,N,PLOW,PHIGH
    5 FORMAT(12X,'PERIOD SPACING INVERSE VARIANCE CALCULATION:'/
     *       12X,I5,' SPACINGS FROM ',F7.2,' TO ',F7.2/)

      AMAX=0.0
      DO I=1,N
        F=FLO+(I-1)*DF
        PSPACE=1.0/F
C
C       NOW CALCULATE O-C FOR EACH PERIOD WHERE O IS THE OBSERVED
C       PERIOD VALUE AND C = PBASE + E*PSPACE
C
        PBASE = 0.0
        DO J=1,NQ
C
C         FIRST CALCULATE THE 'CYCLE' NUMBERS FROM O, THE OBSERVED 
C         PERIODS USING (O-PREF)/PSPACE. NOTE THAT THE CYCLE NUMBERS,
C         E, HAVE TO BE INTEGERS.
C
          E = NINT( (TQ(J)-PREF)/PSPACE )
          CYC(J) = E
C
C         NOW CALCULATE PBASE TO OPTIMISE THE FIT BETWEEN THE OBSERVED
C         PERIODS AND THE ONES FROM THE 'LINEAR EPHEMERIS' 
C         PBASE + E * PSPACE. PBASE IS GIVEN BY THE FORMULA
C             1/N * SUM( O - PSPACE*E )
C
          PBASE = PBASE + (TQ(J)-PSPACE*E)/NQ
        END DO
C
C       NOW CALCULATE O-C
C
        DO J=1,NQ
          OMC(J) = TQ(J) - (PBASE + CYC(J)*PSPACE)
C         CONVERT TO CYCLES
          OMC(J) = OMC(J)/PSPACE
          IF(OMC(J).GT.0.5.OR.OMC(J).LT.-0.5)OMC(J)=OMC(J)-NINT(OMC(J))
        END DO
C
C       NOW CALCULATE THE MEAN AND INVERSE VARIANCE OF THE O-C'S
C
        FF(I)=0.0
        AVEOMC=0.0
        DO J=1,NQ
          AVEOMC=AVEOMC+OMC(J)
        END DO
        AVEOMC=AVEOMC/N
        DO J=1,NQ
          FF(I)=FF(I)+(OMC(J)-AVEOMC)*(OMC(J)-AVEOMC)
        END DO
        FF(I)=(NQ-2)/FF(I)
        IF(AMAX.LT.FF(I))THEN
          AMAX=FF(I)
          FREMAX=F
          BASMAX=PBASE
          AVERES=0.0
          DO J=1,NQ
            PERIOD(J)=TQ(J)
            OMCBES(J)=OMC(J)
            CYCBES(J)=CYC(J)
            AVERES=AVERES+ABS(OMC(J))
          END DO
          AVERES=AVERES/NQ
        END IF
        XNU(I)=F
        KD=I
      END DO
      PRINT 25,1/FREMAX,BASMAX,AMAX,AVERES,
     *         (TQ(J),CYCBES(J),OMCBES(J),J=1,NQ)
   25 FORMAT(12X,'BEST PERIOD SPACING IS            ',F7.2/
     *       12X,'BEST BASE PERIOD IS               ',F7.2/
     *       12X,'CORRESPONDING INVERSE VARIANCE IS ',F7.2/
     *       12X,'AVERAGE ABSOLUTE RESIDUAL IS      ',F7.3/
     *       12X,'PERIOD      CYCLE      O-C'/
     *       (12X,F6.1,4X,F7.2,F9.3))

      IF(MODE.EQ.1)THEN 
C     
C       WRITE THE INVERSE VARIANCE SPECTRUM   OUT TO A FILE 
C     
        READ(*,'(A)')FNAME    
        CLOSE(11) 
        OPEN(11,FILE=FNAME,STATUS='UNKNOWN')
        PRINT 34,FNAME
   34   FORMAT(12X,'*** INVERSE VARIANCE SPECTRUM   BEING DUMPED'/ 
     *         12X,'*** TO FILE ',A25/)   
        WRITE(11,'(A78)')TITLE
        WRITE(11,36)NQ,HELJD(1),KD,XNU(1),XNU(KD),DF
   36   FORMAT('NPTS : ',I6,4X,'HJD : ',F18.7,4X,'NFREQ :',I6/   
     *     'FLOW : ',1PE15.7,4X,'FHIGH : ',1PE15.7,4X,'DF : ',1PE15.7)
        DO 50 K=1,KD     
          WRITE(11,40)K,XNU(K),FF(K)
   40     FORMAT(I10,1P3E20.8)
   50   CONTINUE  
        CLOSE(11) 
      END IF

C     
C     ROUTINE TO PLOT ARRAY TQ AGAINST O-C
C     
      XL=8.0
      YL=4.0
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'
c     CALL FACTOR(2.54) 
c     CALL FRAME(XL+4.0,YL+4.0)          
c     CALL PLOT(1.5,1.5,-3)   
c     CALL SYMBOL(0.,YL+1.5,0.2,TITLE,0.,60)  
      XTITLE='Period (S)'     
c     CALL SCALE(PERIOD,XL,NQ,1)     
      TT=PERIOD(NQ+1)    
      TTT=PERIOD(NQ+2)   
      OMCMAX=-0.5     
      OMCMIN=0.5
      OMCSCL=-1.0/YL  
      YY=OMCMIN 
      YYY=OMCSCL
      YTITLE=' O-C  (Cycles) '    
c     CALL AXIS(0.0,0.0,XTITLE,-10,XL,0.,TT,TTT) 
c     CALL AXIS(0.0,YL,XTITLE,10,XL,0.,TT,TTT)   
c     CALL AXIS(0.,0.,YTITLE,16,YL,90.0,YY,YYY) 
c     CALL AXIS(XL,0.,YTITLE,-16,YL,90.0,YY,YYY)
      DO 100 J=1,NQ  
        X=(PERIOD(J)-TT)/TTT 
        Y=(OMCBES(J)-OMCMIN)/OMCSCL
        IF(Y.LT.0.0)Y=Y+1.0/ABS(OMCSCL) 
        IF(Y.GT.YL)Y=Y-1.0/ABS(OMCSCL)  
c       CALL SYMBOL(X,Y,0.05,1,0.0,-1)    
  100 CONTINUE    
c     CALL PLOT(0.0,0.0,999) 
      RETURN
      END   
