The solar diffusion routine by Thoul et al. is available by
anonymous ftp at the following address:
ftp ftp.sns.ias.edu
cd pub
cd thoul
get routine.f
The routine is based on the calculations developed in the following paper:
Thoul, A. A., Bahcall, J. N., and Loeb, A. 1994, Ap.J. 421, 828.
Here are a few additional explanatory notes about the diffusion routine.
Questions and/or comments are welcome. Please send them to
thoul@astro.ulg.ac.be
-------------------------------------------------------------------------
-------------------------------------------------------------------------
Given M elements of atomic mass numbers A(M), charges Z(M),
and mass fractions X(M), and given the Coulomb logarithms
CL(M,M), the routine returns the diffusion coefficients AP(M),
AT(M), and AX(M,M).
(NOTE: in the routine, the maximum number of elements MMAX has been
set to 20. If the number of elements is larger, this number
should be increased. And NMAX=2*MMAX+2 should also be changed.)
The diffusion function $\xi$, as defined in TBL (Thoul \& al, Ap.J.
421, p.828, 1994), can then be obtained by multiplying these coefficients
by the appropriate gradients:
\xi(M)=AP(M) {d lnp/dr} + AT(M) {d lnT/dr} + \sum_N AX(M,N) {d lnC(N)/dr}
where N is different from 2 (Helium) and M (electrons).
Note that these gradients are dimensionless. The radii are in units of
R_\sun, as defined in TBL.
-------------------------------------------------------------------------
If mass fraction gradients are available instead of concentration fraction
gradients, it is necessary to convert the latter.
To get (dimensionless) concentration gradients from (dimensionless) mass
fraction gradients:
(here, GRADX contains the mass fraction gradients dlnX(N)/dr, and GRADC
contains the concentration fraction gradients dlnC(N)/dr)
TEMP=0.
DO I=1,M-1
TEMP=TEMP+Z(I)*X(I)/A(I)*GRADX(I)
ENDDO
DO J=1,M-1
GRADC(J)=TEMP
ENDDO
TEMP=0.
DO I=1,M-1
TEMP=TEMP+Z(I)*X(I)/A(I)
ENDDO
DO J=1,M-1
GRADC(J)=GRADX(J)-GRADC(J)/TEMP
ENDDO
Alternatively, using the formula given in footnote #2 in TBL, it is possible to
express the diffusion function $\xi$ in terms of the mass fraction gradients
rather than the concentration gradients (NOTE: again, these are dimensionless
gradients):
\xi(M)=AP(M) {d lnp/dr} + AT(M) {d lnT/dr} + \sum_N AX(M,N) {d lnX(N)/dr}
-\sum_N AX(M,N))*(\sum_J Z(J)X(J)/A(J) {d lnX(J)/dr})/(\sum_I Z(I)X(I)/A(I))
where N is different from 2 (Helium) and M (electrons), and the sums over I and J
are over all ionic species.
TEMP1=0.
TEMP2=0.
DO I=1,M-1
TEMP1=TEMP1+Z(I)*X(I)/A(I)
TEMP2=TEMP2+Z(I)*X(I)/A(I)*GRADX(I)
ENDDO
DO I=1,M
SUM1=0.
SUM2=0.
DO J=1,M-1
IF (J.NE.2) THEN
SUM1=SUM1+AX(I,J)*GRADX(J)
SUM2=SUM2+AX(I,J)
ENDIF
ENDDO
XI(I)=AP(I)*GRADP + AT(I)*GRADT + SUM1 - SUM2*TEMP2/TEMP1
ENDDO
To obtain the dimensionless diffusion velocity, simply multiply \xi by
T^{5/2}/\rho where T is in units of 10^7 K, and \rho is in units of 100 g/cm^3.
To obtain the velocity in cgs units, simply multiply the dimensionless
velocity by R_\sun/\tau_0 (see TBL).
-------------------------------------------------------------------------
To obtain the values of the Coulomb logarithms used in TBL,
add the following lines just before calling the diffusion subroutine
in your main program.
NOTES: (1) in the following, RHO and T are the mass density and the
temperature in CGS UNITS!
(2) there is a typo in TBL paper: in eq. 9, defining the
Coulomb logarithms, the quantity in the parenthesis
(4kT\lambda/Z_s Z_t e^2) should be raised to the power 1.2.
c additional variables to be declared:
c intermediate variables:
REAL ZXA,AC,NI,CZ,XIJ,NE,AO,LAMBDAD,LAMBDA,C(M)
c calculate concentrations from mass fractions:
ZXA=0.
DO I=1,M-1
ZXA=ZXA+Z(I)*X(I)/A(I)
ENDDO
DO I=1,M-1
C(I)=X(I)/(A(I)*ZXA)
ENDDO
C(M)=1.
c calculate density of electrons (NE) from mass density (RHO):
AC=0.
DO I=1,M
AC=AC+A(I)*C(I)
ENDDO
NE=RHO/(1.6726E-24*AC)
c calculate interionic distance (AO):
NI=0.
DO I=1,M-1
NI=NI+C(I)*NE
ENDDO
AO=(0.23873/NI)**(1./3.)
c calculate Debye length (LAMBDAD):
CZ=0.
DO I=1,M
CZ=CZ+C(I)*Z(I)**2
ENDDO
LAMBDAD=6.9010*SQRT(T/(NE*CZ))
c calculate LAMBDA to use in Coulomb logarithm:
LAMBDA=MAX(LAMBDAD,AO)
c calculate Coulomb logarithms:
DO I=1,M
DO J=1,M
XIJ=2.3939E3*T*LAMBDA/ABS(Z(I)*Z(J))
CL(I,J)=0.81245*LOG(1.+0.18769*XIJ**1.2)
ENDDO
ENDDO
---------------------------------------------------------------------