      SUBROUTINE ANALYS(DATA,DATIME,NQ,JDH1,TITLE,
     *                  AT,AF,AA,IPART,IND,LL,XNU,GG,FF,ST) 
C     
C     THIS IS AN ADAPTATION OF THE ROUTINE FROM SQUID OF THE SAME NAME  
C     IT IS VIRTUALLY IDENTICAL WITH THE EXCEPTION THAT THE 
C     CALCULATION OF THE AMPLITUDES OF THE OSCILLATION IS DONE BY A     
C     FOURIER TRANSFORM RATHER THAN BY A PERIODOGRAM. THIS IS     
C     EFFECTED BY A CALL TO TRANFQ WHICH IS A ROUTINE TAKEN FROM  
C     DEEMING, T. J. 1975 , ASTR. SP. SCI., 36, 137.  
      DIMENSION DATA(*),AT(*),AF(2,1000),AA(*)
      double precision datime(*),qtime(1500)
      DIMENSION Q(1500),IND(*)
      DIMENSION XNU(*),GG(*),FF(*),ST(*),D(2,256)  
      CHARACTER*80 TITLE
      DOUBLE PRECISION JDH1   
C     
C     SERIES IS PASSED OVER IN DATA AND DATIME FROM THE MAIN PROGRAM    
C     
C     
C     
      PRINT 1,TITLE    
    1 FORMAT(1H0,11X,A/)    
      READ *,PERIOD,NCYC,MPART,NPER,XL,YL,GRID,PP1,PP2  
      NOJUMP=0
      NOPLOT=0
      IF(PERIOD.LT.0.0)THEN
C
C       ENSURE EQUALLY-SPACED DATA ANALYSIS
C
        NOJUMP=1
        PERIOD=-PERIOD
      END IF
      IF(NCYC.LT.0)THEN
C
C       NO PLOTTING REQUIRED
C
        NOPLOT=1
        NCYC=-NCYC
      ELSE
        LOG=99
        OPEN(UNIT=LOG,STATUS='SCRATCH')
      END IF
C     
C     THE VARIABLES:    
C     ==============    
C     
C     NCYC IS NO OF CYCLES OF DATA TO USE IN EACH PERIODOGRAM     
C     MPART IS NO OF SECTIONS TO SPLIT THIS DATA UP INTO,MOVING   
C       ALONG ONE SECTION EACH TIME 
C     NPER IS THE NUMBER OF PERIODS IN EACH PERIODOGRAM     
C     NN IS THE NUMBER OF POINTS USED IN EACH PERIODOGRAM   
C       THIS IS EQUAL TO NCYC*PERIOD WHERE THE PERIOD IS SPECIFIED
C     MM IS THE NUMBER OF POINTS TO SLIDE ALONG EACH TIME   
C     LPART IS THE NUMBER OF PARTITIONS   
C     PERIOD IS THE CENTRAL PERIOD IN THE PERIODOGRAMS AND ALSO GIVES   
C       NN WHERE NN=PERIOD*NCYC     
C     
      DELT=DATIME(2)-DATIME(1)
      IF(NCYC.LE.0) NCYC=40   
      IF(MPART.LE.0) MPART=2  
      IF(NPER.LE.0) NPER=128  
      IF(XL.LE.0.0) XL=5.0    
      IF(YL.LE.0.0) YL=5.0    
      IF(GRID.LE.0) GRID=1.0  
      IF(MPART.GT.10)OVER=FLOAT(MPART)    
      IF(MPART.GT.10)MPART=100.0/(100.0-OVER)   
      IF(OVER.EQ.0.0)THEN     
        IF(MPART.EQ.1)OVER=0.0
        IF(MPART.NE.1)OVER=100.0*(1.0-1.0/MPART)
      END IF
      IF(PP2.EQ.0.0)PP2=99999.9     
      NN=PERIOD*NCYC/DELT     
      MM=NN/MPART 
      KPART=NQ/MM 
      LPART=KPART-(MPART-1)   
      NN=MPART*MM 
      IF(NN.LE.NQ)GO TO 15    
      PRINT 5     
    5 FORMAT(1H0,11X,36HDATA TOO SHORT FOR PROPOSED ANALYSIS/)    
      RETURN
   15 TIMSTP=MM*DELT    
      T=MM*MPART*DELT   
      FRE=1/PERIOD
      DF=1/(T*GRID)     
      DT=DF*PERIOD**2   
      D(1,NPER)=FRE+DF/2
      D(1,NPER+1)=FRE-DF/2    
      NPERM1=NPER-1     
      FNYQ=1.0/(2.0*DELT)     
      DO 10 I=1,NPERM1  
        D(1,NPER-I)=D(1,NPER-I+1)     
        IF(D(1,NPER-I+1)+DF.LE.FNYQ)D(1,NPER-I)=D(1,NPER-I+1)+DF    
        FNEW=D(1,NPER+I)-DF     
        D(1,NPER+I+1)=0.0
        IF(FNEW.GT.0.0)D(1,NPER+I+1)=FNEW   
   10 CONTINUE    
      MAXD=2*NPER 
      IF(D(1,MAXD).GT.0.0)THEN
        PMAX=1/D(1,MAXD)  
      ELSE
        PMAX=99999.0
      END IF
      PMIN=1/D(1,1)     
      IF(NOJUMP.EQ.1)PRINT 20
   20 FORMAT(12X,'ONLY EQUALLY-SPACED SECTIONS ANALYSED'/)
      PRINT  30,PERIOD,NCYC,TIMSTP,PMAX,PMIN,NPER,XL,YL,GRID,MPART
     *          ,OVER   
   30 FORMAT(12X,'PERIOD',F12.4,'   NCYC',I5,'   TIMSTP',F10.2,'   PMAX'
     2,F10.4,'   PMIN',F10.4,//12X,'NPER',I5,'  XL',F6.2,'   YL',F6.2,  
     3'   GRID ',F6.2,'   MPART',I5,'    % OVERLAP ',F5.1/) 
      PRINT 35,DF,DT,FRE,T,NN 
   35 FORMAT(12X,3HDF ,1PE10.4,5H  DT ,0PF6.4,8H   FREQ ,1PE11.4, 
     *       16H  LENGTH OF DATA ,0PF8.1,24H (S)    NUMBER OF POINTS,I5 
     *       /)   
      PRINT 38,PP1,PP2    
   38 FORMAT(12X,'SEARCH FOR MAXIMA CONFINED FROM ',F7.1,' TO ',F7.1/)   
      PRINT 40    
   40 FORMAT(16X,53H CHANNEL      PERIOD     FREQUENCY     AMPLITUDE    
     2 ,9H AV POINT,25H     NI   NF  TIME (SECS)/12X,91(1H=))     
