根据hitran数据库计算气体吸收截面 源码

hitran数据库就不介绍了, 能看到这篇文章的就一定知道是干嘛的. 气体的吸收谱线数据可以使用javahawks抓取, 国内也有些基于hitran数据库查询气体谱线的工具, 但是一般都没有计算的源码, 本文下面的源码是从harvard.edu中转过来的, 所以正确性就不用怀疑了, 附上原文链接先~
hitran数据库计算气体吸收截面

下面是源码, 另存为f文件, 找个fortran编译器编译运行即可:

program cross_sections_from_hitran

implicit real*8 (a - h, o - z)
character*132 input, output
! spectral lines
parameter (maxlines = 60000)
! number of molecules in hitran
parameter (maxmols = 39)
! maximum number of isotopologues per molecule
parameter (maxiso = 8)
! number of points in spectrum
parameter (maxpoints = 600000)

dimension mol (maxlines), iso (maxlines), sigma0 (maxlines), &
  strnth (maxlines), einstein (maxlines), alpha (maxlines), &
  elow (maxlines), coeff (maxlines)
! maximum # species and isotopomers
dimension qpower (maxmols), amu (maxmols, maxiso)
dimension pos (maxpoints), spec (maxpoints), specmod (maxpoints)
dimension voigtx (maxpoints), v (maxpoints)
! constants
pi = 3.14159265358979d0
c = 2.99792458d10
p0 = 1013.25d0
t0 = 296.d0 ! hitran standard
! codata 2002 constants
h = 6.6260693d-27
an = 6.0221415d23
r = 82.057463d0 ! derived
rk = 1.3806505d-16
du = 2.6867773d16 ! derived
c2 = 1.4387752d0

write (*, '(5x, a)') 'enter hitran input file'
read (*, '(a)') input
write (*, '(5x, a)') 'enter cross section output file'
read (*, '(a)') output
open (unit = 1, file = input, status = 'old')
open (unit = 2, file = output, status = 'unknown')

! read logicals and parameters controlling calculations
!
! molnum is the hitran molecule number.
!
! start, step, and npoints define the calculation grid. in order to avoid
! undersampling the grid should be at least as fine as 1/3 of the smallest
! gaussian fwhm of a spectral line, or 0.5550 times the smallest gaussian hw1e.
! this insures that the maximum sampling error is <1.24e-4 of the full-scale
! line-shape. later: add automatic check for undersampling wrt voigt widths
!
! press is the pressure in atmospheres (millibars / 1013.25). temp is the
! temperature in degrees kelvin.
!
! nvoigt is the number of grid points to each side of a spectral line for
! performing the voigt calculation. assuming calculation to <= 1e-4 is desired,
! nvoigt should be the greater of (1) 100 * hwhm_max / step, where hwhm_max is
! the largest lorentzian hwhm of a line and step is the grid spacing;
! (2) 3.035 * hw1e_max / step, where hw1e_max is the largest gaussian hw1e of a
! line.
!
! hw1e is the gaussian slit with at 1/e intensity. hw1e = the full width at
! half-maximum, fwhm, * 0.60056 = fwhm / (2 * sqrt (ln 2)).
!
! nmod gives the option not to write out humongous spectral files by printing
! every nmod^th spectral value
write (*, '(a)') &
  'read molnum, start, step, npoints, press, temp, nvoigt, hw1e, nmod'
read (*, *) molnum, start, step, npoints, press, temp, nvoigt, hw1e, nmod
if (npoints .gt. maxpoints) pause 'npoints .gt. maxpoints'

! initialize cross sections and calculate the spectrum grid.
spec = 0.d0
specmod = 0.d0
do i = 1, npoints
  pos (i) = start + (i - 1) * step
end do

! setup hitran
call hitran_setup (qpower, amu)

! read lines
i = 1
40 read (1, '(i2, i1, f12.6, 2e10.3, 2f5.4, f10.4, f4.2)', end = 50) mol_temp, &
iso_temp, sigma0_temp, strnth_temp, einstein_temp, alpha_temp, dum, elow_temp, &
coeff_temp
write (*, *) mol_temp, molnum
! only count lines for the specified molecule
  if (mol_temp .eq. molnum) then
    mol (i) = mol_temp
    iso (i) = iso_temp
    sigma0 (i) = sigma0_temp
    strnth (i) = strnth_temp
    einstein (i) = einstein_temp
    alpha (i) = alpha_temp
    elow (i) = elow_temp
    coeff (i) = coeff_temp
    i = i + 1
  end if
  go to 40
