      SUBROUTINE VELFIT(LISTA,MFIT,COVAR)
C-------------------------------------------------------------
C DOES LEAST-SQUARES FIT TO TRUE ORBIT,
C THEN FINDS ONE-SIGMA ERROR POINTS FOR EACH FITTED PARAMETER.
C-------------------------------------------------------------
      IMPLICIT UNDEFINED (A-Z)
      INCLUDE 'PARAMETR.INC'
      INCLUDE 'ELLIPSE.INC'

      CHARACTER*1  CHOICE
      CHARACTER*80 INS

      LOGICAL      SINGLSTP

      INTEGER*4    NDATA, I, J, K, KK, ITST, NCA
      INTEGER*4    MFIT, LISTA(7), INITMRQ, NTRY

      REAL*8       ALAMDA, RMOD, THMOD, RAMOD, DECMOD
      REAL*8       OCHISQ, DELCHISQ, NCHISQ, LOOPFACT
      REAL*8       TORB
      REAL*8       A(7), DYDA(7,2)
      REAL*8       COVAR(7,7), ALPHA(7,7), UNCERT(7)
      REAL*8       DA_TRUE, DE_TRUE, DPERIOD, DEPOCH, DINCLIN
      REAL*8       DASC_NODE, DARG_PERI
                   
      PARAMETER (NCA = 7)

C DECIDE WHETHER TO USE THE TRUE ORBITAL ELEMENTS GIVEN IN 
C KEYIN ("GIVEN") OR THE CALCULATED ELEMENTS.  CALCULATED ELEMENTS
C COME FROM CALCTI (ELEMENTS CALCULATED BY THE THIELE-INNES METHOD),
C BUT THEY ARE OVERWRITTEN BY RESULTS OF PREVIOUS INVOCATIONS 
C OF ORBFIT (COMMAND 'L').

