将上面的代码优化，要求可自定义选择半方差函数模型（如球状模型，高斯模型，指数模型）

以下是一个可自定义选择半方差函数模型的优化版本的MATLAB源代码：

function [z] = ordinary_kriging(x,y,z,xx,yy,model)

% 计算距离矩阵
n = length(x);
D = zeros(n,n);
for i = 1:n
    for j = i+1:n
        D(i,j) = sqrt((x(i)-x(j))^2 + (y(i)-y(j))^2);
        D(j,i) = D(i,j);
    end
end

% 计算半方差函数
h = D(:);
z = z(:);
dh = max(h)/30;
distance = 0:dh:max(h);
n_distance = length(distance);
variance = zeros(n_distance,1);
for i=1:n_distance
    ix = find(h >= distance(i)-dh & h < distance(i)+dh);
    if isempty(ix)
        variance(i) = NaN;
    else
        variance(i) = var(z(ix));
    end
end

% 拟合半方差函数
switch model
    case 'spherical'
        [p,s] = polyfit(distance(~isnan(variance)),variance(~isnan(variance)),2);
        C0 = s.normr^2/length(variance);
        a = p(1);
        r = -p(2)/(2*a);
    case 'exponential'
        [p,s] = polyfit(distance(~isnan(variance)),log(variance(~isnan(variance))),1);
        C0 = s.normr^2/length(variance);
        r = -1/p(1);
    case 'gaussian'
        [p,s] = polyfit(distance(~isnan(variance)),log(variance(~isnan(variance))),2);
        C0 = s.normr^2/length(variance);
        r = sqrt(-1/p(2));
end

% 预测插值点的值
n_new = length(xx);
z_new = zeros(n_new,1);
for i = 1:n_new
    d_new = sqrt((x-xx(i)).^2 + (y-yy(i)).^2);
    switch model
        case 'spherical'
            k_new = C0 + a*(1.5*d_new/r - 0.5*(d_new/r).^3).*(d_new<=r) + a*(d_new/r).*(d_new>r);
        case 'exponential'
            k_new = C0 + exp(-d_new/r);
        case 'gaussian'
            k_new = C0 + exp(-(d_new/r).^2);
    end
    k_new = k_new(~isnan(k_new));
    z_new(i) = sum(k_new.*z(~isnan(k_new))) / sum(k_new);
end

% 输出结果
z = z_new;

在此代码中，我们首先计算了所有插值点之间的距离矩阵，并计算了每个距离的半方差函数。然后，我们根据用户选择的模型来拟合半方差函数，并使用插值点的坐标和半方差函数来预测每个插值点的值。最后，我们输出预测结果。请注意，此代码可能需要根据您的数据进行修改。