Actually we are also facing the same problem, instead of (Z,P,k), even if you use (A,B,C,D), as soon as sos will come in syntax with g, your magnitude response will scaled by the attenuation factor. If you are concerned only with plotting magnitude response you can plot as freqz(sos) you will get exact magnitude response without scaling.
The trigonometric functions enter the calculation because of the use the bilinear transform for transforming an analog filter to a discrete-time filter. The bilinear transform warps the frequencies of the analog filter, that's why we have to pre-warp the frequencies of the analog filter in order to obtain the desired cut-off frequencies of the discrete-time ...
In Matlab when designing standard IIR frequency selective filters (Chebyshev, Butterworth, Cauer), the resulting filter order of band pass and band stop filters is always twice the specified order (i.e., $2n$ instead of $n$), because they are obtained by transforming low pass filters of order $n$.
The reason why your maximum gain is not unity has to do ...
The help file for the cheby2 function states that
If Wn is a two-element vector,
Wn = [W1 W2], CHEBY2 returns an
order 2N bandpass filter with
passband W1 < W < W2.
Since you design a bandpass filter, your order will be $2N$ instead of $N$, and the reason for this could possibly be using two filters; like a high-low combination to create your ...