C     WRITE(PR,*) 'ORBFIT:  NPTS = ',NPTS

      WRITE(INS,1001) A_TRUE, VALS(53)
      CALL LOCATE(0,2)
      CALL EGA_RESTORE_DEFAULT
      CALL WRITE_STRING( INS )
      CALL PGCURSE(XR,YR,CHOICE)
      CALL CAPS(CHOICE)
      IF      (CHOICE.EQ.'K') THEN
          A_TRUE = VALS(53)
      ELSE IF (CHOICE.EQ.'O') THEN
          CALL LOCATE(0,2)
          CALL WRITE_STRING('                                         ')
          CALL LOCATE(0,2)
          CALL WRITE_STRING('Value for semimajor axis: ')
          READ(INC,*) A_TRUE
      END IF

      WRITE(INS,1002) E_TRUE, VALS(54)
      CALL LOCATE(0,2)
      CALL EGA_RESTORE_DEFAULT
      CALL WRITE_STRING( INS )
      CALL PGCURSE(XR,YR,CHOICE)
      CALL CAPS(CHOICE)
      IF      (CHOICE.EQ.'K') THEN
          E_TRUE = VALS(54)
      ELSE IF (CHOICE.EQ.'O') THEN
          CALL LOCATE(0,2)
          CALL WRITE_STRING('                                         ')
          CALL LOCATE(0,2)
          CALL WRITE_STRING('Value for ellipticity: ')
          READ(INC,*) E_TRUE
      END IF

      WRITE(INS,1003) PERIOD, VALS(55)
      CALL LOCATE(0,2)
      CALL EGA_RESTORE_DEFAULT
      CALL WRITE_STRING( INS )
      CALL PGCURSE(XR,YR,CHOICE)
      CALL CAPS(CHOICE)
      IF      (CHOICE.EQ.'K') THEN
          PERIOD = VALS(55)
      ELSE IF (CHOICE.EQ.'O') THEN
          CALL LOCATE(0,2)
          CALL WRITE_STRING('                                         ')
          CALL LOCATE(0,2)
          CALL WRITE_STRING('Value for period: ')
          READ(INC,*) PERIOD
      END IF

      WRITE(INS,1004) EPOCH, VALS(49)
      CALL LOCATE(0,2)
      CALL EGA_RESTORE_DEFAULT
      CALL WRITE_STRING( INS )
      CALL PGCURSE(XR,YR,CHOICE)
      CALL CAPS(CHOICE)
      IF      (CHOICE.EQ.'K') THEN
          EPOCH = VALS(49)
      ELSE IF (CHOICE.EQ.'O') THEN
          CALL LOCATE(0,2)
          CALL WRITE_STRING('                                         ')
          CALL LOCATE(0,2)
          CALL WRITE_STRING('Value for epoch: ')
          READ(INC,*) EPOCH
      END IF

      WRITE(INS,1005) INCLIN/RPDEG, VALS(50)
      CALL LOCATE(0,2)
      CALL EGA_RESTORE_DEFAULT
      CALL WRITE_STRING( INS )
      CALL PGCURSE(XR,YR,CHOICE)
      CALL CAPS(CHOICE)
      IF      (CHOICE.EQ.'K') THEN
          INCLIN = VALS(50)*RPDEG
      ELSE IF (CHOICE.EQ.'O') THEN
          CALL LOCATE(0,2)
          CALL WRITE_STRING('                                         ')
          CALL LOCATE(0,2)
          CALL WRITE_STRING('Value for inclination (deg): ')
          READ(INC,*) INCLIN
          INCLIN = INCLIN*RPDEG
      END IF

      WRITE(INS,1006) ASC_NODE/RPDEG, VALS(52)
      CALL LOCATE(0,2)
      CALL EGA_RESTORE_DEFAULT
      CALL WRITE_STRING( INS )
      CALL PGCURSE(XR,YR,CHOICE)
      CALL CAPS(CHOICE)
      IF      (CHOICE.EQ.'K') THEN
          ASC_NODE = VALS(52)*RPDEG
      ELSE IF (CHOICE.EQ.'O') THEN
          CALL LOCATE(0,2)
          CALL WRITE_STRING('                                         ')
          CALL LOCATE(0,2)
          CALL WRITE_STRING('Value for ascending node (deg): ')
          READ(INC,*) ASC_NODE
          ASC_NODE = ASC_NODE*RPDEG
      END IF

      WRITE(INS,1007) ARG_PERI/RPDEG, VALS(51)
      CALL LOCATE(0,2)
      CALL EGA_RESTORE_DEFAULT
      CALL WRITE_STRING( INS )
      CALL PGCURSE(XR,YR,CHOICE)
      CALL CAPS(CHOICE)
      IF      (CHOICE.EQ.'K') THEN
          ARG_PERI = VALS(51)*RPDEG
      ELSE IF (CHOICE.EQ.'O') THEN
          CALL LOCATE(0,2)
          CALL WRITE_STRING('                                         ')
          CALL LOCATE(0,2)
          CALL WRITE_STRING('Value for arg. of periastron (deg): ')
          READ(INC,*) ARG_PERI
          ARG_PERI = ARG_PERI*RPDEG
      END IF

      NDATA           = 2 * (NPTS-DEL_COUNT)

      A(1)            = A_TRUE
      A(2)            = E_TRUE
      A(3)            = PERIOD
      A(4)            = EPOCH
      A(5)            = INCLIN
      A(6)            = ASC_NODE
      A(7)            = ARG_PERI

C GET STARTING CHI SQUARED
      DO I = 1, NPTS
          CALL TRUE2APP(TOBS(I),A,RMOD,THMOD,RAMOD,DECMOD,DYDA,PR)
          R_EST(I) = RMOD
          THETA_EST(I) = THMOD
      END DO
      CALL CHISQR(R_EST,ROBS,THETA_EST,THETAOBS,DRA,DDEC,PHI_ERR,
     1            NPTS,DELETE,OCHISQ,DELRA,DELDEC)
