function lgKappa = opac (X, Z, lgT, lgRHO);
%   opac : Interpolation Of Opacity By Chemical Composition and T and Rho
%	use function HUNT8 replaced with dsearchn
%   use function POLINT8
%	use  BSPLVB8 rewrited

%	Integer i, j, k, X_ind, T_ind, Q_ind, l, m, n
%	Real*8  lgT, lgRHO, lgQ, X, lgKappa, Op(4), DlgKappa, Z
%	Real*8  B_T(4), B_Q(4)
%        Real*8  opr, opt, opx, opz, Opacity0
%        Integer iexp, ier
global T_net Q_net X_net OpaCoef
%  Convert arguments
lgQ = lgRHO + 18.D0 - 3.D0*lgT;
%  Check for limits
if (lgT <= T_net(1) || lgT >= T_net(end) || ...
      lgQ <= Q_net(1) || lgQ >= Q_net(end))
  warning ('parameters outside the net ranges');
  return
end
%  Interpolation of  spline coefficients
T_ind=max(find(T_net-lgT<=0));
Q_ind=max(find(Q_net-lgQ<=0));
%    This strings suppose T_net sorted increasingly and lgT inside T_net
%    the same for Q
B_T = BSPLVB8 (T_net, 4, lgT, T_ind);
B_Q = BSPLVB8 (Q_net, 4, lgQ, Q_ind);

%  Find appropriate X
    X_ind = max(find(X_net-X <= 0));
    X_Num=length(X_net);
    X_ind = min (max( X_ind-1, 1 ), X_Num-3);

%  Compute
   Op=zeros(4,1);
for kk = 1:4,
Op_X=squeeze(OpaCoef(kk - 1 + X_ind,T_ind-3:T_ind,Q_ind-3:Q_ind));
Op(kk)=B_T*Op_X*B_Q';
%    for ki = 1:4,
%	   for kj = 1:4,
%	      kl = ki - 1 + T_ind-3;
%	      km = kj - 1 + Q_ind-3;
%	      kn = kk - 1 + X_ind;
%	      Op(kk) = Op(kk) + B_T(ki)*B_Q(kj)*Op_X(kl,km);
%       end
%    end
end

%  Interpolation between X
p = polyfit(X_net(X_ind:X_ind+3),Op,3);
lgKappa = polyval(p,X);
%
function BIATX = BSPLVB8( T, Jhigh, x, Left )
%  Calculation of non-zero B-splines in x point.
%  Variables:
%    output BIATX(Jhigh) - nonzero values
%    input: T - net of knots of spline (all available)
%           Jhigh - order of spline (polinomial degree plus 1)
%           x - point of calculation
%           Left - index in T for left point in respect of x (T(Left)<=x)
%      Integer ki, kj, jp1
%       Jmax = 20;  - non actual in ML, needed for description of Delta
%      Real*8 DeltaL(Jmax),DeltaR(Jmax), Saved, Term, Zero
    Zero = 0.0d0;

    kj = 1; % index: order of curently calculated spline
    BIATX(1) = 1; % value for spline of first order (constant)
if ( kj >= Jhigh ), return; end; % this check for first order case only 
while ( kj < Jhigh ),
    jp1 = kj + 1;
    DeltaR(kj) = T(Left+kj) - x;
    DeltaL(kj) = x - T(Left+1-kj);
    Saved = Zero;

    for ki = 1:kj,
         Term = BIATX(ki) / ( DeltaR(ki) + DeltaL(jp1-ki) );
         BIATX(ki) = Saved + DeltaR(ki)*Term;
         Saved = DeltaL(jp1-ki)*Term;
    end
      BIATX(jp1) = Saved;
      kj = jp1;
end

