      PROGRAM nbody                                                     00000100
C                                                                       00000200
C                   PROGRAM IN FORTRAN 77                               00000300
C                                                                       00000400
C     ***  HEIDELBERG N-BODY PROGRAM  ***                               00000500
C                            BY J.SCHUBART AND P.STUMPFF                00000600
C                                                                       00000700
C                                                                       00000900
C                                                                       00001000
      DOUBLE PRECISION EM(30),X(30,3),XPUNKT(30,3),XNABLA(30,3), BESCHL 00001100
     1(30,3),A(30,3,12),B(30,3,12),CAI(5,11),H,T,DELTAT,AA,AB,D(30,30), 00001200
     *DIFJUL(30),DISJUL(30),HD,S,GG(11),                                00001300
     *H1,H2,H3,SEVENS                                                   00001400
      DIMENSION IPOUT(30),DIFMAX(30),KOODIZ(30),                        00001500
     *DISMIN(30),KPDIS(30),ILEM(1000)                                   00001600
      COMMON /ABLABL/A,B                                                00001800
      COMMON CAI,EM,X,XPUNKT,BESCHL,XNABLA,H,D                          00001900
C                                                                       00002000
C       INPUT, OUTPUT UNITS     (dat6 for formatted output,             00002100
C                                dat7 for unformatted output)           00002200
C An example of dat5, the input file, is separately available and       00002300
C                                         shows instructions at the end 00002400
C The symbols of the input values are introduced below                  00002500
C                                                                       00002600
      OPEN (5,FILE='dat5',FORM='FORMATTED')
      OPEN (6,FILE='dat6',FORM='FORMATTED')
      OPEN (7,FILE='dat7',FORM='UNFORMATTED')
C                                                                       00002700
 1000 FORMAT(1H )                                                       00002800
 1001 FORMAT(I1,1X,I4)                                                  00002900
 1002 FORMAT(1X,I1,1X,I4)                                               00003000
 1003 FORMAT(I1,1X,I2,I11,2I10)                                         00003100
 1004 FORMAT(1X,I1,1X,I2,I11,2I10)                                      00003200
 1005 FORMAT(I1,1X,I2,3X,1PD23.15,2D15.7)                               00003300
 1006 FORMAT(1X,I1,1X,I2,3X,1PD23.15,1P,2D15.7)                         00003400
 1007 FORMAT(I1,1X,I2,1X,I2,1PD23.15,D30.15)                            00003500
 1008 FORMAT(1X,I1,1X,I2,1X,I2,1PD23.15,1PD30.15,10X,0PF10.1)           00003600
 1009 FORMAT(I1,1X,I2,1X,I2,1PD23.15)                                   00003700
 1010 FORMAT(1X,I1,1X,I2,1X,I2,1PD23.15)                                00003800
 1011 FORMAT(34H  COORDINATES OR MASSES INCOMPLETE)                     00003900
 1012 FORMAT(17H  INITIAL VALUES.)                                      00004000
 1013 FORMAT(3H N=,I6,4X,2HM=,I6,6H -EXP=,I6,' ZA=',I9,' ZN=',I9,       00004100
     *' DELTAZ=',I9,                                                    00004200
     *           6H   KD=,I6,6H   ED=,I6,6H   EG=,I6)                   00004300
 1014 FORMAT(6X,3HH= ,1PD23.15,6X,3HT= ,1PD23.15,6X,9H DELTAT= ,1PD23.1500004400
     *)                                                                 00004500
 1015 FORMAT(14H  BODY NUMBER ,I2,6H  EM= ,1PD23.15)                    00004600
 1018 FORMAT(28H  NO EXECUTION. INPUT WRONG.)                           00004900
 1019 FORMAT(20H    READY       END )                                   00005000
      GG(2)=0.41666666666666667D-1                                      00005900
      GG(1)=0.16666666666666667D0                                       00006000
      GG(3)=0.22222222222222222D-1                                      00006100
      GG(4)=0.14583333333333333D-1                                      00006200
      GG(5)=0.10615079365079365D-1                                      00006300
      GG(6)=0.82258597883597884D-2                                      00006400
      GG(7)=0.66479276895943563D-2                                      00006500
      GG(8)=0.55372299382716049D-2                                      00006600
      GG(9)=0.47180677475816365D-2                                      00006700
      GG(10)=.40919706740019240D-2                                      00006800
      GG(11)=.35997549720267974D-2                                      00006900
      SEVENS = 7.7777777D74                                             00007000
C
C            The following input values are preliminary.
      N = 2                                                             00007400
C            N: number of bodies
      M=5                                                               00007500
      IEXP=14                                                           00007900
C     recommended values for order of integration formula
C                            and accuracy of starting iteration
      IZA=5                                                             00008000
      IZN=100                                                           00008100
      IDELTZ=3                                                          00008200
