Added comments

function p_ul = power_upper_limit(max_power,n_ps,prob)

%POWER_UPPER_LIMIT Upper limit on signal power at given prob for max_power

%measured

%   Detailed explanation goes here

%POWER_UPPER_LIMIT Upper limit on signal at given prob. for measured power

%   This function calculates "the [power(s)] of the signal that, if it were 

%   present in the bin(s) with the lagest power P_max would, with 

%   probability [1-prob], have produced a power greater than P_max in that

%   bin" (Vaughan+1994, Sec. 3.3.). It works by inverting numerically the

%   noncentral chi^2 CDF which describes the probability of measuring P_max

%   based on the number of power spectra and underlying signal power.

%   Inputs:

%   - max_power (scalar or array): measured power(s) to which to associate

%     the underlying signal power upper limit(s) 

%   - n_ps (scalar): number of power spectra summed to obtain max_power

%   - prob (scalar): probability that a signal with P_ul would result in a

%     power <= P_max 

%   Output: 

%   - p_ul (scalar or array): signal power upper limit(s)

    parser=inputParser;

    parser.KeepUnmatched=true;

    addRequired(parser,'prob',@(x) isreal(x) && x<=1);

    % Simple check on 0<prob<=1 

    addRequired(parser,'prob',@(x) isreal(x) && x<=1 && x>0);

    %p_ul = fzero(@(x) ncx2cdf(max_power,2*n_ps,x)-prob,max_power);

    p_ul = arrayfun(@(max_power) fzero(@(x) ncx2cdf(max_power,2*n_ps,x)-prob,max_power),max_power);

end