Orbiting object around local axis of another one with GLM.

I’ve got problem trying to rotate object based on the second one rotation. I read a lot of internet’s posts about translating, rotating, translating, changing order of translations does nothing in this case.
I’ve got a few object (for example 4, representing local coordinates systems) and with the rotation of the one local coord. system by local axis i want next ones to follow this rotation, for example:

  1. Rotate 2nd coo. system by the value of “qz” (done by quaternions) by local x Axis.
  2. Rotate 3rd co. system by the same angle and by the same axis (2nd system’s local x Axis), that it would look like it was orbiting 2nd coo. system.
  3. Rotate 4th co. system like 3rd, but of course, the orbit should be different.

For now, i try to rotate only 2 local coordinates systems.
My code, axisModel is the first local coordinate system:

                        x = sin(qx / 2);
			y = sin(qy / 2);
			z = sin(qz / 2);
			w = cos(qw / 2);
			quat Quaternion = quat(w, x, z, y);
			mat4 quaternionRotationMatrix = toMat4(Quaternion);
			axisModel = axisModel * quaternionRotationMatrix;


now the part to rotate second coordinate system, axisMVPs[currentPart] is second coordinate system:

			tempResultAxis = axisMVPs[currentPart].ModelMatrix;
			xAxis.x = axisModel[0][0];
			xAxis.y = axisModel[0][1];
			xAxis.z = axisModel[0][2];

			yAxis.x = axisModel[1][0];
			yAxis.y = axisModel[1][1];
			yAxis.z = axisModel[1][2];

			zAxis.x = axisModel[2][0];
			zAxis.y = axisModel[2][1];
			zAxis.z = axisModel[2][2];
                        double x = 0, y = 0, z = 0;
			x = -tempResultAxis[3][0];
			y = -tempResultAxis[3][1];
			z = -tempResultAxis[3][2];
			tempResultAxis = glm::translate(tempResultAxis, vec3(-(x + axisModel[3][0]), -(z + axisModel[3][2]), -(y + axisModel[3][1])));
			tempResultAxis = glm::rotate(tempResultAxis, qz, yAxis);
			tempResultAxis = glm::rotate(tempResultAxis, qy, zAxis);
			tempResultAxis = glm::rotate(tempResultAxis, qx, xAxis);
			tempResultAxis = glm::translate(tempResultAxis, vec3((x + axisModel[3][0]), (z + axisModel[3][2]), (y + axisModel[3][1])));

			axisMVPs[currentPart].ModelMatrix = tempResultAxis;

It gives me unpredictable movement of second model, if it is not in the same position of first local coordinate system - if they have the same xyz coordinates and their axes are parallel, they work perfectly. What is wrong with my code?:frowning: