Add line mixing for CO2

pyrad/lbl/hitran/database.py

@@ -1,11 +1,13 @@
+from collections import OrderedDict
 from logging import getLogger
 from sqlite3 import connect
 from urllib.request import urlopen
 
 from numpy import asarray
 
 from .isotopologues import isotopologues
-from .line_parameters import PARAMETERS
+from .line_parameters import linear_molecule_quantum_numbers, linear_molecule_state_vars, \
+                             PARAMETERS
 from .molecules import molecules
 from .spectral_lines import SpectralLines
 from ...utils.database_utilities import ascii_table_records, scrub, SQL_TYPES
@@ -54,15 +56,38 @@ def create_database(self, database):
         with connect(database) as connection:
             cursor = connection.cursor()
             name = scrub(self.molecule)
-            columns = ", ".join(["{} {}".format(x.shortname, SQL_TYPES[x.dtype])
-                                 for x in self.parameters])
-            cursor.execute("CREATE TABLE {} ({})".format(name, columns))
-            value_subst = ", ".join(["?" for _ in self.parameters])
+            shortnames, types = [], []
+            for x in self.parameters:
+                if x.dtype is linear_molecule_quantum_numbers:
+                    q = getattr(self, x.shortname)[0]
+                    types += [SQL_TYPES[i] for i in linear_molecule_state_vars.values()]
+                    state = "upper" if x.api_name == "statep" else "lower"
+                    shortnames += ["{}_{}".format(i, state) for i in
+                                   linear_molecule_state_vars.keys()]
+                else:
+                    types.append(SQL_TYPES[x.dtype])
+                    shortnames.append(x.shortname)
+            print(shortnames)
+            value_subst = ", ".join(["?" for _ in shortnames])
+            header = ", ".join(["{} {}".format(shortname, type_str)
+                                for shortname, type_str in zip(shortnames, types)])
+            cursor.execute("CREATE TABLE {} ({})".format(name, header))
             for i in range(getattr(self, self.parameters[0].shortname).shape[0]):
                 # Extra type cast is needed because sqlite3 cannot handle numpy int
                 # objects as INTEGER values.
-                values = tuple(x.dtype(getattr(self, x.shortname)[i]) for x in self.parameters)
-                cursor.execute("INSERT INTO {} VALUES ({})".format(name, value_subst), values)
+                values = []
+                for x in self.parameters:
+                    if x.dtype is linear_molecule_quantum_numbers:
+                        values += [i(q) for i, q in zip(linear_molecule_state_vars.values(),
+                                                        getattr(self, x.shortname)[i].values())]
+                    else:
+                        try:
+                            values.append(x.dtype(getattr(self, x.shortname)[i]))
+                        except IndexError:
+                            print(x.shortname)
+                            raise
+                cursor.execute("INSERT INTO {} VALUES ({})".format(name, value_subst),
+                               tuple(values))
             connection.commit()
 
     def download_from_web(self, lower_bound=0., upper_bound=10.e6):
@@ -93,9 +118,30 @@ def load_from_database(self, database):
         with connect(database) as connection:
             cursor = connection.cursor()
             name = scrub(self.molecule)
-            columns = ", ".join(["{}".format(x.shortname) for x in self.parameters])
+            shortnames = [x.shortname for x in self.parameters
+                          if x.shortname != "q1" and x.shortname != "q2"]
+            quantum_numbers = len(shortnames) != len(self.parameters)
+            if quantum_numbers:
+                #Quantum numbers needed.
+                all_parameters = self.parameters.copy()
+                self.parameters = [x for x in all_parameters if x.shortname in shortnames]
+            #Read all the non-quantum number parameters out of the database.
+            columns = ", ".join(shortnames)
             cursor.execute("SELECT {} from {}".format(columns, name))
             self.parse_records(cursor.fetchall())
+            if quantum_numbers:
+                #Now read in the quantum numbers and create the ordered dictionaries.
+                for q, state in zip(["q1", "q2"], ["upper", "lower"]):
+                    data = []
+                    shortnames = [x for x in linear_molecule_state_vars.keys()]
+                    types = [x for x in linear_molecule_state_vars.values()]
+                    columns = ", ".join(["{}_{}".format(x, state) for x in shortnames])
+                    cursor.execute("SELECT {} from {}".format(columns, name))
+                    for record in cursor.fetchall():
+                        data.append(OrderedDict((x.lower(), y(z)) for x, y, z in
+                                                zip(shortnames, types, record)))
+                    setattr(self, q, asarray(data))
+                self.parameters = all_parameters
 
     def parse_records(self, records):
         """Parses all database records and stores the data.
@@ -104,16 +150,25 @@ def parse_records(self, records):
             records: A list of dictionaries containing database record values.
         """
         for record in records:
