//
function  mu = cal_mu(X)
//
//  PE
   mu0 = 1.; lambda = 13.47; m = 0.294; a = .381
//   mu0 = 1;  lambda = 2.;    m = .2;    a = 2.;

   absX = abs(X);
   mu = mu0 * (1 + ( lambda*absX )^a)^((m-1)/a);

endfunction
