      SUBROUTINE PDOTP(EPOCH,EPER,ESEC,ETHIRD,CYCLE,OMINC,N,TITLE)     
C     
C     ROUTINE TO SCAN THROUGH A RANGE OF P AND PDOT AND FIND
C     THE BEST VALUES OF EPHOCH + E*EPER * E**2*ESECOND THAT
C     FITS A SET OF TIMINGS
C     
      DOUBLE PRECISION EPOCH,EPER,ESEC,ETHIRD
      REAL*8 QPOCH,QPER,QSEC
      REAL*8 PLOW,PHIGH,DP
      REAL*8 SECLO,SECHI,DSEC
      REAL*8 ASINI,T0,PORB,DELAY(300),ORBPHS(300)
      REAL*8 TIMING(300),CYC(300),OMC(300)
      REAL*8 E,C(300),PREDIC(300)
      REAL*8 TWOPI,ARG,DIFF
      DIMENSION OMINC(1),CYCLE(1),VARIN(105,105)
      DIMENSION OOMICC(105,4),CYCCYC(105,4)
      DIMENSION X(105),Y(105)
      DOUBLE PRECISION PERIOD(105),PDOT(105)
      CHARACTER FNAME*25,TITLE*80,FORMAT*80,HEADER*80,XLABEL*80
      FORMAT=' '
      IDEBUG=0
      TWOPI=2.0D0*3.141592653D0     
      QPOCH=0.0D0 
      QPER=0.0D0  
      QSEC=0.0D0  
      VARBES=-1.0E30
      PRINT 4
    4 FORMAT(12X,'ENTER THE FILE NAME WITH THE TIMINGS')
      READ 5,FNAME
    5 FORMAT(A)
      IF(FNAME.EQ.'FORMAT'.OR.FNAME.EQ.'format')THEN
C
C       READ IN THE DATA FILE INPUT FORMAT
C
        READ 5,FORMAT
C
C       NOW READ THE DATA FILE NAME (AGAIN!)
C
        READ 5,FNAME
      END IF
C
C     NOW READ THE SCAN RANGES IN P AND PDOT
C
      PRINT 6
    6 FORMAT(12X,'ENTER THE FIXED EPOCH OF THE EPHEMERIS')
      READ *,QPOCH
      PRINT 7
    7 FORMAT(12X,'ENTER THE PERIOD SCAN RANGE: PLOW PHIGH DP')
      READ *,PLOW,PHIGH,DP
      PRINT 8
    8 FORMAT(12X,'ENTER THE SECOND TERM SCAN RANGE: SECLO SECHI DSEC')
      READ *,SECLO,SECHI,DSEC
      PRINT 9
    9 FORMAT(12X,'ENTER THE ORBITAL TIME DELAY: ASINI T0 PORB')
      READ *,ASINI,T0,PORB
      PRINT 10,FNAME,ASINI,T0,PORB,QPOCH,PLOW,PHIGH,DP,SECLO,SECHI,DSEC
   10 FORMAT(12X,'DATA TO BE READ FROM FILE ',A/
     *       12X,'ORBITAL TIME DELAY PARAMETERS:'/
     *       12X,'               A SIN I: ',0PF16.6/
     *       12X,'    TIME OF SUP. CONJ.: ',0PF16.6/
     *       12X,'             ORB. PER.: ',0PF16.7/
     *       12X,'EPOCH FIXED AT          ',0PF16.6/
     *       12X,'SCANNING IN PERIOD FROM ',1PE18.10/
     *       12X,'                     TO ',E18.10,' STEP ',E18.10/
     *       12X,'     AND IN SECOND FROM ',1PE18.10/ 
     *       12X,'                     TO ',E18.10,' STEP ',E18.10)
      CLOSE(11)   
      OPEN(11,FILE=FNAME,STATUS='OLD')
      READ(11,*)N
      IF(FORMAT.EQ.' ')THEN
        READ(11,*)(TIMING(J),J=1,N)
      ELSE
        READ(11,FORMAT)(TIMING(J),J=1,N)
      END IF
      REWIND 11   
      CLOSE(11)   
      IF(ABS(QPOCH-999.0D0).LT.0.001)THEN
