function VectorAboutAxis(th, ax, ay, az, px, py, pz)
    local m = ax*ax+ay*ay+az*az;
    if abs(1-m)>1e-6 then
        local r = 1/sqrt(m);
        ax, ay, az = ax*r, ay*r, az*r;
    end

    -- vector projection
    local d = dot(px,py,pz,ax,ay,az);
    local ux, uy, uz = d*ax, d*ay, d*az;
    --   (the vector rejection is simply <px-ux, py-uy, pz-uz>, so we'll use that later)
    -- cross product
    local vx, vy, vz = cross(ax,ay,az,px,py,pz);
    -- sine & cosine of our angle offset
    local c, s = cos(th), sin(th);
    -- return new vector    u + (p-u)*c + v*s
    return      ux  +(px-ux)*c  +vx*s,
            uy  +(py-uy)*c  +vy*s,
            uz  +(pz-uz)*c  +vz*s;
end