C     integration stops after step 100 in this example,
C     output (if any) starts at step 5 and proceeds at an interval of 3 steps
C     IZA must be greater than or equal to M
      KD=0                                                              00008400
C     KD= 1:  output relative to body No 1,
C             otherwise barycentric output
      IED=1                                                             00008500
C        =1: output of coord. and velocities, otherwise no velocities
C                                             by subroutine  DRUCKE
      IEG=0                                                             00008600
C     IEG = 0 , BARYCENTRIC INPUT ***  ELSE INPUT RELATIVE TO BODY 1
C             Note: The program stores barycentric values at final step
      H=3.141592653589793/60.D0                                         00008700
C        step length of integration (in a special unit)
      T=0.0D0                                                           00008800
C        time at start
      DELTAT=3.0D0                                                      00008900
C        interval of time that corresponds to H
      DO 1 I=1,30                                                       00009400
      EM(I)=SEVENS                                                      00009600
      IPOUT(I)= 0                                                       00009700
C      I: Body No          IPOUT(I)=   0        1           2
C      EM: mass              output:  no     to dat6     to dat7
      DO 1 J=1,3                                                        00009800
      X(I,J)=SEVENS                                                     00009900
    1 XPUNKT(I,J)=SEVENS                                                00010000
C      coordinates and velocities     x xdot y ydot z zdot
C                      omit values of body No 1
      GO TO 2                                                           00010100
C
C     EINLESEN  UND WIEDERGABE, HAUPTPROGRAMM.                          00010200
 2000 WRITE(6,1018)                                                     00010300
      GO TO 9000                                                        00010400
C                preliminary or final value of N7 :
    2 N7=0                                                              00010500
C = number of steps IZ with call of subroutine ELEMNT during integration
      WRITE(6,1000)                                                     00011100
    3 READ(5,1001,END=6) I,J                                            00011200
C             values I= 1, 2, 3, 5, 7 are permitted, 8 starts execution.
C      READ(5...) : An example of an input file dat5 with instructions
C                                               is separately available.
      IF(I.GT.0) GO TO 5                                                00011300
    4 READ(5,1001,END=6) I                                              00011400
      IF(I-8)4,2000,4                                                   00011500
    5 WRITE(6,1002)I,J                                                  00011600
C
      GO TO(7,15,17,3,23,3,29,33,3),I                                   00011700
C              go to group I or execute
    6 WRITE(6,1019)                                                     00011900
      GO TO 9000                                                        00012000
C
    7 IF(J.LE.0) GO TO 3                                                00012100
    8 READ (5,1003)J1,J2,J3,J4,J5                                       00012200
      IF(J1.NE.0) GO TO 4                                               00012300
    9 WRITE(6,1004)J1,J2,J3,J4,J5                                       00012400
      GO TO (10,13,14),J2                                               00012500
   10 N=J3                                                              00012600
      M=J4                                                              00012700
      IEXP=J5                                                           00012800
   11 IF((N.LT.2).OR.(N.GT.30)) GO TO 4                                 00012900
   12 J=J-1                                                             00013000
      GO TO 7                                                           00013100
   13 IZA=J3                                                            00013200
      IZN=J4                                                            00013300
      IDELTZ=J5                                                         00013400
      J=J-1                                                             00013500
      GO TO 7                                                           00013600
   14 KD=J3                                                             00013700
      IED=J4                                                            00013800
      IEG=J5                                                            00013900
      J=J-1                                                             00014000
      GO TO 7                                                           00014100
C
   15 READ(5,1005) J1,J2,H1,H2,H3                                       00017500
      IF(J1.NE.0) GO TO4                                                00017600
      WRITE(6,1006) J1,J2,H1,H2,H3                                      00017700
      H = H1                                                            00017900
      T = H2                                                            00018000
      DELTAT = H3                                                       00018100
      GO TO 3                                                           00018300
C
   17 J=J-1                                                             00018700
      IF(J.LT.0) GO TO 3                                                00018800
   18 READ(5,1009)J1,J2,J3,AA                                           00018900
      IF(J1.NE.0) GO TO 4                                               00019000
   19 WRITE(6,1010)J1,J2,J3,AA                                          00019600
      EM(J2)=AA                                                         00019700
      IPOUT(J2)=J3                                                      00019800
      GO TO 17                                                          00020000
C
   23 J = J-1                                                           00022900
      IF(J.LT.0) GO TO 3                                                00023000
   24 READ( 5,1007)J1,J2,J3,AA,AB                                       00023100
      IF(J1.NE.0) GO TO 4                                               00023200
   25 WRITE(6,1008)J1,J2,J3,AA,AB                                       00023300
      X(J2,J3)=AA                                                       00023400
      XPUNKT(J2,J3)=AB                                                  00023500
      GO TO 23                                                          00023600