C
C       SET THE EPOCH TO THE TIMING WHICH IS 
C       NEAREST TO THE AVERAGE OF THE FIRST
C       AND LAST TIMINGS
C
        QPOCH=(TIMING(1)+TIMING(N))/2.0
        DIFF=1.0E30
        DO J=1,N
          IF(ABS(QPOCH-TIMING(J)).LT.DIFF)THEN
            DIFF=ABS(QPOCH-TIMING(J))
            NEAR=J
          END IF
        END DO
        QPOCH=TIMING(NEAR)
      END IF
      EPOCH=QPOCH
      PRINT 15,N
   15 FORMAT(12X,'NO OF DATA POINTS READ: ',I5)
C
C     CALCULATE THE TIMING DELAYS DUE TO THE ORBIT
C
      DO 20 J=1,N
        ORBPHS(J)=(TIMING(J)-T0)/PORB
        IORB=ORBPHS(J)
        ORBPHS(J)=ORBPHS(J)-IORB
        DELAY(J)=ASINI*COS(ORBPHS(J)*TWOPI)
C       SUBTRACT THE DELAY FROM THE OBSERVED TIMING
        TIMING(J)=TIMING(J)-DELAY(J)
   20 CONTINUE
C      OPEN(13,FILE='VAR',STATUS='NEW')
C     
C     DO THE CALCULATIONS
C
C     LOOP OVER PERIOD
C
      I1=0
      QPER=PLOW-DP
   30 QPER=QPER+DP
      IF(QPER.GT.PHIGH)GO TO 150
      I1=I1+1
C
C     LOOP OVER THE SECOND TERM
C
      I2=0
      QSEC=SECLO-DSEC
   40 QSEC=QSEC+DSEC
      IF(QSEC.GT.SECHI)GO TO 30
      I2=I2+1
      IF(I2.EQ.1)PERIOD(I1)=QPER
      IF(I1.EQ.1)PDOT(I2)=2.0*QSEC/QPER
C
C     ASSIGN NEAREST CYCLES TO ALL THE TIMINGS BASED ON CURRENT
C     LINEAR EPHEMERIS AND CORRECT THEM UNTIL CONSISTENT
C
      DO 60 J=1,N
        CYC(J)=(TIMING(J)-QPOCH)/QPER
        ICYC=CYC(J)+0.5
        CYC(J)=ICYC
        LOOP=0
C
C       NOW CALCULATE O-C FOR THE ASSIGNED CYCLES
C
   55   LOOP=LOOP+1
        E=CYC(J)
        C(J)=QPOCH + E*QPER + E*E*QSEC 
        OMC(J)=TIMING(J)-C(J)
C       CONVERT TO CYCLES AND CHECK THAT O-C IN THE RANGE -0.5 TO 0.5
        OMC(J)=OMC(J)/QPER
        IF(OMC(J).GT.0.5.OR.OMC(J).LT.-0.5.AND.LOOP.LT.5)THEN
C         PRINT 56,J,CYC(J),TIMING(J),C(J),OMC(J)
   56     FORMAT(' J=',I4,'  E= ',F12.1,' O=',F12.6,' C= ',F12.6,
     *           ' O-C= ',F12.6)
          CYC(J)=CYC(J)+NINT(OMC(J))
          GO TO 55
        END IF
   60 CONTINUE

C     
C     CALCULATE THE MEAN AND INVERSE VARIANCE OF THE O-C'S
C     
      VARIN(I1,I2)=0.0     
      AVEOMC=0.0
      DO 70 J=1,N
        AVEOMC=AVEOMC+OMC(J)
   70 CONTINUE
      AVEOMC=AVEOMC/N
      DO 80 J=1,N
        VARIN(I1,I2)=VARIN(I1,I2)+(OMC(J)-AVEOMC)*(OMC(J)-AVEOMC)
   80 CONTINUE
      VARIN(I1,I2)=(N-2)/VARIN(I1,I2)
