Is there a way to rotate multiple model matrices on one pivot point in world space?

I am using OpenGL with GLFW and I am trying to rotate multiple objects, each with their respective model matrices on the same pivot point. I have three different models for my rabbit, one body, one front legs, and one back legs, and I am trying to rotate them all in the end on the Y-Axis to face one direction. However, the legs look like they are 'moving' when I am trying to rotate them in place. If I do not rotate them along the Y-Axis, then they are fine. I want to find a way to rotate all three model matrices on one pivot point, so when they are rotated in the end then they do not get rotated on their own origin. Here is my main code:

class rabbit {
public:
    std::vector<colored_model_object*> rabbit_renderers;
    std::vector<glm::mat4> rabbit_models;

    glm::vec3 position = glm::vec3(0.0f);
    float fore_leg_rotation = 0;
    float hind_leg_rotation = 0;

    float whole_rotation = 0;

    int jump_phase = 0;
    float jump_lerp = 0.0f;
    float jump_speed = 0.003f;
    glm::vec3 jump_start, jump_end;

    float y_velocity = 0;
    
    rabbit() {
        rabbit_renderers.push_back(new colored_model_object("res/objects/rabbit/rabbit_body.obj"));
        rabbit_renderers.push_back(new colored_model_object("res/objects/rabbit/rabbit_front_legs.obj"));
        rabbit_renderers.push_back(new colored_model_object("res/objects/rabbit/rabbit_hind_legs.obj"));
    }

    void jump(glm::vec3 direction, float magnitude) {
        if (jump_phase >= 1) return;
        glm::vec3 jump_power = glm::normalize(direction) * magnitude;
        jump_start = position;
        jump_end = position + jump_power;
        jump_phase = 1;
        jump_lerp = 0;
        fore_leg_rotation = 0.0f;
        hind_leg_rotation = 0.0f;
        y_velocity = 1;

        glm::vec2 dir_2d = glm::normalize(glm::vec2(direction.x, direction.z));
        glm::vec2 up = glm::vec2(0, 1);
        float dot = dir_2d.x * up.x + dir_2d.y * up.y;
        float det = dir_2d.x * up.y - dir_2d.y * up.x;
        float body_angle = atan2(det, dot);
        whole_rotation = body_angle;
    }

    void update() {
        position.y += y_velocity;
        y_velocity -= 0.008f;
        if (position.y < 0) position.y = 0;
        if (jump_phase> 0) {
            glm::vec2 jump_start_xz = glm::vec2(jump_start.x, jump_start.z);
            glm::vec2 jump_end_xz = glm::vec2(jump_end.x, jump_end.z);
            glm::vec2 position_xz = glm::mix(jump_start_xz, jump_end_xz, jump_lerp);
            position.x = position_xz.x;
            position.z = position_xz.y;
            jump_lerp += jump_speed;
        }
        if (jump_phase == 1) {
            fore_leg_rotation += 1;
            hind_leg_rotation += 1;
            if (fore_leg_rotation > 45.0f) {
                jump_phase = 2;
            }
        }
        if (jump_phase == 2) {
            fore_leg_rotation -= 0.2f;
            hind_leg_rotation -= 0.2f;
            if (fore_leg_rotation < 0.0f) {
                jump_phase = 3;
            }
        }
        if (jump_lerp > 1) {
            jump_phase = 0;
            jump_lerp = 0;
            fore_leg_rotation = 0.0f;
            hind_leg_rotation = 0.0f;
        }
    }