C
   29 N7=J                                                              00028700
   30 J=J-1                                                             00028800
      IF(J.LT.0) GO TO 3                                                00028900
   31 READ(5,1003)J1,J2,J3                                              00029000
      IF(J1.NE.0) GO TO 4                                               00029100
   32 WRITE(6,1004) J1,J2,J3                                            00029200
      ILEM(J+1)=J3                                                      00029300
      GO TO 30                                                          00029400
C                       Start of Execution
   33 CONTINUE                                                          00029600
      IF(IZA.GE.M) GO TO 35                                             00031800
   34 IZA=M                                                             00031900
      IF(IDELTZ.GE.M) IZA = IDELTZ
   35 DO 36 I=1,N                                                       00032000
      IF(EM(I).EQ.SEVENS) GO TO 39                                      00032100
   36 CONTINUE                                                          00032200
      DO 38 J=1,3                                                       00032300
      DO 38 I=2,N                                                       00032400
      IF(X(I,J).EQ.SEVENS) GO TO 39                                     00032500
   37 IF(XPUNKT(I,J).EQ.SEVENS) GO TO 39                                00032600
   38 CONTINUE                                                          00032700
      GO TO 40                                                          00033600
   39 WRITE(6,1011)                                                     00034000
      GO TO 2                                                           00034100
   40 WRITE(6,1012)                                                     00034200
      WRITE(6,1013)N,M,IEXP,IZA,IZN,IDELTZ,KD,IED,IEG                   00034300
      WRITE(6,1014)H,T,DELTAT                                           00034400
      IZ = 0                                                            00034500
      IYYL = 0                                                          00034600
      DO 42 I = 1,N                                                     00034700
      WRITE(6,1015) I,EM(I)                                             00034800
      DO 41 J = 1,3                                                     00034900
                     WRITE(6,1008) IYYL,I,J,X(I,J),XPUNKT(I,J),T        00035000
   41 CONTINUE                                                          00035200
   42 CONTINUE                                                          00035300
C                                                                       00035400
   44 WRITE(6,1000)                                                     00035500
      DO 440 J=2,N                                                      00035600
      DIFMAX(J)=0.0                                                     00035700
      DIFJUL(J)=0.0                                                     00035800
      KOODIZ(J)=0                                                       00035900
      DISMIN(J)=1.0E+35                                                 00036000
      DISJUL(J)=0.0                                                     00036100
  440 KPDIS(J)=0                                                        00036200
      DO 441 J=1,N                                                      00036300
      DO 441 L=1,N                                                      00036400
  441 D(J,L)=1.0D35                                                     00036500
      CALL KOEFFZ(M)                                                    00036600
      IZ=M                                                              00036700
      IS=1
      CALL ANFITN(N,M,IEXP,IEG)                                         00036800
      DO 45 I=1,M                                                       00037000
   45 T=T+DELTAT                                                        00037300
      GO TO 47                                                          00037400
   46 CALL SCHRTT(N,M)                                                  00037500
      T=T+DELTAT                                                        00037600
      IZ=IZ+1                                                           00037700
      MPMP2=2*M+2                                                       00038500
C                                                                       00038600
      DO 466 J=2,N                                                      00038700
      DO 461 K=1,3                                                      00038800
      ABB=DABS(B(J,K,MPMP2))                                            00038900
      IF(ABB.LE.DIFMAX(J)) GO TO 461                                    00039000
  460 DIFMAX(J)=ABB                                                     00039100
      DIFJUL(J)=T                                                       00039200
      KOODIZ(J)=K                                                       00039300
  461 CONTINUE                                                          00039400
      LE=J-1                                                            00039500
      DO 463 L=1,LE                                                     00039600
      IF(D(J,L).GE.DISMIN(J)) GO TO 463                                 00039700
  462 DISMIN(J)=D(J,L)                                                  00039800
      DISJUL(J)=T                                                       00039900
      KPDIS(J)=100*L+J                                                  00040000
  463 CONTINUE                                                          00040100
      IF(J.GE.N) GO TO 47                                               00040200
  464 LA=J+1                                                            00040300
      DO 466 L=LA,N                                                     00040400
      IF (D(L,J).GE.DISMIN(J)) GO TO 466                                00040500
  465 DISMIN(J)=D(L,J)                                                  00040600
      DISJUL(J)=T                                                       00040700
      KPDIS(J)=100*L+J                                                  00040800
  466 CONTINUE                                                          00040900
   47 IF(IZ-IZA)51,50,48                                                00043000
   48 IF(IS-IDELTZ)49,50,50                                             00043100
   49 IS=IS+1                                                           00043200
      GO TO 51                                                          00043300
   50 IS=1                                                              00043400