50 nlines = i - 1
write (*, *) 'nlines = ', nlines
if (nlines .gt. maxlines) pause 'nlines .gt. maxlines'

! loop over lines to fill out cross section array
do i = 1, nlines
  vg = 4.30140d-7 * sigma0 (i) * dsqrt (temp / amu (mol (i), iso (i)))
  voigta = press * alpha (i) * (t0 / temp)**coeff (i) / vg
  ratio1 = dexp (-elow (i) * c2 / temp) - dexp (-(sigma0 (i) + &
  elow (i)) * c2 / temp)
  ratio2 = dexp (-elow (i) * c2 / t0) - dexp (-(sigma0 (i) + &
  elow (i)) * c2 / t0)
  ratio = ratio1 / ratio2 * ((t0 / temp)**qpower (mol (i)))
  vnorm = ratio * strnth (i) / vg
  nvlo = max (1, ibin (sigma0 (i), pos, npoints) - nvoigt)
  nvhi = min (npoints, ibin (sigma0 (i), pos, npoints) + nvoigt)
  do j = nvlo, nvhi
    voigtx (j) = (pos (j) - sigma0 (i)) / vg
  end do
  call voigt (voigtx, voigta, v, nvlo, nvhi)
  do j = nvlo, nvhi
    spec (j) = spec (j) + vnorm * v (j)
  end do
end do

! convolve with instrument function
call gauss (pos, spec, specmod, npoints, hw1e)

! output results
do i = 1, npoints
  if (nmod .eq. 1 .or. mod (i, nmod) .eq. 1) write (2, '(f10.3, 1p2e13.5)') &
  pos (i), spec (i), specmod (i)
end do

close (unit = 1)
close (unit = 2)
stop
end
!
subroutine hitran_setup (qpower, amu)

implicit real*8 (a - h, o - z)
parameter (maxmols = 39)
parameter (maxiso = 8)

dimension qpower (maxmols), amu (maxmols, maxiso)
qpower = 1.5d0
amu = 0.d0

! hitran numbers:
!   h2o (1)
!   co2 (2)
!    o3 (3)
!   n2o (4)
!    co (5)
!   ch4 (6)
!    o2 (7)
!    no (8)
!   so2 (9)
!   no2 (10)
!   nh3 (11)
!  hno3 (12)
!    oh (13)
!    hf (14)
!   hcl (15)
!   hbr (16)
!    hi (17)
!   clo (18)
!   ocs (19)
!  h2co (20)
!  hocl (21)
!    n2 (22)
!   hcn (23)
! ch3cl (24)
!  h2o2 (25)
!  c2h2 (26)
!  c2h6 (27)
!   ph3 (28)
!  cof2 (29)
!   sf6 (30)
!   h2s (31)
! hcooh (32)
!   ho2 (33)
!     o (34)
!clono2 (35)
!   no+ (36)
!  hobr (37)
!  c2h4 (38)
! ch3oh (39)
! ch3br (40)
! ch3cn (41)
!   cf4 (42)

do i = 1, maxmols
  if (i .eq. 34) qpower (i) = 0.d0
  if (i .eq. 2 .or. i .eq. 4 .or. i .eq. 5 .or. i .eq. 7 .or. i .eq. 8 &
  .or. i .eq. 13 .or. i .eq. 14 .or. i .eq. 15 .or. i .eq. 16 .or. i &
  .eq. 17 .or. i .eq. 18 .or. i .eq. 19 .or. i .eq. 22 .or. i .eq. 23 &
  .or. i .eq. 26 .or. i .eq. 36) qpower (i) = 1.d0
