      SUBROUTINE FNDMED(XX,N,DELMAX)
C THIS SUBROUTINE FINDS THE MEDIAN IF THE N ELEMENTS OF THE ARRAY X.
C IF N IS EVEN THE MEDIAN IS TAKEN TO BE THE MEAN OF THE TWO MIDDLE
C ELEMENTS.
      DIMENSION XX(1)
      INCLUDE 'RFINE.PRM'
      INCLUDE 'E.INC'
      INCLUDE 'FBUF.INC'
C
      IF(N.GT.0)GO TO 100
      WRITE(IOUT,8)
    8 FORMAT('0NULL LIST HAS NO MEDIAN. RUN TERMINATED.')
      STOP
C
C SORT THE LIST
C
  100 CALL QCKSRT(XX,N,FBUF)
C
C CALCULATE NORMAL PROBABILITY PLOT INFORMATION
C
      ILOW=0.05*N
      IHIGH=0.95*N
      NO=IHIGH-ILOW+1
      SUMX=0.0
      SUMX2=0.0
      SUMXY=0.0
      SUMY=0.0
      SUMY2=0.0
      DO 200 I=ILOW,IHIGH
        PROB=FLOAT(2*(N-I)+1)/FLOAT(2*N)
        X=-EXPECT(PROB)
  	Y=XX(I)
	SUMX=SUMX+X
	SUMX2=SUMX2+X*X
	SUMY=SUMY+Y
        SUMY2=SUMY2+Y*Y
	SUMXY=SUMXY+X*Y
200   CONTINUE
      DET=NO*SUMX2-SUMX**2
      A=(SUMX2*SUMY-SUMX*SUMXY)/DET
      B=(-SUMX*SUMY+NO*SUMXY)/DET
      SUMD2=0.0
      SUMX=0.0
      DO 300 I=ILOW,IHIGH
        PROB=FLOAT(2*(N-I)+1)/FLOAT(2*N)
        X=-EXPECT(PROB)
  	Y=XX(I)
	EY=A+B*X
	SUMD2=SUMD2+(EY-Y)**2
300   CONTINUE
      SSE=SUMD2/(NO-2)
      SIGA=SQRT(SSE*SUMX2/DET)
      SIGB=SQRT(SSE*NO/DET)
      R2=(SUMY2-SUMD2)/SUMY2
C
C CALCULATE THE MEDIAN
C
      DO 400 I=1,N
        XX(I)=ABS(XX(I))
400   CONTINUE
      CALL QCKSRT(XX,N,FBUF)
      IMID=N/2+1
      IF(MOD(N,2).EQ.0)THEN
        XMED=0.5*(XX(IMID-1)+XX(IMID))
      ELSE
        XMED=XX(IMID)
      ENDIF
      DELMAX=9.*XMED
C
C PRINT OUT RESULTS
C
      WRITE(IOUT,10)XMED,DELMAX
   10 FORMAT('0Median DELTAF/SIGMA is',F9.4,' New DELMAX is',F9.4)
      WRITE(IOUT,12)NO,B,SIGB,A,SIGA,R2
12    FORMAT(/'0Normal Probability Plot Regression'
     1//' No. Data    Slope  Sigma    Intercept  Sigma    R^2'
     2//I8,4F9.3,F8.3//' If at convergence the SLOPE deviates '
     3,'from 1, make inverse adjustment to "IGNORANCE FACTOR".')
      RETURN
      END
      FUNCTION EXPECT(PROB)
C
C FUNCTION TO CALCULATE EXPECTED PROBABILITY FOR NORMAL PROBABILITY PLOTS.
C
C CODED BY BOB DOWNS, VPI & SU, APRIL 1991
C
      TOL=0.00000001
      TOLX=0.00000001
      XINC=0.100
      X=0.0
      ICNT=0
999   E=ERFCC(X*SQRT(2.0))
      G=.5*(2.0-E)
      IF(ABS(PROB-G).LT.TOL)GOTO 998
      IF(XINC.LT.TOLX)GOTO 998
      IF(PROB.LT.G)THEN
         IF(IFLAG.EQ.1)XINC=XINC/2.0
         X=X-XINC
         IFLAG=-1
      ELSE
         IF(IFLAG.EQ.-1)XINC=XINC/2.0
         X=X+XINC
         IFLAG=1
      ENDIF
      IF(ICNT.GT.500)GOTO 998
      ICNT=ICNT+1
      GOTO 999
