C***********************************************************************
C   SUBROUTINE CHI COMPUTES THE PROBABILITY INTEGRAL TRANSFORMATION OF
C   ANY OBSERVATION OF A CENTRAL CHI-SQUARE VARIABLE.
C
C   CALLING ARGUMENTS:
C   K    THE DEGREES OF FREEDOM
C   Z    THE CHI-SQUARE VARIABLE
C   P    THE PROBABILITY THAT A CENTRAL CHI-SQUARE VARIABLE WITH K
C        DEGREES OF FREEDOM IS LESS THAN THE OBSERVED VALUE OF Z .  P IS
C        RETURNED AS -1. IF EITHER K IS LESS THAN 1 OR Z IS NEGATIVE.
C
C   FOR K AN EVEN INTEGER THE ANSWER CAN BE FOUND IN CLOSED FORM
C   AS THE SUM OF CERTAIN POISSON PROBABILITIES. THIS IS USED ONLY IF
C   K IS LESS THAN OR EQUAL TO 100.
C   FOR K AN ODD INTEGER LESS THAN 100 INTEGRATION BY PARTS IS USED
C   FOLLOWED BY A NUMERICAL EVALUATION OF THE CHI(1) DISTRIBUTION.
C   IF K IS GREATER THAN 100  THE WILSON-HILERTY TRANSFORMATION IS USED
C   TO MAP Z INTO AN APPROXIMATE STANDARD NORMAL VARIABLE(SEE JOHNSON
C   AND KOTZ) THE SQUARE OF THIS NEW VARIABLE IS CHI(1) AND THUS THE
C   DESIRED PROBABILITY MAY BE COMPUTED FROM SUBROUTINE ONE.
C
C     SUBROUTINES CALLED: ONE
C***********************************************************************
      SUBROUTINE CHI (K,Z,P)
      DATA PI /.564189583/
      P=-1.
      IF ((Z.LT.0.).OR.(K.LT.1)) RETURN
      P=1.
C
C   IF Z IS VERY LARGE FOR THE VALUE OF K RETURN WITH P 1.
C
      BIG=(SQRT(FLOAT(K)*2.)*12.)+FLOAT(K)
      IF (Z.GT.BIG) RETURN
      IF (K.GT.100) GO TO 40
      RK=FLOAT(K)*.5
      Z2=Z*.5
C
C   CHECK FOR EVEN DEGREES OF FREEDOM.
C
      IRK=IFIX(RK)
      IF (ABS(FLOAT(IRK)-RK).GT..1) GO TO 20
C
C   EVEN DEGREES OF FREEDOM, COMPUTATION OF P IS EASY
C
      Q=EXP(-Z2)
      P=1.-Q
      IF (IRK.EQ.1) RETURN
      MAX=IRK-1
      DO 10 I=1,MAX
      Q=Q*Z2/FLOAT(I)
   10 P=P-Q
      RETURN
C
C   DEGREES OF FREEDOM ARE ODD. TWO SETS OF CALCULATIONS ARE NEEDED TO
C   GET P.  FIRST FIND CHI(1) AT Z
C
   20 CALL ONE (Z,P)
      IF (K.EQ.1) RETURN
      A=SQRT(Z2)*2.*EXP(-Z2)*PI
      P=P-A
      IF (K.EQ.3) RETURN
      I=IRK-2
   30 A=A*Z2/(RK-FLOAT(I+1))
      P=P-A
      I=I-1
      IF (I.LT.0) RETURN
      GO TO 30
C
C   DEGREES OF FREEDOM GREATER THAN 100. USE THE WILSON-HILFERTY TRANSF.
C   Z IS MAPPED INTO Y WHICH IS APPROXIMATELY A STANDARD NORMAL VARIABLE
C
   40 Y=(((Z/FLOAT(K))**.3333333)-1.+(1./(4.5*FLOAT(K))))*SQRT(4.5*FLOAT
     1 (K))
      SIGN=-1.
      IF (Y.GE.0.) SIGN=+1.
      ZZ=Y*Y
      CALL ONE (ZZ,P)
      P=.5+SIGN*P*.5
      RETURN
      END