Matrix method for calculating source function, mean intensity, and flux in model stellar atmospher

Kurucz, R. L.

English

NASA Technical Reports Server

A MATRIX METHOD FOR CALCULATING- THE SOURCE FUNCTION, 
MEAN INTENSITY, AND FLUX I N  A MODEL ATMOSPHERE 
Robert L. Kurucz 
GPO PRICE $ 
CSFTI PRICE(S) $ 
June 1968 I 
I 
Hard copy (HC) 3 .cm 
Microfiche (MF) * 63” 
1 
ff 653 Julv 65 
Smithsonian Institution 
Astrophysical Observatory 
Cambridge, Massachusetts 021 38 
805-144 
https://ntrs.nasa.gov/search.jsp?R=19680019944 2020-03-12T09:48:21+00:00Z
. 
A MATRIX METHOD FOR CALCULATING THE SOURCE FUNCTION, 
MEAN INTENSITY, AND FLUX IN A MODEL ATMOSPHERE 
ROBERT L .  KURUCZ 
Smiths onian As t r o phy s i ca 1 Ob s e r va t o r y , Cambridge , Ma s s a c hus e t t  s and 
Harvard College Observatory, Cambridge, Massachusetts 
. 
Received 
ABSTRACT 
The integrals for  the mean intensity and flux are expressed  a s  the multi- 
plication of the source-function vector by mat r ix  representations of the 
integral  operators.  
equation for the source function and to an efficient calculation of the mean 
intensity and flux. 
This method leads to a rapid solution of the integral  
' I .  INTRODUCTION 
In a model-atmosphere calculation the source function S ( T  ), the mean v v  
intensity J ( T  ), and the flux H (T  ) must  be evaluated for  each i teration 
at each  frequency v and at each monochromatic optical depth T . The 
mean  intensity and flux a r e  given by the integrals 
v v  v v  
V 
T 
V 00 
0 
T 
V 
T 
V 
1 
where En(x) = J e-xz/zn dz is  the exponential integral  of o r d e r  n. The 
1 
source function is given by the integral  equation 
where a ( T  ) i s  the fraction of opacity due to scattering and B (T  ) is  the 
Planck function. 
the moment integrals J 
exponential integrals must  be evaluated for every value of I T  
will show below that these integrations can be expressed a s  the multiplication 
of the source-function vector by a matrix representation of the integral  
operator  and that this  approach leads to procedures that considerably reduce 
the integration time. 
v v  v v  
The computations involved in the evaluation of S and of 
V 
and H a r e  extremely t ime consuming i f  the 
V V 
- tl . We 
V 
11. DERIVATION O F  THE INTEGRATION MATRIX 
The moment integrals of the specific intensity can be expressed (dropping 
the v subscript)  as 
n- 1 
Mn(T) = 2 sign (t - T) S(t)  E n l t  - T I  dt . ' I  
0 
Then, 
J(T) = M 1 ( ~ )  and H ( T )  = M 2 ( ~ )  . 
We divide the integration range into N subintervals, 
j =  1 T 
j 
2 
We now assume that S can be represented by a parabola in  the interval  
( T j ,  j t l  ). 
3 N 
S(t)  = c tk- Cjki si , 
k =  1 i =  1 
can be, for  example, 
The 'jki 
where the C 
the average of the forward-parabolic and the backward-parabolic interpolation 
coefficients. 
we find that 
a r e  interpolation coefficients. 
jki 
Inserting the expression for  S(t)  into the j th t e r m  of the sum, 
3 N 
T j t l  
Mnl = z 1 sign ( T - dt E It - ~~1 c t k - ' c  CjkiSi . 
j n 
t k =  1 i =  1 j 
This expression can be simplified to  
3 N 
n l  j n l j k  
k =  1 i =  1 
where 
.r 
is a n  integral  that  can be evaluated analytically. To evaluate qnl jk we use 
the indefinite integral  
dx = - E  n t l  (x) 
3 
and integrate by par ts ,  obtaining 
We now collect t e r m s  to  achieve the fo rm of a ma t r ix  operator  on S, 
N 3  N 
Mnl = ‘nljk ‘jki ’i -~ 
j = 1  k = l  i=  1 
N 3 N  =c ‘nljk c s  jki i 
N 
j = 1  
where 
3 N  
3 “ n l  j =c ’nlik C ikj ’ 
k =  1 i =  1 
The matrix= is the desired matr ix  operator n 
M n = Z  S . n 
4 
111. APPLICATIONS 
We define the ma t r i ces  
and 
q j p s c ; '  
-2 
so that 
J = A S  
and 
H = @ S  . 
The integral  equation for S becomes 
S = (I - a) B + a n S  
o r  
(I - a n )  S = (I - a) B , 
where I is the identity matrix and a is now a diagonal matrix. 
solution is 
The formal 
s = (I - a ~ 1 - l  (I - a) B . 
The solufion of the integral  equation is  obtained by d i rec t  inversion, by 
Gauss-  Jordan elimination, o r  by Gauss-Seidel iteration. The timings for 
these methods a r e  approximately proportional, respectively, to  2N , N /2,  
3 
and the number of i terations times N-, where N is the dimension of the 
3 3  
5 
matrix. 
f o r  evaluating S for  the wide range of conditions encountered in model- 
atmosphere calculations. 
W e  have found the Gauss-Seidel method to be the fas tes t  procedure 
The integration mat r ices  A and Q can be pretabulated for  some fixed 
T set ,  where the values of T a r e  chosen to  give accurate  integrations. 
can then evaluate S v ,  J v ,  and H 
s e t  onto the fixed T set ,  integrating, and then interpolating the result ing 
values back onto the T 
fixed T set  must  be t reated separately. However, for  la rge  depths this  
problem is straightforward, since the asymptotic behavior of J 
par t iculary simple; namely, 
We 
by interpolating a and Bv f r o m  the T 
V V V 
se t .  Optical depths l a rge r  than the last value of the 
V 
and HY is 
V 
We i l lustrate  this procedure in  the following manner.  Let the superscr ip t  
in denote the internal fixed T set  and the superscr ipt  ex denote the external  
T set .  Furthermore,  let P be the interpolation operator f r o m  the T se t  
onto the fixed T set and Q be the operator  f rom the fixed T set  back onto the 
T 
can be expressed symbolically as follows: 
V V 
set .  The case where S is prescr ibed and where we wish to  find J and H 
V 
First we interpolate in  
then we integrate 
6 
and finally we interpolate out 
or 
jex = Q A P S ~ ~  
The case where a and B a r e  prescr ibed and where we wish to find S, J, and 
H can be symbolized as follows: 
We interpolate in 
in  ex a = P a  
integrate 
Sin = (I - ain A)-' (I - sin) Bin 
and interpolate out 
Sex = QSin 
7 
If the interpolation operators  are  mat r ices ,  then QAP and Q@P can be con- 
s idered the integration operators Aex and aex on the T ex set, so that 
Jex = flex Sex 
and 
Hex= Qex Sex . 
However, the calculation of these external  opera tors  is slower than if we 
simply interpolate in and out, because two additional mat r ix  multiplications 
are needed. 
The author 's  model-atmosphere program ATLAS uses  the matrix 
methods described above to  calculate S J , and H . These calculations 
a r e  made with 40 values of T~~ and 40 values of T 
routine uses  a weighted average of forward- and backward-parabolic inter-  
polation coefficients. 
by use of a modified Gauss-Seidel i teration that is terminated when the 
maximum relative change i s  less  than 10 This accuracy is sufficient for  
sin , since the accuracy of the parabolic interpolation i s  somewhat l e s s  than 
this  when only 40  points a r e  used. 
V in  . V I  v The interpolation 
The integral equation for  the source function is solved 
-5 . 
Figures  1 and 2 show the convergence propert ies  of this method for  two 
calculations of Sv  in  a model atmosphere with effective temperature  10000" K 
and log g of 2 .  Figures  3 and 4 show the corresponding values of 1 - a 
B V ,  Sva J v ,  and H v .  
method because scattering i s  a large fraction of the opacity a t  this low gravity, 
where hydrogen is mostly ionized. 
exceeded 99 .99  per  cent at the surface were  m o r e  than four Gauss-Seidel 
i terat ions required.  
V' 
This model provides a good test of the Gauss-Seidel 
Only at frequency points where scattering 
8 
The author wishes to thank D r .  Wolfgang Kalkofen for suggesting this 
investigation. Mr. Deane Peterson and Dr. Stephen St rom contributed helpful 
c r i t i c i sm and many useful discussions t o  this work. Support for  the computa- 
tions reported in  this work was provided, i n  part ,  by grant  NGR 22-024-001 
f r o m  the National Aeronautics and Space Administration. 
acknowledges support f rom the S A 0  fellowship program. 
The author a l so  
9 
T =  1000Oo 
LOG g = 2  
X = 911.77A 
L \ \ \  
\ \ \  3 
4 
5 
6 
7 
8 
-2 - 
-3 
-4 
-5 
-6 
-7 
r 
L I I I I 
-6 -5 -4  -3 -2 - I  0 I 
LOG T ' ~  
F i g .  1. - Relative change for each subiteration of the source-function 
calculation at the long-wavelength side of the Lyman discontinuity. 
the scattering fraction is 0. 9999 near the surface and > O .  5 at l a r g e  depths.  
Here 
- 2  
I 
K 
w 
5; - 4  
u 
W 
I 
v) -5 
(u 
I 
2 
0 
-6 
K 
W 
-7 
-0  
T =  1000Oo 
LOG g = 2  
X = 911.77A 
I 
-9 I I I I I I I 
-6 -5 -4 -3 -2 
LOG T, 
- - I  0 I 
Fig. 2. - The input scattering fraction, Planck function, and the resulting 
source  function, mean intensity, and flux corresponding to  the source-function 
calculation shown in Figure 1. 
I 
0 
? 
- I  
- - 2  
z 
m 
I I I I I 1 1 
- 5  
-6 
-7 
-6 -5 
T =  10000" 
LOG g = 2  
-4 -3 -2 -I 
LOG  TI^ 
0 I 
Fig. 3. - Same a s  Figure 1, for the short-wavelength side of the Balmer 
discontinuity where there  is considerably l e s s  scattering, i. e . ,  0. 0 3  at large 
depths. 
i 
0 
- I  
- 
'N 
I 
U 
W 
I- 
# 
- - 2  
I 
- -3 
I 
V 
W 
# 
(u 
I -4 I 
0 
# 
c3 
W -5 
-6 
-7 
- 
T =  1000Oo 
LOG g = 2  
A =  3647.02 A - 
-6 -5 -4  -3 -2  --I 0 I 2 
LOG T, 
Fig. 4. -Same a s  Figure 2, for the short-wavelength side of the Balmer 
discontinuity where there  is considerably l e s s  scattering, i. e. .. 0. 03 at la rge  
depths. 


A matrix method for calculating the source function, mean intensity, and flux in a model atmosphere

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server