C     
C     
C     
      LL=LL+1     
      NI=-MM+1
      LOOP=0
      DO 500 JJ=1,LPART 
        NI=NI+MM    
   45   NF=NI+NN-1  
        IF(NF.GT.NQ)GO TO 500
        NU=0  
        SQ=0  
        DO 50 J=NI,NF     
          IF(NOJUMP.EQ.1.AND.J.GT.NI)THEN
C
C           CHECK EQUALLY-SPACED DATA
C
            IF( (DATIME(J)-DATIME(J-1)) .GT. (DELT*1.1) .OR.
     *          (DATIME(J)-DATIME(J-1)) .LT. (DELT*0.9) )THEN
               NI=J
               GO TO 45
            END IF
          END IF
          Q(J-NI+1)=DATA(J) 
          QTIME(J-NI+1)=DATIME(J)-DATIME(NI)  
          QTIME(J-NI+1)=DATIME(J) 
          IF(Q(J-NI+1).LT.-1.E9) GO TO 50     
          SQ=SQ+Q(J-NI+1)   
          NU=NU+1     
   50   CONTINUE    
        AQ=SQ/NU    
        DO 54 J=1,NN
          Q(J)=Q(J)-AQ    
   54   CONTINUE    
   55   YMX=0.0     
        DO 100 J=1,MAXD   
          XNU(J)=D(1,J)   
  100   CONTINUE    
C     
C       CALCULATE THE FOURIER TRANSFORM AT THE CALCULATED     
C       FREQUENCIES 
C     
        FREMAX=-1.0 
        CALL TRANFQ(Q,QTIME,NN,MAXD,XNU,GG,FF,ST,0,FREMAX)    
