V
varunkumarchundayil
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
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