FUNCTION chi2nc(x, f, theta) RESULT(fn_val)
 
! Code converted using TO_F90 by Alan Miller
! Date: 2000-06-16  Time: 00:34:40

! ALGORITHM AS 275 APPL.STATIST. (1992), VOL.41, NO.2

! Computes the noncentral chi-square distribution function with positive
! real degrees of freedom f and nonnegative noncentrality parameter theta

IMPLICIT NONE
INTEGER, PARAMETER    :: dp = SELECTED_REAL_KIND(12, 60)

REAL (dp), INTENT(IN)  :: x
REAL (dp), INTENT(IN)  :: f
REAL (dp), INTENT(IN)  :: theta
REAL (dp)              :: fn_val

! Local variables

INTEGER    :: ifault
LOGICAL    :: flag
REAL (dp)  :: lam, n, u, v, x2, f2, t, term, bound

! EXTERNAL alogam

INTEGER, PARAMETER    :: itrmax = 50
REAL (dp), PARAMETER  :: errmax = 1.0E-12, zero = 0.0_dp, one = 1.0_dp,  &
                         two = 2.0_dp

fn_val = x
ifault = 2
IF (f <= zero .OR. theta < zero) GO TO 999
ifault = 3
IF (x < zero) GO TO 999
ifault = 0
IF (x == zero) RETURN
lam = theta / two

!       Evaluate the first term

n = one
u = EXP(-lam)
v = u
x2 = x / two
f2 = f / two
t = x2 ** f2 * EXP(-x2) / EXP(lngamma(f2 + one))
term = v * t
fn_val = term

!       Check if (f+2n) is greater than x

flag = .false.
10 IF (f + two * n - x <= zero) GO TO 30

!       Find the error bound and check for convergence

flag = .true.
20 bound = t * x / (f + two * n - x)
IF (bound > errmax .AND. INT(n) <= itrmax) GO TO 30
IF (bound > errmax) ifault = 1
GO TO 999

!       Evaluate the next term of the expansion and then the partial sum

30 u = u * lam / n
v = v + u
t = t * x / (f + two * n)
term = v * t
fn_val = fn_val + term
n = n + one
IF (flag) GO TO 20
GO TO 10

999 WRITE(*, *) 'Error in CHI2NC, IFAULT = ', ifault
RETURN


CONTAINS


FUNCTION lngamma(z) RESULT(lanczos)

!  Uses Lanczos-type approximation to ln(gamma) for z > 0.
!  Reference:
!       Lanczos, C. 'A precision approximation of the gamma
!               function', J. SIAM Numer. Anal., B, 1, 86-96, 1964.
!  Accuracy: About 14 significant digits except for small regions
!            in the vicinity of 1 and 2.

!  Programmer: Alan Miller
!              1 Creswick Street, Brighton, Vic. 3187, Australia
!  Latest revision - 14 October 1996

REAL (dp), INTENT(IN)  :: z
REAL (dp)              :: lanczos

! Local variables

REAL (dp), PARAMETER   :: a(9) = (/   &
                          0.9999999999995183_dp, 676.5203681218835_dp, &
                          -1259.139216722289_dp, 771.3234287757674_dp, &
                          -176.6150291498386_dp, 12.50734324009056_dp, &
                          -0.1385710331296526_dp, 0.9934937113930748D-05, &
                          0.1659470187408462D-06 /),  &
                          lnsqrt2pi = 0.9189385332046727_dp, half = 0.5_dp, &
                          sixpt5 = 6.5_dp, seven = 7.0_dp
INTEGER    :: j
REAL (dp)  :: tmp

IF (z <= zero) THEN
  WRITE(*, *) 'Error: zero or -ve argument for lngamma'
  RETURN
END IF

lanczos = zero
tmp = z + seven
DO j = 9, 2, -1
  lanczos = lanczos + a(j)/tmp
  tmp = tmp - one
END DO
lanczos = lanczos + a(1)
lanczos = LOG(lanczos) + lnsqrt2pi - (z + sixpt5) + (z - half)*LOG(z + sixpt5)
RETURN

END FUNCTION lngamma


END FUNCTION chi2nc