不使用MATLAB自带工具箱，写一个普通克里金插值法的MATLAB源代码

以下是一个基于普通克里金插值方法的MATLAB源代码：

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

% 计算距离矩阵
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

% 拟合半方差函数
p = polyfit(distance(~isnan(variance)),variance(~isnan(variance)),2);

% 预测插值点的值
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);
    k_new = polyval(p,d_new);
    k_new = k_new(~isnan(k_new));
    z_new(i) = sum(k_new.*z(~isnan(k_new))) / sum(k_new);
end

% 输出结果
z = z_new;

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