function [C,EbN0dB] = capa_nakagami(sigma2,m);
%CAPA_NAKAGAMI Numerically evaluates the Shannon limit for m-Nakagami flat fading channel.
%
% CAPA_NAKAGAMI(SIGMA2,M) tries to approximate the theoretical achievable rate
% over a m-Nakagami channel with Gaussian inputs with channel noise variance SIGMA2.
% The parameter M characterizes the fading. Available values are {1,2,3,4,5,6,7,8,9,10,Inf}
%
% [C,EBN0DB] = CAPA_NAKAGAMI(SIGMA2,M) also returns the mimimum signal-to-noise per bit EBN0DB
% in decibel required to achieve the capacity.
%
%
% Example 1 (AWGN):
% SNRdB = 0;m = Inf;[C,EbN0dB] = capa_nakagami(10^(-SNRdB/10),m)
%
% Example 2 (Rayleigh):
% SNRdB = 0.987;m = 1;[C,EbN0dB] = capa_nakagami(10^(-SNRdB/10),m)
% 
% Example 3: (3-Nakagami)
% SNRdB = 0.367;m = 3;[C,EbN0dB] = capa_nakagami(10^(-SNRdB/10),m)

% Initalization
gamma0 = 1/(sigma2);


% Capacity calculation
if ~isfinite(m); % Gaussian case
	C = .5*log2(1+1/sigma2);
else;
switch m;
	case 1;	% Rayleigh fading		
		C = 0.5/log(2)*exp(1/gamma0)*expint(1/gamma0);
	case {2,3,4,5,6,7,8,9,10};
		Str = 'times(times(power(x,m-n-1),log(x)),exp(-1/gamma0*m*x))';
		f = inline(Str,'x','gamma0','m','n');
		

		tmp = 0;
		for n = 0:m-1;
			tmp = tmp+(-1)^n * nchoosek(m-1,n)*quad(f,1,5/sigma2,[],[],gamma0,m,n);
		end;
		C = 0.5/log(2)/gamma(m)*(m/gamma0)^m*exp(m/gamma0)*tmp;

	otherwise; disp('Warning: The parameter m is not 1,2,3,4,5,6,7,8,9,10 or Inf.');C = NaN;

end;
end;

if imag(C)~=0;C = NaN;end;
% EbN0 determination
EbN0dB = -10*log10(2*C*sigma2);