stdgaminv.m

关于高层建筑强风条件下风压系数计算

function x = stdgaminv(p,gam)

% Set parameters for iteration:
abs_tol         = 10^-10;
rel_tol         = 10^-10;
max_iter        = 500;

% 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));