    void update_models() {
        update();

        rabbit_models.clear();
        rabbit_models.push_back(system3d::get_model());
        rabbit_models.push_back(system3d::get_model());
        rabbit_models.push_back(system3d::get_model());
        rabbit_models.at(0) = (system3d::create_model(rabbit_models.at(0),glm::vec3(0.1f, 0.1f, 0.1f), glm::vec3(-100 + position.x, 310.5f + position.y, 0 + position.z), 0.0f, glm::vec3(1, 1, 1)));
        rabbit_models.at(1) = (system3d::create_model(rabbit_models.at(1), glm::vec3(0.1f, 0.1f, 0.1f), glm::vec3(-100 + position.x, 310.5f + position.y, 0 + position.z), fore_leg_rotation, glm::vec3(1, 0, 0), glm::vec3(0, -42.0f, 0)));
        rabbit_models.at(2) = (system3d::create_model(rabbit_models.at(2), glm::vec3(0.1f, 0.1f, 0.1f), glm::vec3(-100 + position.x, 310.5f + position.y, 0 + position.z), hind_leg_rotation, glm::vec3(1, 0, 0), glm::vec3(0, -42.0f, 17.0f)));
        for (int i = 0; i < rabbit_models.size(); i++) {
            rabbit_models.at(i) = glm::rotate(rabbit_models.at(i), glm::radians(whole_rotation), glm::vec3(0, 1, 0));
        }
    }
};

And here is the referenced system3d class:

class system3d {
public:
    static glm::mat4 get_model(glm::vec3 scale_vec, glm::vec3 translate_vec, float rotation, glm::vec3 rotation_vec) {
        glm::mat4 model = glm::mat4(1.0f);
        model = glm::scale(model, scale_vec);
        model = glm::rotate(model, glm::radians(rotation), rotation_vec);
        model = glm::translate(model, translate_vec);

        return model;
    }

    static glm::mat4 get_model(glm::vec3 scale_vec, glm::vec3 translate_vec, float rotation, glm::vec3 rotation_vec, glm::vec3 rotation_pivot) {
        glm::mat4 model = glm::mat4(1.0f);
        model = glm::scale(model, scale_vec);
        model = glm::translate(model, translate_vec);


        model = glm::translate(model, -rotation_pivot);
        model = glm::rotate(model, glm::radians(rotation), rotation_vec);
        model = glm::translate(model, rotation_pivot);


        return model;
    }

    static glm::mat4 create_model(glm::mat4 existing, glm::vec3 scale_vec, glm::vec3 translate_vec, float rotation, glm::vec3 rotation_vec) {
        glm::mat4 model = existing;
        model = glm::scale(model, scale_vec);
        model = glm::rotate(model, glm::radians(rotation), rotation_vec);
        model = glm::translate(model, translate_vec);

        return model;
    }

    static glm::mat4 create_model(glm::mat4 existing, glm::vec3 scale_vec, glm::vec3 translate_vec, float rotation, glm::vec3 rotation_vec, glm::vec3 rotation_pivot) {
        glm::mat4 model = existing;
        model = glm::scale(model, scale_vec);
        model = glm::translate(model, translate_vec);


        model = glm::translate(model, -rotation_pivot);
        model = glm::rotate(model, glm::radians(rotation), rotation_vec);
        model = glm::translate(model, rotation_pivot);


        return model;
    }


    static glm::mat4 get_model() {
        glm::mat4 model = glm::mat4(1.0f);
        return model;
    }
};

I have been trying to find a solution for a few hours now, and I have not yet found one that does not involve Quaternions. (I find quaternions really confusing, and I personally cannot write quaternion code that works.)

Here is a video of the result first when they Y-Axis rotation in 0, then when the y-axis rotation is not 0:

Does anyone have any ideas on what I can do?

Solution

If I understand your question, the issue is:

When the rabbit is in its initial orientation, and you rotate the legs, the legs rotate the way you want.
When you rotate the rabbit, and then rotate the legs, the legs rotate on a different axis than the rabbit (leading to strange behavior such as the legs flying off the body)

If this represents the issue, you need to perform correct order of operations by:

Rotating each of the individual sub-components as if they were in the initial orientation
Rotating the full group by the bulk rotation that you want to apply to the entire rabbit.

Ex: You want to:

Rotate the right hind leg outward (relative to rabbit) - [R1]
Rotate both hind legs back (relative to rabbit fow/back dir) - [R2]
Rotate the entire rabbit to face the camera - [R3]

Rotations applied to each are then:

Right leg gets: R_fin_right_hind_leg = [R3][R2][R1]
Left hind leg gets: R_fin_left_hind_leg = [R3][R2]
Rest of rabbit body gets: R_fin_rabbit_other = [R3]