已知平面法向量（含负nz），如何计算倾向与倾角并修正Matlab代码？

你的代码主要存在三个问题：未归一化法向量、倾向的坐标系转换错误、未处理nz负号的情况，以下是修正后的代码及关键说明：

修正后的代码 #

load('artificial_slope.txt'); % 加载数据
nx = data(:,4);
ny = data(:,5);
nz = data(:,6);

% 1. 归一化法向量：必须转为单位向量，避免acos参数超出[-1,1]范围
norm_vec = sqrt(nx.^2 + ny.^2 + nz.^2);
nx_unit = nx ./ norm_vec;
ny_unit = ny ./ norm_vec;
nz_unit = nz ./ norm_vec;

% 2. 修正法向量方向：确保法向量向上（nz>0）
% 当nz为负时，法向量指向平面下方，反转后不改变平面属性，但能得到正确的下坡倾向
idx = nz_unit < 0;
nx_unit(idx) = -nx_unit(idx);
ny_unit(idx) = -ny_unit(idx);
nz_unit(idx) = -nz_unit(idx);

% 3. 计算倾角（范围0-90度）
dip = acos(nz_unit) * (180/pi);

% 4. 计算倾向（范围0-360度，地质定义：正北为0，顺时针旋转）
% 先计算数学坐标系角度（从东轴逆时针转的角度），再转换为地质方位角
theta = atan2(ny_unit, nx_unit);
theta_deg = theta * (180/pi);
dip_dir = mod(90 - theta_deg, 360);

% 可选：处理水平平面（倾角接近0时，倾向无意义）
idx_flat = dip < 1e-6;
dip_dir(idx_flat) = NaN;

关键问题说明 #

1. 未归一化法向量：如果原始nx/ny/nz不是单位向量，acos(nz_abs)的参数可能超出[-1,1]，导致NaN或错误结果，归一化是必要步骤。
2. nz负号处理：nz为负表示法向量向下，此时计算的倾向会对应上坡方向，反转法向量后才能得到正确的下坡倾向。
3. 倾向坐标系转换：地质上的倾向是从正北顺时针旋转的方位角，而Matlab的atan2(ny, nx)返回的是从东轴逆时针转的角度，必须通过90 - theta_deg转换，再用mod确保结果落在0-360度范围内。

内容的提问来源于stack exchange，提问作者Hrvoje Lukacic