function x = stdgaminv( p, gam )

% set parameters for iteration:
abs_tol = 10^-2;
rel_tol = 10^-2;
max_iter = 10;

% estimate maximum and minimum values of x for interpolation:
x_max = gam+5*sqrt(gam);
iter = 0;
while gammainc(x_max,gam)<max(p)
    iter = iter+1;
    if iter>max_iter
        error(['Maximum specified probability is too high: ' num2str(max(p))]);
    end
    x_max = 1.5*x_max;
end
x_min = gam+sqrt(gam);
iter = 0;
while gammainc(x_min,gam)>min(p)
    iter = iter+1;
    if iter>max_iter
        error(['Minimum specified probability is too low: ' num2str(min(p))]);
    end
    x_min = 0.75*x_min;
end

% assemble a list of values and probabilities for interpolation:
n_check = 1000;
x_check = linspace(x_min,x_max,n_check);
p_check = gammainc(x_check,gam);

[p_check, ind_u, junk] = unique(p_check); % ensure no repititions
x_check = x_check(ind_u);

% interpolate to estimate values corresponding to specified probabilities:
x_est = interp1(p_check, x_check, p);

iter = 0;
x_step = ones(size(x_est)); % initialize vector of step sizes

while any(abs(x_step)>abs_tol) && any(abs(x_step./x_est)>rel_tol)
    iter = iter+1;
    if iter>max_iter
        warning(['Max interations reached in STDGAMINV. Max abs error: ' num2str(max(abs(x_step))) '; Max rel error: ' num2str(max(abs(x_step./x_est))) ]);
        break;
    end

    % compute probabilities associated with estimated values:
    p_est = gammainc(x_est,gam);
    
    % augment the list of values and probabilities for interpolation
    p_check = [p_check(:); p_est(:)];
    x_check = [x_check(:); x_est(:)];
    [p_check, ind_u, junk] = unique(p_check); % ensure no repetitions
    x_check = x_check(ind_u);

    x_interp = interp1(p_check, x_check, p);
    x_step = x_interp-x_est;
    x_est = x_interp;
end
x = reshape(x_est,size(p));