Get the discrete-time poles and zeros from continuous-time poles and zeros How do you implement the following function:

$$[Z^d, P^d, K^d] = \text{fcn} \,(Z^c, P^c, K^c),$$

where $Z^c = [z^c_m, ..., z^c_1]$, $P^c = [p^c_m, ..., p^c_1]$, and $K^c$ are zeros, poles, and gain of the continuous-time TF, and similarly $Z^d$, $P^d$, and $K^d$ are zeros, poles, and gain of the equivalent discrete-time TF with sampling time of $T_s$.

I have tried the following but it did not give me the correct answer:

$Z^d = e^{Z^c \times T_s}$ 

$P^d = e^{P^c \times T_s}$

$K^d = 1-e^{-K^c \times T_s}$    ---> (is the minus sign right?)

I am looking for some conversion formula that can handle all cases (including single/multiple poles or zeros). 

Note: I do not want to use any Matlab built-in functions like c2d(), zp2tf(), or tf().
there are at least 3, likely more, methods.


$z = e^{sT}$
impulse invariant
bilinear transform