C      WRITE(13,90)(QPER-PLOW)/1.0D-10,2.0*QSEC/QPER/1.0D-12,
C     *            VARIN(I1,I2)
   90 FORMAT(3F8.1)
      IF(VARIN(I1,I2).GT.VARBES)THEN
        VARBES=VARIN(I1,I2)
        EPER=QPER
        ESEC=QSEC
        I1BEST=I1
        I2BEST=I2
        DO 100 J=1,N
          OMINC(J)=OMC(J)
          CYCLE(J)=CYC(J)
          PREDIC(J)=C(J)
  100   CONTINUE    
      END IF
C     PRINT 110,QPER,QSEC,AVEOMC,VARIN(I1,I2)
  110 FORMAT(12X,'PERIOD ',1PE14.8,2X,'SECOND ',E10.3,2X,
     *       '<O-C> ',0PF5.2,' 1/VARIANCE ',F7.1)
      GO TO 40

  150 CONTINUE
      XL=8.0
      IF(TITLE.EQ.' ')THEN
        WRITE(HEADER,'(A)')'O-C diagram for best fit'
        TITLE=HEADER
      END IF
      CALL PLOTPH(XL,DUM1,DUM2,DUM3,CYCLE,OMINC,N,TITLE,-2)
      AMAX=-1.0E20
      AMIN=1.0E20
      DO 220 I=1,I1
        DO 210 J=1,I2
          IF(VARIN(I,J).GT.AMAX)AMAX=VARIN(I,J)
          IF(VARIN(I,J).LT.AMIN)AMIN=VARIN(I,J)
  210   CONTINUE
  220 CONTINUE
      IF(I1.EQ.1.AND.I2.GT.1)THEN
C
C       THIS MEANS A SCAN IN PDOT HAS BEEN DONE AT FIXED PERIOD
C
        IF(TITLE.EQ.' ')THEN
          WRITE(HEADER,'(A15,1PE14.8)')'1/Var(O-C) at P= ',QPER
          TITLE=HEADER
        END IF
        DO 225 J=1,I2
          X(J)=PDOT(J)*1.0E15
          Y(J)=VARIN(1,J)
  225   CONTINUE
        WRITE(XLABEL,226)
  226   FORMAT('Pdot x 10**-15')
        CALL PLTVAR(X,Y,I2,TITLE,XLABEL)
      ELSE IF(I1.GT.1.AND.I2.EQ.1)THEN
C
C       THIS MEANS A SCAN IN PERIOD HAS BEEN DONE AT FIXED PDOT
C
        IF(TITLE.EQ.' ')THEN
          WRITE(HEADER,'(A18,1PE10.3)')'1/Var(O-C) at Pdot= ',PDOT(1)
          TITLE=HEADER
        END IF
        DO 230 J=1,I1
          X(J)=(PERIOD(J)-1.91416425D-4)*1.0E11
          Y(J)=VARIN(J,1)
  230   CONTINUE
        WRITE(XLABEL,231)1.91416425D-4
  231   FORMAT('Period x 10**-11 + ',1PE14.8)       
        CALL PLTVAR(X,Y,I1,TITLE,XLABEL)
      ELSE IF(I1.GT.1.AND.I2.GT.1)THEN
        OPEN(12,FILE='VARIMG',STATUS='NEW')
        WRITE(12,*)I2,I1
        WRITE(12,*)10.0*AMIN,10.0*AMAX,10.0*AMIN,10.0*AMAX
        DO 250 I=1,I1
          WRITE(12,240)(VARIN(I,J),J=1,I2)
  240     FORMAT(1P10F8.1)
  250   CONTINUE
        CLOSE(12)
      END IF
      PRINT 260,VARBES,EPER,ESEC,2.0*ESEC/EPER
  260 FORMAT(/12X,'BEST 1/VARIANCE ',F8.1,5X,'PERIOD ',1PE17.10/
     *        12X,'SECOND        ',E10.3, 5X,'PDOT ',E10.3)

      PRINT 300    
  300 FORMAT(12X,'O-C TABLE FOR BEST VARIANCE:'//
     *       12X,'     CYCLE      PH(ORB)     OBS.  TIME      PREDICTED'  
     *       ,3X,'       O-C'/ 
     *       12X,'----------    ---------     ----------     ----------'  
     *       ,3X,'    ------') 
      PRINT 310,(CYCLE(J),ORBPHS(J),TIMING(J),PREDIC(J),OMINC(J),J=1,N)
  310 FORMAT(12X,F10.0,F13.3,F15.5,F15.5,F13.3)  
      IF(I1.EQ.1.AND.I2.EQ.1)RETURN
