Frequency prewarping of a bilinear transform (Tustin transform)
There might be some undocumented features for the c2d function. Namely, if I follow the documented way of specifying pre-warping (using c2dOptions) one would get:


  opts = c2dOptions('Method', 'tustin', 'PrewarpFrequency', b);
  G = c2d(S, T, opts);


When I used this I do get the same result as when I would manually do pre-warping using:


  z = tf('z', h);
  sp = b / tan(b * T / 2) * (z - 1) / (z + 1);
  G = -m * sp^2 / (sp^2 + 2 * a * b * sp + b^2);


However, the notation that you use is not specified in the documentation for as far I could find. It does seem to give a better approximation of the continuous transfer function in this case. However, if I use the following transfer function

$$
P(s) = \frac{s}{s^2 + 0.1\,s + 1},
$$

with a sample time of $T = 1$ your method behaves in an unexpected way as shown in the figure below. In that figure D1 uses pre-warping using c2dOptions, D2 uses your method of pre-warping and D3 uses manual pre-warping (this completely overlaps with D1 in the figure). 



It can be noted that D2 does follow the overal shape of $P(s)$ better, but it does not match it at the specified pre-warp frequency and the phase is much lower (when zoomed out to a bigger frequency spectrum it looks like a delay). On the other hand D1/D3 exactly matches $P(s)$ at the pre-warp frequency as expected, but the overall shape does deviate more the further away you get from the pre-warp frequency.