SpectraRust/tlusty/extracted/butler.f
2026-03-19 14:05:33 +08:00

142 lines
5.1 KiB
Fortran

SUBROUTINE BUTLER (NI,NJ,T,U0,COL,IERR)
C =======================================
C
C Rate coefficients for collisional excitation of hydrogen
C by electrons. Interpolates in Table 3 of Przybilla & Butler
C (2004, ApJ).
C
C
C Input:
C NI Principal quantum number lower level
C NJ "" upper level
C T Temperature
C U0 =h*nu/K/T
C Output:
C COL collisional rate (cm3 s-1)
C IERR error flat (0=ok, 1=T exceeds table range,
C 2=NI higher than 6 or lower than 1
C NJ higher than 7 or lower than 2)
C
INCLUDE 'IMPLIC.FOR'
DIMENSION COLSTR(16,21),TREF(16)
DATA (TREF(I), I=1,16) /
* 2.5d3, 5d3, 7.5d3, 1d4, 1.5d4, 2d4, 2.5d4, 3d4, 4d4, 5d4, 6d4,
* 8d4, 1d5, 1.5d5, 2d5, 2.5d5 /
DATA ((COLSTR(I,J),J=1,21),I=1,16) /
C J=1,21 corresponds to (NI,NJ)={(1,2),(1,3),...,(1,NL),(2,3),...}
C where NL=7 (higher n covered in Table)
C I=1,16 corresponds to T={2.5e3,5e3,7.5e3,1e4,1.5e4,2e4,2.5e4,3e4,
C 4e4,5e4,6e4,8e4,1e5,1.5e5,2e5,2.5e5}
* 6.40d-1, 2.20d-1, 9.93d-2, 4.92d-2, 2.97d-2, 5.03d-2, 2.35d+1,
* 1.07d+1, 5.22d+0, 2.91d+0, 5.25d+0, 1.50d+2, 7.89d+1, 4.13d+1,
* 7.60d+1, 5.90d+2, 2.94d+2, 4.79d+2, 1.93d+3, 1.95d+3, 6.81d+3,
* 6.98d-1, 2.40d-1, 1.02d-1, 5.84d-2, 4.66d-2, 6.72d-2, 2.78d+1,
* 1.15d+1, 5.90d+0, 4.53d+0, 7.26d+0, 1.90d+2, 9.01d+1, 6.11d+1,
* 1.07d+2, 8.17d+2, 4.21d+2, 7.06d+2, 2.91d+3, 3.24d+3, 1.17d+4,
* 7.57d-1, 2.50d-1, 1.10d-1, 7.17d-2, 6.28d-2, 7.86d-2, 3.09d+1,
* 1.23d+1, 6.96d+0, 6.06d+0, 8.47d+0, 2.28d+2, 1.07d+2, 8.21d+1,
* 1.25d+2, 1.07d+3, 5.78d+2, 8.56d+2, 4.00d+3, 4.20d+3, 1.50d+4,
* 8.09d-1, 2.61d-1, 1.22d-1, 8.58d-2, 7.68d-2, 8.74d-2, 3.38d+1,
* 1.34d+1, 8.15d+0, 7.32d+0, 9.27d+0, 2.70d+2, 1.26d+2, 1.01d+2,
* 1.37d+2, 1.35d+3, 7.36d+2, 9.66d+2, 5.04d+3, 4.95d+3, 1.73d+4,
* 8.97d-1, 2.88d-1, 1.51d-1, 1.12d-1, 9.82d-2, 1.00d-1, 4.01d+1,
* 1.62d+1, 1.04d+1, 9.17d+0, 1.03d+1, 3.64d+2, 1.66d+2, 1.31d+2,
* 1.52d+2, 1.93d+3, 1.02d+3, 1.11d+3, 6.81d+3, 6.02d+3, 2.03d+4,
* 9.78d-1, 3.22d-1, 1.80d-1, 1.33d-1, 1.14d-1, 1.10d-1, 4.71d+1,
* 1.90d+1, 1.23d+1, 1.05d+1, 1.08d+1, 4.66d+2, 2.03d+2, 1.54d+2,
* 1.61d+2, 2.47d+3, 1.26d+3, 1.21d+3, 8.20d+3, 6.76d+3, 2.21d+4,
* 1.06d+0, 3.59d-1, 2.06d-1, 1.50d-1, 1.25d-1, 1.16d-1, 5.45d+1,
* 2.18d+1, 1.39d+1, 1.14d+1, 1.12d+1, 5.70d+2, 2.37d+2, 1.72d+2,
* 1.68d+2, 2.96d+3, 1.46d+3, 1.29d+3, 9.29d+3, 7.29d+3, 2.33d+4,
* 1.15d+0, 3.96d-1, 2.28d-1, 1.64d-1, 1.33d-1, 1.21d-1, 6.20d+1,
* 2.44d+1, 1.52d+1, 1.21d+1, 1.14d+1, 6.72d+2, 2.68d+2, 1.86d+2,
* 1.72d+2, 3.40d+3, 1.64d+3, 1.34d+3, 1.02d+4, 7.70d+3, 2.41d+4,
* 1.32d+0, 4.64d-1, 2.66d-1, 1.85d-1, 1.45d-1, 1.27d-1, 7.71d+1,
* 2.89d+1, 1.74d+1, 1.31d+1, 1.17d+1, 8.66d+2, 3.19d+2, 2.08d+2,
* 1.78d+2, 4.14d+3, 1.92d+3, 1.41d+3, 1.15d+4, 8.26d+3, 2.52d+4,
* 1.51d+0, 5.26d-1, 2.95d-1, 2.01d-1, 1.53d-1, 1.31d-1, 9.14d+1,
* 3.27d+1, 1.90d+1, 1.38d+1, 1.18d+1, 1.04d+3, 3.62d+2, 2.24d+2,
* 1.81d+2, 4.75d+3, 2.15d+3, 1.46d+3, 1.26d+4, 8.63d+3, 2.60d+4,
* 1.68d+0, 5.79d-1, 3.18d-1, 2.12d-1, 1.58d-1, 1.34d-1, 1.05d+2,
* 3.60d+1, 2.03d+1, 1.44d+1, 1.19d+1, 1.19d+3, 3.98d+2, 2.36d+2,
* 1.83d+2, 5.25d+3, 2.33d+3, 1.50d+3, 1.34d+4, 8.88d+3, 2.69d+4,
* 2.02d+0, 6.70d-1, 3.55d-1, 2.29d-1, 1.65d-1, 1.35d-1, 1.29d+2,
* 4.14d+1, 2.23d+1, 1.51d+1, 1.19d+1, 1.46d+3, 4.53d+2, 2.53d+2,
* 1.85d+2, 6.08d+3, 2.61d+3, 1.55d+3, 1.49d+4, 9.21d+3, 2.90d+4,
* 2.33d+0, 7.43d-1, 3.83d-1, 2.39d-1, 1.70d-1, 1.37d-1, 1.51d+2,
* 4.56d+1, 2.37d+1, 1.56d+1, 1.20d+1, 1.67d+3, 4.95d+2, 2.65d+2,
* 1.86d+2, 6.76d+3, 2.81d+3, 1.57d+3, 1.63d+4, 9.43d+3, 3.17d+4,
* 2.97d+0, 8.80d-1, 4.30d-1, 2.59d-1, 1.77d-1, 1.39d-1, 1.93d+2,
* 5.31d+1, 2.61d+1, 1.63d+1, 1.19d+1, 2.08d+3, 5.68d+2, 2.83d+2,
* 1.87d+2, 8.08d+3, 3.15d+3, 1.61d+3, 1.97d+4, 9.78d+3, 3.94d+4,
* 3.50d+0, 9.79d-1, 4.63d-1, 2.71d-1, 1.82d-1, 1.39d-1, 2.26d+2,
* 5.83d+1, 2.78d+1, 1.68d+1, 1.19d+1, 2.39d+3, 6.16d+2, 2.94d+2,
* 1.86d+2, 9.13d+3, 3.36d+3, 1.62d+3, 2.27d+4, 1.00d+4, 4.73d+4,
* 3.95d+0, 1.06d+0, 4.88d-1, 2.81d-1, 1.85d-1, 1.40d-1, 2.52d+2,
* 6.23d+1, 2.89d+1, 1.71d+1, 1.19d+1, 2.62d+3, 6.51d+2, 3.02d+2,
* 1.87d+2, 1.00d+4, 3.51d+3, 1.63d+3, 2.54d+4, 1.02d+4, 5.50d+4 /
NL=7
IERR=0
COL=0.0d0
IF (T.LT.2.5d3.OR.T.GE.2.5d5) IERR=1
IF (NI.LT.1.OR.NI.GT.NL-1.OR.NJ.LT.2.OR.NJ.GT.NL) IERR=2
IF (IERR.EQ.0) THEN
J=0
DO I=1,NI-1
J=J+(NL-I)
END DO
DO K=I+1,NJ
J=J+1
END DO
C find out nearest points in TREF
ILOW=1
DO WHILE (T.GE.TREF(ILOW+1))
ILOW=ILOW+1
END DO
IHIG=16
DO WHILE (T.LT.TREF(IHIG-1))
IHIG=IHIG-1
END DO
IF (IHIG.EQ.ILOW) IHIG=IHIG+1
C interpolate linearly (log-log) the collision strength
SL=LOG10(COLSTR(IHIG,J))-LOG10(COLSTR(ILOW,J))
SL=SL/(LOG10(TREF(IHIG))-LOG10(TREF(ILOW)))
OR=LOG10(COLSTR(IHIG,J))-SL*LOG10(TREF(IHIG))
COL=LOG10(T)*SL+OR
COL=10.**COL
C derive the rate
COL=8.631d-6/(2.d0*NI**2)/SQRT(T)*EXP(-U0)*COL
END IF
RETURN
END