C
C     NOW FIND O-C FOR THE NEXT FOUR BEST SOLUTIONS PEAKS
C
      PRINT *,' '
      LOOP=0
  320 LOOP=LOOP+1
C
C     FLATTEN THE PREVIOUS BEST PEAK
C     FIRST: DETERMINEING THE RANGE TO BE FLATTENED IN EACH DIMENSION
C
      DO 330 J=I1BEST+1,I1
        IF(VARIN(J,I2BEST).GT.VARIN(J-1,I2BEST))GO TO 340
  330 CONTINUE
  340 I1HI=J-1
      DO 350 J=I1BEST-1,1,-1
        IF(VARIN(J,I2BEST).GT.VARIN(J+1,I2BEST))GO TO 360
  350 CONTINUE
  360 I1LO=J+1
      DO 370 J=I2BEST+1,I2
        IF(VARIN(I1BEST,J).GT.VARIN(I1BEST,J-1))GO TO 380
  370 CONTINUE
  380 I2HI=J-1
      DO 390 J=I2BEST-1,1,-1
        IF(VARIN(I1BEST,J).GT.VARIN(I1BEST,J+1))GO TO 400
  390 CONTINUE
  400 I2LO=J+1
      PRINT 420,LOOP,I1LO,I1HI,I2LO,I2HI
  420 FORMAT(/12X,'FLATTENING PEAK NUMBER ',I2/
     *        12X,'FROM ',I3,' TO ',I3,' IN I1 AND'/
     *        12X,'FROM ',I3,' TO ',I3,' IN I2')
      DO 440 I=I1LO,I1HI
        DO 440 J=I2LO,I2HI
  440     VARIN(I,J)=0.0
C
C     NOW FIND THE NEXT HIGHEST MAXIMUM
C
      VARBES=AMIN
      DO 460 I=1,I1
        DO 460 J=1,I2
          IF(VARIN(I,J).GT.VARBES)THEN
            VARBES=VARIN(I,J)
            QPER=PLOW+(I-1)*DP
            QSEC=SECLO+(J-1)*DSEC
            I1BEST=I
            I2BEST=J
          END IF
  460 CONTINUE
      PRINT 470,LOOP+1,VARBES,QPER,QSEC,2.0*QSEC/QPER
  470 FORMAT(/12X,I2,' BEST VARIANCE ',F8.1,5X,'PERIOD ',1PE17.10/
     *        12X,'SECOND ',E10.3,'  PDOT ',E10.3/)
C
C     NOW CALCULATE CYCLES AND O-C FOR THIS BEST VALUE
C
      DO 500 J=1,N
        CYC(J)=(TIMING(J)-QPOCH)/QPER
        ICYC=CYC(J)+0.5
        CYC(J)=ICYC
        LP=0

  480   LP=LP+1
        E=CYC(J)
        C(J)=QPOCH + E*QPER + E*E*QSEC 
        OMC(J)=TIMING(J)-C(J)