998   EXPECT=2.0*X
      RETURN
      END
      FUNCTION ERFCC(X)
      Z=ABS(X)
      T=1./(1.+0.5*Z)
      ERFCC=T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(.37409196+
     *    T*(.09678418+T*(-.18628806+T*(.27886807+T*(-1.13520398+
     *    T*(1.48851587+T*(-.82215223+T*.17087277)))))))))
      IF (X.LT.0.) ERFCC=2.-ERFCC
      RETURN
      END
      SUBROUTINE QCKSRT(ARRAY,NUMBER,ISTACK)
C
C QUICK SORT ROUTINE TO SORT DELTA/SIGMA VALUES
C
      DIMENSION ARRAY(1),ISTACK(2,1)
C
C INITIALIZE STACK AND LEFT AND RIGHT POINTERS
C
      NSTACK=0
      IL=1
      IR=NUMBER
      GO TO 200
C
C UNSTACK TOP ENTRY ON LIST - DONE IF STACK EMPTY
C
100   IF(NSTACK.EQ.0)RETURN
      IL=ISTACK(1,NSTACK)
      IR=ISTACK(2,NSTACK)
      NSTACK=NSTACK-1
C
C CHOOSE LEFT ITEM AS PIVOT AND SUBDIVIDE LIST
C
200   I=IL+1
      J=IR
      PIV=ARRAY(IL)
300   IF(ARRAY(I).LE.PIV.AND.I.LE.IR)THEN
        I=I+1
        GO TO 300
      END IF
310   IF(ARRAY(J).GT.PIV)THEN
        J=J-1
        GO TO 310
      END IF
      IF(I.LT.J) THEN
        XT=ARRAY(I)
        ARRAY(I)=ARRAY(J)
        ARRAY(J)=XT
      END IF
      IF(I.LE.J)GO TO 300
      ARRAY(IL)=ARRAY(J)
      ARRAY(J)=PIV
      IF(IL.LT.J-1)THEN
C STACK LEFT SIDE POINTERS
        NSTACK=NSTACK+1
        ISTACK(1,NSTACK)=IL
        ISTACK(2,NSTACK)=J-1
      END IF
      IF(I.LT.IR)THEN
C SET UP RIGHT SIDE FOR NEXT PASS
        IL=I
        GO TO 200
      ELSE
C RIGHT SIDE EMPTY - GO TO UNSTACK A LEFT SIDE
        GO TO 100
      END IF
      END
      SUBROUTINE MEDWT(DELMAX)
C
C THIS SUBROUTINE DETERMINES THE MEDIAN VALUE OF THE FUNCTION
C ABSM((FOBS-FCALC)/SIGMA)
C FROM THE PREVIOUS REFINEMENT CYCLE. THIS VALUE IS MULTIPLIED BY
C 9, AND THE RESULT IS PLACED IN DELMAX.
C
C IT ALSO CALCULATES THE SLOPE AND INTERCEPT OF THE NORMAL PROBABILITY
C PLOT TO CHECK THE CONSISTENCY OF THE WEIGHTS.
C
      INCLUDE 'RFINE.PRM'
      INCLUDE 'LSMAT.INC'
      INCLUDE 'E.INC'
      INCLUDE 'FBUF.INC'
      DIMENSION RMAT(3,3,24),TRANS(3,24)
      EQUIVALENCE (VOL,RMAT(1,1,1),TRANS(1,1))
C
      DELSIG(FOBS,FCALC,SIGMA)=(FOBS-FCALC)/SIGMA
C
      REWIND ISC4
      READ(ISC4)NSYM,VOL,(((RMAT(I,J,K),I=1,3),J=1,3),K=1,NSYM),((TRANS
     1(J,I),J=1,3),I=1,NSYM)
      NREF=0
  300 READ(ISC4)FBUF
      M=1
  310 IF(IFBUF(M).EQ.99)GO TO 330
      IF(IFBUF(M+7).NE.0)GO TO 320
      NREF=NREF+1
      A(MDEGN+NREF)=DELSIG(FBUF(M+3),FBUF(M+4),FBUF(M+8))
  320 M=M+10
      IF(M-500)310,310,300
C     ALL REFLECTIONS READ. NOW SORT THEM.
  330 IX=MDEGN+1
      CALL FNDMED(A(IX),NREF,DELMAX)
      RETURN
      END
