Let's put together a simple mathematical model of an accelerometer - from this, we can work out some calibration options.
Ignoring non-linearity and other nasty effects, the output measurement of an accelerometer is given by:
$$\hat{\mathbf{f}} = \mathbf{M} \mathbf{f} + \mathbf{b}_a + \mathbf{n}_a$$
where \$\hat{\mathbf{f}}\$ is the actual measurement, \$\mathbf{b}_a\$ is the accelerometer bias, \$\mathbf{n}_a\$ is a random noise vector, \$\mathbf{f}\$ is the true specific force (i.e. acceleration) and \$\mathbf{M}\$ is the Scale Factor/Misalignment Matrix.
The individual elements of the SFA matrix are:
$$ \mathbf{M} = \begin{bmatrix} S_x && \gamma_{xy} && \gamma_{xz} \\ \gamma_{yx} && S_{yy} && \gamma_{yz} \\ S_x && \gamma_{zy} && S_{zz} \\ \end{bmatrix} $$
So, each scale factor is represented by an \$S\$ and each cross-axis sensitivity is represented by a \$\gamma\$.
Ideally, if the scale factor is 1 and there is no cross-axis sensitivity, then then the resulting matrix is \$\mathbf{M} = \mathbf{I}\$.
Representing it like this allows us to develop a compensation model. If we happen to know \$\mathbf{M}\$ and \$\mathbf{b}_a\$ and assume \$\mathbf{n}_a\$ to be small (i.e. close to zero), we can make a good estimate of the "true" acceleration from the measurements:
$$
\mathbf{f} = \mathbf{M}^{-1}\left(\hat{\mathbf{f}} - \mathbf{b}_a\right)
$$
The trick is, of course, working out \$\mathbf{M}\$ and \$\mathbf{b}_a\$.
I'll describe a procedure called the six position test, which is an easy and cheap way to calibrate an accelerometer. Step 1 is to mount the accelerometer in a rectangular box with perfectly \$90^\circ\$ sides (or as close as you can get). Place this on a perfectly level surface (or, again, as close as you can get) - you'd be surprised how good you can do this.
At this point, we know what the value should be: gravity on the z-accelerometer:
$$
\mathbf{f}_1 = \begin{bmatrix} 0 \\ 0 \\ g\end{bmatrix}
$$
So, this becomes:
$$
\hat{\mathbf{f}}_1 = \mathbf{M} \mathbf{f}_1 + \mathbf{b}_a + \mathbf{n}_a
$$
Noting that \$\hat{\mathbf{f}}_1\$ will close, but not the same as \$\mathbf{f}_1\$
If we put the box on it's head, the force acting is \$-g\$:
$$
\mathbf{f}_2 = \begin{bmatrix} 0 \\ 0 \\ -g\end{bmatrix}
$$
And when placed on one side:
$$
\mathbf{f}_3 = \begin{bmatrix} 0 \\ -g \\ 0\end{bmatrix}
$$
And so on for the remaining three sides.
Now, let's write out one of the equations longhand:
$$
\hat{\mathbf{f}}_1 = \mathbf{M} \mathbf{f}_1 + \mathbf{b}_a + \mathbf{n}_a
= \begin{bmatrix} S_{xx} f_x + \gamma_{xy} f_y + \gamma_{xz} f_z + b_x \\
\gamma_{yx} f_x + S_{yy} f_y + \gamma_{yz} f_z + b_y \\
\gamma_{xz} f_x + \gamma_{yz} f_y + S_{zz} f_z + b_z \end{bmatrix}
$$
And even longer hand (for the first one):
$$
\hat{\mathbf{f}}_1 = \begin{bmatrix}
f_x S_{xx} + f_y \gamma_{xy} + f_z \gamma_{xz} +
0 \gamma_{yx} + 0 S_{yy} + 0 \gamma_{yz} +
0 \gamma_{xz} + 0 \gamma_{yz} + 0 S_{zz} +
1 b_x + 0 b_y + 0 b_z \\
0 S_{xx} + 0 \gamma_{xy} + 0 \gamma_{xz} +
f_x \gamma_{yx} + f_y S_{yy} + f_z \gamma_{yz} +
0 \gamma_{xz} + 0 \gamma_{yz} + 0 S_{zz} +
0 b_x + 1 b_y + 0 b_z \\
0 S_{xx} + 0 \gamma_{xy} + 0 \gamma_{xz} +
0 \gamma_{yx} + 0 S_{yy} + 0 \gamma_{yz} +
f_x \gamma_{xz} + f_y \gamma_{yz} + f_z S_{zz} +
0 b_x + 0 b_y + 1 b_z
\end{bmatrix}
$$
So we can create a stacked vector of the unknowns
$$
\mathbf{z} = \mathbf{A} \mathbf{\beta}
$$
Where
$$
\mathbf{z} =
\begin{bmatrix}
\hat{\mathbf{f}}_1 \\
\hat{\mathbf{f}}_2 \\
\vdots \\
\hat{\mathbf{f}}_6
\end{bmatrix}
$$
And
$$
\mathbf{\beta} = \begin{bmatrix}
S_{xx} \\
\gamma_{xy} \\
\gamma_{xz} \\
\gamma_{yx} \\
S_{yy} \\
\gamma_{yz} \\
\gamma_{xz} \\
\gamma_{yz} \\
S_{zz} \\
b_x \\
b_y \\
b_z \\
\end{bmatrix}
$$
The design matrix is (for one set of measurements):
$$
\hat{A}_1 = \begin{bmatrix}
f_x && f_y && f_z &&
0 && 0 && 0 &&
0 && 0 && 0 &&
1 && 0 && 0 \\
0 && 0 && 0 &&
f_x && f_y && f_z &&
0 && 0 && 0 &&
0 && 1 && 0 \\
0 && 0 && 0 &&
0 && 0 && 0 &&
f_x && f_y && f_z &&
0 && 0 && 1
\end{bmatrix}
$$
Now, once this is setup, one may solve for \$\mathbf{\beta}\$ (and hence sensitivity and bias) via least squares.
A similar procedure may be performed with a robotic arm if you can precisely control the angles - it simply replies knowing the precise gravity at that angle which, if you know the angle, is easy to calculate.
Somehow, I didn't notice this question at the time it was asked, or I would have answered sooner.
Yes, what you want to do is definitely possible, and there is software out there to help you do it, but unfortunately, some of it is quite expensive.
Basically, you need to use something like a Kalman filter to combine all of your data (raw acceleration, angular rate and magnetic field readings) into a complete model of the sensor that includes both its relative position and its absolute attitude (orientation). (If you also add something like a GPS receiver, you can also get absolute position.)
The filter operates in such a way that it uses the magnetometer to correct for the long-term errors (e.g., offset and scale changes) in the accelerometers and gyros, and then you can use this estimate of the attitude of the sensor to tilt-compensate the magnetic field readings in order to get a magnetic heading. As a bonus, you also get the true "course over ground" of the sensor from the sequence of position estimates, which could differ from the heading if the swimmer is swimming in a cross-current of some sort.
As far as the software for doing this kind of post-processing, one example I'm aware of is called Inertial Explorer, but it's very expensive and is probably way overkill for your application. If you search around, you should be able to find open-source projects that meet your needs.
Best Answer
Assumptions
The transformation
The measured vectors are obtained in IMU (accelerometer) frame. To take a vector resolved in IMU frame, to inertial(?) frame, the transformation as given in the above reference is
$$ \begin{bmatrix} v \end{bmatrix}^I_{3\times1} = \begin{bmatrix} C\psi & -S\psi & 0\\ S\psi & C\psi & 0\\ 0 & 0 & 1\\ \end{bmatrix} \color{red}{ \begin{bmatrix} C\theta & 0 & S\theta\\ 0 & 1 & 0\\ -S\theta & 0 & C\theta\\ \end{bmatrix} \begin{bmatrix} 1 & 0 & 0\\ 0 & C\phi & -S\phi\\ 0 & S\phi & C\phi\\ \end{bmatrix} } \begin{bmatrix} v \end{bmatrix}^{IMU}_{3\times1} $$
The red matrices indicate what I assume is the transformation equation set shown in the question.
Assume the magnetometer data was available in the same frame of reference as the accelerometer. Let that reading be \$[x_M', y_M'z_M']^T\$.
$$ \begin{bmatrix} x_{M2}\\ y_{M2}\\ z_{M2} \end{bmatrix} = \begin{bmatrix} C\theta & S\phi S\theta & C\phi S\theta\\ 0 & C\phi & -S\phi\\ \dots & \dots & \dots \end{bmatrix} \begin{bmatrix} x_{M}'\\ y_{M}'\\ z_{M}' \end{bmatrix} $$
Since the Y and Z axes are inverted for the magnetometer, the above equation changes to $$ \begin{bmatrix} x_{M2}\\ y_{M2}\\ z_{M2} \end{bmatrix} = \begin{bmatrix} C\theta & S\phi S\theta & C\phi S\theta\\ 0 & C\phi & -S\phi\\ \dots & \dots & \dots \end{bmatrix} \begin{bmatrix} x_{M}\\ \color{red}{-}y_{M}\\ \color{red}{-}z_{M} \end{bmatrix} $$
The above is significantly different from your equations.
Sanity check
You have mentioned in the comments that "pitch on y axis". This means that a rotation about pitch should leave the Y component of a vector unchanged (if it was the last operation performed). Equation for
yMin the question doesn't seem to satisfy that logic. Of course, This check is only correct assuming a certain sequence of rotations.Note
I see that your equations seem almost correct if the sequence of rotations to go from inertial frame to body frame is (Yaw, Roll, Pitch). The negation on Y and Z components before applying the equations are still required. So, see if your output becomes correct if you insert
just before the transformation.