C
      IF(IED.NE.1) GOTO 474
      HD=1.D0/H
      M1=2*M+1
      DO 473 K=1,3
      XPUNKT(1,K)=0.0D0
      DO 472 J=2,N
      S=0.0D0
      DO 471 L=1,M1
  471 S=S+GG(L)*B(J,K,L+1)
      XPUNKT(J,K)=HD*XNABLA(J,K)+H*(B(J,K,1)*0.5D0-S)
  472 XPUNKT(1,K)=XPUNKT(1,K)+EM(J)*XPUNKT(J,K)
  473 XPUNKT(1,K)=-XPUNKT(1,K)/EM(1)
  474 CALL DRUCKE(N,M,IZ,T,KD,IED,IPOUT)                                00043500
C
   51 IF(IZ.EQ.IZN) GO TO 53                                            00043800
      IF(N7.LE.0) GO TO 54                                              00043900
   52 IF(IZ.NE.ILEM(N7)) GO TO 54                                       00044000
   53 CALL ELEMNT(N,M,IZ,T,KD)                                          00044100
      WRITE(6,1530)                                                     00044300
      DO 530 J=2,N                                                      00044400
  530 WRITE(6,1531)DIFMAX(J),DIFJUL(J),KOODIZ(J),DISMIN(J),DISJUL(J),   00044500
     *KPDIS(J)                                                          00044600
 1530 FORMAT(21H0   DIFMAX AND DISMIN)                                  00044700
 1531 FORMAT(1H ,1PE15.6,1X,0PF10.1,1X,I1,5X,1PE15.6,1X,0PF10.1,1X,I4)  00044800
      DO 531 J=2,N                                                      00044900
      DIFMAX(J)=0.0                                                     00045000
      DIFJUL(J)=0.0                                                     00045100
      KOODIZ(J)=0                                                       00045200
      DISMIN(J)=1.0E+35                                                 00045300
      DISJUL(J)=0.0                                                     00045400
  531 KPDIS(J)=0                                                        00045500
      N7=N7-1                                                           00045600
C                                                                       00045700
   54 IF(IZ.LT.IZN)  GOTO 46                                            00048700
      GOTO 2                                                            00053900
C
 9000 CONTINUE                                                          00054000
      STOP                                                              00054100
      END                                                               00054200
C
      SUBROUTINE KOEFFZ(M)                                              00095800
      DOUBLE PRECISION CAI(5,11),C(10,5),R(11,11),FAC(11),S,SIG(11),S1,T00095900
     1,T1                                                               00096000
      COMMON /ABLABL/A,B                                                00096100
      DOUBLE PRECISION A(30,3,12),B(30,3,12)                            00096200
      COMMON CAI                                                        00096300
      S=1.D0                                                            00096400
      FAC(1)=1.D0                                                       00096500
      DO 1 I=1,10                                                       00096600
      FAC(I+1)=FAC(I)*S                                                 00096700
    1 S=S+1.D0                                                          00096800
      SIG(1)=1.D0                                                       00096900
      DO 2 I=1,10                                                       00097000
    2 SIG(I+1)=-(SIG(I))                                                00097100
      S=1.D0                                                            00097200
      DO 8 L=1,M                                                        00097300
      L1=2*L                                                            00097400
      L2=L-1                                                            00097500
      S1=S**2                                                           00097600
      DO 7 K=1,L1                                                       00097700
      C(K,L)=0.D0                                                       00097800
      K2=K-2                                                            00097900
      IF(K2) 4,4,3                                                      00098000
    3 C(K,L)=C(K,L)+C(K2,L2)                                            00098100
    4 IF(K+1-L1) 5,6,7                                                  00098200
    5 C(K,L)=C(K,L)-S1*C(K,L2)                                          00098300
      GOTO 7                                                            00098400
    6 C(K,L)=C(K,L)-S1                                                  00098500
    7 CONTINUE                                                          00098600
    8 S=S+1.D0                                                          00098700
      L1=2*M+1                                                          00098800
      L2=L1-1                                                           00098900
      S=1.D0 - S                                                        00099000
      DO 10 L=1,L1                                                      00099100
      R(L1,L)=1.D0                                                      00099200
      DO 9  K=1,L2                                                      00099300
      K2=L1-K                                                           00099400
    9 R(K2,L)=C(K2,M)+S*R(K2+1,L)                                       00099500
   10 S=S+1.D0                                                          00099600
      S=1.D0                                                            00099700
      DO 13 I=1,M                                                       00099800
      DO 12 L=1,L1                                                      00099900
      L2=L1-L+1                                                         00100000
      S1=S**2                                                           00100100
      T=1.D0                                                            00100200
      CAI(I,L)=0.D0                                                     00100300
      DO 11 K=1,L1                                                      00100400
      T1=R(K,L)*S1/T                                                    00100500
      T=T+1.D0                                                          00100600
      CAI(I,L)=CAI(I,L)+T1/T                                            00100700
   11 S1=S1*S                                                           00100800
   12 CAI(I,L)=CAI(I,L)*SIG(L2)/(FAC(L)*FAC(L2))                        00100900
   13 S=S+1.D0                                                          00101000
      RETURN                                                            00101100
      END                                                               00101200
      SUBROUTINE ANFITN(N,M,IEXP,IEG)                                   00101300
