function [kappa, gamma, m, L] = parameters(sigma,mu,thickness,sigma_0)

% kappa, gamma, m, L : model parameters as defined in
% Phys. Rev. B, 78 (2), 024113, July 2008.
% sigma : vector with grb energies, e.g. [0.1 0.2 0.25]
% m : vector with grb mobilities; e.g. [1 1 1]
% thickness : numerical grb thickness, e.g. 1e-7
% sigma_0 :  start_value for the iterative calculation, gamma will equal
% 1.5 for sigma = sigma_0 (it is best to choose sigma_0 so that gamma is between 0.7 and 4 for all elements of sigma)

% example :  [kappa, gamma, m, L] = parameters([0.1 0.2 0.25],[1 1 1],1e-7,0.25) 
% gives kappa =  1.0e-007 * [  0.0739    0.1497    0.1875 ]
% gamma = [0.7412    1.1629    1.4999]
% m =   1.50e+007
% L =  1.0e+007 * [ 1.3538    1.3356    1.3334 ]


%%%%%%%fitting for gamma < 4 (gamma(3:67))
% %%% g_function^2 -- 1/gamma
% 1/gamma =  p1*(g_function^2)^3 + p2*(g_function^2)^2 + p3*(g_function^2)^1 + p4;
%   p1 = -4.8457;
%   p2 = 11.3;
%   p3 = -8.2937;
%   p4 = 2.006;

% %% 1/gamma -- f_interf 
% f_interf = pp1*(1/gamma)^6 + pp2*(1/gamma)^5 + pp3*(1/gamma)^4 + pp4*(1/gamma)^3 + pp5*(1/gamma)^2 + pp6*(1/gamma)^1 + pp7;
%  pp1 = 0.028948
%   pp2 = -0.19627
%   pp3 = 0.53653
%   pp4 = -0.7693
%   pp5 = 0.65764
%   pp6 = -0.45665
%   pp7 = 0.28256
 

  %%for gamma < 4,fitted for gamma(3:67)
  p1 =      -5.288  
  p2 =    -0.09364  
  p3 =       9.965  
  p4 =      -8.183  
  p5 =       2.007  
       
       pp1 =     0.05676  
       pp2 =     -0.2924  
       pp3 =      0.6367  
       pp4 =     -0.7749  
       pp5 =      0.6107  
       pp6 =     -0.4324  
       pp7 =      0.2792  
  
  
m = 6*sigma_0/(thickness);
g_0 = sqrt(2)/3;
sqrt_f0 = sqrt(1/8);
a_0 = sqrt_f0/g_0;
 for k = 1 : length(sigma)
     a_star = a_0;
     a_0 = 0;
     while abs(a_0 - a_star) > 1e-9
         a_0 = a_star;
         kappa_star = sigma(k)*thickness*a_0;
         g = sigma(k)/sqrt(kappa_star*m);
         x = g^2;
         y = p1*x^4 + p2*x^3 + p3*x^2 + p4*x + p5;
         gamma_star = 1/y;
         f_interf =   pp1*y^6 + pp2*y^5 + pp3*y^4 + pp4*y^3 + pp5*y^2 + pp6*y^1 + pp7;
         a_star = sqrt(f_interf)/g;
     end
     kappa(k) = kappa_star; gamma(k) = gamma_star; 
     L(k) = (mu(k)/(thickness))*(1/a_star);
 end

