I tried to execute the following program in Mathematica (the program from paper, A New Formula for the BER of Binary Modulations with Dual-Branch Selection over Generalized-K Composite Fading Channels),

enter code here (*Testing*) (*Declarations*) p = 0.5; q = 1; mm1 = 1; ms1 = 2; mm2 = 1; ms2 = 2; \[CapitalOmega]1 = 1; \[CapitalOmega]2 = 1;  snr = 10^(15/10)(*SNR=0-20dBs*); (*Invoking the implemented EGBMGF module*) B = (1/(2 Gamma[p] Gamma[mm1] Gamma[ms1] Gamma[mm2] Gamma[         ms2])) S[{{{p}, {}}, {{}, {}}}, {{{1}, {}}, {{mm1,         ms1}, {0}}}, {{{1}, {}}, {{mm2, ms2}, {0}}}' {(       mm1 ms1)/(\[CapitalOmega]1 snr q), (       mm2 ms2)/(\[CapitalOmega]2 snr q)}]; (*END*) 

enter code here(Extended Generalized Bivariate Meijer G-Function (EGBMGF))

Clear All; (*Exception*) S::Inconsistentcoeffs = "Inconsistent Coeffs!"; S[{ast_, bst_}, {as_, bs_}, {at_, bt_}, {zs_, zt_}] := Module[{},    (*Gamma product terms with only 000s000 as argument with other \ parameters*)    Pas = Function[u,       Product[Gamma[1 - as[[1, n]] + u], {n, 1, Length[as[[1]]]}]];    Qas = Function[u,       Product[Gamma[as[[2, n]] - u], {n, 1, Length[as[[2]]]}]];    Pbs = Function[u,       Product[Gamma[bs[[1, n]] - u], {n, 1, Length[bs[[1]]]}]];    Qbs = Function[u,       Product[Gamma[1 - bs[[2, n]] + u], {n, 1, Length[bs[[2]]]}]];    Ms = Function[u, Pas[u] Pbs[u]/(Qas[u] Qbs[u])];    (*Gamma product terms with only 't' as argument with other \ parameters*)    Pat = Function[u,       Product[Gamma[1 - at[[1, n]] + u], {n, 1, Length[at[[1]]]}]];    Qat = Function[u,       Product[Gamma[at[[2, n]] - u], {n, 1, Length[at[[2]]]}]];    Pbt = Function[u,       Product[Gamma[bt[[1, n]] - u], {n, 1, Length[bt[[1]]]}]];    Qbt = Function[u,       Product[Gamma[1 - bt[[2, n]] + u], {n, 1, Length[bt[[2]]]}]];    Mt = Function[u, Pat[u] Pbt[u]/(Qat[u] Qbt[u])];     (*Gamma product terms with only 's+    t' as argument with other parameters*)    Past =      Function[u,       Product[Gamma[ast[[1, n]] + u], {n, 1, Length[ast[[1]]]}]];    Qast =      Function[u,       Product[Gamma[1 - ast[[2, n]] - u], {n, 1, Length[ast[[2]]]}]];     Qbst =      Function[u,       Product[Gamma[bst[[2, n]] + u], {n, 1, Length[bst[[2]]]}]];    Mst = Function[u, Past[u]/(Qast[u] Qbst[u])];     (*Contour limiters(Depends on numerator Gamma arguments i.e.it \ must be half of the least valued Gamma arguments)*)    Rs = 1/4;    Rt = 1/4;    (*Assignments and Declarations*)    Zs = zs;    Zt = zt;    W = 50;    (*Final Evaluation*)    Print["Numerical Integration:"];    value =      1/(2 \[Pi]I)^2 NIntegrate[       MT[s, t] Zs^s Zt^t, {s, Rs - IW, Rs + IW}, {t, Rt - IW,         Rt + IW}];     (*Returning back the value*)    Return[value];    ]; (*End of EGBMGF*) 

But I have always (1) as answer (which is not correct)! any help please, I beginner with Mathematica Thank you so much Wafaa