C                                                                       00101400
C     ANFANGSITERATION - STARTING ITERATION                             00101500
C                                                                       00101600
      DOUBLE PRECISION CAI(5,11),EM(30),X(30,3),XPUNKT(30,3),XNABLA(30,300101700
     *),BESCHL(30,3),A(30,3,12),B(30,3,12),H,HQ,FF,FI,EMMA,CAID,H1,H2   00101800
      DIMENSION EMS(30)                                                 00101900
      COMMON /ABLABL/A,B                                                00102000
      COMMON CAI,EM,X,XPUNKT,BESCHL,XNABLA,H                            00102100
C     IEG = 0 , BARYCENTRIC INPUT ***  ELSE INPUT RELATIVE              00102200
C                                                                       00102300
      IF(IEG.EQ.0) GO TO 105                                            00102400
  100 EMMA=0.0D0                                                        00102500
      DO 101 J=1,N                                                      00102600
  101 EMMA=EMMA+EM(J)                                                   00102700
      DO 104 K=1,3                                                      00102800
      X(1,K)=0.0D0                                                      00102900
      XPUNKT(1,K)=0.0D0                                                 00103000
      DO 102 J=2,N                                                      00103100
      X(1,K)=X(1,K)+EM(J)*X(J,K)                                        00103200
  102 XPUNKT(1,K)=XPUNKT(1,K)+EM(J)*XPUNKT(J,K)                         00103300
      X(1,K)=-X(1,K)/EMMA                                               00103400
      XPUNKT(1,K)=-XPUNKT(1,K)/EMMA                                     00103500
      DO 103 J=2,N                                                      00103600
      X(J,K)=X(J,K)+X(1,K)                                              00103700
  103 XPUNKT(J,K)=XPUNKT(J,K)+XPUNKT(1,K)                               00103800
  104 CONTINUE                                                          00103900
  105 M1=M+1                                                            00104000
      CALL WW(N)                                                        00104100
      DO 1 I=1,3                                                        00104200
      DO 1 J=2,N                                                        00104300
      A(J,I,12)=X(J,I)                                                  00104400
    1 A(J,I,M1)=BESCHL(J,I)                                             00104500
      DELTA=10.**(-IEXP)                                                00104600
      DO 3I=2,N                                                         00104700
      EMS(I)=DELTA*(DABS(BESCHL(I,1))+DABS(BESCHL(I,2))+DABS(BESCHL(I,3)00104800
     *))                                                                00104900
      IF(EMS(I).GE.1.0E-28) GO TO 3                                     00105000
    2 EMS(I)=1.0E-28                                                    00105100
    3 CONTINUE                                                          00105200
      MU=0                                                              00105300
      M2=M+M1                                                           00105400
      M3=1+M2                                                           00105500
      M4=M+M                                                            00105600
      HQ=H*H                                                            00105700
      FF=0.0D0                                                          00105800
      DO 31 MM=1,M                                                      00105900
   31 FF=FF-1.0D0                                                       00106000
    4 FI=FF                                                             00106100
      IT=0                                                              00106200
      DO 20 I=1,M2                                                      00106300
      I1=I-M1                                                           00106400
      I2=I1                                                             00106500
      IF(I1) 5,20,6                                                     00106600
    5 I1=-I1                                                            00106700
    6 H1=FI*H                                                           00106800
      IF(MU.GT.0) GO TO 8                                               00106900
    7 H2=0.5D0*H1*H1                                                    00107000
    8 DO 13 J=2,N                                                       00107100
      DO 13 K=1,3                                                       00107200
      IF(MU.EQ.0) GO TO 12                                              00107300
    9 X(J,K)=0.D0                                                       00107400
      DO 11 L=1,M2                                                      00107500
      L1=L                                                              00107600
      IF(I2.GT.0) GO TO 11                                              00107700
   10 L1=M3-L1                                                          00107800
   11 X(J,K)=X(J,K)+CAI(I1,L1)*B(J,K,L)                                 00107900
      X(J,K)=HQ*X(J,K)                                                  00108000
      GO TO 13                                                          00108100
   12 X(J,K)=H2*A(J,K,M1)                                               00108200
   13 X(J,K)=A(J,K,12)+H1*XPUNKT(J,K)+X(J,K)                            00108300
      CALL WW(N)                                                        00108400
      IF(MU.EQ.0) GO TO 18                                              00108500
   14 DO 17 J=2,N                                                       00108600
      IF(IT.GT.0) GO TO 18                                              00108700
   15 DEL=DABS(BESCHL(J,1)-B(J,1,I))+DABS(BESCHL(J,2)-B(J,2,I))+DABS    00108800
     *(BESCHL(J,3)-B(J,3,I))                                            00108900
      IF(DEL.LE.EMS(J)) GO TO 17                                        00109000
   16 IT=1                                                              00109100
   17 CONTINUE                                                          00109200
   18 DO 19 J=2,N                                                       00109300
      DO 19 K=1,3                                                       00109400
   19 A(J,K,I)=BESCHL(J,K)                                              00109500
   20 FI=FI+1.0D0                                                       00109600
      IF(IT.EQ.0) GO TO 23                                              00109700
   21 DO 22 J=2,N                                                       00109800
      DO 22 K=1,3                                                       00109900
      DO 22 L=1,M2                                                      00110000
   22 B(J,K,L)=A(J,K,L)                                                 00110100
      MU=MU+1                                                           00110200
      GO TO 4                                                           00110300
   23 IF(MU.EQ.0) GO TO 21                                              00110400
   24 WRITE(6,25) MU                                                    00110500
   25 FORMAT(40H STARTING ITERATION. NO. OF ITERATIONS =,I3)            00110600
      I2=M-1                                                            00110700
      DO 30 J=2,N                                                       00110800
      DO 30 K=1,3                                                       00110900
      XNABLA(J,K)=0.0D0                                                 00111000
      DO 27 L=1,M2                                                      00111100
      CAID=CAI(M,L)                                                     00111200
      IF(I2.EQ.0) GO TO 27                                              00111300
   26 CAID=CAID-CAI(I2,L)                                               00111400
   27 XNABLA(J,K)=XNABLA(J,K)+CAID*A(J,K,L)                             00111500
      XNABLA(J,K)=H*(XPUNKT(J,K)+H*XNABLA(J,K))                         00111600
      DO 28 I=1,M4                                                      00111700
      L1=M2-I                                                           00111800
      DO 28 L=1,L1                                                      00111900
   28 A(J,K,L)=A(J,K,L+1)-A(J,K,L)                                      00112000
      DO 29 L=1,M2                                                      00112100
      L1=M3-L                                                           00112200
   29 B(J,K,L)=A(J,K,L1)                                                00112300
   30 B(J,K,M3)=0.0D0                                                   00112400
      RETURN                                                            00112500
      END                                                               00112600
      SUBROUTINE WW(N)                                                  00112700
