%
disp(' ');
disp(' noncentral_t_distribution_table.m  ');
disp(' ver 1.0   August 10, 2010 ');
disp(' ');
disp(' by Tom Irvine ');
disp(' ');
%
v=input(' Enter degrees-of-freedom  ');
%
mu=input(' Enter the noncentrality parameter  ');
mu2=mu^2;
%
disp(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    lambda=mu;
    q=lambda/sqrt(2);
%
    gamma_1_2=sqrt(pi);
    gamma_3_2=(1/2)*gamma_1_2;  % 1
    gamma_5_2=(3/2)*gamma_3_2;  % 2
%
    sum=0;
%
    den=gamma_3_2;
%
    for(n=0:1000)
%       
        b=q^(2*n)/den;
%
        if(den<1.0e+100)
        else
            break;
        end
%
        sum=sum+b;
%
        den=den*( 2*n + 3)/2;
%
    end
%
    A=q*exp(-q^2);
%
    P=A*sum;
    P=(1-P)/2;
    if(lambda>6)
        P=0;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%   I = betainc(X,Z,W)  regularized incomplete beta function
%% betainc(.5,(0:10)',3)
%%
%%
%    Cumulative distribution function
%
vd2=v/2;
dx=1;
for(i=1:21)
    x=(i-1)*dx;
    tsum=0;
    for(j=0:110)
        y=x^2/(x^2+v);
        ee=exp(-mu2/2)*(mu2/2)^j;
        pj=(1/factorial(j))*ee;
        qj=mu*ee/(sqrt(2)*gamma(j+1.5));
        j1=j+0.5;
        j2=j+1;
%
        Iy1=betainc(y,j1,vd2);
        Iy2=betainc(y,j2,vd2);    
%
        tsum=tsum+(pj*Iy1 + qj*Iy2);
%
    end
    tsum=tsum*0.5;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    csum=tsum+P;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    out1=sprintf(' x=%8.4g  csum=%8.4g  tsum=%8.4g  P=%8.4g',x,csum,tsum,P);
    disp(out1);
end