D     WRITE(PR,'(A,F8.4)') ' ORBFIT: STARTING CHI^2 = ',OCHISQ

      WRITE(PR,*) 'ORBFIT INITIALIZATION:'
      WRITE(PR,'(A,F8.4)')  '   A_TRUE   =        ', A(1)
      WRITE(PR,'(A,F7.4)')  '   E_TRUE   =         ',A(2)
      WRITE(PR,'(A,F9.4)')  '   PERIOD   =       ',  A(3)
      WRITE(PR,'(A,F12.4)') '   EPOCH    = JD ',     A(4) 
      WRITE(PR,'(A,F9.4)')  '   INCLIN   =       ',  A(5)/RPDEG
      WRITE(PR,'(A,F9.4)')  '   ASC_NODE =       ',  A(6)/RPDEG
      WRITE(PR,'(A,F9.4)')  '   ARG_PERI =       ',  A(7)/RPDEG
      WRITE(PR,*) '  '

      WRITE(INS,1009) A(1),A(2),A(3),A(4)
      CALL LOCATE(0,2)
      CALL WRITE_STRING ( INS )
      WRITE(INS,1012) A(5)/RPDEG,A(6)/RPDEG,A(7)/RPDEG
      CALL LOCATE(0,3)
      CALL WRITE_STRING ( INS )
      CALL LOCATE(0,4)
      CALL WRITE_STRING ('Hit any key                                 ')
      CALL PGCURSE(XR,YR,CHOICE)

      CALL LOCATE(0,2)
      CALL WRITE_STRING
     1 (' PARAMETERS TO FIT ARE:                                      ')
      CALL LOCATE(0,3)
      CALL WRITE_STRING('    1  SEMIMAJOR AXIS                        ')
      CALL LOCATE(0,4)
      CALL WRITE_STRING('    2  ELLIPTICITY')
      CALL LOCATE(0,5)
      CALL WRITE_STRING('    3  PERIOD')
      CALL LOCATE(0,6)
      CALL WRITE_STRING('    4  EPOCH')
      CALL LOCATE(0,7)
      CALL WRITE_STRING('    5  ORBITAL INCLINATION')
      CALL LOCATE(0,8)
      CALL WRITE_STRING('    6  ANGLE OF ASCENDING NODE')
      CALL LOCATE(0,9)
      CALL WRITE_STRING('    7  ARGUMENT OF PERIASTRON')
C      CALL LOCATE(0,10)
C      CALL WRITE_STRING('                              ')
      CALL LOCATE(0,10)
      CALL WRITE_STRING('PARAMETERS TO FIT (e.g., 126 fits '//
     1      'axis, ellipticity, and asc. node): ')
      DO I = 1,7
          LISTA(I)    = 0
      END DO
      MFIT            = 0
      READ(INC,'(7I1)') LISTA
      DO I = 1,7
          IF(LISTA(I).NE.0) MFIT = MFIT+1
      END DO
      WRITE(INS,'(A,I1,A,F7.2)') 
     1      ' FIT ',MFIT,' PARAMETERS; STARTING CHI^2 = ',OCHISQ
      CALL LOCATE(0,11)
      CALL WRITE_STRING ( INS )

      CALL LOCATE(0,12)
      CALL WRITE_STRING('FRACTIONAL CHANGE IN CHISQ FOR CONVERGENCE: ')
      READ(INC,*) DELCHISQ

      SINGLSTP        = .FALSE.

   10 INITMRQ         = -1            ! RETURN HERE FOR REPEAT

      CALL LOCATE(1,13)
      CALL WRITE_STRING('NEWLINENEWLINENEWLINENEWLINENEWLINENEWL    ')
      CALL LOCATE(1,13)
      CALL WRITE_STRING('NEWLINENEWLINENEWLINENEWLINENEWLINENEWL    ')
      CALL LOCATE(1,13)
      CALL WRITE_STRING('INIT. VAL. OF MARQUARDT COEFF. ALAMDA: ')
      READ(INC,*) ALAMDA
   11 CALL LOCATE(0,14)
      CALL WRITE_STRING
     1             ('NEWLINENEWLINENEWLINENEWLINENEWLINENEWLINENEW    ')
      CALL LOCATE(0,14)
      CALL WRITE_STRING
     1             ('NEWLINENEWLINENEWLINENEWLINENEWLINENEWLINENEW    ')
      CALL LOCATE(0,14)
      CALL WRITE_STRING
     1             ('NUMBER TO MULT. ALAMDA BY AS MARQUARDT METHOD')
      CALL LOCATE(0,15)
      CALL WRITE_STRING('NEWLINENEWLINENEWLINENEWLINENEWLINENEWL    ')
      CALL LOCATE(0,15)
      CALL WRITE_STRING('NEWLINENEWLINENEWLINENEWLINENEWLINENEWL    ')
      CALL LOCATE(0,15)
      CALL WRITE_STRING('GOES OVER TO PARABOLIC CHISQ FIT: [<1] ')
      READ(INC,*) LOOPFACT
      IF (LOOPFACT.GE.1.) GO TO 11

      CALL LOCATE(0,16)
      CALL WRITE_STRING('                                           ')
      CALL LOCATE(0,16)
      CALL WRITE_STRING('SINGLE STEP [N]? ')
      READ(INC,*) ANSWER
      CALL CAPS(ANSWER)
      IF (ANSWER.EQ.'Y') THEN
          SINGLSTP    = .TRUE.
      ELSE
          SINGLSTP    = .FALSE.
      END IF

