% 示例用法  
p = [1, 0, 0];  % 原始点  
axis = [0, 0, 1];  % 旋转轴（Z轴），已经是单位向量  
theta = 45;  % 旋转45度（角度） ，代码自动转成弧度计算
p_rotated = rotatePointAroundAxis(p,axis,theta);

代码如下：

function p_rotated = rotatePointAroundAxis(p, axis, theta)  
    % 输入:  
    % p - 原始点坐标，形式为 [x, y, z]  
    % axis - 旋转轴的单位向量，形式为 [ax, ay, az]  
    % theta - 旋转角度（角度）  
    theta=theta/360*2*pi；   %角度换算成弧度
    % 确保axis是单位向量  
    axis_norm = norm(axis);  
    if abs(axis_norm - 1) > eps  
        axis = axis / axis_norm;  
    end  
  
    % 旋转轴向量  
    u = axis;  

    % 旋转轴的一个垂直向量（这里选择(0,0,1)与u的叉积，如果u接近(0,0,1)则选择(1,0,0)）  
    if abs(u(3)) > 0.99  % 避免除以接近零的数  
        v = [1 0 0];  
    else  
        v = cross([0 0 1], u);  
        v = v / norm(v);  % 确保v是单位向量  
    end  
  
    % w = u x v，第三个垂直向量  
    w = cross(u, v);  
  
    % 构建旋转矩阵R  
    costheta = cos(theta);  
    sintheta = sin(theta);  
    R = [u(1)*costheta + v(1)*sintheta*w(1) - v(2)*sintheta*w(2) + v(3)*sintheta*w(3), ...  
         u(2)*costheta + v(1)*sintheta*w(2) + v(2)*sintheta*w(1) - v(3)*sintheta*w(3), ...  
         u(3)*costheta - v(1)*sintheta*w(3) + v(2)*sintheta*w(3) + v(3)*sintheta*w(1);  
           
         u(1)*(-v(2)*w(3) + v(3)*w(2))*sintheta + u(2)*(v(1)*w(3) - v(3)*w(1))*sintheta + u(3)*(-v(1)*w(2) + v(2)*w(1))*sintheta + costheta*v(1), ...  
         u(1)*(v(3)*w(1) - v(1)*w(3))*sintheta + u(2)*(-v(2)*w(1) + v(1)*w(2))*sintheta + u(3)*(v(2)*w(3) - v(3)*w(2))*sintheta + costheta*v(2), ...  
         u(1)*(-v(1)*w(2) + v(2)*w(1))*sintheta + u(2)*(v(3)*w(2) - v(2)*w(3))*sintheta + u(3)*(v(1)*w(3) - v(3)*w(1))*sintheta + costheta*v(3);  
           
         w(1), w(2), w(3)];  
  
    % 应用旋转矩阵  
    p_rotated = R * p';  % 注意MATLAB中矩阵乘法要求列向量  
    p_rotated = p_rotated';  % 如果需要，将结果转置回行向量  
end