C     DIFFERENTIAL EQUATIONS                                            00112800
      DOUBLE PRECISION CAI(5,11),EM(30),X(30,3),XPUNKT(30,3),BESCHL     00112900
     *(30,3),D(30,30),DX(3),PSEUD1(30,3),PSEUD2(30,3,12),PSEUD3(30,3,12)00113000
     *,H,A,EM1                                                          00113100
      COMMON /ABLABL/PSEUD2,PSEUD3                                      00113300
      COMMON CAI,EM,X,XPUNKT,BESCHL,PSEUD1,H,D                          00113400
      DO 2 K=1,3                                                        00113500
      X(1,K)=0.D0                                                       00113600
      DO 1 J=2,N                                                        00113700
    1 X(1,K)=X(1,K)+EM(J)*X(J,K)                                        00113800
    2 X(1,K)=-X(1,K)/EM(1)                                              00113900
      DO 14 J=1,N                                                       00114000
      EM1 = EM(J)                                                       00114200
      DO 3 K=1,3                                                        00114300
    3 BESCHL(J,K)=0.0D0                                                 00114400
      DO 14 I=1,N                                                       00114500
      IF(I-J)4,14,7                                                     00114600
    4 IF(EM(I)) 5,14,5                                                  00114700
    5 A=EM(I)/D(I,J)                                                    00114900
      DO 6 K=1,3                                                        00115000
    6 BESCHL(J,K)=BESCHL(J,K)+A*(X(I,K)-X(J,K))                         00115100
      GO TO 14                                                          00115200
    7 IF(EM(I)) 9,8,9                                                   00115300
    8 IF(EM1) 9,14,9                                                    00115400
    9 DO 10 K=1,3                                                       00115500
   10 DX(K)=X(I,K)-X(J,K)                                               00115600
      A=DX(1)*DX(1)+DX(2)*DX(2)+DX(3)*DX(3)                             00115700
      D(I,J)=DSQRT(A)                                                   00115800
      D(J,I)=A*D(I,J)                                                   00115900
      IF(J-1)11,14,11                                                   00116000
   11 IF(EM(I)) 12,14,12                                                00116100
   12 A=EM(I)/D(J,I)                                                    00116300
      DO 13 K=1,3                                                       00116400
   13 BESCHL(J,K)=BESCHL(J,K)+A*DX(K)                                   00116500
   14 CONTINUE                                                          00116600
      RETURN                                                            00116700
      END                                                               00116800
      SUBROUTINE SCHRTT (N,M)                                           00116900
