C  PROGRAMME TROJANS
      PROGRAM NBODY
C  THIS IS A GENERAL N-BODY PROGRAMME WHERE INTER-BODY FORCES BETWEEN
C  BODIES i AND j ARE OF THE FORM CM(i)*CM(j)*F(Rij,Vij) WHERE F IS A 
C  FUNCTION OF Rij, THE DISTANCE BETWEEN THE BODIES, AND Vij, THE 
C  RELATIVE VELOCITIES OF THE TWO BODIES.  
C  THE STRUCTURE OF THE PROGRAMME IS;
C
C  (I)   THE MAIN PROGRAMME "NBODY"IS WHICH INCLUDES THE RUNGE-KUTTA 
C        ROUTINE WITH AUTOMATIC STEP CONTROL.
C  (II)  SUBROUTINE "START" WHICH ENABLES INPUT OF THE INITIAL BOUNDARY
C        CONDITIONS.
C  (III) SUBROUTINE "ACC" WHICH GIVES THE ACCELERATION OF EACH BODY
C        DUE TO ITS INTERACTIONS WITH ALL OTHER BODIES.
C  (IV)  SUBROUTINE "STORE" WHICH STORES INTERMEDIATE COORDINATES AND
C        VELOCITY COMPONENTS AS THE COMPUTATION PROGRESSES.
C  (V)   SUBROUTINE "OUT" WHICH OUTPUTS THE RESULTS TO DATA FILES.
C                           
C  BY CHANGING THE SUBROUTINES DIFFERENT PROBLEMS MAY BE SOLVED.   
C  THE CM'S CAN BE MASSES OR CHARGES OR BE MADE EQUAL TO UNITY WHILE 
C  THE FORCE LAW CAN BE INVERSE-SQUARE OR ANYTHING ELSE - 
C  e.g. LENNARD-JONES.  SEE COMMENT AT THE BEGINNING OF SUBROUTINE
C  "ACC" FOR THE TYPES OF FORCES OPERATING.
C   
C  THE FOUR-STEP RUNGE-KUTTA ALGORITHM IS USED. THE RESULTS OF TWO 
C  STEPS WITH TIMESTEP H ARE CHECKED AGAINST TAKING ONE STEP WITH 
C  TIMESTEP 2*H.  IF THE DIFFERENCE IS WITHIN THE TOLERANCE THEN THE 
C  TWO STEPS, EACH OF H, ARE ACCEPTED AND THE STEPLENGTH IS DOUBLED FOR
C  THE NEXT STEP.   HOWEVER, IF THE TOLERANCE IS NOT SATISFIED THEN THE
C  STEP IS NOT ACCEPTED AND ONE TRIES AGAIN WITH A HALVED STEPLENGTH.
C  IT IS ADVISABLE, BUT NOT ESSENTIAL, TO START WITH A REASONABLE 
C  STEPLENGTH; THE PROGRAMME QUICKLY FINDS A SUITABLE VALUE.  
C
C  AS PROVIDED THE PROGRAMME HANDLES UP TO 20 BODIES BUT THIS CAN BE 
C  CHANGED FROM 20 TO WHATEVER IS REQUIRED IN THE DIMENSION STATEMENT.
C
      DIMENSION CM(20),X(20,3),V(20,3),DX(20,3,0:4),DV(20,3,0:4),WT(4),
     +XTEMP(2,20,3),VTEMP(2,20,3),XT(20,3),VT(20,3),DELV(20,3)
      COMMON/A/X,V,TOL,H,TOTIME,DELV,XT,VT,NB,IST,TIME,IG,XTEMP,
     +VTEMP,CM
      DATA WT/0.0,0.5,0.5,1.0/
      IST=0
      OPEN(UNIT=9,FILE='LPT1')
C
C  SETTING THE INITIAL BOUNDARY CONDITION CAN BE DONE EITHER 
C  EXPLICITLY AS VALUES OF (X,Y,Z) AND (U,V,W) FOR EACH BODY OR
C  CAN BE COMPUTED.   THIS IS CONTROLLED BY SUBROUTINE "START".
C  OTHER PARAMETERS ARE ALSO SET IN "START" WHICH ALSO INDICATES
C  THE SYSTEM OF UNITS BEING USED.
C
      CALL START       
C
      TIME=0
C  INITIALIZE ARRAYS
      DO 57 I=1,20
      DO 57 J=1,3
      DO 57 K=0,4
        DX(I,J,K)=0
        DV(I,J,K)=0
   57 CONTINUE
