C-------------------------------------------------------------------
C
C   I  N  T  E  G 
C 
C   INTEGRATE EQUATIONS OF MOTION FROM  EPOCH TO D2.
C   THE INTEGRATION STEP IS H DAYS; 
C   D1 AND D2 GIVES THE INTERVAL FOR WHICH THE RESULTS ARE WRITTEN OUT.
C 
C----------------------------------------------------------------------
      SUBROUTINE  INTEG (D1,H,D2,OUTP, INTEGS)
      INTEGER
     .    D1, 
     .    D2, 
     .    INTEGS 
      REAL
     .    H,     
     .    OUTP   
C
      REAL  PI,PI2,CPI,DTORAD,RTODEG,
     .                  GRAVIT,OBLIQ
      COMMON /CONST/    PI,PI2,CPI,DTORAD,RTODEG,
     .                  GRAVIT,OBLIQ
C 
C 
C
      INTEGER MAXOBJ, MAXDAY, MMM, ORDER
      PARAMETER (MAXOBJ=10, MAXDAY=10000, ORDER=7, MMM=3)

      real tmppow(maxobj,maxobj)
      REAL
     .    EPOCH,
     .    SBINC(0:ORDER,0:ORDER), 
     .    DD(0:ORDER),
     .    KK(-MMM:MMM,-MMM:MMM),
     .    MASS(MAXOBJ),
     .    COORD(3,MAXOBJ),
     .    VELOC(3,MAXOBJ),
     .    XDIFF(3,MAXOBJ),       
     .    ACCEL(-MMM:MMM,3,MAXOBJ,0:1),
     .    ADIFF(MAXOBJ,0:ORDER,3),  
     .    DISTNC(MAXOBJ,MAXOBJ)        
      INTEGER SWAP
C
      COMMON /INTTAB/    SBINC,DD,KK,
     .                   COORD,VELOC,MASS,
     .                   XDIFF,ACCEL,ADIFF,DISTNC,SWAP,
     .                   EPOCH
C
C  OUTPUT
C
      REAL
     .    PLANET(3,MAXOBJ),      
     .    OLDPL(2,MAXOBJ)    
      INTEGER
     .    TODAY
      REAL
     .    OUTBUF(MAXDAY,3,MAXOBJ)      
C
      COMMON /OUT/ OUTBUF, PLANET, OLDPL, TODAY
C
C
      INTEGER
     .    I,J,K,NU,P,COUNT,INDX
      REAL
     .    SUM,X0,X1,X2,
     .    CLOCK,ETIME,
     .    OTIME
C
      TODAY=1
      P=MMM
      COUNT=0 
      ETIME=D2+H/4.0E+0
      OTIME=D1
      CLOCK=EPOCH+(MMM+1)*H
C
C  INTEGRATION LOOP BEGINS HERE
C
 10   CONTINUE
         DO 23 NU=0,ORDER-1
            DO 22 K=1,3
               DO 21 I=1,MAXOBJ
                        adiff(i,nu,k) = 0.0
                  DO 20 J=0,NU
                     IF (P-J .GE. -MMM) THEN
                        INDX= P-J
                     ELSE
                        INDX= P-J+ORDER
                     ENDIF 
                     adiff(i,nu,k)=
     +                    adiff(i,nu,k)+sbinc(nu,j)*accel(indx,k,i,swap)
 20               CONTINUE
 21            CONTINUE
 22         CONTINUE
 23      CONTINUE
C
C
         DO 32 K=1,3
            DO 31 I=1,MAXOBJ
               SUM=0.0
               DO 30 J=0,ORDER-1
                  SUM=SUM+DD(J)*ADIFF(i,J,K) 
 30            CONTINUE
               XDIFF(K,I)=XDIFF(K,I)+H*H*SUM
               COORD(K,I)=COORD(K,I)+XDIFF(K,I)
 31         CONTINUE
 32      CONTINUE   
C
C  CHECK IF IT IS TIME TO STORE RESULTS
C
         IF (CLOCK .GE. OTIME) THEN
            CALL DOUT (CLOCK)
            OTIME=OTIME+OUTP  
         ENDIF  
C
         DO 41 I=1,MAXOBJ
            DO 40 J=I+1,MAXOBJ
               X0=COORD(1,I)-COORD(1,J)
               X1=COORD(2,I)-COORD(2,J)
               X2=COORD(3,I)-COORD(3,J)
               DISTNC(I,J)=SQRT(X0*X0+X1*X1+X2*X2)
               DISTNC(J,I)=DISTNC(I,J)
 40         CONTINUE
 41      CONTINUE
C
         do i = 1, maxobj
            do j = i, maxobj
               tmppow(i,j) = distnc(i,j) ** 3
               tmppow(j,i) = tmppow(i,j)
            end do
         end do
C
         IF (P .EQ. MMM) THEN
            P= -MMM
         ELSE
            P= P+1
         ENDIF 
         COUNT=COUNT+1 
C
C 
         DO 52 K=1,3
            DO 51 I=1,MAXOBJ
               SUM=0.0 
               DO 50 J=1,MAXOBJ
                  IF (J .NE. I .AND. MASS(J) .NE. 0.0) THEN
                      tmp = 1.0/tmppow(i,j)
                      SUM=SUM+tmp*(COORD(K,J)-COORD(K,I))* MASS(J)
                  ENDIF
 50            CONTINUE
               ACCEL(P,K,I,SWAP)=GRAVIT*SUM
 51         CONTINUE
 52      CONTINUE
         CLOCK=CLOCK+H
      IF (CLOCK .LE. ETIME) GOTO 10
C
C  READY
C
      INTEGS=COUNT 
      RETURN
      END 