end do
amu (1, 1) = 18.010565  
amu (1, 2) = 20.014811  
amu (1, 3) = 19.014780
amu (1, 4) = 19.016740  
amu (1, 5) = 21.020985  
amu (1, 6) = 20.020956  
amu (2, 1) = 43.989830  
amu (2, 2) = 44.993185  
amu (2, 3) = 45.994076  
amu (2, 4) = 44.994045  
amu (2, 5) = 46.997431  
amu (2, 6) = 45.997400  
amu (2, 7) = 47.998322  
amu (2, 8) = 46.998291  
amu (3, 1) = 47.984745  
amu (3, 2) = 49.988991  
amu (3, 3) = 49.988991  
amu (3, 4) = 48.988960  
amu (3, 5) = 48.988960  
amu (4, 1) = 44.001062  
amu (4, 2) = 44.998096  
amu (4, 3) = 44.998096  
amu (4, 4) = 46.005308  
amu (4, 5) = 45.005278  
amu (5, 1) = 27.994915  
amu (5, 2) = 28.998270  
amu (5, 3) = 29.999161  
amu (5, 4) = 28.999130  
amu (5, 5) = 31.002516  
amu (5, 6) = 30.002485  
amu (6, 1) = 16.031300  
amu (6, 2) = 17.034655  
amu (6, 3) = 17.037475  
amu (7, 1) = 31.989830  
amu (7, 2) = 33.994076  
amu (7, 3) = 32.994045  
amu (8, 1) = 29.997989  
amu (8, 2) = 30.995023  
amu (8, 3) = 32.002234  
amu (9, 1) = 63.961901  
amu (9, 2) = 65.957695  
amu (10, 1) = 45.992904  
amu (11, 1) = 17.026549  
amu (11, 2) = 18.023583  
amu (12, 1) = 62.995644  
amu (13, 1) = 17.002740  
amu (13, 2) = 19.006986  
amu (13, 3) = 18.008915  
amu (14, 1) = 20.006229  
amu (15, 1) = 35.976678  
amu (15, 2) = 37.973729  
amu (16, 1) = 79.926160  
amu (16, 2) = 81.924115  
amu (17, 1) = 127.912297  
amu (18, 1) = 50.963768  
amu (18, 2) = 52.960819  
amu (19, 1) = 59.966986  
amu (19, 2) = 61.962780  
amu (19, 3) = 60.970341  
amu (19, 4) = 60.966371  
amu (19, 5) = 61.971231  
amu (20, 1) = 30.010565  
amu (20, 2) = 31.013920  
amu (20, 3) = 32.014811  
amu (21, 1) = 51.971593  
amu (21, 2) = 53.968644  
amu (22, 1) = 28.006147  
amu (23, 1) = 27.010899  
amu (23, 2) = 28.014254  
amu (23, 3) = 28.007933  
amu (24, 1) = 49.992328  
amu (24, 2) = 51.989379  
amu (25, 1) = 34.005480  
amu (26, 1) = 26.015650  
amu (26, 2) = 27.019005  
amu (27, 1) = 30.046950  
amu (28, 1) = 33.997238  
amu (29, 1) = 65.991722  
amu (30, 1) = 145.962492
amu (31, 1) = 33.987721
amu (31, 2) = 35.983515
amu (31, 3) = 34.987105
amu (32, 1) = 46.005480
amu (33, 1) = 32.997655
amu (34, 1) = 15.994915
amu (35, 1) = 96.956672
amu (35, 2) = 98.953723
amu (26, 1) = 29.997989
amu (37, 1) = 95.921076
amu (37, 2) = 97.919027
amu (38, 1) = 28.031300
amu (38, 2) = 29.034655
amu (39, 1) = 32.026215

return
end
!
subroutine voigt (x, a, v, nvlo, nvhi)

! the following calculated voigt values at all grid values for each point
! subroutine voigt (x, a, v, nx)
! voigt first and second derivatives commented out
! subroutine voigt (x, a, v, dv, d2v, nx)

implicit real*8 (a - h, o - z)
integer capn
real*8 lamda
logical b
parameter (ndim = 600000)
dimension x (ndim), v (ndim)
! dimension x (ndim), v (ndim), dv (ndim), d2v (ndim)
save
data sqrtpi /1.77245385090551d0/, twooverpi /0.63661977236758d0/, &
  fouroverpi /1.27323954473516d0/

! a = 0.
if (a .lt. 1.0d-8) then
  do i = nvlo, nvhi
! do i = 1, nx
    v (i) = dexp (-x (i)**2) / sqrtpi
!   dv (i) = -2.0d0 * x (i) * v (i)
!   d2v (i) = (4.0d0 * x (i) * x (i) - 2.0d0) * v (i)
  end do
! add lorentzian check here, for speed
else
! coefficient for second derivative
! c2v = 4.0d0 * a * a + 2.d0
  sfac = 1.0d0 - a / 4.29d0
  do i = nvlo, nvhi
