获取 MATLAB 中所示的 Tikz/pgfplot 图形（箭头+轮廓）

Question

我有一个相当复杂的 MATLAB 脚本，它将梯度场向量绘制为表面轮廓之间的长度。

我希望能够用 LaTeX 原生绘制此图表。我可以导入箭筒和轮廓数据，但我不知道如何格式化轮廓数据，以便将其作为 MATLABcontour函数的输出提供，文档在这里，在 LaTeX 中。另一种方法是模仿代码的逻辑，让 Tikz 从头开始​​绘制所有内容，但我不知道如何通过这种方式实现它。代码如下

clf

f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1;
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);

Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));

syms x y
dz = gradient(f(x*cos(pi/6)-y*sin(pi/6),y*cos(pi/6)+x*sin(pi/6)),[x,y]);
dz = matlabFunction(dz);

% Create axes
ax = axes('Color', [1 1 1]); 
hold(ax,'on')
NumOfContours = 10;
[M,~] = contour(ax, X,Y,Z,NumOfContours);

i = 1;
Ns = zeros(1,NumOfContours);
counter = 0;
while i<=length(M)
    counter = counter + 1;
    Ns(counter) = M(2,i);
    i = i + M(2,i)+1;
end

X = nan(1,1000); % Pre-allocation
Y = X;

Ms = zeros(1,NumOfContours-1);

startIdx = 2;
overallCounter = 1;
for j=1:NumOfContours-1
    subX = M(1,startIdx:(startIdx + Ns(j) - 1));
    subY = M(2,startIdx:(startIdx + Ns(j) - 1));
    startIdx = startIdx + Ns(j) + 1;

    dX = gradient(subX);
    dY = gradient(subY);

    if j==1
        PrimaryLength = sum(hypot(dX,dY));
        NumVecPerUnit = 4;
        dl = PrimaryLength/(NumVecPerUnit+1);
    end

    Lens = cumsum(hypot(dX,dY));

    counter = 1;
    gatePass = 1;
    while gatePass      
        [~,idx] = min(abs(dl*counter - Lens));
        if idx+1 > length(Lens)
            gatePass = 0;
            Ms(j) = counter-1;
        else
            counter = counter + 1;
            X(overallCounter) = subX(idx);
            Y(overallCounter) = subY(idx);
            overallCounter = overallCounter + 1;
        end
    end
end

X(isnan(X)) = [];
Y(isnan(Y)) = [];

U = zeros(1,length(X));
V = U;

for i=1:length(X)
    d = dz(X(i),Y(i));
    U(i) = d(1);
    V(i) = d(2);
end

% Define the colormap for the quiver arrows.
% cmap can have any number of rows.
cmap = parula(255);
ax.Colormap = cmap; 

scale_size = zeros(1,length(X)); % Pre-allocate scaling array

dU = U./sqrt(U.^2 + V.^2);
dV = V./sqrt(U.^2 + V.^2);

startIdx = Ns(1)+3;
j = 0;
for i=1:length(X)    
    j_gate = i > cumsum([0,Ms]);
    j_switch = find(j_gate == max(j_gate),1,"last");
    if j_switch > j
        j = j_switch;
        x1 = M(1,startIdx:(startIdx + Ns(j+1) - 1));
        y1 = M(2,startIdx:(startIdx + Ns(j+1) - 1));
        startIdx = startIdx + Ns(j+1) + 1;
    end

    x2 = [X(i),X(i)+sqrt(2)*dU(i)];
    y2 = [Y(i),Y(i)+sqrt(2)*dV(i)];

    [x0,y0] = intersections(x1,y1,x2,y2,1);

     false_switch = (isempty(x0) || isempty(y0));
     x0(false_switch) = inf;
     y0(false_switch) = inf;

    scale_size(i) = hypot(X(i)-x0(1),Y(i)-y0(1));
end

dU = scale_size.*dU;
dV = scale_size.*dV;

remove_id = false(1,length(X));

for i=1:length(X)
    if norm([X(i)+dU(i),Y(i)+dV(i)],"inf") > 1
        remove_id(i) = true;
    end
end

X(remove_id) = [];
Y(remove_id) = [];
dU(remove_id) = [];
dV(remove_id) = [];

quiver(ax, X,Y,dU,dV, ...
    'Color', [0.1 0.1 0.1], 'LineWidth',0.3,'MaxHeadSize',2,'AutoScale','off')

% scatter(X,Y)

% Set properties for the main axes
axis equal
xlim(ax, [-1 1])
ylim(ax, [-1 1])

hold off