matlab 画ekman 螺线，三维的，还有随高度有箭头

Ekman 螺线是描述在旋转流中由科氏力引起的粒子运动轨迹。在 MATLAB 中，可以使用 `quiver3` 函数绘制带箭头的三维向量，来表示粒子在不同高度上的速度和方向。为了绘制 Ekman 螺线，可以按照以下步骤进行：

1. 定义旋转流的速度和方向。这可以通过定义一个函数来实现，例如：

function [u,v,w] = ekman_velocity(x,y,z)
    omega = 1; % 旋转速度
    k = 0.1; % 垂直方向的衰减率
    u = -omega*y + k*z;
    v = omega*x;
    w = -k*x;
end

这个函数返回三个向量 `u`、`v` 和 `w`，分别表示在点 `(x,y,z)` 处的速度在 x、y、z 三个方向上的分量。

2. 定义 Ekman 螺线的起始点和终止点。可以选择在 xz 平面上定义一个圆形轨迹作为起始点，例如：

t = linspace(0,2*pi,100); % 参数化圆形轨迹
r = 1; % 圆形轨迹的半径
x0 = r*cos(t);
y0 = zeros(size(t));
z0 = r*sin(t);

起始点和终止点之间的距离可以根据需要进行调整。

3. 在起始点处计算速度向量，并计算在该点处的科氏力方向。可以使用以下代码：

[u0,v0,w0] = ekman_velocity(x0(1),y0(1),z0(1));
u0 = u0/norm([u0,v0,w0]); % 速度向量单位化
v0 = cross([0,1,0], [u0,v0,w0]); % 科氏力方向
v0 = v0/norm(v0); % 科氏力方向单位化

其中，`cross` 函数计算两个向量的叉积，即科氏力方向。

4. 在起始点处绘制速度向量和科氏力箭头。可以使用以下代码：

quiver3(x0(1),y0(1),z0(1),u0,v0(2),w0,'r'); % 速度向量
quiver3(x0(1),y0(1),z0(1),0,v0(2),0,'b'); % 科氏力箭头

其中，`quiver3` 函数的参数分别表示起始点的坐标和箭头的长度和方向。

5. 在每个点处计算速度向量和科氏力箭头，并绘制。可以使用以下代码：

for i = 2:length(x0)
    [u0,v0,w0] = ekman_velocity(x0(i),y0(i),z0(i));
    u0 = u0/norm([u0,v0,w0]);
    v0 = cross([0,1,0], [u0,v0,w0]);
    v0 = v0/norm(v0);
    quiver3(x0(i),y0(i),z0(i),u0,v0(2),w0,'r');
    quiver3(x0(i),y0(i),z0(i),0,v0(2),0,'b');
end

6. 将图形设置为 3D 模式，并设置坐标轴范围和标签。可以使用以下代码：

view(3);
axis equal;
xlabel('x');
ylabel('y');
zlabel('z');

完整的 Ekman 螺线绘制代码如下：

function ekman_spiral()
    t = linspace(0,2*pi,100);
    r = 1;
    x0 = r*cos(t);
    y0 = zeros(size(t));
    z0 = r*sin(t);

    figure;
    hold on;
    for i = 1:length(x0)
        [u0,v0,w0] = ekman_velocity(x0(i),y0(i),z0(i));
        u0 = u0/norm([u0,v0,w0]);
        v0 = cross([0,1,0], [u0,v0,w0]);
        v0 = v0/norm(v0);
        quiver3(x0(i),y0(i),z0(i),u0,v0(2),w0,'r');
        quiver3(x0(i),y0(i),z0(i),0,v0(2),0,'b');
    end
    [u0,v0,w0] = ekman_velocity(x0(1),y0(1),z0(1));
    u0 = u0/norm([u0,v0,w0]);
    v0 = cross([0,1,0], [u0,v0,w0]);
    v0 = v0/norm(v0);
    quiver3(x0(1),y0(1),z0(1),u0,v0(2),w0,'r');
    quiver3(x0(1),y0(1),z0(1),0,v0(2),0,'b');
    view(3);
    axis equal;
    xlabel('x');
    ylabel('y');
    zlabel('z');
end

function [u,v,w] = ekman_velocity(x,y,z)
    omega = 1;
    k = 0.1;
    u = -omega*y + k*z;
    v = omega*x;
    w = -k*x;
end