C     
C     THIS MAIN PROGRAM SIMPLY SETS UP THE ARRAY SIZES AND THEN   
C     CALLS SUBROUTINE MAIN - THE ACTUAL MAIN PROGRAM.
C     
      COMMON /PGPL/PGPLDEV,IPGPL
      CHARACTER PGPLDEV*80
      LOGICAL STATUS

      PARAMETER (MAXFRQ=250000)
      PARAMETER (MAXPT=250000)

      CHARACTER IROUTE*6,INROUT*6
      DATA JUMP/0/
C
      READ 5,INROUT     
    5 FORMAT(A6)  
      CALL UPCASE(INROUT)
      IF(INROUT(1:5).NE.'SETUP')THEN
        JUMP=1    
c       MAXFRQ=1000     
c       MAXPT=5000
      ELSE  
        JUMP=0    
c       READ *,MAXFRQ,MAXPT   
c       IF(MAXFRQ.EQ.0)MAXFRQ=1000  
c       IF(MAXPT.EQ.0)MAXPT=5000    
      END IF
      PRINT 10,MAXFRQ,MAXPT   
   10 FORMAT('0',50X,'PROGRAM EAGLE'/51X,13('-')//   
     *       10X,'PROGRAM SET UP FOR ',I6,' FREQUENCIES AND ',I6, 
     *       ' DATA POINTS'/) 

      CALL MAIN (MAXFRQ, MAXPT, JUMP, INROUT)  
      IF(IPGPL.EQ.1)CALL PGEND
      STOP  
      END   

       SUBROUTINE UPCASE(STRG)
       CHARACTER STRG*(*),ALOW(26)*1,AUP(26)*1
       DATA ALOW/'a','b','c','d','e','f','g','h','i','j','k','l',
     *   'm','n','o','p','q','r','s','t','u','v','w','x','y','z'/
       DATA AUP/'A','B','C','D','E','F','G','H','I','J','K','L',
     *  'M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/
        LSTRG=LEN(STRG)
        DO J=1,LSTRG
           DO K=1,26
           IF(STRG(J:J).EQ.ALOW(K)) STRG(J:J)=AUP(K)
           ENDDO
        ENDDO
        RETURN
        END


      SUBROUTINE MAIN(MAXFRQ,MAXPT,JUMP,IROUTE) 
C     
C     THIS IS THE MAIN PROGRAM IN THE FORM OF A SUBROUTINE.  THE REAL   
C     MAIN PROGRAM CONTAINS CODE FOR DYNAMICALLY DIMENSIONING ARRAY     
C     Q,TQ,FF,GG,XNU.  IF THIS IS NO LONGER NEEDED, SIMPLY TAKE OUT     
C     THE SUBROUTINE STATEMENT, INSERT A SUITABLE DIMENSION STATEMENT   
C     FOR Q,TQ,GG,FF,XNU AND ELIMINATE THE STATEMENT : IF(JUMP.EQ.1)    
C     GO TO 35. FINALLY, MOVE THE PRINT 10 STATEMENT OUT OF SETUP
C     AND DELETE SETUP
C     
      DIMENSION Q(MAXPT),FF(MAXFRQ),GG(MAXFRQ),XNU(MAXFRQ)    
      DOUBLE PRECISION TQ(MAXPT),DF
      DIMENSION SIGMAA(40),SIGMAF(40),SIGMAP(40),SGFDOT(40),ST(MAXFRQ)  
      DIMENSION NX(100),FREQ(40),AMP(40),PHI(40),FDOT(40)     
      DIMENSION SAMTIM(100),LIM(2,100),IQ(10),PHASES(1000),TIME(1000)   
      DIMENSION SIGPH(1000),AMPLTD(1000),SIGAMP(1000),AFREQ(2,1000)     
      DIMENSION IND(1000),E(1000),OMINC(1000),FILPER(26),FILRES(26)     
      DOUBLE PRECISION O(1000),C(1000)
      DOUBLE PRECISION HELJD(100),TMAX,SIGTMX,HJDMAX,HJD1NU
      DOUBLE PRECISION TTT,FRE,AMPLIT,PHASE,OMEG,DARG,DTWOPI,TI
      DOUBLE PRECISION EPOCH,EPER,ESEC,ETHIRD   
      CHARACTER*6 NX,IROUTE,INROUT,TITLE*80,FNAME*25     
      COMMON /PGPL/PGPLDEV,IPGPL
      CHARACTER PGPLDEV*80,PGFILE*80,PGTYPE*80
      INTEGER PGBEG
      COMMON /PLOTFL/PLNAME
      COMMON /ARROWS/NARROW,FARROW(20),HARROW(20)
      CHARACTER PLNAME*12
      common /VER/ iver
      common /ecl/modee,param(30),nparam,ai
      INTEGER TRMLEN
      EXTERNAL TRMLEN
      LOGICAL STATUS
c     POINTER (PQ,Q)
c     POINTER (PTQ,TQ)
c     POINTER (PXNU,XNU)
c     POINTER (PFF,FF)
c     POINTER (PGG,GG)
c     POINTER (PST,ST)
      PLNAME='PLFILE'
      PGPLDEV='/VPS'
      PGTYPE=PGPLDEV
      PGFILE='pgplot.ps'
      IPGPL=0