C
C  WE NOW TAKE TWO STEPS WITH STEP LENGTH H FOLLOWED BY ONE STEP
C  WITH STEP LENGTH 2*H FROM THE SAME STARTING POINT BUT FIRST WE
C  STORE THE ORGINAL SPACE AND VELOCITY COORDINATES AND TIMESTEP.
C
   25 DO 7 IT=1,2
      DO 8 J=1,NB
      DO 8 K=1,3
        XTEMP(IT,J,K)=X(J,K)
        VTEMP(IT,J,K)=V(J,K)
    8 CONTINUE
      HTEMP=H
      DO 10 NOSTEP=1,3-IT      
      DO 11 I=1,4
      DO 12 J=1,NB
      DO 12 K=1,3
        XT(J,K)=XTEMP(IT,J,K)+WT(I)*DX(J,K,I-1)
   12 VT(J,K)=VTEMP(IT,J,K)+WT(I)*DV(J,K,I-1)
C
      CALL ACC
C
      DO 13 J=1,NB
      DO 13 K=1,3
        DV(J,K,I)=IT*HTEMP*DELV(J,K)
        DX(J,K,I)=IT*HTEMP*VT(J,K)
   13 CONTINUE
   11 CONTINUE
      DO 14 J=1,NB
      DO 14 K=1,3
      XTEMP(IT,J,K)=XTEMP(IT,J,K)+(DX(J,K,1)+DX(J,K,4)+2*
     +(DX(J,K,2)+DX(J,K,3)))/6.0 
      VTEMP(IT,J,K)=VTEMP(IT,J,K)+(DV(J,K,1)+DV(J,K,4)+2*
     +(DV(J,K,2)+DV(J,K,3)))/6.0
   14 CONTINUE
   10 CONTINUE  
    7 CONTINUE
C
C  THE ABOVE HAS MADE TWO STEPS OF H AND, FROM THE SAME STARTING POINT,
C  A SINGLE STEP OF 2*H.   THE RESULTS ARE NOW COMPARED
C
      DO 20 J=1,NB
      DO 20 K=1,3
        IF(ABS(XTEMP(1,J,K)-XTEMP(2,J,K)).GT.TOL)THEN
        H=0.5*H
        GOTO 25
        ENDIF                        
   20 CONTINUE
C
C  AT THIS STAGE THE DOUBLE STEP WITH H AGREES WITHIN TOLERANCE WITH
C  THE SINGLE STEP WITH 2*H.   THE TIMESTEP WILL NOW BE TRIED WITH
C  TWICE THE VALUE FOR THE NEXT STEP.   IF IT IS TOO BIG THEN IT WILL
C  BE REDUCED AGAIN. 
C
      H=2*H
      DO 80 J=1,NB
      DO 80 K=1,3
        X(J,K)=XTEMP(1,J,K)
        V(J,K)=VTEMP(1,J,K)
   80 CONTINUE
      TIME=TIME+H
C
      CALL STORE
C
      IF(TIME.GE.TOTIME)GOTO 50
      IF(IG.GT.1000)THEN
      IG=1000
      GOTO 50
      ENDIF
      GOTO 25
C
   50 CALL OUT
C
      STOP
      END


      SUBROUTINE STORE
      DIMENSION CM(20),X(20,3),V(20,3),XSTORE(1000,20,3),
     +XTEMP(2,20,3),VTEMP(2,20,3),DELV(20,3),XT(20,3),VT(20,3)
      COMMON/A/X,V,TOL,H,TOTIME,DELV,XT,VT,NB,IST,TIME,IG,XTEMP,
     +VTEMP,CM
      COMMON/B/NORIG,XSTORE
      DO 21 J=1,NB
      DO 21 K=1,3
        X(J,K)=XTEMP(1,J,K)
        V(J,K)=VTEMP(1,J,K)
   21 CONTINUE
C
C  UP TO 1000 POSITIONS ARE STORED.  THESE ARE TAKEN EVERY 50 STEPS. 
C
      IST=IST+1
      IF((IST/50)*50.NE.IST)GOTO 50
      IG=IST/50
      IF(IG.GT.1000)GOTO 50
      DO 22 J=1,NB
      DO 22 K=1,3 
        XSTORE(IG,J,K)=X(J,K)
   22 CONTINUE
   50 RETURN
      END


      SUBROUTINE START
      DIMENSION CM(20),X(20,3),V(20,3),XSTORE(1000,20,3),
     +XTEMP(2,20,3),VTEMP(2,20,3),DELV(20,3),XT(20,3),VT(20,3)
      COMMON/A/X,V,TOL,H,TOTIME,DELV,XT,VT,NB,IST,TIME,IG,XTEMP,
     +VTEMP,CM
      COMMON/B/NORIG,XSTORE
      OPEN(UNIT=21,FILE='LAST.DAT')
      OPEN(UNIT=22,FILE='FAST.DAT')
