/**
 *
 * tncCdf.java
 *
 * Copyright(c) Visual Numerics, Inc. 2000
 *
 */

package com.imsl.demo.Samples;
import com.imsl.stat.Cdf;
import com.imsl.math.JMath;
import com.imsl.math.Sfun;

// (defun tnc-cdf (T1 DF DELTA &key print (ifault 0) (tol .00001) (count-limit 40))
// "Args:  (T1, DF, DELTA, print)
//      Translated from XLISP-STAT, FROM FORTRAN IN
//      << ALGORITHM AS 243  APPL. STATIST. (1989), VOL.38, NO. 1
//      Cumulative probability at T of the non-central t-distribution
//      with DF degrees of freedom (may be fractional) and non-centrality
//      parameter DELTA. >>  Vectorized."

/**
 *
 * @author     Gregory D. Rodd
 * @created    October 20, 2002
 */
public final class tncCdf  {
    
    final static double R2PI = 1.0 / (Sfun.gamma(1.5) * JMath.sqrt(2.0));
    final static double ITRMAX = 100.1;
    final static double ERRMAX  = 1.0e-06 ;
    final static double ALNRPI =  JMath.log(JMath.sqrt(JMath.PI));
    
    // @return double;
    // @param t  double;
    // @param df double;
    // @param delt double;
    public static double cdf(double t, double df, double delt) {
        boolean NEGDEL = (t < 0.0);
        double TT = (NEGDEL ? -1.0 * t : t);
        double DEL = (NEGDEL ? -1.0 * delt : delt);
        double ERRBD = 0.0;
        int    IFAULT = 2;
        // C
        // C     Initialize twin series (Guenther, J. Statist. Computn. Simuln.
        // C     vol.6, 199, 1978).
        // C
        
        double EN = 1.0;
        double x = t * t / (t * t + df);
        double lmda = DEL * DEL;
        double P = 0.5 * JMath.exp(-0.5 * lmda);
        double Q = R2PI * P * DEL;
        double S = 0.5 - P;
        double A = 0.5;
        double B  = 0.5 * df;
        double RXB = JMath.pow(1.0 - x, B);
        double ALBETA = ALNRPI + Sfun.logGamma(B) - Sfun.logGamma(0.5 + B);
        double XODD  = Cdf.beta(x, 0.5, B);
        double GODD  = 2.0 * RXB * JMath.exp(A * JMath.log(x) - ALBETA);
        double XEVEN = 1.0 - RXB;
        double GEVEN = B * x * RXB;
        double TNC   = P * XODD + Q * XEVEN;
        
        do {
            A = A + 1.0;
            XODD = XODD - GODD;
            XEVEN = XEVEN - GEVEN;
            GODD  = GODD * x * (A + B - 1.0) / A;
            GEVEN = GEVEN * x * (A + B - 0.5) / (A + 0.5);
            P = P * lmda / (2.0 * EN);
            Q = Q * lmda / ((2.0 * EN) + 1);
            S = S - P;
            EN += 1.0;
            TNC = TNC + P * XODD + Q * XEVEN;
            ERRBD =  2.0  * S * (XODD - GODD);
        }  while( (ERRBD > ERRMAX) && (EN < ITRMAX) );
        IFAULT = 1;
        if (EN > ITRMAX) {
            //  TNC = Nan;
            IFAULT = 1;
            return Double.NaN;
        } else {
            IFAULT = 0;
            //  Add on the UPPER tail of the normal
            TNC = TNC + (1.0 - Cdf.normal(DEL));
            return (NEGDEL) ? 1.0 - TNC : TNC;
        }
    }
}