function es=mir_mod(Clay, mv, TCel, roD, freq)

%Clay =  13.0; mg = 0.179; TCel =  20;  roD = 1.4; txt = 'r'; 
%freq=1;

hPlank = 6.626e-34*1e9;
kBoltz = 1.380e-23;

    epsilZer = [240.61 100.55 82.62
                110.62   nan  73.29
                 82.73  64.48 42.25 
                 10.55   nan  6.12];               
               
    betavZer =  [-0.23  -0.25  0.04
                  0.01   nan   0.18
                 -0.98  -0.57 -2.10
                 -4.70   nan  -2.30]*1.e-3;
                
    delHsR   = [-338.65 -274.09  1758.89
                 415.87    nan   2447.41
                 -61.35 -217.89  1530.89
                  93.22    nan   3456.40];
                
    delSsR   = [ -10.46  -8.65  1.64
                  -9.02   nan   4.49
                  -9.45  -8.28  0.72
                 -10.23   nan   8.39];            
                 
    epsilInf = [20.98;  4.31; 17.06;  4.48]; 
    betavInf = [-1.10; -0.12; -0.29; -0.34]*1.e-3;
   TCelseInf = [20.00; 20.00; -15.0; -15.0];    
    
             
  TCelseZer = [ 20   20   20
                 20   nan  20
                -20  -20  -20
                -20   nan -20];
            
  
    roBUI = [1;1;0.917];

es=0;
mg=mv/roD;

       if TCel >= 0
               FPQ = log((epsilZer(1:2,:)-1)./(epsilZer(1:2,:)+2));
            Tse0PQ =     TCelseZer(1:2,:);
          betav0PQ =      betavZer(1:2,:);
          
              delH =        delHsR(1:2,:);
              delS =        delSsR(1:2,:);              
       else
               FPQ = log((epsilZer(3:4,:)-1)./(epsilZer(3:4,:)+2));
            Tse0PQ =     TCelseZer(3:4,:);
          betav0PQ =      betavZer(3:4,:);
                    
              delH =        delHsR(3:4,:);
              delS =        delSsR(3:4,:);          
       end
       
        tauPQ = hPlank/kBoltz/(TCel+273.16)*exp(delH/(TCel+273.16)-delS);       
       
         E0   = exp(FPQ - betav0PQ.*(TCel - Tse0PQ));        
       eps0PQ = (1 + 2*E0)./(1 - E0); 
     
       if TCel >= 0
                 FPQ = log((epsilInf(1:2)-1)./(epsilInf(1:2)+2));
            TseInfPQ =     TCelseInf(1:2);
          betavInfPQ =      betavInf(1:2); 
       else
                 FPQ = log((epsilInf(3:4)-1)./(epsilInf(3:4)+2));
            TseInfPQ =     TCelseInf(3:4);
          betavInfPQ =      betavInf(3:4);          
       end 

         EInf   = exp(FPQ - betavInfPQ.*(TCel - TseInfPQ));        
       epsInfPQ = (1 + 2*EInf)./(1 - EInf);  % 2 x 1 vector              

        if TCel >= 0
            sigm = ([-32.90; -23.8; nan] + [5.80; 5.20;     nan]*Clay) + ...
                   ([ -1.10; -1.20; nan] + [0.15; 0.16;     nan]*Clay)*(TCel - 20);
             mg1 =  7.9 + 3.4*Clay;
        else
            sigm = ([ -8.90; nan; -2.30] + [1.40;  nan;  0.3500]*Clay) + ...
                   ([ -0.59; nan;  0.41] + [0.11;  nan; -9.2e-4]*Clay)*(TCel + 15);
             mg1 = (6.6 + 1.6*Clay) + (-7.1 + 2.3*Clay)*exp(TCel/6.9);               
        end

        sigm = sigm*1.e-3;
         mg1 =  mg1*1.e-3;
       

  
  
       epsP = nan(3,1);   % complex _eps from (9)-(10) of the ref.
       
   epsP(1) =   (eps0PQ(1,1) - eps0PQ(1,2))/(1 - 1i*2*pi*freq*tauPQ(1,1)) ...
            +  (eps0PQ(1,2) - eps0PQ(1,3))/(1 - 1i*2*pi*freq*tauPQ(1,2)) ...
            +  (eps0PQ(1,3) - epsInfPQ(1))/(1 - 1i*2*pi*freq*tauPQ(1,3)) ...
            + epsInfPQ(1);
   
        if TCel >= 0, j = 2; else j = 3; end
        
   epsP(j) =   (eps0PQ(2,1) - eps0PQ(2,3))/(1 - 1i*2*pi*freq*tauPQ(2,1)) ...
            +  (eps0PQ(2,3) - epsInfPQ(2))/(1 - 1i*2*pi*freq*tauPQ(2,3)) ...
            + epsInfPQ(2);
        

        if    mg <= mg1

          
           t1 = 0.4 + (sqrt(epsP(1)) - 1)/roBUI(1)*mg;
  
           t  = 1 + roD*t1 ;
           
           t2 = mg *sigm(1)/roBUI(1);
                                            
           es = t*t + 1i*roD*t2/(2*pi*freq*8.854e-3)  ;
                
 
        else
 
      
           t1 = 0.4 + (sqrt(epsP(1)) - 1)/roBUI(1)* mg1 ...
                    + (sqrt(epsP(j)) - 1)/roBUI(j)*(mg - mg1);
                
           t  = 1 + roD*t1 ;
           
           t2 = mg1*sigm(1)/roBUI(1) + (mg - mg1)*sigm(j)/roBUI(j);           
                                            
           es = t*t + 1i*roD*t2/(2*pi*freq*8.854e-3);                
    
 
        end
 

  
  end