D     WRITE(PR,*) 'ORBFIT:  GOING TO MRQMIN FIRST TIME'
      CALL MRQMIN(NDATA,A,LISTA,MFIT,COVAR,ALPHA,
     1            NCA,ALAMDA,LOOPFACT,INITMRQ)
D     WRITE(PR,*) 'ORBFIT:  BACK FROM MRQMIN'

      K               = 1
      NTRY            = 1
      ITST            = 0

   20 K               = K+1     ! LOOP TO HERE IF NO CONVERGENCE     
      OCHISQ          = CHISQ

      IF (SINGLSTP .AND. NTRY.GT.1) THEN
          CALL LOCATE(0,16)
          CALL WRITE_STRING
     1            ('                                           ')
          CALL LOCATE(0,16)
          CALL WRITE_STRING('SINGLE STEP [N]? ')
          READ(INC,*) ANSWER
          CALL CAPS(ANSWER)
          IF (ANSWER.EQ.'Y') THEN
              SINGLSTP    = .TRUE.
          ELSE
              SINGLSTP    = .FALSE.
          END IF
      END IF

      KK              = MOD(K,6)
      CALL LOCATE(0,17+KK)
      WRITE(INS,1013) CHISQ 
      CALL WRITE_STRING ( INS )
C      CALL LOCATE(0,24)
C      WRITE(INS,1015) NTRY
C      CALL WRITE_STRING ( INS )

C CALL MRQMIN AND OBSERVE CHANGE IN CHISQ

D     WRITE(PR,'(A,2(G11.4,1X))') ' ORBFIT:  CHISQ, OCHISQ = ',
D    1      CHISQ,OCHISQ
D     WRITE(PR,*) 'ORBFIT:  GOING TO MRQMIN AFTER STMT 20'

      CALL MRQMIN(NDATA,A,LISTA,MFIT,COVAR,ALPHA,
     1            NCA,ALAMDA,LOOPFACT,INITMRQ)
      NTRY            = NTRY + 1

D     WRITE(PR,'(A,I3)') ' ORBFIT:  BACK FROM MRQMIN; NTRY = ',NTRY
D     WRITE(PR,'(A,F8.4,2X,F8.4)')
D    1      ' ORBFIT BEFORE "IF":  CHISQ, OCHISQ = ',CHISQ,OCHISQ
D     WRITE(PR,'(A,G10.4,2X,G10.4)')
D    1      ' ORBFIT:  ABS(OCHISQ/CHISQ-1.), DELCHISQ = ',
D    2      ABS(OCHISQ/CHISQ-1.), DELCHISQ

      IF(SINGLSTP) THEN
          WRITE(INS,1010) A(1),A(2),A(3),A(4)
          CALL LOCATE(0,2)
          CALL WRITE_STRING ( INS )
          WRITE(INS,1012) A(5)/RPDEG, A(6)/RPDEG, A(7)/RPDEG
          CALL LOCATE(0,3)
          CALL WRITE_STRING ( INS )
          CALL LOCATE(0,4)
          CALL WRITE_STRING('Hit any key to continue fitting')
          CALL PGCURSE(XR,YR,CHOICE)
      END IF

C      IF ((CHISQ/OCHISQ-1.).GT.DELCHISQ/10.) THEN
C          ITST        = 0
      IF (CHISQ.GE.OCHISQ) THEN
          ITST        = 0
      ELSE
          IF (ABS(OCHISQ/CHISQ-1.).LT.DELCHISQ) THEN
              ITST    = ITST + 1
          END IF
          OCHISQ      = CHISQ
      END IF

