ComputeIgCDF(X,Lambda,Mu)

function [F] = ComputeIgCDF(X,Lambda,Mu)
%
%  See O'Reilly and Rueda paper referenced in main program for explanation of this subroutine, pgs. 387-388.
%
%
q = max([Lambda/Mu,0]);
F = zeros(1,length(X));
%
%
%
for ii = 1:length(X)

    x = X(1,ii);

    if x>0
        z = x/Mu;
        c = sqrt(q/z);
        %
        a = c*(z-1);
        A = normcdf(a,0,1);
        %
        b= -c*(z+1);
        B = normcdf(b,0,1);
        %
        F(1,ii) = A + exp(2*q)*B;
    else
        F(1,ii) = 0;
    end

end