c $1BFTC MAIN M94/2,X_?
      COMMON /COM/ V, F, PI, RAD, VI
      DOUBLE PRECISION V(16,7),F(80,4),PI,RAD,VI(16)
      COMMON /COM/ ITRAT,KOUNT
      COMMON /COM/ NBODYI, NBODY
      COMMON /COM/ IPRT(4), IPRTI(4)
      COMMON /COM/ KL, IPG, LINCT, LINPGE
      COMMON /COM/ NBLOCK
      open(5,file='prg.data')
    1 CALL INPUT
    2 CALL EPHEM
      write(*,'(''AP EPHEM'')')
    3 CALL VMASS
      write(*,*)'AP VMASS',itrat
      IF(ITRAT .EQ. 1) GO TO 4
      IF(ITRAT .EQ. 2) GO TO 7
C INITIALIZATION OF VIRTUAL MASS-DEPENDENT VALUES
      V(7,7) = V(6,1)
      write(*,*)'A v(7,7) ',V(7,7),' v(6,1) ',V(6,1)
      DO 600 J=2,4
          V(10,J+3)=V(8,J)
          V(9,J)=V(10,J)
  600     V(11,J)=V(12,J)
      V(9,1)=V(10,1)
      V(8,5)=1.
      V(5,1)=V(6,1)
      GO TO 9
    4 ITRAT=2
C VIRTUAL MASS AVERAGE MAGNITUDE AND VELOCITY
    5 V(7,7)=.5*V(5,1)+.5*V(6,1)
      write(*,*)'B  v(7,7) ',V(7,7),' v(6,1) ',V(6,1)
      DO 390 J=2,4
          V(10,J+3)=(V(6,J)-V(5,J))/V(7,6)
  390     V(11,J)=V(3,J)-V(10,J+3)
    9 write(*,'(''AV VECTOR'')')
      CALL VECTOR
      write(*,'(''AP VECTOR'')')
      IF(KOUNT .LT. 0) GO TO 10
      IF(ITRAT .EQ. 1) GO TO 2
      IF(ITRAT .EQ. 2) GO TO 3
C VIRTUAL MASS AVERAGE ACCELERATIONS
    7 V(8,6)=(V(6,1)-V(5,1)-V(7,1)*V(7,6))/V(8,5)
      DO 340 J=2,4
  340     V(10,J+3)=(V(6,J)-V(5,J)-V(7,J)*V(7,6))/V(8,5)