D     WRITE(PR,'(A,F8.4,2X,F8.4,2X,I1)')
D    1      ' ORBFIT AFTER "IF":  CHISQ, OCHISQ,ITST = ',
D    2        CHISQ,OCHISQ,ITST

      IF (.NOT.SINGLSTP .AND. MOD(NTRY,10).EQ.0) THEN
          CALL LOCATE(0,23)
          WRITE(INS,1016) NTRY
          CALL WRITE_STRING ( INS )
          READ(INC,*) CHOICE
C          CALL EGA_RESTORE_DEFAULT
          IF ((.NOT.(CHOICE.EQ.'N'.OR.CHOICE.EQ.'n'))
     1        .AND.(ITST.LT.2)) THEN
              GO TO 20
          ELSE
              GO TO 25
          END IF
      END IF

      IF (ITST.LT.2) GO TO 20

C WHEN CHISQ HAS CONVERGED, CALL MRQMIN AGAIN WITH INITMRQ = 1
C TO GET COVARIANCE AND CURVATURE MATRICES

   25 INITMRQ         = 1 

D     WRITE(PR,*) 'ORBFIT:  GOING TO MRQMIN LAST TIME'
      CALL MRQMIN(NDATA,A,LISTA,MFIT,COVAR,ALPHA,
     1            NCA,ALAMDA,LOOPFACT,INITMRQ)
D     WRITE(PR,*) 'ORBFIT:  BACK FROM MRQMIN'

      CALL LOCATE(0,17+KK)
      CALL WRITE_STRING ('                             ')
      CALL LOCATE(0,17+KK)
      CALL WRITE_STRING ('                REPEAT? [N]')
      READ(INC,*) COMND
      CALL EGA_RESTORE_DEFAULT
      IF (COMND.EQ.'Y'.OR.COMND.EQ.'y') GO TO 10

      ABEST(1)        = A(1)
      ABEST(2)        = A(2)
      ABEST(3)        = A(3)
      ABEST(4)        = A(4)
      ABEST(5)        = A(5)
      ABEST(6)        = A(6)
      ABEST(7)        = A(7)

C PUT THE RESULTS OUT IN A WAY THAT CAN BE READ

30    A_TRUE          = ABEST(1)
      E_TRUE          = ABEST(2)
      PERIOD          = ABEST(3)
      EPOCH           = ABEST(4)
      INCLIN          = ABEST(5)
      ASC_NODE        = ABEST(6)
      ARG_PERI        = ABEST(7)

      DA_TRUE         = COVAR(1,1)
      DE_TRUE         = COVAR(2,2)
      DPERIOD         = COVAR(3,3)
      DEPOCH          = COVAR(4,4)
      DINCLIN         = COVAR(5,5)
      DASC_NODE       = COVAR(6,6)
      DARG_PERI       = COVAR(7,7)

      WRITE(PR,*) 'ORBFIT: '

      WRITE(PR,'(1X,A,F9.4,A,F8.5)')
     1 'SEMIMAJOR AXIS  = ',A_TRUE,' +/- ',DA_TRUE

      WRITE(PR,'(1X,A,F7.4,A,F10.8)')
     1 'ELLIPTICITY     = ',E_TRUE,' +/- ',DE_TRUE

      WRITE(PR,'(1X,A,F10.5,A,F9.6)')
     1 'PERIOD          = ',PERIOD,' +/- ',DPERIOD

      WRITE(PR,'(1X,A,F12.4,A,F8.5)')
     1 'EPOCH           = JD ',EPOCH ,' +/- ',DEPOCH 

      WRITE(PR,'(1X,A,F9.4,A,F8.5)')
     1 'INCLINATION     = ',INCLIN/RPDEG,' +/- ',DINCLIN/RPDEG
                          
      WRITE(PR,'(1X,A,F9.4,A,F8.5)')
     1 'ASCENDING NODE  = ',ASC_NODE/RPDEG,' +/- ',DASC_NODE/RPDEG

      WRITE(PR,'(1X,A,F9.4,A,F8.5)')
     1 'ARG. PERIASTRON = ',ARG_PERI/RPDEG,' +/- ',DARG_PERI/RPDEG

      WRITE(PR,'(/1X,A,F10.5,A,I3,A)')
     1 'CHI SQUARED = ',CHISQ,' FOR ',NDATA-MFIT,' DEGREES OF FREEDOM'

      WRITE(PR,*) 'COVARIANCE MATRIX:'
      DO J=1,7
          WRITE(PR,1014) (COVAR(J,K), K=1,7)
      END DO


