%
disp(' ');
disp(' integrate_t_distribution.m  ');
disp(' ver 1.0   August 10, 2010 ');
disp(' ');
disp(' by Tom Irvine ');
disp(' ');
disp(' Student t-distribution ');
disp(' ');
%
gamma=input(' Enter area (gamma)  ');
%
v=input(' Enter degrees-of-freedom  ');
%
iflag=0;
b_before=0.00001;
Q_before = quad(@t_distribution,0,b_before,[],[],v);
b=300;
Q = quad(@t_distribution,0,b,[],[],v);
dx=b-b_before;
%
for(i=1:100)
% 
%%%    c2=(gamma-Q_before)/(Q-Q_before);
%%%    c1=1-c2;
%%%    b_int= c1*b_before + c2*b;
%
    b_int=(b+b_before)/2.;   % works better
%    
    Q_int = quad(@t_distribution,-1000,b_int,[],[],v);
%
%%    out1=sprintf(' %8.4g  %8.4g  %8.4g ',b_before,b_int,b);
%%    out2=sprintf(' %8.4g  %8.4g  %8.4g \n',Q_before,Q_int,Q);
%%    disp(out1);
%%    disp(out2);
%
    if(Q_int==gamma)
        break;
    end
%
    if(Q_before < gamma && gamma < Q_int)
       Q=Q_int;
       b=b_int;
       dx=b-b_before;
    end
    if(Q_int < gamma && gamma < Q)
       Q_before=Q_int;
       b_before=b_int;
       dx=b-b_before;
    end
%    
end
%
out1=sprintf('\n x=%10.5g  area=%10.5g  (one-sided)',b_int,Q_int);
disp(out1);    