C TEST FOR STOPPING CONDITIONS
  400 IF (V(2,5) .GT. V(2,1)) GO TO 401
      CALL SPACE (3)
      WRITE (6,4000)
 4000 FORMAT (//53H STOPPING CONDITION--EXCEEDED MAXIMUM TRAJECTORY TIME
     $)
      GO TO 995
  401 IF (F(3,4) .GT. V(2,6)) GO TO 402
      CALL SPACE (3)
      WRITE (6,4010)
 4010 FORMAT (//35H STOPPING CONDITION--IMPACTED EARTH)
      GO TO 995
  402 IF (F(7,4) .GT. V(2,7)) GO TO 403
      CALL SPACE (3)
      WRITE (6,4020)
 4020 FORMAT (//35H STOPPING CONDITION--IMPACTED MOON )
      GO TO 995
  403 CONTINUE
      IF(KOUNT .EQ. 0) GO TO 11
      KOUNT = 0
   10 CALL PRINT
      write(*,*)'AP PRINT',kount
      IF(KOUNT .LT. 0) GO TO 1
   11 CALL ESTMT
      write(*,*)'AP ESTMT'
      GO TO 5
  995 KOUNT=-1
      GO TO 10
      END

      SUBROUTINE INPUT
      COMMON /COM/ V, F, PI, RAD, VI
      DOUBLE PRECISION V(16,7),F(80,4),PI,RAD,VI(16)
      COMMON /COM/ ITRAT,KOUNT
      COMMON /COM/ NBODYI, NBODY
      COMMON /COM/ IPRT(4), IPRTI(4)
      COMMON /COM/ KL, IPG, LINCT, LINPGE
      COMMON /COM/ NBLOCK
      character*6 icard
      COMMON /GCDIN/ ICARD(14)
      
      character*6 aprt,aaprt
      DIMENSION INCHK(6),FMT(2),CRDTYP(2,6),APRT(4),AAPRT(2,4),IX(3)
      DIMENSION NBLK(4)
      DOUBLE PRECISION WD(16)
      EQUIVALENCE (WRD,WD(1))
      DATA INCHK /0,1,1,1,0,1 /
      DATA NBLK /10,11,9,12 /
      DATA FMT(1) /12H(A24,4D24,0)/
      DATA CRDTYP (1,1) /12HPRINT       /
      DATA CRDTYP (1,2) /12HPOSITION    /
      DATA CRDTYP (1,3) /12HVELOCITY    /
      DATA CRDTYP (1,4) /12HEFEMERIS    /
      DATA CRDTYP (1,5) /12HACCURACY    /
      DATA CRDTYP (1,6) /12HSTOP        /
      DATA AAPRT /6H      ,6H      ,
     $            6H      ,6H   1  ,
     $            6H      ,6H   2  ,
     $            6H      ,6H   3   /
      DATA INERR /0/
      DATA PRINT    / 6HPRINT  /
      DATA POSITI   / 6HPOSITI /
      DATA VELOCI   / 6HVELOCI /
      DATA EFEMER   / 6HEFEMER /
      DATA ACCURA   / 6HACCURA /
      DATA STOP     / 6HSTOP   /
      DATA PROBLE   / 6HPROBLE /
      DATA IENTRY /0/
      IPG=0
      KL=-77777
    1 DO 10 I=1,16
         DO 10 J=1,7
   10        V(I,J)=0.
      DO 20 I=1,80
          DO 20 J=1,4
   20        F(I,J)=0.
      V( 3,5)=VI(1)
      V( 1,1)=VI(2)
      V( 1,2)=VI(3)
      V( 1,3)=VI(4)
      V( 1,4)=VI(5)
      V( 3,2)=VI(6)
      V( 3,3)=VI(7)
      V( 3,4)=VI(8)
      V( 3,1)=VI(9)
      V( 1,5)=VI(10)
      V( 1,6)=VI(11)
      V( 1,7)=VI(12)
      V( 4,5)=VI(13)
      V( 2,5)=VI(14)
      V( 2,6)=VI(15)
      V( 2,7)=VI(16)
      DO 30 I=1,4
   30     IPRT(I)=IPRTI(I)
  100 CALL DINPT (5,NOUT,FMT,WD)
      IF(WRD.EQ.PRINT ) GO TO 110
      IF(WRD.EQ.POSITI) GO TO 120
      IF(WRD.EQ.VELOCI) GO TO 130
      IF(WRD.EQ.EFEMER) GO TO 140 
      IF(WRD.EQ.ACCURA) GO TO 150
      IF(WRD.EQ.STOP  ) GO TO 160
      IF(WRD.EQ.PROBLE) GO TO 170
      GO TO 500
  110 CONTINUE
      V(3,5)=WD(2)
      BACKSPACE 5
      READ (5,111)IX
  111 FORMAT (41X3I1)
      DO 113 I=2,4
  113     IPRT(I)=0
      DO 112 I=1,3
          ISUB=IX(I)
          IF (ISUB.EQ.0) GO TO 112
          IPRT(ISUB+1)=1
  112 CONTINUE
      INCHK(1)=0
      GO TO 100
  120 CONTINUE
      V(1,1)=WD(2)
      V(1,2)=WD(3)
      V(1,3)=WD(4)
      V(1,4)=WD(5)
      INCHK(2)=0
      GO TO 100
  130 CONTINUE
      V(3,2)=WD(2)
      V(3,3)=WD(3)
      V(3,4)=WD(4)
      INCHK(3)=0
      GO TO 100
  140 CONTINUE
      V(3,1)=WD(2)
      V(1,5)=WD(3)
      V(1,6)=WD(4)
      V(1,7)=WD(5)
      INCHK(4)=0
      GO TO 100
  150 CONTINUE
      V(4,5)=WD(2)
      INCHK(5)=0
      GO TO 100
  160 CONTINUE
      V(2,5)=WD(2)
      V(2,6)=WD(3)
      V(2,7)=WD(4)
      INCHK(6)=0
      GO TO 100
  170 CONTINUE
      KLT=WD(2)
      GO TO 600
  500 IF (INERR.NE.1) CALL NEWPGE
      CALL SPACE (5)
      WRITE (6,501)ICARD
  501 FORMAT (/56H AN ERROR HAS BEEN DETECTED READING THE FOLLOWING CARD
     $--/10X,13A6,A2/18H EXECUTION DELETED/1H )
      INERR=1
      GO TO 100
  600 IF (IENTRY.NE.0) GO TO 650
      DO 601 I=1,6
  601     IF (INCHK(I).NE.0) GO TO 602
      GO TO 649
  602 CALL NEWPGE
      WRITE (6,604)
  604 FORMAT (/91H THE FOLLOWING REQUIRED CARD TYPES WERE NOT INPUT ON T
     $HE FIRST PROBLEM -- EXECUTION DELETED)
      DO 605 I=1,6
          IF (INCHK(I).EQ.0) GO TO 605
          WRITE (6,606) (CRDTYP(J,I),J=1,2)
  606 FORMAT (10X,2A6)
  605 CONTINUE
      INERR=1
  649 IENTRY=1
      DO 610 I=1,4
  610     IPRTI(I)=IPRT(I)
      VI( 1)=V(3,5)
      VI( 2)=V(1,1)
      VI( 3)=V(1,2)
      VI( 4)=V(1,3)
      VI( 5)=V(1,4)
      VI( 6)=V(3,2)
      VI( 7)=V(3,3)
      VI( 8)=V(3,4)
      VI( 9)=V(3,1)
      VI(10)=V(1,5)
      VI(11)=v(1,6)
      VI(12)=V(1,7)
      VI(13)=V(4,5)
      VI(14)=V(2,5)
      VI(15)=V(2,6)
      VI(16)=V(2,7)
  650 IF (INERR.NE.0) GO TO 100
      NBODY=NBODYI
      NBLOCK=0
      DO 660 I=1,4
  660     IF (IPRT(I).NE.0) NBLOCK=NBLOCK+NBLK(I)
      CALL NEWPGE
      KL=KLT
      DO 651 I=1,4
          J=IPRT(I)+1
  651     APRT(I)=AAPRT(J,I)
      WRITE (6,652) V(3,5),
     $              APRT,V(1,1),(V(1,J),J=2,4),(V(3,J),J=2,4),
     $              V(3,1),(V(1,J),J=5,7),V(4,5),(V(2,J),J=5,7),KL
  652 FORMAT (12H PRINT      1PD20.11,6X4A6/
     $        12H POSITION   4D20.11/
     $        12H VELOCITY   3D20.11/
     $        12H EFEMERIS   4D20.11/
     $        12H ACCURACY    D20.11/
     $        12H STOP       3D20.11/
     $        12H PROBLEM    I5)
      CALL NEWPGE
C DIMENSIONAL CONVERSION FACTORS
      V(1,5)=V(1,5)/RAD
      V(4,7)=V(1,5)*V(1,6)
      V(5,5)=V(4,7)*V(1,6)
      V(5,6)=V(5,5)*V(1,5)
      V(6,5)=V(5,6)*V(1,6)
      V(6,6)=V(6,5)*V(1,5)
C EPHEMERIS DATA
      V(1,7)=-V(1,7)
      V(5,7)=1.+V(1,7)
      F(1,4)=V(5,7)
      F(5,4)=-V(1,7)
      NBODY = 2
      NBODY=4*NBODY-3
C NONDIMENSIONALIZATION
      DO 50 J=2,3
          V(J,5)=V(J,5)*V(1,5)
   50     V(2,J+4)=V(2,J+4)/V(1,6)
      V(2,1)=V(1,1)*V(1,5)
      V(4,1)=V(3,1)*V(1,5)
      DO 51 J=2,4
          V(2,J)=V(1,J)/V(1,6)
   51     V(4,J)=V(3,J)/V(4,7)
C INITIALIZATION OF MISCELANEOUS VALUES
      V(3,6)=DEXP(1.13756474179255 + .509713741462307*DLOG(V(4,5))
     $          +.14560181279278D-2 * DLOG(V(4,5))**2 )
      ITRAT = 3
      KOUNT = 1
      V(13,3)=V(2,1)+V(3,5)
      RETURN
      END

      SUBROUTINE EPHEM
      COMMON /COM/ V, F, PI, RAD, VI
      DOUBLE PRECISION V(16,7),F(80,4),PI,RAD,VI(16)
      COMMON /COM/ ITRAT,KOUNT
      COMMON /COM/ NBODYI, NBODY
      COMMON /COM/ IPRT(4), IPRTI(4)
      COMMON /COM/ KL, IPG, LINCT, LINPGE
      COMMON /COM/ NBLOCK
      V(14,7)=DSIN(V(4,1))
      V(14,6)=DCOS(V(4,1))
c      DO 101 I=1,5,4
      DO 101 I=1,NBODY,4
         DO 100 J=1,2
  100         F(I,J)=V(I,7)*V(14,J+5)
         F(I+1,1)=-F(I,2)
  101    F(I+1,2)=F(I,1)
      RETURN
      END

      SUBROUTINE VMASS
      COMMON /COM/ V, F, PI, RAD, VI
      DOUBLE PRECISION V(16,7),F(80,4),PI,RAD,VI(16)
      COMMON /COM/ ITRAT,KOUNT
      COMMON /COM/ NBODYI, NBODY
      COMMON /COM/ IPRT(4), IPRTI(4)
      COMMON /COM/ KL, IPG, LINCT, LINPGE
      COMMON /COM/ NBLOCK
C VIRTUAL. MASS POSITION AND MAGNITUDE
      V(12,5)=0.
      DO 201 I=1,NBODY,4
          DO 200 J=1,3
  200         F(I+2,J)=V(2,J+1)-F(I,J)
          F(I+2,4)=DSQRT(F(I+2,1)**2+F(I+2,2)**2+F(I+2,3)**2)
          F(I+3,4)=F(I,4)/F(I+2,4)**3
  201     V(12,5)=V(12,5)+F(I+3,4)
      DO 203 J=1,3
          V(6,J+1)=0.
          DO 202 I=1,NBODY,4
  202          V(6,J+1)=V(6,J+1)+F(I+3,4)*F(I,J)
          V(6,J+1)=V(6,J+1)/V(12,5)
  203     V(10,J+1)=V(2,J+1)-V(6,J+1)
      V(10,1)=DSQRT(V(10,2)**2+V(10,3)**2+V(10,4)**2)
      V(6,1)=V(10,1)**3*V(12,5)
      write(*,*)'VMASS A V(6,1) ',v(6,1)     
C VIRTUAL MASS VELOCATY AND MAGNITUDE RATE
      V(12,6)=0.
      DO 301 I=1,NBODY,4
          DO 300 J=1,3
  300         F(I+3,J)=V(4,J+1)-F(I+1,J)
      F(I+1,4)=3.*(F(I+2,1)*F(I+3,1)+F(I+2,2)*F(I+3,2)+F(I+2,3)*F(I+3,3)
     1)/F(I+2,4)**2
  301      V(12,6)=V(12,6)-F(I+1,4)*F(I+3,4)
      DO 303 J=1,3
          V(8,J+1)=0.
          DO 302 I=1,NBODY,4
  302         V(8,J+1)=v(8,J+1)+F(I+3,4)*(F(I+1,J)-F(I,J)*F(I+1,4))
          V(8,J+1)=(V(8,J+1)-V(6,J+1)*V(12,6))/V(12,5)
  303     V(12,J+1)=V(4,J+1)-V(8,J+1)
      write(*,*)'vmass',v(10,1),v(12,6),v(10,1)**v(12,6)
      
      V(8,1)=V(6,1)*(3.*(V(10,2)*V(12,2)+v(10,3)*v(12,3)+v(10,4)*V(12,4)
     1)/V(10,1)**V(12,6)/V(12,5))
      DO 901 I=1,NBODY,4
          write(*,*)'body',1+i/4
          write(*,*)'  xyzµ ',(f(i  ,j),j=1,4)
          write(*,*)'  xyzµ''',(f(i+1,j),j=1,4)
          write(*,*)' Sxyzµ ',(f(i+2,j),j=1,4)
          write(*,*)' Sxyzµ''',(f(i+3,j),j=1,4)
  901 continue       
      RETURN
      END

      SUBROUTINE VECTOR
      COMMON /COM/ V, F, PI, RAD, VI
      DOUBLE PRECISION V(16,7),F(80,4),PI,RAD,VI(16)
      COMMON /COM/ ITRAT,KOUNT
      COMMON /COM/ NBODYI, NBODY
      COMMON /COM/ IPRT(4), IPRTI(4)
      COMMON /COM/ KL, IPG, LINCT, LINPGE
      COMMON /COM/ NBLOCK
C VECTOR ORBITAL ELEMENTS
      DO 401 J=2,4
          DO 400 I=9,11,2
  400         V(I,J+3)=V(I,J)
  401     V(16,J)=V(9,J+1)*V(11,J+2)-V(9,J+2)*V(11,J+1)
      DO 403 J=2,4
          DO 402 I=11,16,5
  402         V(I,J+3)=V(I,J)
  403 V(14,J)=-V(9,J)/V(9,1)-(V(16,J+1)*V(11,J+2)-V(16,J+2)*V(11,J+1))/V
     1(7,7)
      DO 404 I=14,16,2
          V(I,5)=V(I,2)**2+V(I,3)**2+V(I,4)**2
  404     V(I,1)=DSQRT(V(I,5) )
      V(13,5)=1.-V(14,5)
      V(13,6)=DSQRT(DABS(V(13,5)))
      V(13,7)=V(16,5)/V(7,7)
      write(*,*)'C v(7,7) ',V(7,7)
      write(*,*)V(13,7),V(16,5)
      IF (ITRAT.EQ.3) RETURN
C SPACECRAFT FINAL POSITION AND VELOCITY
      V(6,7)=V(7,6)+V(7,6)*V(7,6)*V(6,7)
      DO 410 J=2,4
  410     V(9,J+3)=V(9,J)+V(6,7)*V(11,J)
      V(8,7)=V(13,7)/(V(14,2)*V(9,5)+V(14,3)*V(9,6)+V(14,4)*V(9,7)+DSQRT
     1(V(9,5)**2+V(9,6)**2+V(9,7)**2))
      DO 411 J=2,4
  411     V(10,J)=V(8,7)*V(9,J+3)
      V(10,1)=DSQRT(V(10,2)**2+V(10,3)**2+V(10,4)**2)
      DO 414 J=2,4
  414     V(12,J)=V(14,J)+V(10,J)/V(10,1)
      DO 413 J=2,4
          DO 412 I=12,16,4
  412         V(I,J+3)=V(I,J)
          V(12,J)=(V(16,J+1)*V(12,J+2)-V(16,J+2)*V(12,J+1))/V(13,7)
          V(2,J)=V(10,J)+V(6,J)
  413     V(4,J)=V(12,J)+V(10,J+3)
C KEPLERIAN TIME OF FLIGHT
      IF(V(14,1) .NE. 0.) GO TO 501
  500 M=9
      GO TO 502
  501 M=14
  502 NN=16-M
      DO 504 J=2,4
          DO 503 I=M,16,NN
  503         V(I,J+3)=V(I,J)
  504 V(11,J+3)=(V(16,J+1)*V(M,J+2)-V(M,J+1)*V(16,J+2))/(V(16,1)*V(M,1)) 
      DO 517 I=9,10
          IF(V(13,5) .EQ. 0.) GO TO 506
  505     V(I+6,5)=V(13,7)/V(13,6)
          write(*,*)v(i+6,5),v(13,7),v(13,6)
          GO TO 507
  506     V(I+6,5)=2./(V(I,1)-V(13,7))
  507   V(I+6,5)=(V(11,5)*V(I,2)+V(11,6)*V(I,3)+V(11,7)*V(I,4))/V(I+6,5)
          write(*,*)v(i+6,5)
          IF(V(13,5)) 508,509,510
  508         V(I+6,6)=DLOG(V(I+6,5)+DSQRT(V(I+6,5)**2+1.))
              GO TO 517

  509         V(I+6,6)=V(14,1)**2*(V(13,7)*V(I+6,5))**3/3.
              GO TO 517 

  510         IF(DABS(V(I+6,5)) .LT. 1.) GO TO 524
                  IF (DABS(V(I+6,5)) .LE. 1.0001) GO TO 597
                      KOUNT=-1
                      CALL SPACE(2)
                      WRITE (6,596)
  596 FORMAT (/27H UNACCEPTABLE ERROR IN ATAN)
                      RETURN
  597             CALL SPACE (2)
                  WRITE (6,598)
  598 FORMAT (/27H   ACCEPTABLE ERROR IN ATAN)
                  V(I+6,6)=DSIGN(PI/2.,V(I+6,5))
                  GO TO 523
  524         V(I+6,6)=DATAN(V(I+6,5)/DSQRT(1.-V(I+6,5)**2))
C ?  
  523         IF(V(14,1) .GT. 0.) GO TO 512
  511            V(I+6,7)=V(I,2)*V(9,2)+V(I,3)*V(9,3)+V(I,4)*V(9,4)
                 GO TO 513
  512         V(I+6,7)=V(13,7)/V(13,5)-V(I,1)
  513         IF(V(I+6,7) .GE. 0.) GO TO 517
                 IF(V(I+6,5) .GE. 0.) GO TO 516
                     V(I+6,6)=-PI-V(I+6,6)
                     GO TO 517
  516            V(I+6,6)=PI-V(I+6,6)
  517 CONTINUE 
      IF(V(13,5) .EQ. 0.) GO TO 519
          V(8,7)=V(7,7)**2/V(16,1)**3*V(13,5)*V(13,6)
          GO TO 520
  519 V(8,7)=V(16,1)*V(14,1)**2/2.
  520 V(7,5)=V(16,6)-V(15,6)+V(14,1)*(V(15,5)-V(16,5))
      IF(V(7,5) .GE. 0.) GO TO 522
          IF (v(8,7).LT.0.) GO TO 522
              V(7,5)=V(7,5)+2.*PI
  522 V(7,5)=V(7,5)/V(8,7)
C UPDATE ESTIMATE OF_CURVE
      IF (V(7,5).EQ.0.) GO TO 525
          V(6,7)=(V(6,7)-V(7,5) )/V(7,5)**2 
  525 RETURN
      END

      SUBROUTINE PRINT
      COMMON /COM/ V, F, PI, RAD, VI
      DOUBLE PRECISION V(16,7),F(80,4),PI,RAD,VI(16)
      COMMON /COM/ ITRAT,KOUNT
      COMMON /COM/ NBODYI, NBODY
      COMMON /COM/ IPRT(4), IPRTI(4)
      COMMON /COM/ KL, IPG, LINCT, LINPGE
      COMMON /COM/ NBLOCK
      DOUBLE PRECISION TEMP(16),RV,VV,RS,VS
      DO 351 J=2,4
          DO 350 I=1,5,4
              V(I,J)=V(1,6)*V(I+1,J)
              V(I+2,J)=V(4,7)*V(I+3,J)
              F(I,J-1)=V(1,6)*F(I,J-1)
  350         F(I+1,J-1)=V(4,7)*F(I+1,J-1)
  351 CONTINUE
      DO 352 J=1,4
  352     V(15,J)=V(5,5)*V(16,J)
      V(1,1)=V(2,1)/V(1,5)
      V(3,1)=V(4,1)/V(1,5)
      V(5,1)=V(6,1)*V(6,5)
      V(7,1)=V(8,1)*V(6,6)
      write(*,*)v(7,1),V(8,1),V(6,6)
      RV=DSQRT(V(5,2)**2+V(5,3)**2+V(5,4)**2)
      VV=DSQRT(V(7,2)**2+V(7,3)**2+V(7,4)**2)
      RS=DSQRT(V(1,2)**2+V(1,3)**2+V(1,4)**2)
      VS=DSQRT(V(3,2)**2+V(3,3)**2+V(3,4)**2)
  410 CALL SPACE (NBLOCK)
      WRITE (6,411) (V(1,I),I=1,4),RS,(V(3,I),I=2,4),VS
  411 FORMAT (/////20H TRAJECTORY TIME =  1PD20.11//
     $             40H SPACECRAFT INERTTIAL TRAJECTORY        /
     $             40H     POSITION . . . . . . . . . . . . . 4D20.11/
     $             40H     VELOCITY . . . . . . . . . . . . . 4D20.11)
  420 IF (IPRT(2).EQ.0) GO TO 430
      TEMP( 1) = F(1,1)
      TEMP( 2) = F(1,2)
      TEMP(3 ) = F(1,3)
      TEMP(4) =DSQRT(TEMP(1)**2+TEMP(2)**2+TEMP(3)**2)
      TEMP( 5) = F(2,1)
      TEMP(6 ) = F(2,2)
      TEMP( 7) = F(2,3)
      TEMP(8) =DSQRT(TEMP(5)**2+TEMP(6)**2+TEMP(7)**2)
      TEMP( 9) = F(5,1)
      TEMP(10) = F(5,2)
      TEMP(11) = F(5,3)
      TEMP(12)=DSQRT(TEMP(9)**2+TEMP(10)**2+TEMP(11)**2)
      TEMP(13) = F(6,1)
      TEMP(14) = F(6,2)
      TEMP(15) = F(6,3)
      TEMP(16)=DSQRT(TEMP(13)**2+TEMP(14)**2+TEMP(15)**2)
      WRITE (6,421) V(3,1),(TEMP(I),I=1,16)
  421 FORMAT (/// 20H EPHEMERIS TIME   =  1PD20.11//
     $            40H EPHEMERIS DATA                         /
     $            40H     POSITION OF EARTH  . . . . . . . . 4D20.11/
     $            40H     VELOCITY OF EARTH  . . . . . . . . 4D20.11//
     $            40H     POSITION OF MOON   . . . . . . . . 4D20.11/
     $            40H     VELOCITY OF MOON   . . . . . . . . 4D20.11)
  430 IF (IPRT(3).EQ.0) GO TO 440
          TEMP ( 1) = V(1,2)-F(1,1)
          TEMP ( 2) = V(1,3)-F(1,2)
          TEMP ( 3) = V(1,4)-F(1,3)
          TEMP( 4)=DSQRT(TEMP(1)**2+TEMP(2)**2+TEMP(3)**2)
          TEMP ( 5) = V(3,2)-F(2,1)
          TEMP ( 6) = V(3,3)-F(2,2)
          TEMP ( 7) = V(3,4)-F(2,3)
          TEMP( 8)=DSQRT(TEMP(5)**2+TEMP(6)**2+TEMP(7)**2)
          TEMP (9 ) = V(1,2)-F(5,1)
          TEMP (10) = V(1,3)-F(5,2)
          TEMP (11) = V(1,4)-F(5,3)
          TEMP(12)=DSQRT(TEMP(9)**2+TEMP(10)**2+TEMP(11)**2)
          TEMP (13) = V(3,2)-F(6,1)
          TEMP (14) = V(3,3)-F(6,2)
          TEMP (15) = V(3,4)-F(6,3)
          TEMP(16)=DSQRT(TEMP(13)**2+TEMP(14)**2+TEMP(15)**2)
          WRITE (6,431) (TEMP(I),I=1,16)
  431 FORMAT(///   40H SPACECRAFT RELATIVE TRAJECTORIES       /
     $             40H     POSITION REL. TO EARTH . . . . . . 4D20.11/
     $             40H     VELOCITY REL. TO EARTH . . . . . . 4D20.11//
     $             40H     POSITION REL. TO MOON  . . . . . . 4D20.11/
     $             40H     VELOCITY REL. TO MOON  . . . . . . 4D20.11)
  440 IF (IPRT(4).EQ.0) RETURN
      TEMP(1)=V( 9,2)*V(1,6)
      TEMP(2)=V( 9,3)*V(1,6)
      TEMP(3)=V( 9,4)*V(1,6)
      TEMP(4)=V( 9,1)*V(1,6)
      TEMP(5)=V(11,2)*V(4,7)
      TEMP(6)=V(11,3)*V(4,7)
      TEMP(7)=V(11,4)*V(4,7)
      TEMP(8)=DSQRT(TEMP(5)**2+TEMP(6)**2+TEMP(7)**2)
      WRITE (6,441) (V(5,I),I=2,4),RV,(V(7,I),I=2,4),VV,
     $     (TEMP(I),I=1,8),
     $              (V(15,I),I=2,4),V(15,1),(V(14,I),I=2,4),V(14,1),
     $              V(5,1),(7,1)
  441 FORMAT(//  / 20H VIRTUAL MASS DATA  /
     $             40H     VIRTUAL MASS POSITION  . . . . . . 4(1PD20.11
     $)/           40H     VIRTUAL MASS VELOCITY  . . . . . . 4D20.11/
     $             40H     SPACECRAFT POS. REL. TO V.M. . . . 4D20.11/
     $             40H     SPACECRAFT VEL. REL. TO V.M. . . . 4D20.11/
     $             40H     KEPLER (ANG. NON.) VECTOR. . . . . 4D20.11/
     $             40H     ECCENTRICITY VECTOR  . . . . . . . 4D20.11/
     $             24H     V.M. MAGN.      =  D20.11/
     $             24H     V.M. MAGN. RATE =  D20.11)
      RETURN
      END

      SUBROUTINE ESTMT
      COMMON /COM/ V, F, PI, RAD, VI
      DOUBLE PRECISION V(16,7),F(80,4),PI,RAD,VI(16)
      COMMON /COM/ ITRAT,KOUNT
      COMMON /COM/ NBODYI, NBODY
      COMMON /COM/ IPRT(4), IPRTI(4)
      COMMON /COM/ KL, IPG, LINCT, LINPGE
      COMMON /COM/ NBLOCK
C INDEX VALUES IN V ARRAY 
      DO 361 I=1,9,2
          DO 360 J=1,4
  360         V(I,J)=V(I+1,J)
  361 CONTINUE
      ITRAT=1
C ESTABLISH COMPUTING TIME INCREMENT
  395 V(13,4)=V(7,6)*V(16,5)/(V(15,5)-V(16,5))
      IF (V(7,6).EQ. 0.) V(13,4)=-1.
      V(11,1)=DSQRT(V(11,2)**2+V(11,3)**2+V(11,4)**2)
      V(7,6)=V(3,6)*V(9,1)/V(11,1)
      IF(V(13,4) .LT. 0.) GO TO 394
          IF(V(13,4) .GT. 1.1*V(7,6)) GO TO 394
              V(7,6)=V(13,4)
              CALL SPACE (5)
              WRITE (6,6000)
 6000 FORMAT (//2X,19HMAJOR AXIS CROSSING//1H )
              GO TO 400
  394 IF(V(1,1)+1.1*V(7,6) .LT. V(13,3)) GO TO 378
  390     V(7,6)=V(13,3)-V(1,1)
          V(13,3)=V(13,3)+V(3,5)
  400     KOUNT = 1
C INCREMENT TIMES
  378 DO 379 I=1,3,2
  379     V(I+1,1)=V(I,1)+V(7,6)
      V(8,5)=V(7,6)**2
      IF (V(2,1).GE.V(13,3)) V(13,3)=V(13,3)+V(3,5)
C ESTIMATE VIRTUAL MASS FINAL POSITION AND MAGNITUDE
      V(6,1)=V(5,1)+V(7,1)*V(7,6)+V(8,6)*V(8,5)
      write(*,*)'ESTMT A v(6,1) ',V(6,1)
      
      DO 380 J=2,4
  380     V(6,J)=V(5,J)+V(7,J)*V(7,6)+V(10,J+3)*V(8,5)
      RETURN
      END

      BLOCK DATA
      COMMON /COM/ V, F, PI, RAD, VI
      DOUBLE PRECISION V(16,7),F(80,4),PI,RAD,VI(16)
      COMMON /COM/ ITRAT,KOUNT
      COMMON /COM/ NBODYI, NBODY
      COMMON /COM/ IPRT(4), IPRTI(4)
      COMMON /COM/ KL, IPG, LINCT, LINPGE
      COMMON /COM/ NBLOCK
      DATA PI   /3.141592653589793 /
      DATA RAD  /57.29577951308232 /
      DATA NBODYI  /     2         /
      DATA IPRTI   / 1,0,0,0       /
      DATA LINPGE/60/
      DATA VI / 5.D0, 11*0.D0,1.D-7, 3*0.D0 /
      END

      SUBROUTINE DINPT (X,XX,XXX,WD)
      DOUBLE PRECISION WD(5)
      character*6 icard
      COMMON/GCDIN/ ICARD(14)
      READ (5,1) ICARD
    1 FORMAT (13A6,A2)
c      write(*,'(13A6,''icard '',A2)')icard
      BACKSPACE 5
      READ (5,2) WD
    2 FORMAT (A6,2X,4D18.0)
c      write(*,'(A6,2X,4D18.5)')wd
      RETURN
      END

      SUBROUTINE SPACE(LINES)
      COMMON /COM/ V, F, PI, RAD, VI
      DOUBLE PRECISION V(16,7),F(80,4),PI,RAD,VI(16)
      COMMON /COM/ ITRAT,KOUNT
      COMMON /COM/ NBODYI, NBODY
      COMMON /COM/ IPRT(4), IPRTI(4)
      COMMON /COM/ KL, IPG, LINCT, LINPGE
      COMMON /COM/ NBLOCK
      IF (LINPGE.LT.(LINCT+LINES) ) CALL NEWPGE
      LINCT=LINCT+LINES 
      RETURN
      END

      SUBROUTINE NEWPGE
      COMMON /COM/ V, F, PI, RAD, VI
      DOUBLE PRECISION V(16,7),F(80,4),PI,RAD,VI(16)
      COMMON /COM/ ITRAT,KOUNT
      COMMON /COM/ NBODYI, NBODY
      COMMON /COM/ IPRT(4), IPRTI(4)
      COMMON /COM/ KL, IPG, LINCT, LINPGE
      COMMON /COM/ NBLOCK
      IPG=IPG+1
      WRITE (6,1)
    1 FORMAT                                                       (120H
     $1   V I R T U A L    M A S S    P R O G R A M    F O R    C O M P
     $U T I N G    S P A C E    T R A J E C T O R I E S     )
C
C        WHEN KL = -7777, ONLY TITLE AND PAGE NUMBER ARE GIVEN, AS
C      THIS SIGNALS INPUT DIAGNOSTICS ARE TO BE GIVEN, OR INPUT DATA
C      IS TO UE LISTED,
C
      IF (KL.EQ.-77777) GO TO 10
          WRITE (6,2) KL,IPG
    2 FORMAT (90X8HPROBLEM I5,6X5HPAGE I4///
     $                                  40X80H          X - COMP.
     $   Y - COMP,            Z - COMP,          RESULTANT   )
          LINCT=6
          RETURN
   10 WRITE (6,3) IPG
    3 FORMAT (109X5HPAGE I4)
      LINCT=2
      RETURN
      END