C PREPARE TO DISPLAY THE RESULTANT ELLIPSE ON THE SCREEN

      TORB             = EPOCH

C      WRITE(PR,'(A)')
C     1  ' I    TORB-EPOCH     RMOD THMOD    RA_MOD   DEC_MOD'

      IF(ORBPTS.GT.ELLPTS) ORBPTS = ELLPTS      ! Arrays are dimensioned 
                                                ! with ELLPTS
      DO I = 1, ORBPTS
          TORB         = TORB + PERIOD/(ORBPTS)
          CALL TRUE2APP(TORB,A,RMOD,THMOD,RAMOD,DECMOD,DYDA,PR)
          RA_MOD(I)    = RAMOD
          DEC_MOD(I)   = DECMOD
C          WRITE(PR,'(1X,I3,2X,F11.3,2X,F6.2,2X,F6.2,2X,F6.2,2X,F6.2)')
C     1      I,TORB,RMOD,THMOD,RA_MOD(I),DEC_MOD(I)
      END DO

      WRITE(INS,1011) A(1),A(2),A(3),A(4)
      CALL LOCATE(0,2)
      CALL WRITE_STRING ( INS )
      WRITE(INS,1012) A(5)/RPDEG,A(6)/RPDEG,A(7)/RPDEG
      CALL LOCATE(0,3)
      CALL WRITE_STRING ( INS )
      CALL LOCATE(0,4)
      CALL WRITE_STRING ('Hit any key                                 ')
      CALL PGCURSE(XR,YR,CHOICE)

      TRUFIT           = .TRUE.

      RETURN

 1001 FORMAT(1X,'SEMIMAJ: [C]alculated (',F9.4,'), [K]EYIN (',
     1          F9.4,') or [O]ther value? [C]  ')      ! 78 characters
 1002 FORMAT(1X,'ELLIPTICITY: [C]alc. (',F9.4,'), [K]EYIN (',
     1          F9.4,') or [O]ther value? [C]   ')     ! 78 characters
 1003 FORMAT(1X,'PERIOD: [C]alculated (',F9.4,'), [K]EYIN (',
     1          F9.4,') or [O]ther value? [C]   ')     ! 78 characters
 1004 FORMAT(1X,'EPOCH: [C]alculated (',F11.3,'), [K]EYIN (',
     1          F11.3,') or [O]ther value? [C]')       ! 78 characters
 1005 FORMAT(1X,'INCLIN: [C]alculated (',F7.2,'), [K]EYIN (',
     1          F7.2,') or [O]ther value? [C]       ') ! 78 characters
 1006 FORMAT(1X,'ASC. NODE: [C]alculated (',F7.2,'), [K]EYIN (',
     1          F7.2,') or [O]ther value? [C]    ')    ! 78 characters
 1007 FORMAT(1X,'ARG. PERI: [C]alculated (',F7.2,'), [K]EYIN (',
     1          F7.2,') or [O]ther value? [C]    ')    ! 78 characters
 1008 FORMAT(1X,'PERIOD: (K)EYIN = ',F10.5,' or (D)atafile = ',
     1       F10.5,' VALUE? [C]     ')
 1009 FORMAT(1X,'INPUT ELEMENTS: ',F8.4,1X,F7.4,1X,F9.4,1X,F12.4)
 1010 FORMAT(1X,'CURRENT ELEMENTS: ',F8.4,1X,F7.4,1X,F9.4,1X,F12.4)
 1011 FORMAT(1X,'FINAL ELEMENTS: ',F8.4,1X,F7.4,1X,F9.4,1X,F12.4)
 1012 FORMAT(1X,'                ',F9.4,1X,F9.4,1X,F9.4)
 1013 FORMAT(1X,'ORBFIT: CHISQ = ',F11.3)
 1014 FORMAT(1X,7(E9.3,2X))
 1015 FORMAT(1X,'NUMBER OF TRIALS: ',I4)
 1016 FORMAT(1X,'NUMBER OF TRIALS: ',I4,'; KEEP TRYING [Y]?')

      END

