function [ber]=berfadingamc(snr,M,L,upper,lower,m)
 %All Rights Reserved.
 %NOTICE:  All information contained herein is, and remains
 %the property of Chen Dong and Southampton University.
 
% It is the code of the appendix of http://eprints.soton.ac.uk/348686
% This program extends berfading.
% snr is Eb/No. snr is signal to noise ratio.
% M is the value of MQAM. M=2 for BPSK.
% L is diversity order. (single selective diversity order)
% upper and lower are the integral thresholds.
% m is parameter of Nakagami-m channel == diversity order in berfading(snr,'qam',M,diversity order)
% e.g. for Rayleigh channel, qpsk:berfadingamc(snr,4,1,100,0,1)
% if m is >=5, there is a slight unaccurate
if M==2
    if upper==lower
        ber=erfc(sqrt(upper*10^(snr/10)))/2;
    elseif upper~=lower
    %%%%%%%denominator
    temp_gammainc1=0;temp_gammainc2=0;
    for m_1=0:m-1
    temp_gammainc1=temp_gammainc1+(upper*m)^m_1/factorial(m_1);
    temp_gammainc2=temp_gammainc2+(lower*m)^m_1/factorial(m_1);
    end
    temp_gammainc1=temp_gammainc1*exp(-upper*m);
    temp_gammainc2=temp_gammainc2*exp(-lower*m);
    temp_gammaincl=0;
    for l=0:L
    temp_gammaincl=temp_gammaincl+nchoosek(L,l)*(-1)^l*(temp_gammainc1^l-temp_gammainc2^l);
    end
    denominator1=temp_gammaincl;
    %%%%%%%denominator
    ber=int_gamma(upper,lower,m,L,2,snr)/denominator1;
    end
elseif  M>2
    if upper==lower&(M<=64)
    if M==4
         ber=erfc(sqrt(2*upper*10^(snr/10)/2))/2;
    elseif M==16
         ber=3/4*erfc(sqrt(4*upper*10^(snr/10)/10))/2+...
             2/4*erfc(3*sqrt(4*upper*10^(snr/10)/10))/2-...
             1/4*erfc(5*sqrt(4*upper*10^(snr/10)/10))/2;
    elseif M==64
         ber=7/12*erfc(sqrt(6*upper*10^(snr/10)/42))/2+...
             6/12*erfc(3*sqrt(6*upper*10^(snr/10)/42))/2-...
             1/12*erfc(5*sqrt(6*upper*10^(snr/10)/42))/2+...
             1/12*erfc(9*sqrt(6*upper*10^(snr/10)/42))/2-...
             1/12*erfc(13*sqrt(6*upper*10^(snr/10)/42))/2;
    end
    elseif upper~=lower|(upper==lower&M>64)
    if upper==lower
        upper=lower+0.0000000001;
    end

m_amc=m;
M_factor=3/(M-1);ber_QAM_matrix=0;
for i=1:log2(sqrt(M))
    length_ber_QAM_matrix=length(ber_QAM_matrix);
    for j2=1:length_ber_QAM_matrix
    for j1=length_ber_QAM_matrix+1:length_ber_QAM_matrix*2
        ber_QAM_matrix(j2,j1)=ber_QAM_matrix(j2,length_ber_QAM_matrix*2+1-j1)+1;
    end
    end
    for j2=1:length_ber_QAM_matrix
    for j1=1:length_ber_QAM_matrix
        ber_QAM_matrix(j2+length_ber_QAM_matrix,j1+length_ber_QAM_matrix)=ber_QAM_matrix(j2,j1);
    end
    end
    for j2=1:length_ber_QAM_matrix
    for j1=1+length_ber_QAM_matrix:2*length_ber_QAM_matrix
        ber_QAM_matrix(j2+length_ber_QAM_matrix,j1-length_ber_QAM_matrix)=ber_QAM_matrix(j2,j1);
    end
    end
        
