I implement a torpedo object that chases a rotating planet. Specifically, it will turn toward the planet each update. Initially my implement was:

void move() {
        vector3<float> to_target = target - get_position();
        to_target.normalize();
        position += (to_target * speed);
}

which works perfectly for torpedo that is a solid sphere. Now my torpedo is actually a model, which has a forward vector, so using this method looks odd because it doesn't actually turn toward but jump toward. So I revised it a bit to get,

double get_rotation_angle(vector3<float> u, vector3<float> v) const {
    u.normalize();
    v.normalize();
    double cosine_theta = u.dot(v);
    // domain of arccosine is [-1, 1]
    if (cosine_theta > 1) {
        cosine_theta = 1;
    }
    if (cosine_theta < -1) {
        cosine_theta = -1;
    }
    return math3d::to_degree(acos(cosine_theta));
}

vector3<float> get_rotation_axis(vector3<float> u, vector3<float> v) const {
    u.normalize();
    v.normalize();
    // fix linear case
    if (u == v || u == -v) {
        v[0] += 0.05;
        v[1] += 0.0;
        v[2] += 0.05;
        v.normalize();
    }
    vector3<float> axis = u.cross(v);
    return axis.normal();
}

void turn_to_face() {
    vector3<float> to_target = (target - position);
    vector3<float> axis = get_rotation_axis(get_forward(), to_target);
    double angle = get_rotation_angle(get_forward(), to_target);
    double distance = math3d::distance(position, target); 
    gl_matrix_mode(GL_MODELVIEW);
        gl_push_matrix(); {
        gl_load_identity();
        gl_translate_f(position.get_x(), position.get_y(), position.get_z());
        gl_rotate_f(angle, axis.get_x(), axis.get_y(), axis.get_z());
        gl_get_float_v(GL_MODELVIEW_MATRIX, OM);
    } gl_pop_matrix();
    move();
}

void move() {
    vector3<float> to_target = target - get_position();
    to_target.normalize();
    position += (get_forward() * speed);
}

The logic is simple, I find the rotation axis by cross product, the angle to rotate by dot product, then turn toward the target position each update. Unfortunately, it looks extremely odds since the rotation happens too fast that it always turns back and forth. The forward vector for torpedo is from the ModelView matrix, the third column A:

        MODELVIEW MATRIX
--------------------------------------------------
R           U           A           T
--------------------------------------------------
1           0           0           0
0           1           0           0
0           0           1           0
0           0           0           1
-------------------------------------------------- 

Any suggestion or idea would be greatly appreciated.