%
disp(' ');
disp(' integrate_normal_distribution.m  ');
disp(' ver 1.1   August 11, 2010 ');
disp(' ');
disp(' by Tom Irvine ');
disp(' ');
%
disp(' Select analysis ');
disp('  1=calculate probability for given lambda value ');
disp('  2=calculate lambda value for given probability ');
ia=input(' ');
%
if(ia==1)
%
     lambda=input(' Enter limit value lambda  ');
%
     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;
%
     out1=sprintf('\n      P[ -lambda <= Z <=lambda ] = %15.9g ',P);
     disp(out1);
%
     out1=sprintf('\n              P[ abs(Z)>lambda ] = %15.9g ',1-P);
     disp(out1);
%
     out1=sprintf('\n   P[ - infinity <= Z <=lambda ] = %15.9g ',(1-P)/2 + P);
     disp(out1);
%
else
%
    zprob=input(' Enter probability  for + and - lambda  ');
%
    lambda_before=-7.;
    P_before=0.;
    lambda=7;
    P=1.;
    dx=lambda-lambda_before;
%
    for(i=1:100)
%
        lambda_int=(lambda+lambda_before)/2.;   % works better than secant
%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
        sum=0;
        den=(1/2)*sqrt(pi);
        q=lambda_int/sqrt(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
%
        P_int=q*exp(-q^2)*sum;
%
%%        out1=sprintf('\n  P_before=%8.4g   P_int=%8.4g   P=%8.4g   ',P_before,P_int,P);
%%        disp(out1);
%%        out1=sprintf('  lambda_before=%8.4g   lambda_int%8.4g   lambda=%8.4g   ',lambda_before,lambda_int,lambda);
%%        disp(out1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
        if(P_int==zprob)
            break;
        end
%
        if(P_before < zprob && zprob < P_int)
           P=P_int;
           lambda=lambda_int;
        end
        if(P_int < zprob && zprob < P)
           P_before=P_int;
           lambda_before=lambda_int;
        end
        dx=lambda-lambda_before;
%    
    end
    out1=sprintf('\n lambda=%8.4g ',lambda_int);
    disp(out1);
end