C
C  Dynamic memory allocation
C
c     PQ=MALLOC(SIZEOF(Q(1))*MAXPT)
c     PTQ=MALLOC(SIZEOF(TQ(1))*MAXPT)
c     PXNU=MALLOC(SIZEOF(XNU(1))*MAXFRQ) 
c     PFF=MALLOC(SIZEOF(FF(1))*MAXFRQ)
c     PGG=MALLOC(SIZEOF(GG(1))*MAXFRQ)
c     PST=MALLOC(SIZEOF(ST(1))*MAXFRQ)
C
      iver=0
      NX(1)='DISKRQ'    
      NX(2)='TIMES '    
      NX(3)='LSTSQ '    
      NX(4)='TITLE '    
      NX(5)='BRIDGE'    
      NX(6)='QUIST '    
      NX(7)='STOP  '    
      NX(8)='STPLOT'    
      NX(9)='PREPQ '    
      NX(10)='TRANFQ'   
      NX(11)='PLOTP '   
      NX(12)='PLOTW '   
      NX(13)='PRINT '   
      NX(14)='PREWHI'   
      NX(15)='ANALYS'   
      NX(16)='FOLDPH'   
      NX(17)='FOLDAM'   
      NX(18)='FILLTA'   
      NX(19)='DUMPTA'   
      NX(20)='PGPLOT'   
      NX(21)='PRINTQ'   
      NX(22)='PRINTT'   
      NX(23)='      '   
      NX(24)='TRANFA'   
      NX(25)='PREWHA'   
      NX(26)='PUTINQ'   
      NX(27)='PLOTQ '   
      NX(28)='DELETE'   
      NX(29)='FIT   '   
      NX(30)='WRITEQ'   
      NX(31)='PREPPQ'   
      NX(32)='OMINC '   
      NX(33)='PHASES'   
      NX(34)='NONLIN'   
      NX(35)='XCORR'   
      NX(36)='      '   
      NX(37)='FILQ  '   
      NX(38)='DTREND'   
      NX(39)='PLOG'     
      NX(40)='PSQLOG'   
      NX(41)='GENWRT'   
      NX(42)='GENRED'   
      NX(43)='DATGEN'   
      NX(44)='      '   
      NX(45)='      '   
      NX(46)='      '   
      NX(47)='      '   
      NX(48)='TIMBAS'   
      NX(49)='      '   
      NX(50)='      '   
      NX(51)='SCALE '   
      NX(52)='PSQAMP'   
      NX(53)='NORMAL'   
      NX(54)='PHWRIT'   
      NX(55)='BIN   '   
      NX(56)='HTOEJD'   
      NX(57)='EPHEM '   
      NX(58)='PLAMP '   
      NX(59)='ZEROQ '   
      NX(60)='NLLSFD'   
      NX(61)='PHREAD'   
      NX(62)='      '   
      NX(63)='MAG   '   
      NX(64)='LSQRV '   
      NX(65)='LSQRV1'   
      NX(66)='PLFILE'
      NX(67)='SORTIM'
      NX(68)='ADDQ'
      NX(69)='PLOTAM'
      NX(70)='PLOTPH'
      NX(71)='PLTOMC'
      NX(72)='DTRANF'
      NX(73)='AVTRAN'
      NX(74)='PLTRNS'
      NX(75)='ARROWS'
      NX(76)='DIV2CH'
      NX(77)='LOGQ'
      NX(78)='SQRTQ'
      NX(79)='AUTOCP'
      NX(80)='TQSCAL'
      NX(81)='TQINVR'
      NX(82)='CCDIN '
      NX(83)='PDOTP '
      NX(84)='SUBQ  '
      NX(85)='EQUALP'
      NX(86)='QEDIN '
      NX(87)='LINCOR'
      NX(88)='SETHJD'
      NX(89)='SORTTQ'
      NX(90)='RESIDL'
      NX(91)='ECLFIT'
      NX(92)='BIN1'
      NX(93)='ECLFT1'
      NX(94)='ECLFT2'

      TITLE=' '
      IFIRST=0                  
      NSHIFT=0    
      LL=0  
      NORUN=0
      NQ=0
      NRNOLD=0
      NQOLD=0     
      ITPHAS=0    
      PI=3.141592653    
      TWOPI=2.0*PI
      DTWOPI=2.0D0*3.141592654D0    
      DO J=1,20
        HARROW(J)=0.0
      END DO

      IF(JUMP.EQ.1)GO TO 35   
   20 INROUT='      '   
      READ(*,30,END=32)INROUT 
   30 FORMAT(A6)  
   32 IF(INROUT.EQ.'      ')GO TO 20
      CALL UPCASE(INROUT)
      IROUTE=INROUT
   35 PRINT 40,IROUTE   
   40 FORMAT(//10X,A6/10X,'------'/)
      DO 50 I=1,94
        IF(IROUTE.EQ.NX(I))GO TO 45 
        GO TO 50  
   45   GO TO (100,200,300,400,500,600,700,800,900,1000,    
     *         1100,1200,1300,1400,1500,1600,1700,1800,     
     *         1900,2000,2100,2200,2300,2400,2500,2600,     
     *         2700,2800,2900,3000,3100,3200,3300,3400,     
     *         3500,3600,3700,3800,3900,4000,4100,4200,     
     *         4300,4400,4500,4600,4700,4800,4900,5000,     
     *         5100,5200,5300,5400,5500,5600,5700,5800,     
     *         5900,6000,6100,6200,6300,6400,6500,6600,
     *         6700,6800,6900,7000,7100,7200,7300,7400,
     *         7500,7600,7700,7800,7900,8000,8100,8200,
     *         8300,8400,8500,8600,8700,8800,8900,9000,
     *         9100,9200,9300,9400),I  
   50 CONTINUE    
      PRINT 60    
   60 FORMAT(10X,'UNRECOGNISED COMMAND') 
      GO TO 700
C     
  100 NQOLD=NQ
      NRNOLD=NORUN
      CALL DISKRQ(Q,NQ,LIM,SAMTIM,HELJD,NORUN,TITLE,NQOLD,NRNOLD,
     *            ITPHAS) 
      IF(NORUN.EQ.1)NSHIFT=0  
      CALL TIMES(TQ,NQ,NORUN,SAMTIM,LIM,HELJD,NQOLD,NRNOLD)  
      GO TO 20    
C     
  200 PRINT *,'*** NOW DEFUNCT ***' 
      GO TO 20    
C     
  300 CONTINUE    
      READ *,NI,NF,MD,NFREQ,(FREQ(J),J=1,NFREQ)    
      IF(NI.LE.0.OR.NI.GE.NQ)NI=1   
      IF(NF.LE.0.OR.NF.GT.NQ)NF=NQ  
      IF(NFREQ.EQ.0)THEN
        IF(FREMAX.GT.0.0)THEN 
          NFREQ=1 
          FREQ(1)=FREMAX
        ELSE
          PRINT *,'LEAST SQUARES FITTING FREQUENCY IS ZERO' 
          GO TO 700
        END IF    
      END IF
      PRINT 310,NFREQ,NI,NF   
      IF(MD.EQ.1)PRINT 320    
  310 FORMAT(12X,'LEAST SQUARES FITTING OF ',I3,' FREQUENCIES',   
     *       ' TO POINTS ',I6,' TO ',I6)  
  320 FORMAT(12X,'THE DATA WILL BE REPLACED BY THE RESIDUALS')    
      CALL LSTSQ(Q,TQ,NI,NF,NFREQ,FREQ,AMP,PHI,SIGMAA,SIGMAP,AFIT,
     *           MD,YMEAN,YSTDEV,RMEAN,RSTDEV)     
      TSPAN=TQ(NF)-TQ(NI)     
      NPT=NF-NI+1 
      PRINT 330,YMEAN,YSTDEV,RMEAN,RSTDEV 
  330 FORMAT(12X,'DATA HAS MEAN ',1PE16.6,' AND STANDARD DEVIATION ',   
     *       E13.6/12X,'RESIDUALS : MEAN ',E13.6,' AND STANDARD', 
     *       ' DEVIATION ',E13.6)   
      PRINT 340   
  340 FORMAT(/12X,'RESULTS : ',2X,'N',8X,'FREQUENCY',10X,'SIGMA',10X,   
     *       'AMPLITUDE',10X,'SIGMA',10X,'PHASE',10X,'SIGMA'/12X,9('='),
     *       1X,100('-'))     
      DO 360 I=1,NFREQ  
        SIGF=RSTDEV/(TSPAN*AMP(I)) * SQRT(6./(NPT*PI*PI))   
        PRINT 350,I,FREQ(I),SIGF,AMP(I),SIGMAA(I),PHI(I),SIGMAP(I)
  350   FORMAT(22X,I3,1PE17.6,E15.6,E19.6,E15.6,E15.6,E15.6)
  360 CONTINUE    
      TMAX=-PHI(1)/(TWOPI*FREQ(1))  
      SIGTMX=SIGMAP(1)/(TWOPI*FREQ(1))    
      HJDMAX=HELJD(1)+TMAX/86400.0D0
      PRINT 370,AFIT,TMAX,SIGTMX,HJDMAX,SIGTMX/86400.0D0
  370 FORMAT(22X,100('-')/22X,'CONSTANT TERM ',1PE13.6/     
     *       22X,'MAXIMUM OF F1 ',1PE13.6,5X,'ERROR ',E13.6/
     *       36X,0PF13.6,11X,F13.6//)     
      GO TO 20    
C     
  400 READ 401,TITLE
  401 FORMAT(A)
      PRINT 402,TITLE
  402 FORMAT(12X,A)  
      GO TO 20    
C     
  500 CALL BRIDGE(Q,NQ) 
      GO TO 20    
C     
  600 CALL QUIST(TQ,Q,NQ,NORUN,SAMTIM,LIM)
      GO TO 20    
C     
  700 IF(IPGPL.EQ.1)CALL PGEND
      STOP  
C
  800 READ 801,PLNAME
  801 FORMAT(A)
      PRINT 802,PLNAME
  802 FORMAT(12X,'ALL PLOTTING TO GO TO FILE ',A) 
      GO TO 20    
C     
  900 READ *,NPOLY,NORM    
      CALL PREP(Q,TQ,1,NQ,SAMTIM(1),NPOLY,NORM)    
      GO TO 20    
C     
 1000 READ *,F1,F2,DF,MODE 
      IF(DF.EQ.0.0)DF=1.0/(5.0D0*(TQ(NQ)-TQ(1)))
      FN=(F2-F1)/DF+1.  
      KD=FN 
      IF(FN-KD.NE.0.0)KD=KD+1 
      IF(KD.GT.MAXFRQ)THEN
        PRINT *,'NUMBER OF FREQUENCIES IS   ',KD
        PRINT *,'MAXIMUM ALLOWED NUNMBER IS ',MAXFRQ
        PRINT *,'RERUN THE PROGRAM WITH SETUP PARAMETERS REVISED'
        GO TO 20
      END IF
      DO 1020 J=1,KD    
        XNU(J)=F1+(J-1)*DF    
        IF(XNU(J).GT.F2)XNU(J)=F2   
 1020 CONTINUE    
      PRINT 1030,KD,XNU(1),XNU(KD),DF,1./MAX(XNU(1),1.0E-10)
     *          ,1./XNU(KD),MODE   
 1030 FORMAT(12X,'POWER SPECTRUM CALCULATED FOR ',I6,' FREQUENCIES ',   
     *       'FROM ',1PE13.6,' TO ',1PE13.6,' IN STEPS OF ',
     *       1PE13.6/42X,'CORRESPONDING PERIODS : ',0PG13.6,4X,   
     *       G13.6,5X,'PROGRAM IN MODE ',I2)    
      IF(MODE.GE.2)THEN 
        READ(*,'(A)')FNAME    
        CLOSE(11) 
        OPEN(11,FILE=FNAME,STATUS='NEW')   
        IF(MODE.EQ.2.OR.MODE.EQ.3)PRINT 1012,FNAME
 1012   FORMAT(12X,'*** SINE AND COSINE SUMS BEING DUMPED'/ 
     *         12X,'*** TO FILE ',A25/)   
        IF(MODE.EQ.4)PRINT 1014,FNAME
 1014   FORMAT(12X,'*** AMPLITUDE SPECTRUM   BEING DUMPED'/ 
     *         12X,'*** TO FILE ',A25/)   
        IF(MODE.EQ.5)PRINT 1015,FNAME
 1015   FORMAT(12X,'*** SPECTRAL WINDOW      BEING DUMPED'/ 
     *         12X,'*** TO FILE ',A25/)   
        WRITE(11,401)TITLE
        IF(MODE.LT.5)THEN
          WRITE(11,1016)NQ,HELJD(1),KD,XNU(1),XNU(KD),DF
 1016     FORMAT('NPTS : ',I6,4X,'HJD : ',F18.7,4X,'NFREQ :',I6/   
     *     'FLOW : ',1PE15.7,4X,'FHIGH : ',1PE15.7,4X,'DF : ',1PE15.7)
        ELSE
          WRITE(11,1016)NQ,HELJD(1),2*KD-1,-XNU(KD),XNU(KD),DF
        END IF
      END IF
      CALL TRANFQ(Q,TQ,NQ,KD,XNU,GG,FF,ST,MODE,FREMAX)
      GO TO 20    
C     
 1100 CALL PLOTPW(XNU,FF,GG,KD,0,TITLE,IROUTE)    
      GO TO 20    
C     
 1200 CALL PLOTPW(XNU,FF,GG,KD,1,TITLE,IROUTE)    
      GO TO 20    
C     
 1300 CONTINUE    
      CALL PRINT(XNU,FF,GG,KD)
      GO TO 20    
C     
 1400 READ *,NI,NF,IPOLY,AFIT,BFIT,CFIT,
     *       NFREQ,(FREQ(J),AMP(J),PHI(J),J=1,INFREQ)
      IF(NI.LE.0.OR.NI.GE.NQ)NI=1   
      IF(NF.LE.0.OR.NF.GT.NQ)NF=NQ  
C     
C     PREWHITEN THE DATA
C     
      DO I=NI,NF
        TI=TQ(I)    
        QI=Q(I)     
        Q0=0.0
        DO J=1,NFREQ  
          DARG=DMOD(DTWOPI*FREQ(J)*TI+PHI(J),DTWOPI)  
          CI=DCOS(DARG)    
          Q0=Q0+AMP(J)*CI       
        END DO      
        DO J=1,IPOLY  
          IF(J.EQ.1)Q0=Q0+AFIT
          IF(J.EQ.2)Q0=Q0+BFIT*TI
          IF(J.EQ.3)Q0=Q0+CFIT*TI*TI
        END DO          
        Q(I)=Q(I)-Q0
      END DO
      GO TO 20    
C     
 1500 CALL ANALYS(Q,TQ,NQ,HELJD(1),TITLE,TIME,AFREQ,AMPLTD,NSHIFT,  
     *            IND,LL,XNU,GG,FF,ST)    
      GO TO 20    
C     
 1600 CALL FOLDPH(TIME,AFREQ,PHASES,SIGPH,OMINC,NSHIFT,HELJD(1),IND,E
     *            ,TITLE)    
      GO TO 20    
C     
 1700 CALL FOLDAM(TIME,AFREQ,AMPLTD,SIGAMP,NSHIFT,HELJD(1),IND,E,    
     *            TITLE)     
      GO TO 20    
C     
 1800 CALL IO(TIME,AFREQ,AMPLTD,PHASES,SIGAMP,SIGPH,NSHIFT,HELJD(1),IND,
     *        E,OMINC,0)
      GO TO 20    
C     
 1900 CALL IO(TIME,AFREQ,AMPLTD,PHASES,SIGAMP,SIGPH,NSHIFT,HELJD(1),IND,
     *        E,OMINC,1)
      GO TO 20    
C     
 2000 CONTINUE                                  
      READ 801,PGPLDEV
      PGTYPE=PGPLDEV
      PGFILE=' '           
      DO J=80,1,-1
        IF(PGPLDEV(J:J).EQ.'/')THEN
          PGTYPE=PGPLDEV(J:)
          IF(J.GT.1)PGFILE=PGPLDEV(1:J-1)
          GO TO 2001
        END IF
      END DO
 2001 CALL UPCASE(PGTYPE)
      PGPLDEV=PGTYPE
      IF(PGFILE.NE.' ')PGPLDEV=PGFILE(1:TRMLEN(PGFILE))//PGPLDEV
      PRINT 2002,PGPLDEV,PGFILE,PGTYPE
 2002 FORMAT(12X,'ALL PGPLOT PLOTTING ON DEVICE ',A30/
     *       12X,'                       FILE   ',A30/
     *       12X,'                       TYPE   ',A30)
      GO TO 20    
C     
 2100 READ *,KS,KE
      IF(KS.EQ.0)KS=1   
      IF(KE.EQ.0)KE=NQ  
      PRINT 2110,KS,KE  
 2110 FORMAT(12X,13HARRAY Q FROM ,I6,5H  TO ,I6,4H   :/)    
      PRINT 2120,(Q(K),K=KS,KE)     
 2120 FORMAT(12X,10G12.4)     
      GO TO 20    
C     
 2200 READ *,KS,KE
      IF(KS.EQ.0)KS=1   
      IF(KE.EQ.0)KE=NQ  
      PRINT 2210,KS,KE  
 2210 FORMAT(12X,16HTIME ARRAY FROM ,I6,5H  TO ,I6,4H   :/) 
      IF(ITPHAS.EQ.0)PRINT 2220,(TQ(K),K=KS,KE) 
 2220 FORMAT(12X,10F12.1)     
      IF(ITPHAS.EQ.1)PRINT 2230,(TQ(K),K=KS,KE) 
 2230 FORMAT(12X,10F12.4)     
      GO TO 20    
C     
 2300 CONTINUE          
      GO TO 20    
C     
 2400 CONTINUE    
      CALL TRANFA(TIME,AMPLTD,NSHIFT,HELJD,XNU,FF,GG,ST,KD,FREMAX)
      GO TO 20    
C     
 2500 SUM=0.0     
      DO 2510 J=1,NSHIFT
        SUM=SUM+AMPLTD(J)     
 2510 CONTINUE    
      AMEAN=SUM/NSHIFT  
      DO 2520 J=1,NSHIFT
        AMPLTD(J)=AMPLTD(J)-AMEAN   
 2520 CONTINUE    
      READ *,FRE,AMPLIT,PHASE    
      PRINT 2530,FRE,AMPLIT,PHASE   
 2530 FORMAT(12X,30HPRE-WHITENING OF AMPS: FREQ = ,1PE10.4, 
     *       11H  AMPLIT = ,E10.3,10H  PHASE = ,E10.3)
      OMEG=DTWOPI*FRE   
      DO 2550 J=1,NSHIFT
        DARG=OMEG*TIME(J)+PHASE     
        DARG=DMOD(DARG,DTWOPI)
        AMPLTD(J)=AMPLTD(J)-(AMPLIT*DCOS(DARG)) 
 2550 CONTINUE    
      GO TO 20    
C     
 2600 CALL PUTINQ(Q,TQ,NQ,NORUN,SAMTIM,LIM,HELJD)     
      GO TO 20    
C     
 2700 IF(NFREQ.EQ.0)THEN
        IF(FREMAX.GT.0.0)THEN 
          NFREQ=1 
          FREQ(1)=FREMAX
        END IF    
      END IF
      CALL PLOTQ(Q,TQ,NQ,NORUN,LIM,HELJD,NFREQ,FREQ,AMP,PHI,AFIT,
     *           BFIT,CFIT,TITLE,ITPHAS,O,C,NSHIFT)
      GO TO 20    
C     
 2800 CALL DELETE(Q,TQ,NQ,NORUN,SAMTIM,LIM,HELJD)     
      GO TO 20    
C     
 2900 CALL FIT(TIME,AFREQ,AMPLTD,SIGAMP,HELJD(1),NSHIFT,TITLE)   
      GO TO 20    
C     
 3000 READ 801,FNAME  
      LUWRIT=20
      IHJD=HELJD(1)     
      RHJD=HELJD(1)-IHJD
      IF(RHJD.LT.0.0000001)RHJD=0.0 
      PRINT 3005,FNAME,TITLE,IHJD,RHJD,SAMTIM(1),NQ  
 3005 FORMAT(12X,'OUTPUT OF ARRAY Q TO FILE ',A/     
     *       12X,'TITLE ',A/12X,'H.J. DAY =',I10/12X,'FRACTION '
     *       ,'OF DAY = ',F8.7/12X,'SAMPLING TIME = ',F6.1/ 
     *       12X,'NO OF POINTS = ',I6)    
      OPEN(LUWRIT,FILE=FNAME,STATUS='NEW')
      WRITE(LUWRIT,3010)TITLE,IHJD,RHJD,SAMTIM(1),NQ  
 3010 FORMAT(A36,I10,F8.7,10X,F4.1,2X,2I6)
      NC=(NQ+9)/10
      NCM=NC-1    
      DO 3050 J=0,NCM   
        DO 3020 I=1,10  
          IQ(I)=Q(I+J*10)     
 3020   CONTINUE  
        WRITE(LUWRIT,3040)(IQ(KK),KK=1,10),J    
 3040   FORMAT(10I7,I10)
 3050 CONTINUE    
      CLOSE(LUWRIT)
      GO TO 20    
C     
 3100 READ *,NPOLY,NORM    
      IF(NORM.EQ.999)GO TO 3180     
      DO 3150 J=1,1000  
        IF(J.GT.NORUN)GO TO 20
        NI=LIM(1,J)     
        NF=LIM(2,J)     
        DELT=SAMTIM(J)  
        CALL PREP(Q,TQ,NI,NF,DELT,NPOLY,NORM)
 3150 CONTINUE    
      PRINT 3175  
 3175 FORMAT(12X,23HBAD RETURN FROM PREPPQ )    
      GO TO 20    
C     
C     THIS IS THE SECTION TO DO PREP-ING IN SECTIONS SPECIFIED    
C     BY THE USER AND NOT SPECIFIED BY THE RUNS 
C     
 3180 READ *,NORM    
      READ *,NI,N,ISHIFT,NSHIFT  
      IF(NI.LE.0)NI=1   
      IF(N.LE.0.OR.N.GT.NQ)THEN     
        PRINT *,'BAD VALUE FOR NUMBER OF POINTS TO BE FITTED'     
        GO TO 20  
      END IF
      IF(ISHIFT.EQ.0)ISHIFT=N 
      IF(NSHIFT.EQ.0)NSHIFT=(NQ/ISHIFT)-(N/ISHIFT)+1  
      DO 3190 I=1,NSHIFT
        NII=NI+(I-1)*ISHIFT   
        NFI=NII+N-1     
        IF(NFI.GT.NQ)GO TO 3190     
        IF(NFI.GT.NQ)NFI=NQ   
        CALL PREP(Q,TQ,NII,NFI,DELT,NPOLY,NORM)    
 3190 CONTINUE    
      GO TO 20    
C     
 3200 CONTINUE    
      CALL SOMINC(Q,TQ,NQ,HELJD,TITLE,TIME,AFREQ,AMPLTD,PHASES,  
     *          SIGAMP,SIGPH,NSHIFT,IND,E,O,C,OMINC,LL,PHCORR,
     *          FILPER,FILRES)
      GO TO 20    
C     
 3300 CONTINUE    
      CALL SPHASE(Q,TQ,NQ,HELJD,TITLE,TIME,AFREQ,AMPLTD,PHASES,  
     *            SIGAMP,SIGPH,IND,NSHIFT,LL,PHCORR,FILPER,FILRES)
      GO TO 20    
C     
 3400 READ *,NI,NF,MD,IPOLY,INFREQ,(FREQ(J),AMP(J),PHI(J),J=1,INFREQ)
      IF(NI.LE.0.OR.NI.GE.NQ)NI=1   
      IF(NF.LE.0.OR.NF.GT.NQ)NF=NQ  
      IF(INFREQ.GT.0)THEN     
        NFREQ=INFREQ    
      ELSE IF(IPOLY.EQ.0.AND.NFREQ.EQ.0)THEN    
        IF(FREMAX.GT.0.0)THEN 
          NFREQ=1 
          FREQ(1)=FREMAX
          AMP(1)=0.0    
          PHI(1)=0.0    
        ELSE
          PRINT *,'NO LEAST SQUARES FITTING FREQUENCY'
          GO TO 700
        END IF    
      END IF
      CALL NONLIN(Q,TQ,NI,NF,IPOLY,AFIT,BFIT,CFIT,SIGA,SIGB,SIGC, 
     *            NFREQ,FREQ,AMP,PHI,SIGMAF,SIGMAA,SIGMAP,IPRINT,MD)    
      GO TO 20    
C     
C
 3500 IF(NFREQ.EQ.0)THEN
        IF(FREMAX.GT.0.0)THEN 
          NFREQ=1 
          FREQ(1)=FREMAX
        END IF    
      END IF
      CALL XCORR(Q,TQ,NQ,NORUN,LIM,HELJD,NFREQ,FREQ,AMP,PHI,AFIT,
     *           BFIT,CFIT,TITLE,ITPHAS,O,C,NSHIFT)
      GO TO 20    

C     
 3600 CONTINUE                                              
      GO TO 20    
C     
 3700 DELT=SAMTIM(1)    
      CALL FILQ(TQ,Q,NQ,DELT,PHCORR,FILPER,FILRES)    
      GO TO 20    
C     
 3800 CALL DTREND(Q,NQ) 
      GO TO 20    
C     
 3900 CONTINUE    
 4000 DO 4010 J=1,KD    
        IF(IROUTE.EQ.'PSQLOG')FF(J)=FF(J)*FF(J) 
        FF(J)=ALOG10(FF(J))   
 4010 CONTINUE    
      CALL PLOTPW(XNU,FF,GG,KD,2,TITLE,IROUTE)   
      GO TO 20    
C     
 4100 CALL GENWRT(Q,TQ,NQ,HELJD)    
      GO TO 20    
C     
 4200 CALL GENRED(Q,TQ,NQ,NORUN,SAMTIM,LIM,HELJD,ITPHAS)     
      IF(NORUN.EQ.1)NSHIFT=0
      GO TO 20    
C     
 4300 CALL DATGEN(Q,TQ,NQ,NORUN,SAMTIM,LIM,HELJD)     
      GO TO 20    
C     
 4400 CONTINUE
      GO TO 20    
C     
 4500 CONTINUE                      
      GO TO 20    
C     
 4600 CONTINUE                                              
      GO TO 20    
C     
 4700 CONTINUE                                              
      GO TO 20    
C     
 4800 CONTINUE    
      CALL TIMBAS(TQ,NQ,HELJD,LIM,NORUN,EPOCH,EPER,ESEC,FREMAX,ITPHAS,
     *            HJDMAX)  
      GO TO 20    
C     
 4900 CONTINUE                                        
      GO TO 20    
C     
 5000 CONTINUE                                        
      GO TO 20    
C     
 5100 READ *,SCALEF     
      PRINT 5105,SCALEF 
 5105 FORMAT(12X,'DATA TO BE SCALED BY FACTOR ',F15.5)
      DO 5110 J=1,NQ    
        Q(J)=Q(J)*SCALEF
 5110 CONTINUE    
      GO TO 20    
C     
 5200 DO 5210 J=1,KD    
 5210 FF(J)=FF(J)*FF(J) 
      CALL PLOTPW(XNU,FF,GG,KD,0,TITLE,IROUTE)   
      GO TO 20    
C     
 5300 DO 5350 I=1,NORUN 
        READ *,FACTOR   
        PRINT 5310,I,LIM(1,I),LIM(2,I),FACTOR   
 5310   FORMAT(12X,'RUN ',I2,' FROM ',I6,' TO ',I6,' NORMALIZED', 
     *         ' BY FACTOR ',F15.5) 
        DO 5320 J=LIM(1,I),LIM(2,I) 
          Q(J)=Q(J)/FACTOR    
 5320   CONTINUE  
 5350 CONTINUE    
      GO TO 20    
C     
 5400 CALL PHWRIT(Q,TQ,NQ,LIM,SAMTIM,HELJD,NORUN)     
      GO TO 20    
C     
 5500 IF(ITPHAS.EQ.0)THEN     
        PRINT 5510
 5510   FORMAT(12X,'I CANT BIN BECAUSE THE TIME ARRAY HAS'/ 
     *         12X,'NOT BEEN TRAMSFORMED TO PHASE BY E.G.'/ 
     *         12X,'INSTRUCTION TIMBAS')  
        PRINT *,'BIN ERROR'
        GO TO 700
      END IF
      CALL BIN(Q,TQ,NQ) 
      GO TO 20    
C     
 5600 CALL HTOEJD(Q,TQ,NQ,LIM,SAMTIM,HELJD,NORUN)     
      GO TO 20    
C     
 5700 CONTINUE    
      CALL EPHEM(EPOCH,EPER,ESEC,ETHIRD,E,OMINC,NSHIFT,TITLE)    
      GO TO 20    
C     
 5800 CONTINUE    
      CALL PLAMP(TIME,AFREQ,AMPLTD,NSHIFT,HELJD(1),IND)
      GO TO 20    
C     
 5900 CONTINUE    
      DO 5910 J=1,MAXPT 
        Q(J)=0.0  
        TQ(J)=0.0 
 5910 CONTINUE    
      NORUN=0     
      NQ=0
      NQOLD=0
      NRNOLD=0  
      TITLE=' '
      GO TO 20    
C     
 6000 CONTINUE    
      READ *,NI,NF,MD,INFREQ,(FREQ(J),AMP(J),PHI(J),J=1,INFREQ)
      IF(NI.LE.0.OR.NI.GE.NQ)NI=1   
      IF(NF.LE.0.OR.NF.GT.NQ)NF=NQ  
      IF(INFREQ.GT.0)THEN     
        NFREQ=INFREQ    
      ELSE IF(IPOLY.EQ.0.AND.NFREQ.EQ.0)THEN    
        IF(FREMAX.GT.0.0)THEN 
          NFREQ=1 
          FREQ(1)=FREMAX
          AMP(1)=0.0    
          PHI(1)=0.0    
        ELSE
          PRINT *,'NO LEAST SQUARES FITTING FREQUENCY'
          GO TO 700
        END IF    
      END IF
      CALL NLLSFD(Q,TQ,NI,NF,NFREQ,FREQ,AMP,PHI,FDOT,SIGMAF,SIGMAA,     
     *            SIGMAP,SGFDOT,IPRINT,MD)
      GO TO 20    
C     
 6100 CONTINUE    
      CALL PHREAD(Q,TQ,NQ,LIM,SAMTIM,HELJD,NORUN,ITPHAS)    
      GO TO 20    
C     
 6200 CONTINUE    
      GO TO 20    
C     
 6300 CONTINUE    
      READ *,ZEROPT     
      IF(ZEROPT.GE.0.0)THEN
        DO 6310 J=1,NQ    
 6310   Q(J)=-2.5*ALOG10(Q(J))+ZEROPT 
        PRINT 6320,ZEROPT 
 6320   FORMAT(12X,'DATA CONVERTED TO MAG SCALE WITH ZERO PT ',F8.3)
      ELSE
        DO 6330 J=1,NQ    
 6330   Q(J)=10.0**(0.4*(-ZEROPT-Q(J)))
        PRINT 6340,ZEROPT 
 6340   FORMAT(12X,'DATA CONVERTED FROM A MAG SCALE WITH ZERO PT ',
     *         F8.3)
      END IF
      GO TO 20    
C     
 6400 CONTINUE    
      READ *,NI,NF,MD,FREQ(1)    
      IF(NI.LE.0.OR.NI.GE.NQ)NI=1   
      IF(NF.LE.0.OR.NF.GT.NQ)NF=NQ  
      NFREQ=1     
      IF(FREQ(1).EQ.0)THEN    
        IF(FREMAX.GT.0.0)THEN 
          NFREQ=1 
          FREQ(1)=FREMAX
        ELSE
          PRINT *,'LEAST SQUARES FITTING FREQUENCY IS ZERO' 
          GO TO 700
        END IF    
      END IF
      PRINT 310,NFREQ,NI,NF   
      IF(MD.EQ.1)PRINT 320    
      CALL LSQRV(Q,TQ,NI,NF,NFREQ,FREQ,AMP,PHI,SIGMAA,SIGMAP,AFIT,
     *           SIGA,MD,YMEAN,YSTDEV,RMEAN,RSTDEV)  
      NPT=NF-NI+1 
      PRINT 330,YMEAN,YSTDEV,RMEAN,RSTDEV 
      PRINT 6440  
 6440 FORMAT(/12X,'RESULTS : ',8X,'FREQUENCY',10X,'AMPLITUDE',10X,
     *       'SIGMA',10X,'PHASE',10X,'SIGMA',5X,'GAMMA',5X,'SIGMA'/     
     *       12X,9('='),1X,100('-'))
      PRINT 6450,FREQ(1),AMP(1),SIGMAA(1),PHI(1),SIGMAP(1),AFIT,SIGA     
 6450 FORMAT(22X,1PE17.6,E19.6,E15.6,E15.6,E15.6,0PF10.5,0PF10.5) 
      TMAX=-PHI(1)/(TWOPI*FREQ(1))  
      SIGTMX=SIGMAP(1)/(TWOPI*FREQ(1))    
      HJDMAX=HELJD(1)+TMAX/86400.0D0
      PRINT 6470,TMAX,SIGTMX,HJDMAX,SIGTMX/86400.0D0  
 6470 FORMAT(22X,100('-')/    
     *       22X,'MAXIMUM OF F1 ',1PE13.6,5X,'ERROR ',E13.6/
     *       36X,0PF13.6,11X,F13.6//)     
      GO TO 20    
C     
 6500 CONTINUE    
      READ *,NI,NF,MD,FREQ(1)    
      IF(NI.LE.0.OR.NI.GE.NQ)NI=1   
      IF(NF.LE.0.OR.NF.GT.NQ)NF=NQ  
      NFREQ=1     
      IF(FREQ(1).EQ.0)THEN    
        IF(FREMAX.GT.0.0)THEN 
          NFREQ=1 
          FREQ(1)=FREMAX
        ELSE
          PRINT *,'LEAST SQUARES FITTING FREQUENCY IS ZERO' 
          GO TO 700
        END IF    
      END IF
      PRINT 310,NFREQ,NI,NF   
      IF(MD.EQ.1)PRINT 320    
      CALL LSQRV1(Q,TQ,NI,NF,NFREQ,FREQ,AMP,PHI,SIGMAA,SIGMAP,AFIT,
     *           SIGA,MD,YMEAN,YSTDEV,RMEAN,RSTDEV)  
      NPT=NF-NI+1 
      PRINT 330,YMEAN,YSTDEV,RMEAN,RSTDEV 
      PRINT 6440  
      PRINT 6450,FREQ(1),AMP(1),SIGMAA(1),PHI(1),SIGMAP(1),AFIT,SIGA     
      TMAX=-PHI(1)/(TWOPI*FREQ(1))  
      SIGTMX=SIGMAP(1)/(TWOPI*FREQ(1))    
      HJDMAX=HELJD(1)+TMAX/86400.0D0
      PRINT 6470,TMAX,SIGTMX,HJDMAX,SIGTMX/86400.0D0  
      GO TO 20    
C
 6600 GO TO 800
C     
 6700 CALL SORTIM(Q,TQ,NQ,NORUN,SAMTIM,LIM,HELJD,0)
      GO TO 20
C
 6800 CONTINUE
      DO 6850 I=1,NORUN 
        READ *,ADDQ
        AVNEW=0.0
        DO 6820 J=LIM(1,I),LIM(2,I) 
          Q(J)=Q(J)+ADDQ      
          AVNEW=AVNEW+Q(J)
 6820   CONTINUE  
        AVNEW=AVNEW/(LIM(2,I)-LIM(1,I)+1)
        PRINT 6810,I,LIM(1,I),LIM(2,I),ADDQ,AVNEW
 6810   FORMAT(12X,'RUN ',I2,' FROM ',I6,' TO ',I6,' INCREMENTED',
     *         ' BY ',G12.3,' : NEW AVERAGE IS ',G12.3) 
 6850 CONTINUE    
      GO TO 20
C
 6900 CONTINUE
      XL=8.0
      IEB=-1
      CALL PLOTAM(XL,TIME,AMPLTD,SIGAMP,NSHIFT,TITLE,IEB)
      GO TO 20
C
 7000 CONTINUE
      XL=8.0
      IEB=1
      CALL PLOTPH(XL,TIME,PHASES,SIGPH,E,OMINC,NSHIFT,TITLE,IEB)
      GO TO 20
C
 7100 CONTINUE
      XL=8.0
      IEB=-1
      CALL PLOTPH(XL,TIME,PHASES,SIGPH,E,OMINC,NSHIFT,TITLE,IEB)
      GO TO 20
C
 7200 CONTINUE
C
C     DISCRETE FOURIER TRANSFORM WITH DOUBLE PRECISION FREQUENCY ARRAY
C
      IF(IFIRST.EQ.0)THEN
C
C       ALLOCATE DYNAMIC MEMORY FOR THE DOUBLE PRECISION FREQUENCY
C       ARRAY
C
        IFIRST=1

C  Vax statements; replaced by followind code
C        STATUS = LIB$GET_VM( MAXFRQ*8, ADXNU )
C        IF(.NOT.STATUS) STOP 'Error in dynamic mem allocation for DXNU'
        

c        PRINT 7220
c 7220   FORMAT(12X,'DOUBLE PRECISION FREQUENCIES MEMORY ALLOCATED')
      END IF
      CALL DTRANF(Q,TQ,NQ,KD,XNU,MXFRQ,GG,FF,ST,MODE,FREMAX,FNAME,
     *            TITLE,HELJD)
      GO TO 20
C
 7300 CONTINUE
      CALL AVTRAN(KD,XNU,FF,HELJD)
      GO TO 20
C
 7400 CONTINUE
      CALL PLTRNS(XNU,FF,KD,TITLE)    
      GO TO 20    
C


 7500 CONTINUE
      READ *,NARROW,(FARROW(J),J=1,IABS(NARROW))
      IF(NARROW.LT.0.0)READ *,(HARROW(J),J=1,IABS(NARROW))
      NARROW=IABS(NARROW)
      GO TO 20
C
 7600 CONTINUE
      CALL DIV2CH(Q,NQ,LIM,SAMTIM,HELJD,NORUN,TITLE) 
      GO TO 20
C     
 7700 CONTINUE    
      READ *,ZEROPT     
      DO J=1,NQ    
        Q(J)=ALOG10(Q(J))+ZEROPT 
      END DO
      PRINT 7720,ZEROPT 
 7720 FORMAT(12X,'DATA CONVERTED TO LOG SCALE WITH ZERO PT ',F8.3)
      GO TO 20    
C
 7800 CONTINUE
      DO J=1,NQ    
        Q(J)=SQRT(Q(J)) 
      END DO
      PRINT 7820
 7820 FORMAT(12X,'DATA CONVERTED TO SQUARE ROOT OF ITSELF ')
      GO TO 20
C
 7900 CONTINUE
      CALL AUTOCP(Q,TQ,NQ,KD,XNU,GG,FF,ST,MODE,FREMAX)
      GO TO 20
C
 8000 CONTINUE
      READ *,SCALET
      PRINT 8005,SCALET 
 8005 FORMAT(12X,'TIME ARRAY TO BE SCALED BY FACTOR ',1PE12.4)
      DO J=1,NQ    
        TQ(J)=TQ(J)*SCALET
      END DO
      HELJD(1)=HELJD(1)*SCALET 
      GO TO 20
C
 8100 CONTINUE
      PRINT 8105
 8105 FORMAT(12X,'TIME ARRAY TO BE REPLACED BY ITS INVERSE')
      DO J=1,NQ    
        IF(TQ(J).NE.0.0)THEN
          TQ(J)=1.0/TQ(J)
        ELSE
          TQ(J)=1.0E30
        END IF
      END DO
C
C     NOW INVERT THE ORDERING SO THE SMALLEST TIME IS NEAR
C     THE START OF THE ARRAY
C
      NQQ=NQ/2
      DO J=1,NQQ
        STOREQ=Q(J)
        STORTQ=TQ(J)
        Q(J)=Q(NQ+1-J)
        TQ(J)=TQ(NQ+1-J)
        Q(NQ+1-J)=STOREQ
        TQ(NQ+1-J)=STORTQ
      END DO
      GO TO 20
C
 8200 CONTINUE
      CALL CCDIN(Q,TQ,NQ,NORUN,SAMTIM,LIM,HELJD)     
      GO TO 20
C
 8300 CONTINUE
      CALL PDOTP(EPOCH,EPER,ESEC,ETHIRD,E,OMINC,NSHIFT,TITLE)    
      GO TO 20
C
 8400 CONTINUE
      READ *,SUBQ
      PRINT 8405,SUBQ 
 8405 FORMAT(12X,1PE12.4,' TO BE SUBTRACTED FROM DATA ARRAY Q')
      DO J=1,NQ    
        Q(J)=Q(J)-SUBQ
      END DO
      GO TO 20
C
 8500 CONTINUE
      CALL EQUALP(Q,TQ,NQ,KD,XNU,FF,FREMAX,TITLE,HELJD)
      GO TO 20
C     
 8600 CALL QEDIN(Q,TQ,NQ,NORUN,SAMTIM,LIM,HELJD,TITLE)     
      IF(NORUN.EQ.1)NSHIFT=0
      GO TO 20    
C     
 8700 CONTINUE    
      CALL LINCOR(E,OMINC,NSHIFT,TITLE)    
      GO TO 20    
C     
 8800 CONTINUE
      READ *,HJD1NU
      PRINT 8805,HJD1NU 
 8805 FORMAT(12X,'NEW HJD1: ',F18.7)
      DO J=NORUN,1,-1
        HELJD(J) = HELJD(J) + (HJD1NU-HELJD(1))
      END DO
      GO TO 20
C     
 8900 CONTINUE    
      CALL SORTIM(Q,TQ,NQ,NORUN,SAMTIM,LIM,HELJD,1)
      GO TO 20    
C     
 9000 READ *,F1,F2,DF,NFREQ  
C
C     KEEP THE LOWEST FREQUENCY AWAY FROM ZERO AS THINGS
C     GO A BIT FUNNY THERE
C
      IF(F1.EQ.0.0)F1=F2/20.0

      IF(DF.EQ.0.0)DF=1.0/(5.0D0*(TQ(NQ)-TQ(1)))
      FN=(F2-F1)/DF+1.  
      KD=FN 
      IF(FN-KD.NE.0.0)KD=KD+1 
      IF(KD.GT.MAXFRQ)THEN
        PRINT *,'NUMBER OF FREQUENCIES IS   ',KD
        PRINT *,'MAXIMUM ALLOWED NUMBER IS ',MAXFRQ
        PRINT *,'RERUN THE PROGRAM WITH SETUP PARAMETERS REVISED'
        GO TO 20
      END IF
      DO J=1,KD    
        XNU(J)=F1+(J-1)*DF    
        IF(XNU(J).GT.F2)XNU(J)=F2   
      END DO      
      PRINT 9030,KD,XNU(1),XNU(KD),DF,1./MAX(XNU(1),1.0E-10)
     *          ,1./XNU(KD),NFREQ 
 9030 FORMAT(12X,'RESIDUALGRAM CALCULATED FOR ',I6,' FREQUENCIES ',   
     *       'FROM ',1PE13.6,' TO ',1PE13.6,' IN STEPS OF ',
     *       1PE13.6/42X,'CORRESPONDING PERIODS : ',0PG13.6,4X,   
     *       G13.6,5X,'NUMBER OF FOURIER TERMS ',I2)
      IF(NFREQ.LT.0)THEN
        READ(*,'(A)')FNAME    
        CLOSE(11) 
        OPEN(11,FILE=FNAME,STATUS='NEW')   
        PRINT 9031,FNAME
 9031   FORMAT(12X,'*** RESIDUAL SPECTRUM BEING DUMPED'/ 
     *         12X,'*** TO FILE ',A25/)   
        WRITE(11,401)TITLE
        WRITE(11,1016)NQ,HELJD(1),KD,XNU(1),XNU(KD),DF
      END IF
      CALL RESIDL(Q,TQ,NQ,KD,XNU,FF,FREMAX,NFREQ, 
     *            FREQ,AMP,PHI,SIGMAA,SIGMAP)
      GO TO 20    
C     
 9100 CONTINUE    
      READ *,NI,NF,MD
      read *,epoch,eper
      IF(NI.LE.0.OR.NI.GE.NQ)NI=1   
      IF(NF.LE.0.OR.NF.GT.NQ)NF=NQ  
      NPARM=3
c
c     find rough value for mid eclipse
c
      sum1=0.0
      sum2=0.0
      do j=1,10
        sum1=sum1+q(ni+j-1)/10.0
        sum2=sum2+q(nf+1-j)/10.0
      end do
      qmid=0.5*(sum1+sum2)
      do j=ni+1,nf-1
        qav=(q(j-1)+q(j)+q(j+1))/3.0
        if(sum2.gt.sum1.and.qav.gt.qmid)then
          CFIT=tq(j-1)
          go to 9110
        end if
        if(sum2.lt.sum1.and.qav.lt.qmid)then
          CFIT=tq(j-1)
          go to 9110
        end if
      end do
 9110 AFIT=MAX(sum1,sum2)
      BFIT=ABS(sum2-sum1)
      if(sum2.le.sum1)modee=1
      if(sum2.gt.sum1)modee=-1
      CALL ECLFIT(Q,TQ,NI,NF,NPARM,AFIT,BFIT,CFIT,SIGA,SIGB,SIGC, 
     *            IPRINT,MD,MODEE)    
      open(11,file='times.dat',status='unknown',access='append')
      ttt=heljd(1)+cfit/86400.0d0
      cycl=nint((ttt-epoch)/eper)
      write(11,*)ttt,cycl
      close(11)
      GO TO 20    
C     
 9200 CALL BIN1(Q,TQ,NQ,NORUN,SAMTIM,LIM,HELJD)     
      GO TO 20    
C     
 9300 CONTINUE    
      READ *,NI,NF,MD
      read *,nparam,ai 
      IF(NI.LE.0.OR.NI.GE.NQ)NI=1   
      IF(NF.LE.0.OR.NF.GT.NQ)NF=NQ  
      CALL ECLFT1(Q,TQ,NI,NF,IPRINT,MD,TITLE)
      GO TO 20    
C     
 9400 CONTINUE    
      READ *,NI,NF,MD
      read *,nparam,ai 
      IF(NI.LE.0.OR.NI.GE.NQ)NI=1   
      IF(NF.LE.0.OR.NF.GT.NQ)NF=NQ  
      CALL ECLFT2(Q,TQ,NI,NF,IPRINT,MD,TITLE)
      GO TO 20    
C     
      END   
