Документ взят из кэша поисковой машины. Адрес оригинального документа : http://crydee.sai.msu.ru/ftproot/users/vab/opacity/opac.m Дата изменения: Thu Mar 24 22:52:34 2005 Дата индексирования: Tue Oct 2 15:34:08 2012 Кодировка: Поисковые слова: http astrokuban.info astrokuban

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