C     
C       NOW TRANSFER THE RESULTS BACK TO ARRAYS THAT ARE USED 
C       IN THE SQUID VERSION OF ANALYS. ALSO CALCULATE THE MAXIMUM  
C       AMPLITUDE IN THIS SPECTRUM    
C     
        DO 200 J=1,MAXD   
          IF(XNU(J).GT.0.0)THEN
            PP=1.0/XNU(J)   
          ELSE
            PP=99999.9
          END IF
          D(2,J)=FF(J)    
          IF(FF(J).LT.YMX.OR.PP.LT.PP1.OR.PP.GT.PP2)GO TO 200 
          IMX=J     
          YMX=FF(J) 
  200   CONTINUE    
  440   IF(D(1,IMX).GT.0.0)THEN
          PERIOD=1/D(1,IMX) 
        ELSE
          PERIOD=99999.9
        END IF
        LOOP=LOOP+1
        IPART=IPART+1     
        AT(IPART)=(DATIME(NI)+DATIME(NF))/2.
        PRINT 450,IPART,IMX,PERIOD,D(1,IMX),D(2,IMX),AQ,NI,NF,AT(IPART)   
  450   FORMAT(12X,I3,I9,F12.5,1P3E14.6,I7,I5,0PF13.1)     
        AF(1,IPART)=D(1,IMX)    
        AF(2,IPART)=PERIOD
        AA(IPART)=D(2,IMX)
        IND(IPART)=LL     
        IF(NOPLOT.EQ.0)WRITE(LOG,475)JJ,(D(2,KKK),KKK=1,256)  
  475   FORMAT(I3,(16E8.3))     
  500 CONTINUE    
      LPART=LOOP
      PRINT 600   
  600 FORMAT(12X,91(1H=))     
      PRINT 700,JDH1    
  700 FORMAT(12X,25HTIMES RELATIVE TO HJD OF ,F18.7/) 
      PRINT 800,LPART,IPART
  800 FORMAT(12x,'LPART ',I5,'    IPART ',I5)
      IF(NOPLOT.EQ.1)RETURN
      ENDFILE LOG 
      REWIND LOG  
      CALL PLOTTE(LOG,MAXD,AA,D,IPART,LPART,XL,YL,TITLE)
      CLOSE(LOG)  
      RETURN
      END
C     
C     
C     
      SUBROUTINE PLOTTE(LOG,MAXD,AA,D,IPART,LPART,XL,YL,TITLE) 
      common /plotfl/plname
      DIMENSION D(2,1),AA(1),SKY(256)
      character plname*12,TITLE*36 
c     call setcc
c     call plfile(plname)
c     call plots(ierr,idummy,28)
c     if(ierr.ne.0)stop 'Error on opening call to Spooltek'
c     CALL FACTOR(2.54) 
      QQQ=(XL+2.0) 
      PPP=(YL+2.0) 
c     CALL FRAME(QQQ,PPP)  
c     CALL PLOT(1.0,0.75,-3)  
      YMAX=0
      DO 10 J=1,LPART 
        READ(LOG,45)JJJ,(D(2,KKK),KKK=1,256)
        DO 5 KKK=1,256
          IF(D(2,KKK).GT.YMAX) YMAX=D(2,KKK)
    5   CONTINUE
   10 CONTINUE    
      REWIND LOG  
      RANGE=D(1,1)-D(1,MAXD)  
      DH=YL/(LPART+1)   
      H=DH*2.0    
      DFDX=RANGE/XL     
      FREMIN=D(1,MAXD)  
c     CALL FAXIS(0.0,0.0,'FREQUENCY ( H )',-15,XL,FREMIN*1000.,   
c    *                                            DFDX*1000.)     
c     CALL SYMBOL(XL/2.0+4.0*0.14-0.07,-0.68,0.10,'M  Z',0.0,4)   
      DO 40 J=1,MAXD    
   40   SKY(J)=0    
      YH=-DH
      DO 100 JJ=1,LPART 
        READ(LOG,45)JJJ,(D(2,KKK),KKK=1,256)
   45   FORMAT(I3,(16E8.3))     
        IF(JJJ.NE.JJ)PRINT 46,JJJ,JJ  
   46   FORMAT(12X,43HPERIODOGRAMS DONT MATCH:CHECK PROGRAM JJJ= ,  
     *         I5,4H JJ=,I5)    
        YH=YH+DH    
        MJ=MOD(JJ,2)
        DO 60 J=1,MAXD    
          I=J*(1-2*MJ)+(MAXD+1)*MJ
          Y=H*D(2,I)/YMAX   
          Y=Y+YH
          X=(D(1,I)-D(1,MAXD))/DFDX     
          IF(Y.LT.SKY(I)) GO TO 50
          SKY(I)=Y    
c         IF(J.EQ.1)CALL PLOT(X,Y,3)    
c         IF(J.GT.1)CALL PLOT(X,Y,2)    
          GO TO 60    
   50     CONTINUE    
          X=(D(1,I+1-2*MJ)-D(1,MAXD))/DFDX    
          Y=H*D(2,I+1-2*MJ)/YMAX  
          Y=Y+YH
c         CALL PLOT (X,Y,3) 
   60   CONTINUE    
  100 CONTINUE    
c     CALL PLOT(0.0,0.0,3)    
c     CALL PLOT(0.0,YL,2)     
      H1=YL+0.88  
c     CALL SYMBOL(0.0,H1,.21,TITLE,0.0,36)     
c     CALL AXI(XL,DFDX,FREMIN,YL)   
c     CALL PLOT(XL,0.0,2)     
c     CALL PLOT(0.0,0.0,999) 
      RETURN
      END   