! do i = 1, nx
    absx = dabs (x (i))
    if ((a .lt. 4.29d0) .and. (absx .lt. 5.33d0)) then
      s = sfac * dsqrt (1.d0 - (x (i) / 5.33d0)**2)
      h = 1.6d0 * s
      h2 = 2.0d0 * h
      capn = 6.d0 + 23.0d0 * s
      lamda = h2**capn
      nu = 9.0d0 + 21.0d0 * s
          else
      h = 0.0d0
      capn = 0
      nu = 8
    end if
    b = (h .eq. 0.0d0) .or. (lamda .eq. 0.0d0)
    r1 = 0.0d0
    r2 = 0.0d0
    s1 = 0.0d0
    s2 = 0.0d0
    n = nu
    do in = 1, nu + 1
      np1 = n + 1
      t1 = a + h + dfloat (np1) * r1
      t2 = absx - dfloat (np1) * r2
      c = .5d0 / (t1 * t1 + t2 * t2)
      r1 = c * t1
      r2 = c * t2
      if ((h .gt. 0.0d0) .and. (n .le. capn)) then
        t1 = lamda + s1
        s1 = r1 * t1 - r2 * s2
        s2 = r2 * t1 + r1 * s2
        lamda = lamda / h2
      end if
      n = n - 1
    end do
    if (b) then
      v (i) = twooverpi * r1
!     dv (i) = fouroverpi * (a * r2 - absx * r1)
    else
      v (i) = twooverpi * s1
!     dv (i) = fouroverpi * (a * s2 - absx * s1)
    end if
!   dv (i) = -dsign (dv (i), x (i))
!   d2v (i) = fouroverpi * a - (c2v + 4.d0 * x (i) * x (i)) * &
!   v (i) - 4.d0 * x (i) * dv (i)
  end do
end if

return
end
!
integer function ibin (target, array, nentries)

! binary search in an array of real numbers in increasing order.
! returned is the number of the last entry which is less than target, or
! 0 if not within array. (this was written to find values enclosing
! target for a linear interpolation scheme.) 4/9/84 john lavagnino;
! adapted from jon bentley, cacm february 1984, vol. 27, no. 2, p. 94.

implicit real*8 (a - h, o - z)
parameter (maxpoints = 600000)
dimension array (maxpoints)
integer upper

lower = 0
upper = nentries + 1
do while (lower + 1 .ne. upper)
  middle = (lower + upper) / 2
  if (array (middle) .lt. target) then
    lower = middle
  else
    upper = middle
  end if
end do

! at this point, either array (lower) <= target <= array (upper), or
! lower = 0, or upper = nentries + 1 (initial values).
if (lower .gt. 0 .and. upper .ne. nentries + 1) then
  ibin = lower
else
  ibin = 0
end if

end function ibin
!
subroutine gauss (pos, spec, specmod, npoints, hw1e)

! convolves input spectrum with gaussian slit function of specified hw1e
! (half-width at 1/e intensity). Assumes input spectrum has constant
! spacing (will need to modify or replace with nm version when and if
! needed).

implicit real*8 (a - h, o - z)
parameter (maxpts = 600000)
parameter (maxslit = 10000)
dimension pos (maxpts), spec (maxpts), specmod (maxpts), slit (maxpts)

write (*, *) 'hw1e = ', hw1e
if (hw1e .eq. 0.d0) then
  write (*, *) ' no gaussian convolution applied.'
  do i = 1, npoints
    specmod (i) = spec (i)
  end do
  return
end if

emult = - 1.d0 / (hw1e**2)
delpos = pos (2) - pos (1)

! apply slit function convolution
write (*, *) 'applying slit function convolution'
! calculate slit function values out to 0.001 times x0 value, normalize
! so that sum = 1.

slitsum = 1.d0
slit0 = 1.d0
do i = 1, maxslit
  slit (i) = exp (emult * (delpos * i)**2)
  slitsum = slitsum + 2.d0 * slit (i)
  if (slit (i) / slit0 .le. 0.001) then
    nslit = i
    write (*, *) ' nslit = ', nslit
    go to 40
  end if
end do
40 if (i .eq. maxslit) write (*, *) 'maxslit is too small'

slit0 = slit0 / slitsum
do i = 1, nslit
  slit (i) = slit (i) / slitsum
end do

! convolve spectrum. don't reflect at endpoints (for now).
do i = 1, npoints
  specmod (i) = slit0 * spec (i)
  do j = 1, nslit
    nlo = i - j
    nhi = i + j
    if (nlo .ge. 1) specmod (i) = specmod (i) + slit (j) * spec (nlo)
    if (nhi .le. npoints) specmod (i) = specmod (i) + slit (j) * &
    spec (nhi)
  end do
end do

return
end