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

101 lines
3.9 KiB
Fortran

subroutine h2minus(t,anh2,ane,fr,oph2m)
C =======================================
C
C H- free-free opacity
C
C K L Bell 1980 J. Phys. B: At. Mol. Phys. 13 1859, Table 1
C The first column is theta=5040/T(K)
C The first row are names for each row corresponding to lambda (angstroms)
C The last row for 10.0 is linearly extrapolated
C The units of everything else is 10^26 cm4/dyn-1
C
INCLUDE 'IMPLIC.FOR'
INCLUDE 'BASICS.FOR'
dimension FFthet(9),FFlamb(18),FFkapp(18,9)
data FFthet / 0.5, 0.8, 1.0, 1.2, 1.6, 2.0,
* 2.8, 3.6, 10.0 /
data nthet /9/
data FFlamb /151883., 113913., 91130., 60753.,
* 45565., 36452., 30377., 22783.,
* 18226., 15188., 11391., 9113., 7594.,
* 6509., 5696., 5063., 4142., 3505./
data nlamb /18/
data FFkapp /
* 7.16e+01,4.03e+01,2.58e+01,1.15e+01,6.47e+00,
* 4.15e+00,2.89e+00,1.63e+00,1.05e+00,7.36e-01,
* 4.20e-01,2.73e-01,1.92e-01,1.43e-01,1.10e-01,
* 8.70e-02,5.84e-02,4.17e-02,9.23e+01,5.20e+01,
* 3.33e+01,1.48e+01,8.37e+00,5.38e+00,3.76e+00,
* 2.14e+00,1.39e+00,9.75e-01,5.64e-01,3.71e-01,
* 2.64e-01,1.98e-01,1.54e-01,1.24e-01,8.43e-02,
* 6.10e-02,1.01e+02,5.70e+01,3.65e+01,1.63e+01,
* 9.20e+00,5.92e+00,4.14e+00,2.36e+00,1.54e+00,
* 1.09e+00,6.35e-01,4.22e-01,3.03e-01,2.30e-01,
* 1.80e-01,1.46e-01,1.01e-01,7.34e-02,1.08e+02,
* 6.08e+01,3.90e+01,1.74e+01,9.84e+00,6.35e+00,
* 4.44e+00,2.55e+00,1.66e+00,1.18e+00,6.97e-01,
* 4.67e-01,3.39e-01,2.59e-01,2.06e-01,1.67e-01,
* 1.17e-01,8.59e-02,1.18e+02,6.65e+01,4.27e+01,
* 1.91e+01,1.08e+01,6.99e+00,4.91e+00,2.84e+00,
* 1.87e+00,1.34e+00,8.06e-01,5.52e-01,4.08e-01,
* 3.17e-01,2.55e-01,2.10e-01,1.49e-01,1.11e-01,
* 1.26e+02,7.08e+01,4.54e+01,2.04e+01,1.16e+01,
* 7.50e+00,5.28e+00,3.07e+00,2.04e+00,1.48e+00,
* 9.09e-01,6.33e-01,4.76e-01,3.75e-01,3.05e-01,
* 2.53e-01,1.82e-01,1.37e-01,1.38e+02,7.76e+01,
* 4.98e+01,2.24e+01,1.28e+01,8.32e+00,5.90e+00,
* 3.49e+00,2.36e+00,1.74e+00,1.11e+00,7.97e-01,
* 6.13e-01,4.92e-01,4.06e-01,3.39e-01,2.49e-01,
* 1.87e-01,1.47e+02,8.30e+01,5.33e+01,2.40e+01,
* 1.38e+01,9.02e+00,6.44e+00,3.90e+00,2.68e+00,
* 2.01e+00,1.32e+00,9.63e-01,7.51e-01,6.09e-01,
* 5.07e-01,4.27e-01,3.16e-01,2.40e-01,2.19e+02,
* 1.26e+02,8.13e+01,3.68e+01,2.18e+01,1.46e+01,
* 1.08e+01,7.18e+00,5.24e+00,4.17e+00,3.00e+00,
* 2.29e+00,1.86e+00,1.55e+00,1.32e+00,1.13e+00,
* 8.52e-01,6.64e-01/
c locate position in temperature array
theta=5040./t
call locate(FFthet,nthet,theta,j,nthet)
if (j.eq.0) then
write(*,*)
write(*,'(a,f6.0,a)')
* 'Error: requested temperature is outside the ranges'
write(*,'(a)') 'h2minus:Stop'
write(*,*)
stop
end if
flamb=2.997925e18/fr
c locate position in wavelength array
call locate(FFlamb,nlamb,flamb,i,nlamb)
c linearly interpolate in frequency and temperature
if (j.eq.nthet) then
c hold values constant if off high temperature end of table
y1=FFkapp(i,j)
y2=FFkapp(i+1,j)
tt=(flamb-FFlamb(i))/(FFlamb(i+1)-FFlamb(i))
Fkappa=(1.-tt)*y1 + tt*y2
else if (i.eq.0 .or. i.eq.nlines) then
c set values to 0 if off frequency table
Fkappa=0.0
else
c interpolate linearly within table
y1=FFkapp(i,j)
y2=FFkapp(i+1,j)
y3=FFkapp(i+1,j+1)
y4=FFkapp(i,j+1)
tt=(flamb-FFlamb(i))/(FFlamb(i+1)-FFlamb(i))
uu=(theta-FFthet(j))/(FFthet(j+1)-FFthet(j))
Fkappa=(1.-tt)*(1.-uu)*y1 + tt*(1.-uu)*y2 + tt*uu*y3 +
* (1.-tt)*uu*y4
end if
pe=ane*BOLK*t
oph2m= anh2 * 1.0E-26 *pe * Fkappa
return
end