function ret = angularDistributionFunction(this, theta)
    %This function calculate the angle distribution of atoms coming out
    %from a single channel.
    
    KnudsenNumber = this.MeanFreePath/this.NozzleLength;

    alpha=0.5 - (3*this.Beta^2)^-1 * ...
        ((1 - (2*this.Beta^3) + ((2*this.Beta^2) - 1) * sqrt(1+this.Beta^2)) / ...
        (sqrt(1+this.Beta^2) - (this.Beta^2 * asinh((this.Beta^2)^-1))));
    
    eta0 = alpha;
    eta1 = 1 - alpha;
    
    delta = (eta0./sqrt(2*KnudsenNumber*(eta1-eta0)))./sqrt(cos(theta));
    
    F = 2/sqrt(pi)* (1-eta1)/eta0 * delta.* exp( -(delta*eta1/eta0).^2 );
    q = this.Beta^-1 * tan(theta);
    R = acos(q) - (q .* sqrt(1-q.^2));
    
    if abs(q) >= 1 
        t = linspace(0,1,10000);
        S = sum(sqrt(1-t.^2).* ( erf(delta.*(1 + (t.*(eta1-eta0)./(q.*eta0)) ))-erf(delta)))*(t(2)-t(1));
        if S == 0 || isnan(S)
            ret = eta0*cos(theta);
        else
            ret = eta0*cos(theta)+ 2/sqrt(pi)*eta0*cos(theta) * (exp(delta.^2)/delta) * S;
        end
    else
        t = linspace(0,q,10000);
        S = sum(sqrt(1-t.^2).* ( erf(delta.*(1 + (t.*(eta1-eta0)./(q.*eta0)) ))-erf(delta)))*(t(2)-t(1));
        if isnan(S)
            S=0;
        end
        ret = 2/sqrt(pi)*eta0*cos(theta)*(exp(delta.^2)/delta) * (R./2*(erf(delta*eta1/eta0)-erf(delta)+F)+S)+eta0*cos(theta);
    end
end