C       CONVERT TO CYCLES AND CHECK THAT O-C IN THE RANGE -0.5 TO 0.5
        OMC(J)=OMC(J)/QPER
        IF(OMC(J).GT.0.5.OR.OMC(J).LT.-0.5.AND.LP.LT.5)THEN
C         PRINT 56,J,CYC(J),TIMING(J),C(J),OMC(J)
          CYC(J)=CYC(J)+NINT(OMC(J))
          GO TO 480
        END IF
        CYCCYC(J,LOOP)=CYC(J)
        OOMICC(J,LOOP)=OMC(J)
        IF(LOOP.LE.2)
     *  PRINT 310,CYC(J),ORBPHS(J),TIMING(J),C(J),OMC(J)
  500 CONTINUE
      XL=8.0
      WRITE(HEADER,'(A16,I1,A9)')'O-C diagram for ',LOOP+1,' best fit'
      TITLE=HEADER
      IF(I1.GT.1.AND.I2.GT.1..AND.LOOP.LE.2)
     *CALL PLOTPH(XL,DUM1,DUM2,DUM3,CYCCYC(1,LOOP),OOMICC(1,LOOP),N,
     *            TITLE,-2)
      IF(LOOP.LE.3)GO TO 320

      RETURN
      END   
      SUBROUTINE PLTVAR(X,FF,KD,TITLE,XLABEL)
      COMMON /PLOTFL/PLNAME
      CHARACTER PLNAME*12
      DIMENSION X(1),FF(1),DUM(4)   
      CHARACTER TITLE*80,XLABEL*80
      INTEGER TRMLEN
      EXTERNAL TRMLEN

      LENTIT=TRMLEN(TITLE)
      LXLAB=0
      DO 6 J=80,1,-1
        IF(XLABEL(J:J).NE.' '.AND.LXLAB.EQ.0)LXLAB=J
    6 CONTINUE
C
C     INITIALISE THE PLOTTING
C
c     CALL SETCC
c     CALL PLFILE(PLNAME)
c     CALL PLOTS(IERR,IDUMMY,28)
      IF(IERR.NE.0) STOP 'ERROR ON OPENING SPOOLTEK CALL'
      XLEN=6.0 
      YLEN=5.0
      AMAXFF=0.0
      AMINFF=0.0
      XHI=-1.0E9   
      XLO=1.0E9 
      DO 10 J=1,KD
        IF(XHI.LT.X(J))XHI=X(J)   
        IF(XLO.GT.X(J))XLO=X(J)   
   10 CONTINUE    
      DELTAV=(XHI-XLO)/XLEN     
      X(KD+1)=XLO   
      X(KD+2)=DELTAV  
c     CALL FACTOR(2.54) 
c     CALL FRAME(XLEN+3.5,YLEN+2.5)
c     CALL PLOT(2.,1.0,-3)   
      DO 20 J=1,KD    
        IF(FF(J).GT.AMAXFF)AMAXFF=FF(J) 
        IF(FF(J).LT.AMINFF)AMINFF=FF(J)
   20 CONTINUE  
      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)   
c     IF(TITLE.NE.' ')CALL SYMBOL(0.0,YLEN+0.84,0.14,TITLE,0.0,LENTIT)  
c     CALL FAXIS(0.,0.,XLABEL,-LXLAB,XLEN,XLO,DELTAV)    
c     CALL AAXIS(0.,0.,'1/Var(O-C)',10,YLEN,AMINFF,SCALEF)  
c     CALL PLOT(0.,YLEN,3)
c     CALL PLOT(XLEN,YLEN,2)
c     CALL AAXIS(XLEN,0.,' ',-1,YLEN,AMINFF,SCALEF) 
      DO   30 J=1,KD    
        XXX=(X(J)-XLO)/DELTAV   
        YYY=(FF(J)-AMINFF)/SCALEF
c       IF(J.EQ.1)CALL PLOT(XXX,YYY,3)    
c       CALL PLOT(XXX,YYY,2)  
   30 CONTINUE    
c     CALL PLOT(0.0,0.0,999) 
      RETURN
      END   
