      SUBROUTINE MRQCOFK(X,Y1,Y2,SIG1,SIG2,NDATA,A,MA,LISTA,MFIT,ALPHA,
     *                  BETA,CHISQ,DO_1,DO_2)

C FROM NUMERICAL RECIPES; ADAPTED FOR RADIAL VELOCITY CURVE FITTING
C  4 SEPT 1991        J.T. ARMSTRONG
C  14 JULY 1994       J.T. ARMSTRONG
C                       Adapted for radial velocity curve fitting
C                       using K1 and K2 rather than a1 sin i and
C                       a2 sin i

      IMPLICIT UNDEFINED (A-Z)
      INCLUDE 'VELPARM.INC'

      REAL*8 X(MAXVDATA),Y1(MAXVDATA),Y2(MAXVDATA)
      REAL*8 SIG1(MAXVDATA),SIG2(MAXVDATA)
      REAL*8 A(MMAX),DVELDA(MMAX,2),ALPHA(MMAX,MMAX),BETA(MMAX)
      REAL*8 SIG2I,CHISQ,Y1MOD,Y2MOD,DYDA(MMAX,2)
      REAL*8 DY,WT

      INTEGER*4 LISTA(MMAX)
      INTEGER*4 J,K,I,MA,MFIT,NDATA

      LOGICAL DO_1,DO_2

      DO J=1,MFIT
        DO K=1,J
          ALPHA(J,K)=0.
        END DO
        BETA(J)=0.
      END DO
      CHISQ=0.

D     OPEN(UNIT=14,FILE='TEST.OUT',STATUS='NEW',FORM='FORMATTED')
D     WRITE(14,'(A,A)') 
D    *         'star i j k   sig2i     dy    dy/da(j)    ',
D    *         'wt     d(beta) dy/da(k) d(alpha)'
C# ## ## ## #.##e+## #.##e+## #.##e+## #.##e+## #.##e+## #.##e+## #.##e+##
  101 FORMAT(1X,A,1X,I2,'       ',2(G8.2,1X))
  102 FORMAT(1X,'     ',I2,'                      ',3(G8.2,1X))
  103 FORMAT(1X,'        ',I2,
     *     '                                              ',2(G8.2,1X))

      DO I=1,NDATA
        CALL TRUK2VEL(X(I),A,Y1MOD,Y2MOD,DYDA)

        IF(DO_1 .AND. Y1(I).LT.999.) THEN
          SIG2I=1./(SIG1(I)*SIG1(I))
          DY=Y1(I)-Y1MOD
D         WRITE(14,101) '1',I,SIG2I,DY
          DO J=1,MFIT
            WT=DYDA(LISTA(J),1)*SIG2I
            BETA(J)=BETA(J)+DY*WT
D           WRITE(14,102) LISTA(J),DYDA(LISTA(J),1),WT,DY*WT
            DO K=1,J
              ALPHA(J,K)=ALPHA(J,K)+WT*DYDA(LISTA(K),1)
D             WRITE(14,103) LISTA(K),DYDA(LISTA(K),1),
D    *                      WT*DYDA(LISTA(K),1)
            END DO
          END DO
          CHISQ=CHISQ+DY*DY*SIG2I
        END IF

        IF(DO_2.AND.Y2(I).LT.999.) THEN
          SIG2I=1./(SIG2(I)*SIG2(I))
          DY=Y2(I)-Y2MOD
D         WRITE(14,101) '2',I,SIG2I,DY
          DO J=1,MFIT
            WT=DYDA(LISTA(J),2)*SIG2I
            BETA(J)=BETA(J)+DY*WT
D           WRITE(14,102) LISTA(J),DYDA(LISTA(J),2),WT,DY*WT
            DO K=1,J
              ALPHA(J,K)=ALPHA(J,K)+WT*DYDA(LISTA(K),2)
D             WRITE(14,103) LISTA(K),DYDA(LISTA(K),2),
D    *                      WT*DYDA(LISTA(K),2)
            END DO
          END DO
          CHISQ=CHISQ+DY*DY*SIG2I
        END IF

      END DO

D     WRITE(OUTC,*) ' MRQCOF:  ALPHA = '
D     DO J=1,MFIT
D       WRITE(OUTC,'(1X,2(G12.5,1X))') 
D    *                (ALPHA(J,K),K=1,MFIT)
D     END DO
D     WRITE(14,*) ' MRQCOF:  ALPHA = '
D     DO J=1,MFIT
D       WRITE(14,'(1X,2(G12.5,1X))') 
D    *                (ALPHA(J,K),K=1,MFIT)
D     END DO
D     CLOSE(14)

      DO J=2,MFIT
        DO K=1,J-1
          ALPHA(K,J)=ALPHA(J,K)
        END DO
      END DO
      RETURN
      END

