- Email: [email protected]
Radiative transfer effects in Cepheid atmospheres
J. Quant. Spectrosc. Radiat. Transfer. Vol. 11, pp. 647-653. Pergamon Press 1971. Printed in Great Britain
R A D I A T I V E T R A N S F E R EFFECTS IN CEPHEID ATMOSPHERES* CECIL G. DAVIS, JR. Los Alamos Scientific Laboratory, Los Alamos, N.M. Abstract--The transfer equation, without retardation and considering scattering as absorption, is transformed by the moments expansion method for use in studying the dynamics of Cepheid atmospheres. The 0th and 1st moment equations are solved implicitly, using a variable Eddington factor (f) and 13 frequencygroups. Thef is calculated each time step by a characteristic ray plane transfer method. Results in a "bump" model (3.52M0, T©rf = 4585°K) of R. CHRISTYare presented and are compared to an observed 10.9 period Cepheid V.X. Per.") 1. INTRODUCTION THE USE of the variable Eddington equations of radiative transfer in a model of a Cepheid atmosphere should result in a considerable improvement in understanding over the usual diffusion models t2,3,4) without too large an increase in computing time. The necessity of calculating the Eddington factor to obtain the exact transfer solution is realized by using an explicit plane transfer method which adds only 30 per cent to the computing time. The variable Eddington equations and the calculation o f f are given in Section 2. A short description of the model used and the effects due to the transfer equations is described in Section 3. A detail discussion of the model and comparisons to a diffusion model will be given elsewhere. In Section 4 a general discussion of the effects in the Cepheid atmosphere due to the use of the correct transfer equations is given along with a comparison to an observed Cepheid, V.X. Per. 2. THE TRANSFER EQUATIONS The equations of radiative transfer used, and the moments methods for obtaining their implicit solution, start from the basic transfer equation: 1 tgl v
(1--[A
c Ot -~ - - r
2) ~I v
1
aI v
~+P~r
v
f
+(~a+tr~)Iv = SV+a~
P(/2,p')P(ff)d/2'
(1)
-1
where I v is the intensity at r, in the solid angle f~ at a time t, i.e. IV(r, ~, t), S v is the source function, P(#,/2') is the scattering phase function, which can be carried along in the method to be described even to a 1st order C o m p t o n (and inverse Compton) approximation, ts~ In our case, since scattering is small, it will be considered only as absorption. The wave * Done under the auspices of the U.S. Atomic Energy Commission. 647
CECIL G. DAVIS. JR.
648
nature of the radiative transfer equation is contained in the retardation term (l/c aZ”at). This can also be carried along easily, but since the radiation field is quasi-static, i.e. l/c aZvat is small compared to the other terms, it will be dropped and the resulting transfer equation to be solved is:
arv -+pg
bpz
( Iap r
=
S’-a&I’
(2)
where a; = a’ +a and ai is absorption corrected for induced emission by ai = a,(1 -e -hv’kT).‘The”approach used in solving this equation is the moments expansion method where, we have : M, = 2n J’_ 1 1: dp or 1
M,, = 2~
s
I, dp = E,lc
-1
(3)
E, is the radiation energy density, Fr the net radiation flux, and pr the radiation pressure tensor, higher order expansion introduces combinations of these terms. We consider only the 0th and 1st moment terms of the transfer equation (l), namely V.F? = -caL(E,-T) and
(4)
Here L.T.E. (Local Thermodynamic Equilibrium) has been assumed in the use of the source function S’ = a’B’(T) (Kirchhoff’s Law). In Eddington’s approximation the equations are closed to 1st order by assuming that, V . Fr = *VE,, but, in our case, we will set V . F, = f VE, and assume a correct value off can be found. The transport equation (1) is then correctly approximated by the 0th and 1st moment equations if f can be calculated separately. The final equations are : 0th
1 a(r*P) = ca;(T-Ej r’ar
(5)
and 1st
(6)
and with the energy equation, ai(cE’ - 4nB’( T)) dv 0
(7)
Radiative transfer effects in Cepheid atmospheres
649
give a complete set for the frequency dependent radiation transfer solution if the f f values can be found. The use of the superscript (v) for the frequency dependence now replaces the radiation subscript (r). These equations are like the usual diffusion equations and can be made tri-diagonal; i.e. set in the form zaizl - 1Lq -~n+1"k-BiE~+ l + C i +
r;. + 1 -L T D i = O,
where the constants will not be given here, then; by backward substitution - {6~ one solves for E "÷1 and from (6) we find F "+1. I f f v is set equal to ½ then"
c dE ~ 3a v Or
F v --
and if E,(J'~ E ~ dr) is in equilibrium at the material temperature (T) i.e., equal to aT'*, we have the usual diffusion flux:
F,= f
ac OT 4
(8)
The method we use for obtaining f f is based on a plane geometry solution to the transport equation by FREEMAN and DAVIS. (7) Using double-gaussian expansion on angle and an integration by parts of (1) we have;
I
OB~'
[aBVI_SBV I
ri=[s-~-~},+l_&l, &l,+,
]e_a/2
(9)
+[1L,_/ ~ 0B~
where z = aAx and A is the zone optical thickness (zi-zi-1). The second term in (9) is an approximation to 02BV/O~2. By gaussian integration then we have K
E~
=
E A,I• 1
and K
P* = ~ A~,p2 I~v
(10)
1
then f v = Ev/pv where A, and # are the gaussian weights and angles (# = cos 0). The solution for K = 3 and 6 was compared to a complete spherical transfer solution using a direct integration scheme {7} and a Y line (impact ray) for each zone. Even with K = 3 the fluxes from the two methods agreed to within a few percent. The method of coupling in frequency groups to solve (5)-(7) is due to FREEMAN et al. °~ and uses an auxiliary energy equation, i.e.,
D(bj~) o,
4aT3b~[ -
,-I
I p,
~3E.~ DV]
4aT 3
(11)
650
CECILG. DAVaS,JR.
where Em and Pm are the material energy and pressure and the equation is linearized by setting ~ = T 4. The bSs are the interval of the normalized Planck function: •i+l
b~
lIB(T)
t ~ 2hv3/c2(e - h v / k T - 1)- 1 dr. Vi
The effect of (11) is to estimate a temperature at n + 1 from the difference in emission and absorption at n. These implicit frequency grouped solutions are then summed to give E,, F,, and P,. In the actual coupling of the radiation field to the hydrodynamics, a test is made on the optical thickness in each zone by the frequency group and if the zone is thick we use V. P , the diffusion form, but if the zone is thin we use (7) which limits correctly to the Planck emission term if a Planck mean is used, i.e., E~=
tr~B " dv vi
B" dv B v(r°~ = ~p(To)BV(To).
(12)
0
At present we use a Rosseland mean tr• instead of % (Note: trp can be 10 to 100 times trR in some cases) in (12). The Rosseland mean is, lla~R = f.~;+'lltrvdB'l~'r dv S~ ~BV/dz dv
(13)
3. TRANSFER EFFECTS IN A CEPHEID ATMOSPHERE The details of the Cepheid model used and the comparisons to the diffusion model will be presented elsewhere. The diffusion methods of calculating the dynamics of pulsating stars are described in detail by Cox et al. ~2) and CHRISTY. (a) The calculational approach used in this paper is similar except now use is made of the monochromatic absorption coefficients in the transfer equation instead of the grey Rosseland means (also supplied by Cox and STEWART(a)) in the diffusion equation. The calculation time increases linearly with the number of frequency groups used. The results of a 13 group transfer calculation, requiring approximately one hour C.D.C. 6600 computer time for each period, is compared to a grey diffusion calculation, which takes approximately 20 min for each period, in Fig. 1. The results of the transfer calculation of the Christy model agree in general with the diffusion results. A shifting of the velocity bump at the outer edge of the atmosphere is the largest difference in the hydrodynamics. The velocity during a period at the depth where the iron lines may be formed does not show this shift. The velocity bump is in phase with the luminosity bump. In static atmosphere calculations the constant flux approximation is used. The effect of this approximation on the emergent spectrum also appears to be in good agreement with the dynamic results in Cepheids. (9) Apparently there are no appreciable differences between transfer results and static results in this case. To test the question of using snapshots of diffusion models for spectral results the transfer spectrum is compared to snapshots of a diffusion model near light minimum and light maximum and with the Planck distributions in Fig. 2. These comparisons indicate a general agreement with the true spectrum at light minimum. The result near light maximum though shows a considerable
Radiative transfer effects in Cepheid atmospheres
651
u.-~ •
-
,
/7i\
1:5 groups f = l / 5 Grey d i f f u s i o n
....
/
/
,/ 'i_,
l/i\
13 grouptransport Grey diffusion
--"
~1
/'
/7
'
',,.2' 05
,~
li I
/
,'7 it
'
".-,, IX x,.,/
I0
Phase,
*
//~
/4
/,
..,= .>
f-
1.0
0.5
t/Tr o
Phase,
1.5
t/~ o
FIG. 1. The diffusion calculation is compared to the transfer results of luminosity and velocity versus phase, with the variable Eddington factor ( f ) set to ~ and as calculated by a plane characteristic ray method.
increase in flux above the region where the Balmer jump occurs (3.4 eV). To understand this effect a model using 19 frequency groups with increased resolution around the Balmer edge was run. The result, as shown in Fig. 2, is to apparently remove the increased flux above the Balmer limit at light maximum. Actually the increase in flux above the Planck distribution remains, it is shifted in phase and occurs only in a narrow range of phase (0.14).2); around light maximum (not shown in Fig. 2c). 4. D I S C U S S I O N
AND
COMPARISON
TO V.X. P E R
The effects of radiative transfer in the atmosphere of a Cepheid have been calculated using the variable Eddington method with an exact determination of f. The importance of an improvement in angular resolution over the diffusion method appears slight (obs. (a) Transport
(c)
(b)
Sr~pshot
,--~-",,
~
i
"--.,.Ligh t maximum
"
maximum
\-<
\\Lig
v,
',,,
• igh,
minimum 2
4
6
Photon energy, eV
8
10
l
2
i
4
\~, t
6
\ "
J
8
10
Photon energy, eV
FIG. 2. The transfer spectrum, at light m i n i m u m and light m a x i m u m is
2
4
6
Photon energy,
8 eV
compared to snapshots of
diffusion model results (a) and Planck distributions (b). A resolved Balmer jump (19 frequency group) spectra is given in (c).
652
CECILG. DAVIS,JR.
Fig. 1). There does seem to be a significant effect due to the frequency discretization around light maximum but only over a narrow range in phase. The overall effects of using the correct transfer equations, in the model described, appear to be small. These effects will have to be re-evaluated in stars of higher luminosity to mass ratio such as R. V. Tauri or W. Virginis stars. A comparison of the Christy model to observations made in the UBV of MITCHELL e t al. ~1~ for V.X. Per are made in Fi~ 3. These results have not been corrected for line blanketing effects or adjusted for the spectral band response. The light amplitudes differ considerably but the location of the bump for a 10 day Cephcid is in about the right phase relation using a mass that is approximately ½ the evolutionary mass, the mass suggested by Christy. Higher resolution spectral scans will make these comparisons more meaningful.
VX Per
I0 d 8 9
86 88
I §F-,
90
,q/
92 9.4 96 98
-
\
:
V
i 0~i~ \
J / - 0 - \ b---.o-\ /--
~o. ~ ~ .
m 14
enl
~
OIi
~
1.2:t 141 0
"~'X3"A3"x2"/ I
02
1
04
I
06
I
0.8
I
1.0
I
12
14
Phase
FIG. 3. Comparison to
V.X. Per. The
points are the calculated values where V was adjusted to centre the visible luminosity on the observed.
Acknowledgements---The author would like to acknowledge the considerable help of Mrs. J. BENDTin setting up and running the models as well as doing much of the programming and the many discussions with A. N. Cox,
J. P. Cox and D. KING.
Radiative transfer effects in Cepheid atmospheres
653
REFERENCES 1. 2. 3. 4. 5. 6.
R. MITCHELL, B. IRIARTE, D. STEIMETZand H. JOI-INSON, Bol. Obs. Tonantzintla y Tacubaya 3, 153 (1964). J. Cox, A. Cox, K. OLSON, D. KING and D. EILERS, Ap J. 144, 1038 (1966). R. CHRISTY,Rev. rood. Phys. 36, 555 (1964). R. STOBEIE,Mon. Not. R. astr. Soc. 144, 461 (1969). B. FREEMAN,et al., 3SIR-29 vol I (Systems Science & Software, San Diego, Calif.) (1968). R. RICHTMYERand K. MORTON, Difference Methods for Initial Value Problems, Interscience Pub. Inc., N.Y., p. 193 (1957). 7. B. FREEMANand C. DAVIS, JR., A.F.W.L.-TR6-143, vols. III and IV (1966). 8. A. Cox, J. STEWARTand D. EILERS, Ap J. Suppl. 94 (1965). 9. C. KELLERand P. MUTSCHLECHNER,Ap J. 161, 217 (1970).
R A D I A T I V E T R A N S F E R EFFECTS IN CEPHEID ATMOSPHERES* CECIL G. DAVIS, JR. Los Alamos Scientific Laboratory, Los Alamos, N.M. Abstract--The transfer equation, without retardation and considering scattering as absorption, is transformed by the moments expansion method for use in studying the dynamics of Cepheid atmospheres. The 0th and 1st moment equations are solved implicitly, using a variable Eddington factor (f) and 13 frequencygroups. Thef is calculated each time step by a characteristic ray plane transfer method. Results in a "bump" model (3.52M0, T©rf = 4585°K) of R. CHRISTYare presented and are compared to an observed 10.9 period Cepheid V.X. Per.") 1. INTRODUCTION THE USE of the variable Eddington equations of radiative transfer in a model of a Cepheid atmosphere should result in a considerable improvement in understanding over the usual diffusion models t2,3,4) without too large an increase in computing time. The necessity of calculating the Eddington factor to obtain the exact transfer solution is realized by using an explicit plane transfer method which adds only 30 per cent to the computing time. The variable Eddington equations and the calculation o f f are given in Section 2. A short description of the model used and the effects due to the transfer equations is described in Section 3. A detail discussion of the model and comparisons to a diffusion model will be given elsewhere. In Section 4 a general discussion of the effects in the Cepheid atmosphere due to the use of the correct transfer equations is given along with a comparison to an observed Cepheid, V.X. Per. 2. THE TRANSFER EQUATIONS The equations of radiative transfer used, and the moments methods for obtaining their implicit solution, start from the basic transfer equation: 1 tgl v
(1--[A
c Ot -~ - - r
2) ~I v
1
aI v
~+P~r
v
f
+(~a+tr~)Iv = SV+a~
P(/2,p')P(ff)d/2'
(1)
-1
where I v is the intensity at r, in the solid angle f~ at a time t, i.e. IV(r, ~, t), S v is the source function, P(#,/2') is the scattering phase function, which can be carried along in the method to be described even to a 1st order C o m p t o n (and inverse Compton) approximation, ts~ In our case, since scattering is small, it will be considered only as absorption. The wave * Done under the auspices of the U.S. Atomic Energy Commission. 647
CECIL G. DAVIS. JR.
648
nature of the radiative transfer equation is contained in the retardation term (l/c aZ”at). This can also be carried along easily, but since the radiation field is quasi-static, i.e. l/c aZvat is small compared to the other terms, it will be dropped and the resulting transfer equation to be solved is:
arv -+pg
bpz
( Iap r
=
S’-a&I’
(2)
where a; = a’ +a and ai is absorption corrected for induced emission by ai = a,(1 -e -hv’kT).‘The”approach used in solving this equation is the moments expansion method where, we have : M, = 2n J’_ 1 1: dp or 1
M,, = 2~
s
I, dp = E,lc
-1
(3)
E, is the radiation energy density, Fr the net radiation flux, and pr the radiation pressure tensor, higher order expansion introduces combinations of these terms. We consider only the 0th and 1st moment terms of the transfer equation (l), namely V.F? = -caL(E,-T) and
(4)
Here L.T.E. (Local Thermodynamic Equilibrium) has been assumed in the use of the source function S’ = a’B’(T) (Kirchhoff’s Law). In Eddington’s approximation the equations are closed to 1st order by assuming that, V . Fr = *VE,, but, in our case, we will set V . F, = f VE, and assume a correct value off can be found. The transport equation (1) is then correctly approximated by the 0th and 1st moment equations if f can be calculated separately. The final equations are : 0th
1 a(r*P) = ca;(T-Ej r’ar
(5)
and 1st
(6)
and with the energy equation, ai(cE’ - 4nB’( T)) dv 0
(7)
Radiative transfer effects in Cepheid atmospheres
649
give a complete set for the frequency dependent radiation transfer solution if the f f values can be found. The use of the superscript (v) for the frequency dependence now replaces the radiation subscript (r). These equations are like the usual diffusion equations and can be made tri-diagonal; i.e. set in the form zaizl - 1Lq -~n+1"k-BiE~+ l + C i +
r;. + 1 -L T D i = O,
where the constants will not be given here, then; by backward substitution - {6~ one solves for E "÷1 and from (6) we find F "+1. I f f v is set equal to ½ then"
c dE ~ 3a v Or
F v --
and if E,(J'~ E ~ dr) is in equilibrium at the material temperature (T) i.e., equal to aT'*, we have the usual diffusion flux:
F,= f
ac OT 4
(8)
The method we use for obtaining f f is based on a plane geometry solution to the transport equation by FREEMAN and DAVIS. (7) Using double-gaussian expansion on angle and an integration by parts of (1) we have;
I
OB~'
[aBVI_SBV I
ri=[s-~-~},+l_&l, &l,+,
]e_a/2
(9)
+[1L,_/ ~ 0B~
where z = aAx and A is the zone optical thickness (zi-zi-1). The second term in (9) is an approximation to 02BV/O~2. By gaussian integration then we have K
E~
=
E A,I• 1
and K
P* = ~ A~,p2 I~v
(10)
1
then f v = Ev/pv where A, and # are the gaussian weights and angles (# = cos 0). The solution for K = 3 and 6 was compared to a complete spherical transfer solution using a direct integration scheme {7} and a Y line (impact ray) for each zone. Even with K = 3 the fluxes from the two methods agreed to within a few percent. The method of coupling in frequency groups to solve (5)-(7) is due to FREEMAN et al. °~ and uses an auxiliary energy equation, i.e.,
D(bj~) o,
4aT3b~[ -
,-I
I p,
~3E.~ DV]
4aT 3
(11)
650
CECILG. DAVaS,JR.
where Em and Pm are the material energy and pressure and the equation is linearized by setting ~ = T 4. The bSs are the interval of the normalized Planck function: •i+l
b~
lIB(T)
t ~ 2hv3/c2(e - h v / k T - 1)- 1 dr. Vi
The effect of (11) is to estimate a temperature at n + 1 from the difference in emission and absorption at n. These implicit frequency grouped solutions are then summed to give E,, F,, and P,. In the actual coupling of the radiation field to the hydrodynamics, a test is made on the optical thickness in each zone by the frequency group and if the zone is thick we use V. P , the diffusion form, but if the zone is thin we use (7) which limits correctly to the Planck emission term if a Planck mean is used, i.e., E~=
tr~B " dv vi
B" dv B v(r°~ = ~p(To)BV(To).
(12)
0
At present we use a Rosseland mean tr• instead of % (Note: trp can be 10 to 100 times trR in some cases) in (12). The Rosseland mean is, lla~R = f.~;+'lltrvdB'l~'r dv S~ ~BV/dz dv
(13)
3. TRANSFER EFFECTS IN A CEPHEID ATMOSPHERE The details of the Cepheid model used and the comparisons to the diffusion model will be presented elsewhere. The diffusion methods of calculating the dynamics of pulsating stars are described in detail by Cox et al. ~2) and CHRISTY. (a) The calculational approach used in this paper is similar except now use is made of the monochromatic absorption coefficients in the transfer equation instead of the grey Rosseland means (also supplied by Cox and STEWART(a)) in the diffusion equation. The calculation time increases linearly with the number of frequency groups used. The results of a 13 group transfer calculation, requiring approximately one hour C.D.C. 6600 computer time for each period, is compared to a grey diffusion calculation, which takes approximately 20 min for each period, in Fig. 1. The results of the transfer calculation of the Christy model agree in general with the diffusion results. A shifting of the velocity bump at the outer edge of the atmosphere is the largest difference in the hydrodynamics. The velocity during a period at the depth where the iron lines may be formed does not show this shift. The velocity bump is in phase with the luminosity bump. In static atmosphere calculations the constant flux approximation is used. The effect of this approximation on the emergent spectrum also appears to be in good agreement with the dynamic results in Cepheids. (9) Apparently there are no appreciable differences between transfer results and static results in this case. To test the question of using snapshots of diffusion models for spectral results the transfer spectrum is compared to snapshots of a diffusion model near light minimum and light maximum and with the Planck distributions in Fig. 2. These comparisons indicate a general agreement with the true spectrum at light minimum. The result near light maximum though shows a considerable
Radiative transfer effects in Cepheid atmospheres
651
u.-~ •
-
,
/7i\
1:5 groups f = l / 5 Grey d i f f u s i o n
....
/
/
,/ 'i_,
l/i\
13 grouptransport Grey diffusion
--"
~1
/'
/7
'
',,.2' 05
,~
li I
/
,'7 it
'
".-,, IX x,.,/
I0
Phase,
*
//~
/4
/,
..,= .>
f-
1.0
0.5
t/Tr o
Phase,
1.5
t/~ o
FIG. 1. The diffusion calculation is compared to the transfer results of luminosity and velocity versus phase, with the variable Eddington factor ( f ) set to ~ and as calculated by a plane characteristic ray method.
increase in flux above the region where the Balmer jump occurs (3.4 eV). To understand this effect a model using 19 frequency groups with increased resolution around the Balmer edge was run. The result, as shown in Fig. 2, is to apparently remove the increased flux above the Balmer limit at light maximum. Actually the increase in flux above the Planck distribution remains, it is shifted in phase and occurs only in a narrow range of phase (0.14).2); around light maximum (not shown in Fig. 2c). 4. D I S C U S S I O N
AND
COMPARISON
TO V.X. P E R
The effects of radiative transfer in the atmosphere of a Cepheid have been calculated using the variable Eddington method with an exact determination of f. The importance of an improvement in angular resolution over the diffusion method appears slight (obs. (a) Transport
(c)
(b)
Sr~pshot
,--~-",,
~
i
"--.,.Ligh t maximum
"
maximum
\-<
\\Lig
v,
',,,
• igh,
minimum 2
4
6
Photon energy, eV
8
10
l
2
i
4
\~, t
6
\ "
J
8
10
Photon energy, eV
FIG. 2. The transfer spectrum, at light m i n i m u m and light m a x i m u m is
2
4
6
Photon energy,
8 eV
compared to snapshots of
diffusion model results (a) and Planck distributions (b). A resolved Balmer jump (19 frequency group) spectra is given in (c).
652
CECILG. DAVIS,JR.
Fig. 1). There does seem to be a significant effect due to the frequency discretization around light maximum but only over a narrow range in phase. The overall effects of using the correct transfer equations, in the model described, appear to be small. These effects will have to be re-evaluated in stars of higher luminosity to mass ratio such as R. V. Tauri or W. Virginis stars. A comparison of the Christy model to observations made in the UBV of MITCHELL e t al. ~1~ for V.X. Per are made in Fi~ 3. These results have not been corrected for line blanketing effects or adjusted for the spectral band response. The light amplitudes differ considerably but the location of the bump for a 10 day Cephcid is in about the right phase relation using a mass that is approximately ½ the evolutionary mass, the mass suggested by Christy. Higher resolution spectral scans will make these comparisons more meaningful.
VX Per
I0 d 8 9
86 88
I §F-,
90
,q/
92 9.4 96 98
-
\
:
V
i 0~i~ \
J / - 0 - \ b---.o-\ /--
~o. ~ ~ .
m 14
enl
~
OIi
~
1.2:t 141 0
"~'X3"A3"x2"/ I
02
1
04
I
06
I
0.8
I
1.0
I
12
14
Phase
FIG. 3. Comparison to
V.X. Per. The
points are the calculated values where V was adjusted to centre the visible luminosity on the observed.
Acknowledgements---The author would like to acknowledge the considerable help of Mrs. J. BENDTin setting up and running the models as well as doing much of the programming and the many discussions with A. N. Cox,
J. P. Cox and D. KING.
Radiative transfer effects in Cepheid atmospheres
653
REFERENCES 1. 2. 3. 4. 5. 6.
R. MITCHELL, B. IRIARTE, D. STEIMETZand H. JOI-INSON, Bol. Obs. Tonantzintla y Tacubaya 3, 153 (1964). J. Cox, A. Cox, K. OLSON, D. KING and D. EILERS, Ap J. 144, 1038 (1966). R. CHRISTY,Rev. rood. Phys. 36, 555 (1964). R. STOBEIE,Mon. Not. R. astr. Soc. 144, 461 (1969). B. FREEMAN,et al., 3SIR-29 vol I (Systems Science & Software, San Diego, Calif.) (1968). R. RICHTMYERand K. MORTON, Difference Methods for Initial Value Problems, Interscience Pub. Inc., N.Y., p. 193 (1957). 7. B. FREEMANand C. DAVIS, JR., A.F.W.L.-TR6-143, vols. III and IV (1966). 8. A. Cox, J. STEWARTand D. EILERS, Ap J. Suppl. 94 (1965). 9. C. KELLERand P. MUTSCHLECHNER,Ap J. 161, 217 (1970).