%
disp(' ');
disp(' noncentral_t_distribution.m  ');
disp(' ver 1.0   August 10, 2010 ');
disp(' ');
disp(' by Tom Irvine ');
disp(' ');
%
v=input(' Enter degrees-of-freedom  ');
vd2=v/2;
%
mu=input(' Enter the noncentrality parameter  ');
mu2=mu^2;
%
disp(' ');
disp(' Select analysis ');
disp('  1=calculate probability for given t value ');
disp('  2=calculate t value for given probability ');
ia=input(' ');
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(ia==1)
%
    t=input(' Enter t value ');
%
%    Cumulative distribution function
%
    x=t;
    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(' probability =%8.4g  from -infinity to t ',csum);
    disp(out1);
else
%
    zprob=input(' Enter probability  for -infinity to +lambda  ');
%
    lambda_before=-20.;
    Pr_before=0.;
    lambda=20;
    Pr=1.;
    dx=lambda-lambda_before;
%
    for(i=1:30)
%
        lambda_int=(lambda+lambda_before)/2.;   % works better than secant
%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
        x=lambda_int;
        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;
%
       Pr_int=tsum+P;
%
%%        out1=sprintf('\n  Pr_before=%8.4g   Pr_int=%8.4g   Pr=%8.4g   ',Pr_before,Pr_int,Pr);
%%        disp(out1);
%%        out1=sprintf('  lambda_before=%8.4g   lambda_int%8.4g   lambda=%8.4g   ',lambda_before,lambda_int,lambda);
%%        disp(out1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
        if(Pr_int==zprob)
            break;
        end
%
        if(Pr_before < zprob && zprob < Pr_int)
           Pr=Pr_int;
           lambda=lambda_int;
        end
        if(Pr_int < zprob && zprob < Pr)
           Pr_before=Pr_int;
           lambda_before=lambda_int;
        end
        dx=lambda-lambda_before;
%    
    end
    out1=sprintf('\n lambda=%8.4g ',lambda_int);
    disp(out1);    
end