-            try:
-                data = [x.dtype(y) for x, y in zip(self.parameters, record)]
-            except ValueError as e:
-                if "#" in str(e):
-                    warning("bad data value in database record:\n{}".format(record))
-                    continue
-                else:
-                    raise
+            data = []
+            for x, y in zip(self.parameters, record):
+                try:
+                    data.append(x.dtype(y))
+                except ValueError as e:
+                    if "#" in str(e):
+                        if x.shortname == "n_self":
+                            data.append(None)
+                        else:
+                            warning("bad data value in database record:\n{}".format(record))
+                            continue
+                    else:
+                        raise
+            if len(data) < len(self.parameters):
+                continue
             for x, datum in zip(self.parameters, data):
                 self.__dict__[x.shortname].append(datum)
+            if "n_self" in self.__dict__ and self.n_self[-1] is None:
+                self.n_self[-1] = self.n_air[-1]
         for x in self.parameters:
             setattr(self, x.shortname, asarray(getattr(self, x.shortname)))
         info("Found data for {} lines for {}.".format(getattr(self, self.parameters[0].shortname).shape[0],

pyrad/lbl/hitran/line_mixing.f90

@@ -0,0 +1,147 @@
+!Calculate the relaxation matrix.
+subroutine create_relaxation_matrix(nlines, temperature, iso, li, lf, gamma_hwhm, &
+                                    population, dk0, ji, jf, &
+                                    b0pp, b0pq, b0pr, &
+                                    b0qp, b0qq, b0qr, &
+                                    b0rp, b0rq, b0rr, &
+                                    w0pp, w0pq, w0pr, &
+                                    w0qp, w0qq, w0qr, &
+                                    w0rp, w0rq, w0rr, w)
+  use, intrinsic :: iso_fortran_env, only: int32, real64
+  implicit none
+
+  integer(kind=int32), intent(in) :: nlines
+  real(kind=real64), intent(in) :: temperature
+  integer(kind=int32), intent(in) :: iso, lf, li
+  integer(kind=int32), dimension(nlines), intent(in), target :: jf, ji
+  real(kind=real64), dimension(nlines), intent(in) :: dk0, gamma_hwhm, population
+  real(kind=real64), dimension(0:9,0:9,0:130,0:130), intent(in) :: b0pp, b0pq, b0pr, &
+                                                                   b0qp, b0qq, b0qr, &
+                                                                   b0rp, b0rq, b0rr, &
+                                                                   w0pp, w0pq, w0pr, &
+                                                                   w0qp, w0qq, w0qr, &
+                                                                   w0rp, w0rq, w0rr
+  real(kind=real64), dimension(nlines,nlines), intent(inout) :: w
+
+  integer(kind=int32) :: i, j, lf_, li_
+  integer(kind=int32), dimension(:), pointer :: jf_, ji_
+  logical :: iso_flag
+  real(kind=real64) :: b0, logt, sumlw, sumup, w0, ycal
+  real(kind=real64), parameter :: t0 = 296._real64
+
+  w(:,:) = 0._real64
+  logt = log(t0/temperature)
+  iso_flag = iso .gt. 2 .and. iso .ne. 7 .and. iso .ne. 10
+
+  !Define the transition so it is "downward".
+  li_ = min(li, lf)
+  lf_ = max(li, lf)
+  if (li .le. lf) then
+    ji_ => ji
+    jf_ => jf
+  else
+    ji_ => jf
+    jf_ => ji
+  endif
+
+  do i = 1, nlines
+    do j = 1, nlines
+      if (ji_(j) .gt. ji_(i)) cycle
+      if (iso_flag .and. mod(abs(ji(i) - ji(j)), 2) .ne. 0) cycle
+      if (ji_(i) .gt. jf_(i) .and. ji_(j) .gt. jf_(j)) then
+        w0 = w0pp(li_,lf_,ji_(i),ji_(j))
+        b0 = b0pp(li_,lf_,ji_(i),ji_(j))
+      elseif(ji_(i) .gt. jf_(i) .and. ji_(j) .eq. jf_(j))then
+        w0 = w0pq(li_,lf_,ji_(i),ji_(j))
+        b0 = b0pq(li_,lf_,ji_(i),ji_(j))
+      elseif(ji_(i) .gt. jf_(i) .and. ji_(j) .lt. jf_(j))then
+        w0 = w0pr(li_,lf_,ji_(i),ji_(j))
+        b0 = b0pr(li_,lf_,ji_(i),ji_(j))
+      elseif(ji_(i) .lt. jf_(i) .and. ji_(j) .gt. jf_(j))then
+        w0 = w0rp(li_,lf_,ji_(i),ji_(j))
+        b0 = b0rp(li_,lf_,ji_(i),ji_(j))
+      elseif(ji_(i) .lt. jf_(i) .and. ji_(j) .eq. jf_(j))then
+        w0 = w0rq(li_,lf_,ji_(i),ji_(j))
+        b0 = b0rq(li_,lf_,ji_(i),ji_(j))
+      elseif(ji_(i) .lt. jf_(i) .and. ji_(j) .lt. jf_(j))then
+        w0 = w0rr(li_,lf_,ji_(i),ji_(j))
+        b0 = b0rr(li_,lf_,ji_(i),ji_(j))
+      elseif(ji_(i) .eq. jf_(i) .and. ji_(j) .gt. jf_(j))then
+        w0 = w0qp(li_,lf_,ji_(i),ji_(j))
+        b0 = b0qp(li_,lf_,ji_(i),ji_(j))
+      elseif(ji_(i) .eq. jf_(i) .and. ji_(j) .eq. jf_(j))then
+        w0 = w0qq(li_,lf_,ji_(i),ji_(j))
+        b0 = b0qq(li_,lf_,ji_(i),ji_(j))
+      elseif(ji_(i) .eq. jf_(i) .and. ji_(j) .lt. jf_(j))then
+        w0 = w0qr(li_,lf_,ji_(i),ji_(j))
+        b0 = b0qr(li_,lf_,ji_(i),ji_(j))
+      endif
+      ycal = exp(w0 - b0*logt)
+      w(j,i) = ycal
+      w(i,j) = ycal*population(i)/population(j)
+    enddo
+  enddo
+
+  do i = 1, nlines
+    do j = 1, nlines
+      if (i .ne. j) then
+        w(i,j) = -abs(w(i,j))
+      endif
+    enddo
+  enddo
+
+  do i = 1, nlines
+    w(i,i) = gamma_hwhm(i)
+  enddo
+
+  do i = 1, nlines
+    sumlw = 0._real64
+    sumup = 0._real64
+    do j = 1, nlines
+      if (iso_flag .and. mod(abs(ji(i) - ji(j)), 2) .ne. 0) cycle
+      if (j .gt. i) then
+        sumlw = sumlw + abs(dk0(j))*w(j,i)
+      else
+        sumup = sumup + abs(dk0(j))*w(j,i)
+      endif
+    enddo
+    do j = i + 1, nlines
+      if (sumlw .eq. 0._real64) then
+        w(j,i) = 0._real64
+        w(i,j) = 0._real64
+      else
+        w(j,i) = w(j,i)*(-sumup/sumlw)
+        w(i,j) = w(j,i)*population(i)/population(j)
+      endif
+    enddo
+  enddo
+end subroutine create_relaxation_matrix
+
+
+!Calculate the first-order line-mixing coefficients.
+subroutine calculate_coefficients(nlines, iso, dipole, ji, v, w, y)
+  use, intrinsic :: iso_fortran_env, only: int32, real64
+  implicit none
+
+  integer(kind=int32), intent(in) :: iso, nlines
+  integer(kind=int32), dimension(nlines), intent(in) :: ji
+  real(kind=real64), dimension(nlines), intent(in) :: dipole, v
+  real(kind=real64), dimension(nlines,nlines), intent(in) :: w
+  real(kind=real64), dimension(nlines), intent(inout) :: y
+
+  logical :: iso_flag
+  integer(kind=int32) :: i, j
+  real(kind=real64) :: dv
+
+  iso_flag = iso .gt. 2 .and. iso .ne. 7 .and. iso .ne. 10
+  do i = 1, nlines
+    y(i) = 0._real64
+    do j = 1, nlines
+      if (j .eq. i) cycle
+      if (iso_flag .and. mod(abs(ji(i) - ji(j)), 2) .ne. 0) cycle
+      dv = v(i) - v(j)
+      if (abs(dv) .lt. 1.d-4) dv = 1.d-4
+      y(i) = y(i) + 2._real64*abs(dipole(j))/abs(dipole(i))*(1._real64/dv)*w(j,i)
+    enddo
+  enddo
+end subroutine calculate_coefficients