SUBROUTINE FPCFIT
*-- Author :  R. Henderson
      SUBROUTINE FPCFIT(IP,Y,W,SLOPE,ZERO,COVSLZ,CHISQ,PBCHI)
**: FPCFIT.......SM. Bug fix. Protect against small Chisq for PROB.                                           
**----------------------------------------------------------------------                                      
C-------------------------------------------------------------                                                
C                                                                                                             
C---  Fit four points y(i) with weights w(i) at positions z(i)                                                
C---  Returns SLOPE (dy/dz) ZERO (y(0))                                                                       
C---          COVSLZ(2,2) their correlation matrix                                                            
C---  and     CHISQ and probability from chisquare pbchi                                                      
C                                                                                                             
C-------------------------------------------------------------                                                
C---                                                                                                          
*KEEP,FPLGEO.                                                                                                 
C---                                                                                                          
      COMMON /FPLGEO/   ZPLAN(36)   , TP(9)   , YP(26)    , PLANE(3,9), 
     1                 RMAX    , RMIN    , YSTART    , YSPACE    ,      
     2                 X0      , Y0      , PZSTRU (8), STAGER   ,       
     3                 RESOL   , ACUT    , CTP(9)    , STP(9)           
C---                                                                                                          
*KEND.                                                                                                        
C---                                                                                                          
      DOUBLE PRECISION Y(4),Z(4)                                        
      REAL W(4),COVSLZ(2,2)                                             
C                                                                                                             
C---  Section of code varies with w(i)                                                                        
C                                                                                                             
      Z(1) = 0.0                                                        
      Z(2) = ZPLAN( IP + 1) - ZPLAN( IP )                               
      Z(3) = ZPLAN( IP + 2) - ZPLAN( IP )                               
      Z(4) = ZPLAN( IP + 3) - ZPLAN( IP )                               
      ZSUM = 0.0                                                        
      ZSUM2 = 0.0                                                       
      WSUM = 0.0                                                        
      DO 20 I=1,4                                                       
      ZSUM = ZSUM + Z(I)*W(I)                                           
      ZSUM2 = ZSUM2 + Z(I)*Z(I)*W(I)                                    
      WSUM = WSUM + W(I)                                                
   20 CONTINUE                                                          
      DET = WSUM * ZSUM2  - (ZSUM * ZSUM)                               
C                                                                                                             
C---  return slope and constant unphysical if det is 0                                                        
C                                                                                                             
      IF( DET .EQ. 0.0 )THEN                                            
         SLOPE = 0.0                                                    
         CONST = 2000.0                                                 
      RETURN                                                            
      ENDIF                                                             
C                                                                                                             
C---  calculate error matrix                                                                                  
C                                                                                                             
      COVSLZ(1,1) = WSUM/DET                                            
      COVSLZ(2,2) = ZSUM2/DET                                           
      COVSLZ(1,2) = -ZSUM/DET                                           
      COVSLZ(2,1) = COVSLZ(1,2)                                         
C                                                                                                             
C---  Initialization per fit                                                                                  
C                                                                                                             
      SLOPE = 0.0                                                       
      CONST = 0.0                                                       
      YSUM = 0.0                                                        
      YZSUM = 0.0                                                       
C                                                                                                             
C---  Calculate required sums                                                                                 
C                                                                                                             
      NDF = 0                                                           
      DO 25 I=1,4                                                       
      IF(W(I) .NE. 0.0) NDF = NDF + 1                                   
      YSUM = YSUM + Y(I) * W(I)                                         
      YZSUM = YZSUM + Y(I) * Z(I) * W(I)                                
   25 CONTINUE                                                          
C                                                                                                             
C---  Calculate slopes and data zeros                                                                         
C                                                                                                             
      SLOPE = ( WSUM * YZSUM  -  ZSUM * YSUM ) /DET                     
      ZERO  = ( ZSUM2 * YSUM - ZSUM * YZSUM )  /DET                     
C                                                                                                             
C---  Calculate chisquare                                                                                     
C                                                                                                             
      CHISQ = 0.0                                                       
      DO 23 I=1,4                                                       
      CHISQ =   CHISQ   +                                               
     1         ( Y(I)   -                                               
     2           SLOPE*Z(I) - ZERO )**2   *  W(I)                       
   23 CONTINUE                                                          
C                                                                                                             
C---  Calculate probability from chisquare                                                                    
C                                                                                                             
      NDF = NDF - 2                                                     
*     Fix for v. small Chisq...                                                                               
      IF(CHISQ .LT. 0.001) THEN                                         
       PBCHI = 0.9999999                                                
      ELSE                                                              
       PBCHI = PROB( ABS(CHISQ),NDF )                                   
      ENDIF                                                             
      END                                                               
*