C
C  THE PROGRAMME AS PROVIDED IS FOR AN INVERSE-SQUARE LAW AND
C  USES UNITS FOR WHICH THE UNIT MASS IS THAT OF THE SUN, THE
C  THE UNIT OF DISTANCE IS THE ASTRONOMICAL UNIT (MEAN SUN-EARTH 
C  DISTANCE),THE UNIT OF TIME IS THE YEAR AND THE GRAVITATIONAL 
C  CONSTANT IS 4*PI**2.  
C
      WRITE(6,'('' INPUT THE NUMBER OF BODIES'')')
      READ(5,*)NB
      WRITE(6,'('' INPUT THE VALUES OF CM IN SOLAR-MASS UNITS. '')')
      WRITE(6,'('' THE TROJAN ASTEROID MASSES CAN BE PUT AS ZERO'')')
      DO 1 I=1,NB
        WRITE(6,500)I
  500   FORMAT(25H READ IN THE VALUE OF CM[,I3, 1H])    
        READ(5,*)CM(I)
    1 CONTINUE
      WRITE(6,'('' INPUT THE INITIAL TIMESTEP [years]'')')
      READ(5,*)H
      WRITE(6,'('' INPUT TOTAL TIME FOR THE SIMULATION [years]'')')
      READ(5,*)TOTIME
C
C  THE PROGRAMME ASKS THE USER TO SPECIFY A TOLERANCE, THE MAXIMUM 
C  ABSOLUTE ERROR THAT CAN BE TOLERATED IN ANY POSITIONAL COORDINATE
C  (X, Y OR Z).  IF THIS IS SET TOO LOW THEN THE PROGRAMME CAN BECOME
C  VERY SLOW.  FOR COMPUTATIONS INVOLVING PLANETS A TOLERANCE OF 1.0E-6
C  (c. 150 KM) IS USUALLY SATISFACTORY.  
C
      WRITE(6,'('' INPUT THE TOLERANCE '')')
      WRITE(6,'('' SEE COMMENT ABOVE THIS STATEMENT IN LISTING'')')
      READ(5,*)TOL
      WRITE(6,'('' THE CALCULATION CAN BE DONE RELATIVE TO AN '')')
      WRITE(6,'('' ARBITRARY ORIGIN OR WITH RESPECT TO ONE OF '')')
      WRITE(6,'('' THE BODIES AS ORIGIN.   INPUT ZERO FOR AN '')')
      WRITE(6,'('' ARBITRARY ORIGIN OR THE NUMBER OF THE BODY.'')')
      WRITE(6,'('' IF A BODY IS CHOSEN AS ORIGIN THEN ALL ITS'')')
      WRITE(6,'('' POSITIONAL AND VELOCITY VALUES ARE SET TO ZERO'')') 
      READ(5,*)NORIG
      DO 31 J=1,NB
        WRITE(6,100)J
  100   FORMAT(23H INPUT [X,Y,Z] FOR BODY,I3)      
        READ(5,*)X(J,1),X(J,2),X(J,3)
        WRITE(6,200)J
  200   FORMAT(32H INPUT [XDOT,YDOT,ZDOT] FOR BODY,I3)
        READ(5,*)V(J,1),V(J,2),V(J,3)
   31 CONTINUE
      RETURN
      END




      SUBROUTINE ACC
      DIMENSION CM(20),X(20,3),V(20,3),XSTORE(1000,20,3),R(3),
     +XTEMP(2,20,3),VTEMP(2,20,3),DELV(20,3),XT(20,3),VT(20,3),DD(3)
      COMMON/A/X,V,TOL,H,TOTIME,DELV,XT,VT,NB,IST,TIME,IG,XTEMP,
     +VTEMP,CM
      COMMON/B/NORIG,XSTORE
