*-- Author :  R. Henderson
      SUBROUTINE FPRCHI(PSEG,PERSEG,RSEG,RERSEG,Z,PROBCH)
**: FPRCHI.......SM. Bug fix. Protect against small Chisq for PROB.                                           
**----------------------------------------------------------------------                                      
C                                                                                                             
C---  Routine calculates the probability of PSEG and RSEG                                                     
C---  being the same vector at Z given that PERSEG is the covariance                                          
C---  matrix of PSEG both calculated at at Z=0                                                                
C---  (these are the quantities supplied by FPFYUV)                                                           
C---  and RERSEG is the covariance matrix of RSEG at Z                                                        
C                                                                                                             
C---  The probability fron chisquare is returned as PROB                                                      
C                                                                                                             
      DOUBLE PRECISION COVAR                                            
C---                                                                                                          
      DIMENSION PSEG(4) , PERSEG(4,4)                                   
      DIMENSION RSEG(4) , RERSEG(4,4)                                   
      DIMENSION COVAR(4,4) , DERIV(4,4) , WORK(16) , WM(4,4) , WV(4)    
C                                                                                                             
C---  Construct the derivative matrix                                                                         
C                                                                                                             
      CALL VZERO(DERIV,16)                                                                             
      DO 10 I = 1,4                                                     
      DERIV(I,I) = 1.0                                                  
   10 CONTINUE                                                          
      DERIV(1,3) = Z                                                    
      DERIV(2,4) = Z                                                    
C                                                                                                             
C---  Calculate X and Y at Z                                                                                  
C                                                                                                             
      PSEG(1) = PSEG(1)  +  Z * PSEG(3)                                 
      PSEG(2) = PSEG(2)  +  Z * PSEG(4)                                 
C                                                                                                             
C---  Propagate PERSEG FROM Z=0 TO Z                                                                          
C                                                                                                             
      CALL VZERO(WM,16)                                                                                
      DO 24 I = 1,4                                                     
        DO 25 J = 1,4                                                   
         DO 26 K = 1,4                                                  
         WM(I,J) = WM(I,J) + PERSEG(I,K) * DERIV(J,K)                   
   26    CONTINUE                                                       
   25   CONTINUE                                                        
   24 CONTINUE                                                          
C---                                                                                                          
      CALL VZERO(COVAR,32)                                                                             
      DO 54 M = 1,4                                                     
        DO 55 K = 1,4                                                   
         DO 56 I = 1,4                                                  
         COVAR(M,K) = COVAR(M,K) + DERIV(M,I) * WM(I,K)                 
   56    CONTINUE                                                       
   55   CONTINUE                                                        
   54 CONTINUE                                                          
C---                                                                                                          
C      DO 111 I = 1,4                                                                                         
C      WRITE(6,1001)(COVAR(I,J),J=1,4)                                                                        
C 1001 FORMAT(1X,'COVAR',4e15.6)                                                                              
C  111 CONTINUE                                                                                               
C      WRITE(6,*)'------------------------------'                                                             
C                                                                                                             
C---  Add on the radial covariance matrix AT Z = Z                                                            
C                                                                                                             
      DO 30 I = 1,4                                                     
      DO 31 J = 1,4                                                     
      COVAR(I,J) = COVAR(I,J) + RERSEG(I,J)                             
   31 CONTINUE                                                          
   30 CONTINUE                                                          
C                                                                                                             
C---  Calculate weight matrix                                                                                 
C                                                                                                             
      CALL DINV(4,COVAR,4,WORK,IFAIL)                                                                  
        IF(IFAIL .NE. 0)THEN                                            
          WRITE(6,*)' Matrix failed to invert in FPRCHI'                
          PROBCH = 0.0                                                  
          RETURN                                                        
        ENDIF                                                           
C---                                                                                                          
C      DO 131 I = 1,4                                                                                         
C      WRITE(6,1003)(COVAR(I,J),J=1,4)                                                                        
C 1003 FORMAT(1X,'COVAR',4e15.6)                                                                              
C  131 CONTINUE                                                                                               
C      WRITE(6,*)'------------------------------'                                                             
C                                                                                                             
C---  Now calculate chisquare                                                                                 
C                                                                                                             
      CHISQ = 0.0                                                       
      DO 40 I = 1,4                                                     
      WV(I) = 0.0                                                       
      DO 41 J = 1,4                                                     
      WV(I) = WV(I) +                                                   
     1        COVAR(I,J) * (PSEG(J) - RSEG(J))                          
   41 CONTINUE                                                          
      CHISQ = CHISQ + (PSEG(I) - RSEG(I)) * WV(I)                       
   40 CONTINUE                                                          
C                                                                                                             
C---  Find probability                                                                                        
C                                                                                                             
      IF(CHISQ .LT. 0.001) THEN                                         
       PROBCH = 0.99999                                                 
      ELSE                                                              
       PROBCH = PROB(CHISQ,4)                                           
      ENDIF                                                             
C      WRITE(6,*)PROBCH                                                                                       
C      WRITE(6,*)'------------------------------'                                                             
      END                                                               
*