end
for i=1:length(ber_QAM_matrix)
    for j=i+1:length(ber_QAM_matrix)
        positive_Q(i,j)=j-i;
        positive_Q(length(ber_QAM_matrix)+1-i,length(ber_QAM_matrix)+1-j)=j-i;
    end
end
% %%%%%%%%%find the ber_matrix
% positive_Q=(positive_Q.*2-1).*sign(positive_Q);
% negative_Q=positive_Q+2;negative_Q=negative_Q-eye(sqrt(M));
% negative_Q(:,1)=0;negative_Q(:,length(ber_QAM_matrix))=0;
% for i=1:sqrt(M)
%     for j=1:sqrt(M)
%         if i~=j
%             if j~=1&j~=sqrt(M)
%             ber_matrix(i,j)=int_gamma(100,0,1,1,positive_Q(i,j)^2*M_factor*log2(M),snr)-int_gamma(100,0,1,1,negative_Q(i,j)^2*M_factor*log2(M),snr);
%             else
%             ber_matrix(i,j)=int_gamma(100,0,1,1,positive_Q(i,j)^2*M_factor*log2(M),snr);
%             end
%         end
%     end
%     ber_matrix(i,i)=1-sum(ber_matrix(i,:));
% end
% %%%%%%%%%%%%
negative_Q=positive_Q+1;
negative_Q(:,1)=0;negative_Q(:,length(ber_QAM_matrix))=0;
for i=1:length(ber_QAM_matrix)
Q_vector(i)=sum(sum((positive_Q==i).*ber_QAM_matrix))-sum(sum((negative_Q==i).*ber_QAM_matrix));
end

Q_vector=Q_vector/(0.5*log2(M))/(0.5*sqrt(M))/2;
temp=0;
for i=1:length(ber_QAM_matrix)
temp=temp+int_gamma(upper,lower,m_amc,L,(i*2-1)^2*M_factor*log2(M),snr)*Q_vector(i);
end
%%%%%%%denominator
  temp_gammainc1=0;temp_gammainc2=0;
    for m_1=0:m-1
    temp_gammainc1=temp_gammainc1+(upper*m)^m_1/factorial(m_1);
    temp_gammainc2=temp_gammainc2+(lower*m)^m_1/factorial(m_1);
    end
    temp_gammainc1=temp_gammainc1*exp(-upper*m);
    temp_gammainc2=temp_gammainc2*exp(-lower*m);
    temp_gammaincl=0;
    for l=0:L
    temp_gammaincl=temp_gammaincl+nchoosek(L,l)*(-1)^l*(temp_gammainc1^l-temp_gammainc2^l);
    end
    denominator1=temp_gammaincl;
%%%%%%%denominator
ber=temp/denominator1;
    end%upper~=lower
end
end

%int_gamma(20,0,1,1,1/5*4,0)(upper,lower,m_amc,L,B,snr)*3/4+int_gamma(20,0,1,1,9/5*4,0)*2/4-int_gamma(20,0,1,1,25/5*4,0)/4

% simualtion of int_T1^T2 nakachannel. m is 1.2
% T1=0.32;T2=1.85;
% seg=0.0001;
% result=0;
% snr=0;
% snr_EbN0=10^(snr/10);
% m=1.2;
% for gamma_ins=T1:seg:T2
%     result=result+gamma_ins^(m-1)*exp(-m*gamma_ins)*erfc(sqrt(2*snr_EbN0*gamma_ins)/sqrt(2))/2;
% end
% result=result*m^m/gamma(m)*seg/(gammainc(T2*m,m)-gammainc(T1*m,m))
% 
%%%%%%%%%%%%%%%%%%int_0^pi/2 formula d\phi.
% result=0;
% for angle_ins=0:seg:pi/2
%     result=result+1/pi/(m+1/sin(angle_ins)^2)^m*(gammainc(T2*(m+1/sin(angle_ins)^2),m)-gammainc(T1*(m+1/sin(angle_ins)^2),m));
% end
% result=result*m^m*seg/(gammainc(T2*m,m)-gammainc(T1*m,m))