C                                                                       00117000
C     EXTRAPOLATIONSSCHRITT,ORDNUNG=2*M+1 - STEP OF INTEGRATION         00117100
C                                                                       00117200
      DOUBLE PRECISION CAI(5,11),EM(30),X(30,3),XPUNKT(30,3),BESCHL(30,300117300
     *),XNABLA(30,3),A(30,3,12),B(30,3,12),DD(11),H,HQ,S                00117400
      COMMON /ABLABL/A,B                                                00117500
      COMMON CAI,EM,X,XPUNKT,BESCHL,XNABLA,H                            00117600
      DD(1)=0.0D0                                                       00117700
      DD(2)=0.83333333333333333D-1                                      00117800
      DD(3)=0.83333333333333333D-1                                      00117900
      DD(4)=0.79166666666666667D-1                                      00118000
      DD(5)=0.75D-1                                                     00118100
      DD(6)=0.71345899470899471D-1                                      00118200
      DD(7)=0.68204365079365079D-1                                      00118300
      DD(8)=0.65495756172839506D-1                                      00118400
      DD(9)=0.63140432098765432D-1                                      00118500
      DD(10)=.61072649861712362D-1                                      00118600
      DD(11)=.59240564123376623D-1                                      00118700
      HQ=H*H                                                            00118800
      M1=2*M+1                                                          00118900
      M2=M1+1                                                           00119000
      DO 2 J=2,N                                                        00119100
      DO 2 K=1,3                                                        00119200
      S=0.0D0                                                           00119300
      DO 1 L=2,M1                                                       00119400
    1 S=S+DD(L)*B(J,K,L+1)                                              00119500
      XNABLA(J,K)=XNABLA(J,K)+HQ*(B(J,K,1)+S)                           00119600
    2 X(J,K)=X(J,K)+XNABLA(J,K)                                         00119700
      CALL WW(N)                                                        00119800
      DO 4 J=2,N                                                        00119900
      DO 4 K=1,3                                                        00120000
      A(J,K,1)=BESCHL(J,K)                                              00120100
      DO 3 L=1,M1                                                       00120200
    3 A(J,K,L+1)=A(J,K,L)-B(J,K,L)                                      00120300
      DO 4 I=1,M2                                                       00120400
    4 B(J,K,I)=A(J,K,I)                                                 00120500
      RETURN                                                            00120600
      END                                                               00120700
      SUBROUTINE DRUCKE (N,M, IZ,T,KD,IED,IPOUT)                        00120800
C                                                                       00120900
C     IPOUT=  0 NOTHING,  1 PRINT, 2 OUTPUT UNFORMATTED, 3 BOTH         00121000
C                                                                       00121100
C     KD= 1:  RELATIVE, otherwise BARYCENTRIC, IED=1: WITH VELOCITIES   00121200
C                                                                       00121300
      DOUBLE PRECISION CAI(5,11),EM(30),X(30,3),XPUNKT(30,3),BEX( 60,3),00121400
     *A(30,3,12),B(30,3,12),S(3),SV(3),T,H                              00121500
      DIMENSION IPOUT(30)                                               00121600
      COMMON /ABLABL/A,B                                                00121700
      COMMON CAI,EM,X,XPUNKT,BEX,H                                      00121800
 1000 FORMAT(4HBODY,I3,3H T=,1PD14.7,7H  X,Y,Z,3D14.6)                  00122200
 1001 FORMAT(34H0 0.1*H*H*MAX(LAST USED DIFF.)  = ,1PE10.2)             00122500
 1002 FORMAT(4HBODY,I3,3H T=,1PD14.7,7H  X,Y,Z,3D14.6,8H  VELOC.,3D14.6)
      M1=2*M+2                                                          00122600
      IP=0                                                              00122700
      DO 2 J=1,N                                                        00122800
      I=IPOUT(J)+1                                                      00122900
      GO TO (2,1,2,1),I                                                 00123000
    1 IP=1                                                              00123100
    2 CONTINUE                                                          00123200
    4 DO 18 J=1,N                                                       00123700
      IF(J.GE.2) GO TO 6                                                00123800
    5 IF(KD.EQ.1) GO TO 18                                              00123900
    6 I=IPOUT(J)                                                        00124000
      IF(I.EQ.0) GO TO 18                                               00124100
    7 DO 9 K=1,3                                                        00124200
      S(K)=X(J,K)                                                       00124300
      IF(IED.EQ.1) SV(K)=XPUNKT(J,K)
      IF(KD.NE.1) GO TO 9                                               00124400
    8 S(K)=S(K)-X(1,K)                                                  00124500
      IF(IED.EQ.1) SV(K)=SV(K)-XPUNKT(1,K)
    9 CONTINUE                                                          00124600
      IF(I.EQ.2) GO TO 14                                               00124700
      IF(IED.EQ.1) GOTO 12
   11 WRITE(6,1000) J,T,S                                               00125300
      GOTO 14
   12 WRITE(6,1002) J,T,S,SV
   14 IF(I.EQ.1) GO TO 18                                               00125600
      IF(IED.EQ.1) GOTO 16
   15 WRITE(7) J,T,S                                                    00125800
      GOTO 18
   16 WRITE(7) J,T,S,SV
   18 CONTINUE                                                          00126100
      IF(IP.EQ.0)  GOTO 21                                              00126200
   19 BM=0.0                                                            00126300
      DO 20 J=2,N                                                       00126400
      DO 20 K=1,3                                                       00126500
      BN=DABS(B(J,K,M1))                                                00126600
   20 BM=AMAX1(BM,BN)                                                   00126700
      BM  =BM*0.1*H*H                                                   00126800
      WRITE(6,1001)BM                                                   00126900
   21 RETURN                                                            00127000
      END                                                               00127100
      SUBROUTINE ELEMNT (N,M,IZ,T,KD)                                   00127200