C
C  THE PROGRAMME AS PROVIDED IS FOR AN INVERSE-
C  SQUARE LAW AND USES UNITS FOR WHICH THE UNIT MASS IS THAT OF THE SUN,
C  THE UNIT OF DISTANCE IS THE ASTRONOMICAL UNIT (MEAN SUN-EARTH DISTANCE),
C  THE UNIT OF TIME IS THE YEAR AND THE GRAVITATIONAL CONSTANT IS 4*PI**2.  
C  HOWEVER, THE USER MAY MODIFY THE SUBROUTINE "ACC" TO CHANGE TO ANY OTHER
C  FORCE LAW AND/OR ANY OTHER SYSTEM OF UNITS.
C
C  SET THE VALUE OF G IN ASTRONOMICAL UNITS
      PI=4.0*ATAN(1.0)
      G=4*PI*PI
      DO 1 J=1,NB
      DO 1 K=1,3
        DELV(J,K)=0
    1 CONTINUE
C  THE FOLLOWING PAIR OF DO LOOPS FINDS INTERACTIONS FOR ALL PAIRS
C  OF BODIES
      DO 2 J=1,NB-1
      DO 2 L=J+1,NB
      DO 3 K=1,3
        R(K)=XT(J,K)-XT(L,K)
    3 CONTINUE
      RRR=(R(1)**2+R(2)**2+R(3)**2)**1.5
      DO 4 K=1,3
C  THE NEXT TWO STATEMENTS GIVE THE CONTRIBUTIONS TO THE THREE
C  COMPONENTS OF ACCELERATION DUE TO BODY J ON BODY L AND DUE TO
C  BODY L ON BODY J.  
        DELV(J,K)=DELV(J,K)-G*CM(L)*R(K)/RRR
        DELV(L,K)=DELV(L,K)+G*CM(J)*R(K)/RRR
    4 CONTINUE
    2 CONTINUE
C  IF ONE OF THE BODIES IS TO BE THE ORIGIN THEN ITS ACCELERATION IS
C  SUBTRACTED FROM THAT OF ALL OTHER BODIES.
      IF(NORIG.EQ.0)GOTO 10
      DD(1)=DELV(NORIG,1)
      DD(2)=DELV(NORIG,2)
      DD(3)=DELV(NORIG,3)
      DO 6 J=1,NB
      DO 6 K=1,3
        DELV(J,K)=DELV(J,K)-DD(K)
    6 CONTINUE
   10 RETURN
      END


      SUBROUTINE OUT
      DIMENSION CM(20),X(20,3),V(20,3),XSTORE(1000,20,3),
     +XTEMP(2,20,3),VTEMP(2,20,3),DELV(20,3),XT(20,3),VT(20,3)
      COMMON/A/X,V,TOL,H,TOTIME,DELV,XT,VT,NB,IST,TIME,IG,XTEMP,
     +VTEMP,CM
      COMMON/B/NORIG,XSTORE 
C     (X,Y) VALUES FOR THE LEADING ASTEROID ARE PLACED IN DATA FILE 
C     LAST.DAT AND FOR THE FOLLOWING ASTEROID IN DATA FILE FAST.DAT.
C
C  FOR THE TROJAN ASTEROID PROBLEM THE JUPITER RADIUS VECTOR IS
C  ROTATED TO PUT IT ON THE Y AXIS.   THE POSITIONS OF THE ASTEROIDS
C  RELATIVE TO JUPITER ARE PLOTTED.
      LAST=MIN(IG,1000)
      DO 60 I=1,LAST
C  THETA IS THE ANGLE BETWEEN THE X AXIS AND THE JUPITER RADIUS VECTOR
        THETA=ATAN2(XSTORE(I,2,2),XSTORE(I,2,1))
        XSTORE(I,2,2)=SQRT(XSTORE(I,2,2)**2+XSTORE(I,2,1)**2)
        XSTORE(I,2,1)=0
C  NOW THE ASTEROID RADIUS VECTORS ARE ROTATED BY PI/2-THETA
      DO 61 J=1,2
        AA=XSTORE(I,J+2,1)*SIN(THETA)-XSTORE(I,J+2,2)*COS(THETA)
        BB=XSTORE(I,J+2,1)*COS(THETA)+XSTORE(I,J+2,2)*SIN(THETA)
        XSTORE(I,J+2,1)=AA
        XSTORE(I,J+2,2)=BB
   61 CONTINUE
   60 CONTINUE
C  THE MODIFIED POSITIONS ARE NOW OUTPUT TO DATA FILES.
      DO 63 J=3,4
        N=18+J
        REWIND N
      DO 64 I=1,LAST
        WRITE(N,*)XSTORE(I,J,1),XSTORE(I,J,2)
   64 CONTINUE
   63 CONTINUE
      RETURN
      END            

