function C_water=water_select(H,SNR,Nt,Nr,Lt);
%**************************************************************************%

%**************************************************************************
%clear all
%clc
%format long

%SNR= 20;

%Nr=6;
%Nt=6;
%Lt=4;
%H=wgn(Nr,Nt,0,'complex');
 
        % [VIN,EIG]=eig(H*H'); %extract eigenvalues
         [U,S,V]=svd(H);
         eigen=diag(S);
         M=size(find(eigen),1);
         gamma=zeros(M,1);
         snr=10^(SNR/10);
         r=M;%set rank = number of antennas (full rank)
         p=1;%initial count
     for i=1:r
 %        mu=getmu(r,SNR,T,p,M);%determine mu value (see eqn. (1.30))
          temp2=0;
          for k=1:r-p+1
            temp1=1/(eigen(k)^2);
            temp2=temp2+temp1;
            temp1=0;
          end
          if r-p+1~=0
              mu=(M/(r-p+1))*(1+(1/snr)*temp2);
          else 
              mu=0;
          end
        gamma(i,1)=mu-(M/snr)*(1/(eigen(i)^2)); %calculate gamma (see eqn (1.31))
         if gamma(i,1)<0
             gamma(i,1)=0;% if gamma < 0, set it to zero i.e. discard it
             p=p+1;% increment count
             mu=0;%clear register
         else
            mu=0; % if gamma >0 store it and clear register
         end
     end
   
    if(M<=Nt)
        gam=[gamma;zeros(Nt-M,1)];
    else
        gam=gamma(1:Nt,1);
    end
    EIGNEW=diag(gam);
%solve eqn (1.26)using the determinant form
    K=V*EIGNEW*V';
    kw=diag(K);
    [kw_opt,kw_index]=sort(kw);
    hw_opt(1:Lt,1)=kw_index(Nt-Lt+1:Nt,1);
    
    hw_optr(1:Lt,1)=kw_opt(Nt-Lt+1:Nt,1);
    H4=[];
    C_water=0;
    for ki=1:length(hw_opt)
       H4=[H4 H(:,hw_opt(ki))];
      % C_water=C_water+log2(hw_optr(ki)*/(10^(SNR/10)*Lt));
    end
   C_water=log2(det(eye(Nr,Nr)+(10^(SNR/10)).*H4*H4'./Lt));
  %I+(SNR/M)*diag(gamma).*diag(eigen);
   %  a=det(a);
	% y(K)=log2(a);