C                                                                       00127300
C     COORDINATES AND VELOCITIES, KD=1 RELATIVE, otherwise BARYCENTRIC
C                                                                       00127500
      DOUBLE PRECISION CAI(5,11),EM(30),X(30,3),XPUNKT(30,3),BESCHL(30, 00127600
     *3),XNABLA(30,3),A(30,3,12),B(30,3,12),GG(11),V(6),H,HD,S,T        00127700
      COMMON /ABLABL/A,B                                                00127800
      COMMON CAI,EM,X,XPUNKT,BESCHL,XNABLA,H                            00127900
      GG(1)=0.16666666666666667D0                                       00128000
      GG(2)=0.41666666666666667D-1                                      00128100
      GG(3)=0.22222222222222222D-1                                      00128200
      GG(4)=0.14583333333333333D-1                                      00128300
      GG(5)=0.10615079365079365D-1                                      00128400
      GG(6)=0.82258597883597884D-2                                      00128500
      GG(7)=0.66479276895943563D-2                                      00128600
      GG(8)=0.55372299382716049D-2                                      00128700
      GG(9)=0.47180677475816365D-2                                      00128800
      GG(10)=.40919706740019240D-2                                      00128900
      GG(11)=.35997549720267974D-2                                      00129000
      HD=1.D0/H                                                         00129100
      M1=2*M+1                                                          00129200
      M2=M1+1                                                           00129300
      J1=0                                                              00129400
  999 FORMAT(1H0/1H )                                                   00129500
 1000 FORMAT(51H     COORDINATES AND VELOCITIES FOR STEP NUMBER     ,I9)00129600
 1001 FORMAT(9H0     T = ,1PD17.9)                                      00129700
 1004 FORMAT(32H00.01*H*MAX(LAST USED DIFF.)  = ,1PE10.2  )             00129900
 1005 FORMAT(14H0 BODY NUMBER ,I2,6H  EM= ,1PD23.15)                    00130000
 1006 FORMAT(9H COORD.= ,1PD23.15,6X,9HVELOCITY=,1PD23.15)              00130100
      DO 3 K=1,3                                                        00130300
      XPUNKT(1,K)=0.0D0                                                 00130400
      DO 2 J=2,N                                                        00130500
      S=0.0D0                                                           00130600
      DO 1 L=1,M1                                                       00130700
    1 S=S+GG(L)*B(J,K,L+1)                                              00130800
      XPUNKT(J,K)=HD*XNABLA(J,K)+H*(B(J,K,1)*0.5D0-S)                   00130900
    2 XPUNKT(1,K)=XPUNKT(1,K)+EM(J)*XPUNKT(J,K)                         00131000
    3 XPUNKT(1,K)=-XPUNKT(1,K)/EM(1)                                    00131100
      WRITE(6,999)                                                      00131300
      WRITE(6,1000)IZ                                                   00131400
      WRITE(6,1001)T                                                    00131500
C
      DO 19 J=1,N                                                       00132400
      IF(J-2)6,7,7                                                      00132500
    6 IF(KD.EQ.1) GO TO 19                                              00132600
    7 DO 9 K=1,3                                                        00132700
      V(K)=X(J,K)                                                       00132800
      V(K+3)=XPUNKT(J,K)                                                00132900
      IF(KD.NE.1) GO TO 9                                               00133000
    8 V(K)=V(K)-X(1,K)                                                  00133100
      V(K+3)=V(K+3)-XPUNKT(1,K)                                         00133200
    9 CONTINUE                                                          00133300
      WRITE(6,1005)J,EM(J)                                              00133500
      DO 12 K=1,3                                                       00134100
      WRITE(6,1006)V(K),V(K+3)                                          00134300
   12 CONTINUE                                                          00134400
   19 CONTINUE                                                          00135500
C
      BM=0.                                                             00135600
      DO 20 J=2,N                                                       00135700
      DO 20 K=1,3                                                       00135800
      BN=DABS(B(J,K,M2))                                                00135900
   20 BM=AMAX1(BM,BN)                                                   00136000
      BM=BM*0.01*H                                                      00136100
      WRITE(6,1004)BM                                                   00136300
      RETURN                                                            00136400
      END                                                               00136500
