!**************************************************************************
!*                  CHI2 AND INVERSE CHI2 FUNCTIONS                       *
!* ---------------------------------------------------------------------- *
!*  Reference:                                                            *
!*                                                                        *
!*  http://www.fourmilab.ch/rpkp/experiments/analysis/chiCalc.html        *
!*                                                                        *
!* The following JavaScript functions for calculating chi-square probabi- *
!* lities and critical values were adapted by John Walker from C imple-   *
!* mentations written by Gary Perlman of Wang Institute, Tyngsboro.       *
!* Both the original C code and the JavaScript edition are in the public  *
!* domain.                                                                *
!* ---------------------------------------------------------------------- *
!* SAMPLE RUN:                                                            *
!*                                                                        *
!*    CHI2 Law:                                                           *
!*    ========                                                            *
!*  Value of variable: 100                                                *
!*  Degree of freedom: 100                                                *
!*                                                                        *
!*  Probability = 0.481191684527957                                       *
!*                                                                        *
!*  Verify Inverse Law:                                                   *
!*    X = 99.9999999178790                                                *
!*                                                                        *
!*                                   F90 Version By J-P Moreau, Paris.    *
!*                                          (www.jpmoreau.fr)             *
!**************************************************************************
PROGRAM TEST_CHI2

   real*8  P,X,Y, chi2, invchi2
   integer df

   print *,' '
   print *,'  CHI2 Law:'
   print *,'  ======== '
   write(*,10,advance='no'); read *, df
   write(*,20,advance='no'); read *, X
   P = chi2(X, df)
   print *,' '
   print *,' Probability = ', P
   print *,' Verify Inverse Law:'
   Y = invchi2(P, df)
   print *,'   X = ', Y
   print *,' '
   stop

10 format('  Degree of freedom: ')
20 format('  Value of variable: ')

END


real*8 Function poz (z)
real*8 z
!--------------------------------------------------------------------------
! POZ  --  probability of standard normal z value
!
! Adapted from a polynomial approximation in:
! Ibbetson D, Algorithm 209
! Collected Algorithms of the CACM 1963 p. 616
!
! Note:
! This routine has six digits accuracy, so it is only useful for absolute
! z values < 6.  For z values >=  6.0, poz() returns 0.0.                 
!--------------------------------------------------------------------------
    real*8 Y, X, w, zmax
    zmax = 6.d0
    if (z == 0.0) then
      X = 0.d0
    else
      Y = 0.5 * dabs(z)
      if (Y >= zmax * 0.5d0) then
        X = 1.d0
      else if (Y < 1.d0) then
        w = Y * Y
        X = ((((((((0.000124818987*w -0.001075204047)*w + 0.005198775019)*w -0.019198292004)*w + 0.059054035642)*w &
             -0.151968751364)*w + 0.319152932694) * w - 0.5319230073)*w + 0.7978845605929999)*Y*2.d0
      else
        Y = Y - 2.d0
        X = (((((((((((((-0.000045255659*Y + 0.00015252929)*Y -0.000019538132)*Y - 6.769049860000001E-04)*Y  &
             +0.001390604284)*Y-0.00079462082)*Y -0.002034254874)*Y + 0.006549791214)*Y - 0.010557625006)*Y  &
             +0.011630447319)*Y-9.279453341000001E-03)*Y + 0.005353579108)*Y - 0.002141268741)*Y             &
             +0.000535310849)*Y + 0.999936657524;
      end if
    end if
    if (z > 0.d0) then
      poz = (X + 1.d0) * 0.5d0
    else
      poz = (1.d0 - X) * 0.5d0
    end if
    return    
    end


real*8 Function chi2(X, df)
real*8 X
integer df
!--------------------------------------------------------- 
! Adapted From:
! Hill, I. D. and Pike, M. C. Algorithm 299
! Collected Algorithms for the CACM 1967 p. 243
! Updated for rounding errors based on remark in
! ACM TOMS June 1985, page 185   
!---------------------------------------------------------
    integer even   ! parity of df
    real*8  a,bigx,c,e,lnpi, ipi,s,Y,z, PI, poz
	PI = 4.d0*datan(1.d0)
    bigx = 200.d0
    lnpi = dlog(dsqrt(PI))
    ipi = 1.d0 / dlog(PI)
    if (X <= 0.d0.or.df < 1)  then
	  chi2 = 1.d0
	  return
    end if
    a = 0.5d0 * X
    if (2*(df/2) == df) then
      even = 1
    else
      even = 0
    end if
    if (df > 1)  Y = dexp(-a)
    if (even==1) then
      s = Y
    else
      s = 2.d0 * poz(-dsqrt((X)))
    end if
    if (df > 2) then
        X = 0.5d0 * (df - 1.d0)
        if (even==1) then
          z = 1.d0
        else
          z = 0.5d0
        end if
        if (a > bigx) then
            if (even==1) then
              e = 0.d0
            else
              e = lnpi
            end if
            c = log(a)
            do while (z <= X)
              e = log(z) + e
              s = s + dexp(c * z - a - e)
              z = z + 1.d0
            end do
			chi2 = s
            return
        else
            if (even==1) then
              e = 1.d0
            else
              e = ipi / dsqrt(a)
            end if
            c = 0.d0
            do while (z <= X)
              e = e * (a / z)
              c = c + e
              z = z + 1.d0
            end do
			chi2 = c*Y + s
            return
        end if
    else
      chi2=s
    end if
    return
    end

real*8 Function invchi2 (P, idf)
    real*8 P
    real*8 chimax,chisave,chisqval,epsilon,minchisq,maxchisq,chi2
    epsilon = 1.d-6
    chimax  = 1.d6
    minchisq = epsilon
    maxchisq = chimax
    if (P <= 0.d0) then
      invchi2 = maxchisq
	  return
    else
      if (P >= 1.d0)  then
	    invchi2 = 0.d0
		return
      end if
    end if
    chisqval = (1.d0*idf) / dsqrt(P)
    do while (maxchisq - minchisq > epsilon)
	  chisave=chisqval !chisqval is modified by call to chi2 
      if (chi2(chisqval, idf) < P) then
		maxchisq = chisave
      else
        minchisq = chisave
      end if
      chisqval = (maxchisq + minchisq)/2.d0
    end do
	invchi2 = chisqval
    return
	end

! end of file chi2.f90