      SUBROUTINE MRQCOF(NDATA,A,LISTA,MFIT,ALPHA,BETA,NALP)
C ALSO FROM NUMERICAL RECIPES, P. 527
C Adapted for fitting binary orbits from nightly separation and 
C position angle by J.T. Armstrong.
C     25 Apr 1991     J.T. Armstrong  Revised for proper use of
C                                     error ellipse

      IMPLICIT UNDEFINED (A-Z)
      INCLUDE 'PARAMETR.INC'
      INCLUDE 'ELLIPSE.INC'

      INTEGER*4 LISTA(7),MFIT,NALP,NDATA,I,J,K

      REAL*8    SIG2RA,SIG2D,WTRA,WTD
      REAL*8    RMOD,THMOD,DR,RDTHETA,CHIRA,CHIDEC
      REAL*8    RAMOD,DECMOD
      REAL*8    A(7)
      REAL*8    DRDA,DTHDA,THEST_PHI,THOBS_PHI,SINO_PHI,SINE_PHI
      REAL*8    COSO_PHI, COSE_PHI
      REAL*8    COVAR(7,7), ALPHA(7,7), BETA(7), DYDA(7,2)


D     WRITE(PR,*) 'MRQCOF:  ENTERING'

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

      CHISQ=0.

D     WRITE(PR,'(A)') ' MRQCOF:  ELEMENTS TO TRUE2APP:'
D     WRITE(PR,'(A,F14.10,A,F13.10,A,F15.10,A,F18.10 / 
D    1      A,F15.10,A,F15.10,A,F15.10)')
D    1      '  A: ',A(1),'  E: ',A(2),'  P: ',A(3),'  T: ',A(4),
D    2      '  I: ',A(5)/RPDEG,' ASC: ',A(6)/RPDEG,'  ARG: ',
D    3      A(7)/RPDEG
C     WRITE(PR,'(A)')
C    1  '   T(JD)     E      V    R_TRUE  R_EST  TH_EST  DRHO/Detc.'

C     WRITE(PR,'(A,14(6X,I1,6X))') ' I', (LISTA(J),LISTA(J), J=1,MFIT)
      DO I=1,NPTS
          IF(.NOT.DELETE(I)) THEN

              CALL TRUE2APP(TOBS(I),A,RMOD,THMOD,RAMOD,DECMOD,DYDA,PR)
              R_EST(I)=RMOD
              THETA_EST(I)=THMOD
C             WRITE(PR,'(1X,F8.3,A,1X,F8.4,1X,F7.3,14(1X,G11.3))')
C    1            TOBS(I)-2440000.,' IN MRQCOF AFTER TRUE2APP',
C    2            R_EST(I), THETA_EST(I),
C    3            DYDA(1,1),DYDA(2,1),DYDA(3,1),DYDA(4,1),DYDA(5,1),
C    4            DYDA(6,1),DYDA(7,1),DYDA(1,2),DYDA(2,2),DYDA(3,2),
C    5            DYDA(4,2),DYDA(5,2),DYDA(6,2),DYDA(7,2)
C             DO J=1,MFIT
C                 WRITE(PR,'(A,I1,2(2X,G11.4))') 
C    1                ' LISTA(J), DYDA(LISTA(J),1&2): ',
C    2                LISTA(J),DYDA(LISTA(J),1),DYDA(LISTA(J),2)
C             END DO
              SIG2RA=1./(DRA(I)*DRA(I))
              SIG2D=1./(DDEC(I)*DDEC(I))
              THOBS_PHI=THETAOBS(I)-PHI_ERR(I)+PI/2.
              THEST_PHI=THETA_EST(I)-PHI_ERR(I)+PI/2.
              SINO_PHI=SIN(THOBS_PHI)
              SINE_PHI=SIN(THEST_PHI)
              COSO_PHI=COS(THOBS_PHI)
              COSE_PHI=COS(THEST_PHI)
              DELRA(I) =ROBS(I)*SINO_PHI - R_EST(I)*SINE_PHI
              DELDEC(I)=ROBS(I)*COSO_PHI - R_EST(I)*COSE_PHI
C=======================================================================
C
C This code uses the variable PHI_ERR to describe the position angle
C of the error ellipse for a measured separation and position angle.
C
C Rotate derivatives to angle of error ellipse

              DO J=1,MFIT
                  DRDA=DYDA(LISTA(J),1)
C                 WRITE(PR,'(A,I1,A,G11.4)') ' J=',J,' ; DRDA=',DRDA
C                 ----Deriv. of R w.r.t. an orbital element
                  DTHDA=DYDA(LISTA(J),2)
C                 WRITE(PR,'(A,I1,A,G11.4)') ' J=',J,' ; DTHDA=',DTHDA
C                 ----Deriv. of Theta w.r.t. an orbital element
C                 WRITE(PR,'(A,I1,A,3(2X,G11.4))') ' J=',J,
C    1                ' ; THEST_PHI, SIN, COS=',THEST_PHI/RPDEG,
C    2                SINE_PHI,COSE_PHI
C                 WRITE(PR,'(A,I1,A,G11.4)') ' J=',J,' ; R_EST=',
C    1                R_EST(I)
                  DYDA(LISTA(J),1) = DRDA * SINE_PHI
     1                     + DTHDA * R_EST(I) * COSE_PHI
                  DYDA(LISTA(J),2) = DRDA * COSE_PHI
     1                     - DTHDA * R_EST(I) * SINE_PHI
                  WTRA=DYDA(LISTA(J),1)*SIG2RA
                  WTD =DYDA(LISTA(J),2)*SIG2D
                  DO K=1,J
                      ALPHA(J,K)=ALPHA(J,K)+WTRA*DYDA(LISTA(K),1)
     1                      +WTD*DYDA(LISTA(K),2)
                  END DO
                  BETA(J)=BETA(J)+DELRA(I)*WTRA+DELDEC(I)*WTD
              END DO
C             WRITE(PR,'(I2,14(2X,G11.4))') 
C    1            I,(DYDA(LISTA(J),1),DYDA(LISTA(J),2), J=1,MFIT)
          END IF
      END DO  
D     WRITE(PR,*) 'MRQCOF:  -0.5*GRAD(CHI^2)'
D     WRITE(PR,'(7(2X,G15.10))') (BETA(J), J=1,MFIT)
C     WRITE(PR,*) 'MRQCOF GOING TO CHISQR; NPTS = ',NPTS
      CALL CHISQR(R_EST,ROBS,THETA_EST,THETAOBS,DRA,DDEC,
     1            PHI_ERR,NPTS,DELETE,CHISQ,DELRA,DELDEC)
D     WRITE(PR,'(A,G15.10)') 
D    1    ' MRQCOF:  BACK FROM CHISQR; CHISQ = ',CHISQ
      DO J=2,MFIT
          DO K=1,J-1
              ALPHA(K,J)=ALPHA(J,K)
          END DO  
      END DO  

      RETURN
      END

