      SUBROUTINE BIN(Q,TQ,NQ) 
C     
C     ROUTINE TO BIN DATA.  THIS ROUTINE ASSUMES THE ARRAY TQ ALREADY   
C     CONTAINS TIMES IN UNITS OF PHASE AS CARRIED OUT BY 'TIMBAS' 
C     
      DIMENSION Q(NQ),BINS(50000),NBINS(50000)     
      DOUBLE PRECISION TQ(NQ)
      IDIV=0
      READ *,BIN1,BIN2,NB     
      IF(NB.LT.0)THEN   
        IDIV=1    
        NB=-NB    
      END IF
      PRINT 10,BIN1,BIN2,NB   
   10 FORMAT(12X,'DATA TO BE BINNED INTO THE PHASE INTERVAL'/     
     *       12X,F7.4,' TO ',F7.4,' USING ',I4,' BINS')     
      IF(IDIV.GT.0)PRINT 12   
   12 FORMAT(12X,'DATA WILL NOT BE NORMALIZED BY NUMBER IN BIN')  
      IF(NB.GT.50000)THEN
        PRINT 15  
   15   FORMAT(12X,'TOO MANY BINS REQUESTED - MAX IS 50000') 
        STOP
      END IF
      DO 20 J=1,50000    
        BINS(J)=0.0     
        NBINS(J)=0.0    
   20 CONTINUE    
      DO 30 J=1,NQ
        IPHASE=TQ(J)    
        PHASE=TQ(J)-IPHASE    
   25   IF(PHASE.LT.0.0)THEN  
          PHASE=PHASE+1.0     
          GO TO 25
        END IF    
        IF(BIN1.GE.0.0)THEN   
          IF(BIN2.LT.1.0)THEN
            IF(PHASE.LT.BIN1.OR.PHASE.GT.BIN2)GO TO 30  
            IB=(PHASE-BIN1)/(BIN2-BIN1)*NB+0.99999
          END IF
          IF(BIN2.GE.1.0)THEN
            IF(PHASE.GT.BIN2-1.0.AND.PHASE.LT.BIN1)GO TO 30
            IF(PHASE.LT.BIN1)PHASE=PHASE+1.0
            IB=(PHASE-BIN1)/(BIN2-BIN1)*NB+0.99999
          END IF
        ELSE
          IF(PHASE-1.0.LT.BIN1.AND.PHASE.GT.BIN2)GO TO 30   
          IF(PHASE.LE.BIN2)IB=(PHASE-BIN1)/(BIN2-BIN1)*NB+0.99999 
          IF(PHASE-1.0.GE.BIN1)IB=(PHASE-1.0-BIN1)/(BIN2-BIN1)*   
     *                                             NB+0.99999     
        END IF    
        IB=MAX0(IB,1)   
        IB=MIN0(IB,NB)  
        BINS(IB)=BINS(IB)+Q(J)
        NBINS(IB)=NBINS(IB)+1 
   30 CONTINUE    
      NQ=0  
      DO 50 J=1,NB
        Q(J)=BINS(J)    
        TQ(J)=(J-0.5)/NB*(BIN2-BIN1) + BIN1     
        IF(IDIV.EQ.0.AND.NBINS(J).GT.0)Q(J)=Q(J)/NBINS(J)   
        IF(ABS(ABS(BIN2-BIN1)-1.0).LT.0.0001)THEN     
          Q(J+NB)=Q(J)  
          TQ(J+NB)=TQ(J)+1.0  
          NQ=NB   
        END IF    
   50 CONTINUE    
      PRINT 60,(NBINS(J),J=1,NB)    
   60 FORMAT(12X,20I6)  
      NQ=NQ+NB    
      RETURN
      END   
