diff --git a/pyrad/lbl/hitran/database.py b/pyrad/lbl/hitran/database.py index b28d9cb..6bceafa 100644 --- a/pyrad/lbl/hitran/database.py +++ b/pyrad/lbl/hitran/database.py @@ -1,3 +1,4 @@ +from collections import OrderedDict from logging import getLogger from sqlite3 import connect from urllib.request import urlopen @@ -5,7 +6,8 @@ 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], diff --git a/pyrad/lbl/hitran/line_mixing.f90 b/pyrad/lbl/hitran/line_mixing.f90 new file mode 100644 index 0000000..a1a80dd --- /dev/null +++ b/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 diff --git a/pyrad/lbl/hitran/line_mixing.py b/pyrad/lbl/hitran/line_mixing.py new file mode 100644 index 0000000..e0d4cd8 --- /dev/null +++ b/pyrad/lbl/hitran/line_mixing.py @@ -0,0 +1,319 @@ +from logging import getLogger + +from netCDF4 import Dataset +from numpy import copy, exp, identity, int32, log, nonzero, ones, power, sqrt, where, zeros + +#from .line_mixing_fortran import calculate_coefficients, create_relaxation_matrix +from .line_parameters import c2, linear_molecule_state_vars, reference_temperature +from ..tips import TIPS_REFERENCE_TEMPERATURE +from .wigner3j import wigner_3j + + +info = getLogger(__name__).info + + +class LineMixing(object): + """Line mixing container. Adapted from doi: 10.1016/j.jqsrt.2004.11.011. + + Attributes: + b0pp: Fitted relaxation matrix temperature dependence powers for p-p transitions. + b0pq: Fitted relaxation matrix temperature dependence powers for p-q transitions. + b0pr: Fitted relaxation matrix temperature dependence powers for p-r transitions. + b0qp: Fitted relaxation matrix temperature dependence powers for q-p transitions. + b0qq: Fitted relaxation matrix temperature dependence powers for q-q transitions. + b0qr: Fitted relaxation matrix temperature dependence powers for q-r transitions. + b0rp: Fitted relaxation matrix temperature dependence powers for r-p transitions. + b0rq: Fitted relaxation matrix temperature dependence powers for r-q transitions. + b0rr: Fitted relaxation matrix temperature dependence powers for r-r transitions. + w0pp: Fitted relaxation matrix values for p-p transitions. + w0pq: Fitted relaxation matrix values for p-q transitions. + w0pr: Fitted relaxation matrix values for p-r transitions. + w0qp: Fitted relaxation matrix values for q-p transitions. + w0qq: Fitted relaxation matrix values for q-q transitions. + w0qr: Fitted relaxation matrix values for q-r transitions. + w0rp: Fitted relaxation matrix values for r-p transitions. + w0rq: Fitted relaxation matrix values for r-q transitions. + w0rr: Fitted relaxation matrix values for r-r transitions. + """ + def __init__(self, path): + """Reads the (w0, b0) fitted parameters taken from doi: 10.1016/j.jqsrt.2014.09.017. + + Args: + path - Path to netcdf file. + """ + info("Reading line-mixing relaxation matrix data from {}.".format(path)) + with Dataset(path, "r") as dataset: + for branch in ["pp", "pq", "pr", "qp", "qq", "qr", "rp", "rq", "rr"]: + setattr(self, "w0{}".format(branch), + copy(dataset.variables["w0{}".format(branch)], order="F")) + setattr(self, "b0{}".format(branch), + copy(dataset.variables["b0{}".format(branch)], order="F")) + + @staticmethod + def calculate_coefficients(relaxation_matrix, dipole, line_center, I_off, mask): + """Calculates first-order line-mixing coefficients using equation 6 of + doi: 10.1016/j.jqsrt.2004.11.011. + + Args: + relaxation_matrix: Relaxation matrix. + line_center: Line center wavenumber [cm-1]. + I_off: 2D matrix with 0s on diagnal and 1s off-diagonal. + mask: Matrix used to mask out certain transitions. + + Returns: + First-order line-mixing coefficients. + """ + y = zeros(line_center.size) + dipole_ratio = dipole/dipole.reshape((dipole.size, 1)) + x = line_center.reshape((line_center.size, 1)) - line_center + dw = where(abs(x) < 1.e-4, 1.e-4, x) + x = mask*I_off*dipole_ratio/dw + for i in range(line_center.size): + y[i] = 2.*sum(x[i, :]*relaxation_matrix[:, i]) + return y + + def create_relaxation_matrix(self, l2_i, l2_f, j_i, j_f, population, temperature, + gamma, dk0, mask, I, I_off): + """Calculates the relaxation matrix as in doi: 10.1016/j.jqsrt.2014.09.017. + + Args: + l2_i: Initial state l2 angular momentum. + l2_f: Final state l2 angular momentum. + j_i: Array of initial state rotational quantum numbers. + j_f: Array of final state rotational quantum_numbers. + population: Array of lower level populations. + temperature: Temperature [K]. + gamma: Array of pressure-broadened halfwidth coefficients [cm-1 atm-1]. + dk0: Rigid rotor dipole matrix elements. + mask: Matrix used to mask out specific transitions. + I: Identity matrix. + I_off: Matrix with zeros on-diagonal and ones off-diagonal. + + Returns: + Relaxation matrix. + """ + tlog = log(reference_temperature/temperature) + w = zeros((j_i.size, j_i.size)) + if l2_i <= l2_f: + a, b, j1, j2 = l2_i, l2_f, j_i[:], j_f[:] + else: + a, b, j1, j2 = l2_f, l2_i, j_f[:], j_i[:] + #xmask = where((j1.reshape(j1.size, 1) > j1), 0., 1.)*mask + xmask = where((j1 > j1.reshape(j1.size, 1)), 0., 1.)*mask + for i in range(j_i.size): + branch1 = branch(j1[i], j2[i]) + w0, b0 = zeros(j_i.size), zeros(j_i.size) + indices = nonzero(xmask[i, :])[0] + for j in indices: + branch2 = branch(j1[j], j2[j]) + w0[j] = getattr(self, "w0{}{}".format(branch1, branch2))[a, b, j1[i], j1[j]] + b0[j] = getattr(self, "b0{}{}".format(branch1, branch2))[a, b, j1[i], j1[j]] + ycal = exp(w0[indices] - b0[indices]*tlog) + w[indices, i] = ycal[:] + w[i, indices] = ycal[:]*population[i]/population[indices] + + # Set diagonal terms to the line-broadening parameters. + w = -1.*abs(w)*I_off + gamma*I + + for i in range(gamma.size): + sum_lw = sum(mask[i, i+1:]*abs(dk0[i+1:])*w[i+1:, i]) + sum_up = sum(mask[i, :i+1]*abs(dk0[:i+1])*w[:i+1, i]) + if sum_lw == 0.: + w[i+1:, i] = 0. + w[i, i+1:] = 0. + else: + w[i+1:, i] = -1.*w[i+1:, i]*sum_up/sum_lw + w[i, i+1:] = w[i+1:, i]*population[i]/population[i+1:] + return w + + def first_order_coefficients(self, spectral_lines, temperature, pressure, mixing_ratio): + """Calculates 1st order line-mixing coefficients. + + Args: + spectral_lines: SpectralLines object. + temperature: Temperature [K]. + pressure: Pressure [atm]. + mixing_ratio: Mixing ratio. + """ + bands = create_bands(spectral_lines) + y = zeros(spectral_lines.v.size) + info("Calculating line-mixing coefficents ({} bands)".format(len(bands.keys()))) + for i, (nums, indices) in enumerate(bands.items()): + #Determine the initial and final l values. + li = int32(spectral_lines.q2[indices[0]]["l2"]) + lf = int32(spectral_lines.q1[indices[0]]["l2"]) + if li > 8 or lf > 8 or abs(li - lf) > 1: + continue + + #Create array of initial and final j values. + ji, jf = zeros(len(indices), dtype=int), zeros(len(indices), dtype=int32) + for j, (qns2, qns1) in enumerate(zip(spectral_lines.q2[indices], + spectral_lines.q1[indices])): + ji[j] = qns2["j"] + jf[j] = qns1["j"] + + #Adjust the lower level populations for temperature. + n = temperature_adjust_population(spectral_lines.iso[indices[0]], + spectral_lines.population[indices], + spectral_lines.en[indices], + temperature, + spectral_lines.q.total_partition_function) + + #Adjust the doppler halfwidth for temperature. + gamma = temperature_adjust_halfwidth(spectral_lines.n_air[indices], + spectral_lines.n_self[indices], + spectral_lines.gamma_air[indices], + spectral_lines.gamma_self[indices], + mixing_ratio, + temperature) + + #Sort parameters in decreasing order of intensity. + s = spectral_lines.v[indices]*n[:]*power(spectral_lines.dipole[indices], 2) + sorted_indices = s.argsort() + ji_s = ji[sorted_indices[::-1]] + jf_s = jf[sorted_indices[::-1]] + s_s = s[sorted_indices[::-1]] + v_s = spectral_lines.v[indices][sorted_indices[::-1]] + gamma_s = gamma[sorted_indices[::-1]] + dair_s = spectral_lines.d_air[indices][sorted_indices[::-1]] + n_s = n[sorted_indices[::-1]] + dk0_s = spectral_lines.dk0[indices][sorted_indices[::-1]] + dipole_s = spectral_lines.dipole[indices][sorted_indices[::-1]] + + #Create some useful constants. + I = identity(gamma_s.size) + unity = ones((gamma_s.size, gamma_s.size)) + I_off = unity - I + x = mask(spectral_lines.iso[indices[0]], ji_s) + + #Create relaxation matrix. + w = self.create_relaxation_matrix(li, lf, ji_s, jf_s, n_s, temperature, + gamma_s, dk0_s, x, I, I_off) +# w = zeros((ji_s.size, ji_s.size), order="F") +# create_relaxation_matrix(ji_s.size, temperature, spectral_lines.iso[indices[0]], +# li, lf, gamma_s, n_s, dk0_s, ji_s, jf_s, +# self.b0pp, self.b0pq, self.b0pr, +# self.b0qp, self.b0qq, self.b0qr, +# self.b0rp, self.b0rq, self.b0rr, +# self.w0pp, self.w0pq, self.w0pr, +# self.w0qp, self.w0qq, self.w0qr, +# self.w0rp, self.w0rq, self.w0rr, w) + + #Calculate first order line-mixing coefficients. + yt_s = self.calculate_coefficients(w, dipole_s, v_s, I_off, x) + yt = zeros(ji_s.size) +# yt_s, yt = zeros(ji_s.size), zeros(ji_s.size) +# calculate_coefficients(ji_s.size, spectral_lines.iso[indices[0]], +# dipole_s, ji_s, v_s, w, yt_s) + yt[sorted_indices[::-1]] = yt_s[:] + y[indices] = yt[:] + return y + + +def temperature_adjust_population(iso, population, en, temperature, partition_function): + """Adjusts level populations to input temperature. + + Args: + iso: Isotopologue id. + population: Upper level population. + en: Lower state enegery [cm-1]. + temperature: Temperature [K]. + partition_function: Total partition function routine. + + Returns: + Temperature-corrected population. + """ + return population*partition_function(TIPS_REFERENCE_TEMPERATURE, iso) * \ + exp(c2*en*(1./temperature - 1./TIPS_REFERENCE_TEMPERATURE)) / \ + partition_function(temperature, iso) + + +def temperature_adjust_halfwidth(n_air, n_self, gamma_air, gamma_self, mixing_ratio, + temperature): + """Adjusts input halfwidths to input temperature. + + Args: + n_air: Air-broadened halfwidth temperature dependence power. + n_self: Self-broadened halfwidth temperature dependence power. + gamma_air: Air-broadened halfwidth at reference temperature [atm-1 cm-1]. + gamma_self: Self-broadened halfwidth at reference temperature [atm-1 cm-1]. + mixing_ratio: Molecular mixing ratio. + temperature: Temperature [K]. + + Returns: + Temperature-corrected halfwidth [atm-1 cm-1]. + """ + return (1. - mixing_ratio)*gamma_air*power(reference_temperature/temperature, n_air) + \ + mixing_ratio*gamma_self*power(reference_temperature/temperature, n_self) + + +def create_bands(spectral_lines): + """Groups the lines in the same vibrational bands. + + Args: + spectral_lines: SpectralLines object. + + Returns: + Dictionary of lists of indices. + """ + bands = {} + for i, (upper, lower) in enumerate(zip(spectral_lines.q1, spectral_lines.q2)): + key = tuple([spectral_lines.iso[i]] + + [upper[x.lower()] for x in linear_molecule_state_vars.keys() if x != "J"] + + [lower[x.lower()] for x in linear_molecule_state_vars.keys() if x != "J"]) + if key in bands: + bands[key].append(i) + else: + bands[key] = [i] + return bands + + +def branch(ji, jf): + """Returns the branch designation string for the transition. + + Args: + ji: Initial state rotational quantum number. + jf: Final state rotational quantum number. + + Returns: + Character describing transition. + """ + if ji > jf: + return "p" + elif ji == jf: + return "q" + else: + return "r" + + +def mask(iso, ji): + """Creates a mask for transitions. + + Args: + iso: Isotopologue id. + ji: Initial state rotational quantum number. + + Returns: + Array of zeros and ones used to mask out transitions. + """ + if iso > 2 and iso != 7 and iso != 10: + return where((ji.reshape((ji.size, 1)) - ji) % 2 != 0, 0., 1.) + else: + return ones((ji.size, ji.size)) + + +def rigid_rotor_dipole_matrix_element(ji, jf, l2i, l2f): + """Computes a rigid rotor dipole matrix element from equation 3 in + doi: 10.1016/j.jqsrt.2004.11.011 + + Args: + ji: Initial state rotational quantum number. + jf: Final state rotational quantum number. + l2i: Initial state l2 angular momentum value. + l2f: Final state l2 angular momentum value. + + Returns: + Rigid rotor dipole matrix element. + """ + return power(-1., jf + l2f + 1)*sqrt(2*jf + 1) * \ + wigner_3j(ji, 1, jf, l2i, l2f - l2i, -l2f) diff --git a/pyrad/lbl/hitran/line_mixing_fortran.pyf b/pyrad/lbl/hitran/line_mixing_fortran.pyf new file mode 100644 index 0000000..a6b2535 --- /dev/null +++ b/pyrad/lbl/hitran/line_mixing_fortran.pyf @@ -0,0 +1,50 @@ +! -*- f90 -*- +! Note: the context of this file is case sensitive. + +python module line_mixing_fortran ! in + interface ! in :line_mixing_fortran + 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) ! in :line_mixing_fortran:line_mixing.f90 + integer(kind=4) intent(in) :: nlines + real(kind=8) intent(in) :: temperature + integer(kind=4) intent(in) :: iso + integer(kind=4) intent(in) :: li + integer(kind=4) intent(in) :: lf + real(kind=8) dimension(nlines),intent(in),depend(nlines) :: gamma_hwhm + real(kind=8) dimension(nlines),intent(in),depend(nlines) :: population + real(kind=8) dimension(nlines),intent(in),depend(nlines) :: dk0 + integer(kind=4), target,dimension(nlines),intent(in),depend(nlines) :: ji + integer(kind=4), target,dimension(nlines),intent(in),depend(nlines) :: jf + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: b0pp + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: b0pq + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: b0pr + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: b0qp + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: b0qq + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: b0qr + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: b0rp + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: b0rq + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: b0rr + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: w0pp + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: w0pq + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: w0pr + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: w0qp + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: w0qq + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: w0qr + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: w0rp + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: w0rq + real(kind=8) dimension(0:9,0:9,0:130,0:130),intent(in) :: w0rr + real(kind=8) dimension(nlines,nlines),intent(inout),depend(nlines,nlines) :: w + end subroutine create_relaxation_matrix + subroutine calculate_coefficients(nlines,iso,dipole,ji,v,w,y) ! in :line_mixing_fortran:line_mixing.f90 + integer(kind=4) intent(in) :: nlines + integer(kind=4) intent(in) :: iso + real(kind=8) dimension(nlines),intent(in),depend(nlines) :: dipole + integer(kind=4) dimension(nlines),intent(in),depend(nlines) :: ji + real(kind=8) dimension(nlines),intent(in),depend(nlines) :: v + real(kind=8) dimension(nlines,nlines),intent(in),depend(nlines,nlines) :: w + real(kind=8) dimension(nlines),intent(inout),depend(nlines) :: y + end subroutine calculate_coefficients + end interface +end python module line_mixing_fortran + +! This file was auto-generated with f2py (version:2). +! See http://cens.ioc.ee/projects/f2py2e/ diff --git a/pyrad/lbl/hitran/line_parameters.py b/pyrad/lbl/hitran/line_parameters.py index d355fde..096f42d 100644 --- a/pyrad/lbl/hitran/line_parameters.py +++ b/pyrad/lbl/hitran/line_parameters.py @@ -4,13 +4,37 @@ from numpy import float64 +c2 = -1.4387768795689562 # (hc/k) [K cm]. +reference_temperature = 296. # [K]. +linear_molecule_state_vars = OrderedDict([("ElecStateLabel", str), ("v1", int), ("v2", int), + ("l2", int), ("v3", int), ("J", int), ("r", int), + ("kronigParity", str)]) + + +class LinearMoleculeQuantumNumbersError(BaseException): + """Signals an error while parsing linear molecule quantum numbers.""" + pass + + def linear_molecule_quantum_numbers(value): - m = match(r"ElecStateLabel=X;v1=([0-9]+);v2=([0-9]+);l2=([0-9]+);v3=([0-9]+);J=([0-9]+);", - value) + """Parses quantum numbers from HITRAN http responses. + + Args: + value: String containing Hitran state information. + + Returns: + An ordered dictionary of quantum numbers. + """ + #Build the regular expression. + capture = [] + for name, type in linear_molecule_state_vars.items(): + group = "([0-9]+)" if type is int else "([A-Za-z]+)" + capture.append("{}={}".format(name, group)) + m = match(r"{}".format(".*?;".join(capture)), value) if not m: - raise ValueError("invalid quantum numbers in {}".format(value)) - return OrderedDict((x, int(y)) for x, y in zip(["v1", "v2", "l2", "v3", "j"], - m.groups())) + raise LinearMoleculeQuantumNumbersError("invalid quantum numbers in {}".format(value)) + return OrderedDict((x.lower(), y(z)) for (x, y), z in zip(linear_molecule_state_vars.items(), + m.groups())) HitranParameter = namedtuple("HitranParameter", ["api_name", "longname", "shortname", diff --git a/pyrad/lbl/hitran/lorentz.py b/pyrad/lbl/hitran/lorentz.py index 1f7485b..49a9c46 100644 --- a/pyrad/lbl/hitran/lorentz.py +++ b/pyrad/lbl/hitran/lorentz.py @@ -10,7 +10,7 @@ class Lorentz(object): """ def __init__(self): - self.parameters = ["gamma_air", "gamma_self", "n_air"] + self.parameters = ["gamma_air", "gamma_self", "n_air", "n_self"] self.halfwidth = None def update(self, spectral_lines, temperature, pressure, partial_pressure): @@ -23,7 +23,8 @@ def update(self, spectral_lines, temperature, pressure, partial_pressure): partial_pressure: Partial pressure [atm]. """ self.halfwidth = pressure_broadened_halfwidth(pressure, partial_pressure, temperature, - spectral_lines.n_air, spectral_lines.gamma_air, + spectral_lines.n_air, spectral_lines.n_self, + spectral_lines.gamma_air, spectral_lines.gamma_self) def profile(self, spectral_lines, v, index): @@ -41,22 +42,23 @@ def profile(self, spectral_lines, v, index): def pressure_broadened_halfwidth(pressure, partial_pressure, temperature, - n, gamma_air, gamma_self): + n_air, n_self, gamma_air, gamma_self): """Calculates pressure-broadened line halfwidth. Args: pressure: Pressure [atm]. partial_pressure: Partial pressure [atm]. temperature: Temperature [K] - n: Air-broadened temperature dependence powers. + n_air: Air-broadened temperature dependence powers. + n_self: Self-broadened temperature dependence powers. gamma_air: Air-broadened halfwidth [cm-1 atm-1]. gamma_self: Self-broadened halfwidth [cm-1 atm-1]. Returns: Pressure-broadened line halfwidth [cm-1]. """ - return power((296./temperature), n) * \ - (gamma_air*(pressure - partial_pressure) + gamma_self*partial_pressure) + return power((296./temperature), n_air)*(gamma_air*(pressure - partial_pressure)) + \ + power((296./temperature), n_self)*gamma_self*partial_pressure def lorentz_profile(dv, halfwidth): diff --git a/pyrad/lbl/hitran/spectral_lines.py b/pyrad/lbl/hitran/spectral_lines.py index d5dd37c..403efc2 100644 --- a/pyrad/lbl/hitran/spectral_lines.py +++ b/pyrad/lbl/hitran/spectral_lines.py @@ -2,25 +2,27 @@ from numpy import asarray, copy, exp, searchsorted, sqrt, zeros +from .line_mixing import rigid_rotor_dipole_matrix_element +from .line_parameters import c2, reference_temperature from ..tips import TIPS_REFERENCE_TEMPERATURE -c2 = -1.4387768795689562 # (hc/k) [K cm]. -reference_temperature = 296. # [K]. - - class SpectralLines(object): """Hitran spectral line parameters for a molecule. Attributes: delta: Numpy array of air-broadended pressure shifts [cm-1 atm-1] (lines). + dipole: Numpy array of transition moments [cm] (lines). + dk0: Numpy array of rigid rotor dipole matrix elements (lines). en: Numpy array of transition lower state energies [cm-1] (lines). + g: Numpy array of lower state statistical weights (lines). gamma_air: Numpy array of air-broadened halfwidths [cm-1 atm-1] (lines). gamma_self: Numpy array of self-broadened halfwidths [cm-1 atm-1] (lines). id: HITRAN molecule id. iso: Numpy array of HITRAN isotopologue ids (lines). mass: Numpy array of isotopologue masses [g] (lines). n: Numpy array of air-broadened temperature dependence powers (lines). + population: Numpy array of lower state populations (lines). s: Numpy array of line strengths [cm] (lines). q: TotalPartitionFunction object. v: Numpy array of transition wavenumbers [cm-1] (lines). @@ -47,7 +49,6 @@ def __init__(self, database, total_partition_function): # Get the mass of the isotopologues. self.mass = asarray([float(database.isotopologues[x-1].mass) for x in self.iso]) - self.line_profile = database.line_profile self.q = total_partition_function @@ -55,6 +56,17 @@ def __init__(self, database, total_partition_function): self.s[:] *= self.temperature_correct_line_strength(self.q, TIPS_REFERENCE_TEMPERATURE, self.iso, self.en, self.v) + if self.line_profile.line_mixing is not None: + self.population = self.g[:]*exp(c2*self.en[:]/reference_temperature) / \ + self.q.total_partition_function(TIPS_REFERENCE_TEMPERATURE, self.iso[:]) + self.dipole = sqrt(self.s[:]/(self.g[:]*self.v[:])) + self.dk0 = zeros(self.q2.size) + for i in range(self.q2.size): + self.dk0[i] = rigid_rotor_dipole_matrix_element(self.q2[i]["j"], + self.q1[i]["j"], + self.q2[i]["l2"], + self.q1[i]["l2"]) + def absorption_coefficient(self, temperature, pressure, partial_pressure, wavenumber, cut_off=25.): """Calculates the absorption coefficient. @@ -70,10 +82,10 @@ def absorption_coefficient(self, temperature, pressure, partial_pressure, wavenu Numpy array of absorption coefficients [cm2] (wavenumber). """ lines = shallow_copy(self) - lines.s = self.correct_line_strengths(temperature) - lines.v = lines.pressure_shift_transition_wavenumbers(pressure) profile = shallow_copy(lines.line_profile) profile.update(lines, temperature, pressure, partial_pressure) + lines.s = lines.correct_line_strengths(temperature) + lines.v = lines.pressure_shift_transition_wavenumbers(pressure) k = zeros(wavenumber.size) for i in range(lines.s.size): left = searchsorted(wavenumber, lines.v[i] - cut_off, side="left") diff --git a/pyrad/lbl/hitran/voigt.py b/pyrad/lbl/hitran/voigt.py index 747ae8f..923891d 100644 --- a/pyrad/lbl/hitran/voigt.py +++ b/pyrad/lbl/hitran/voigt.py @@ -2,6 +2,7 @@ from scipy.special import wofz from .doppler import doppler_broadened_halfwidth +from .line_mixing import LineMixing from .lorentz import pressure_broadened_halfwidth @@ -10,14 +11,22 @@ class Voigt(object): Attributes: doppler_halfwidth: Doppler-broadened halfwidth [cm-1]. + line_mixing: LineMixing object. parameters: List of HITRAN parameter names. pressure_halfwidth: Pressure-broadened halfwidth [cm-1]. """ - def __init__(self): - self.parameters = ["center", "gamma_air", "gamma_self", "n_air"] + def __init__(self, relaxation_matrix_path=None): + self.parameters = ["center", "gamma_air", "gamma_self", "n_air", "n_self"] self.doppler_halfwidth = None self.pressure_halfwidth = None + if relaxation_matrix_path is not None: + self.parameters += ["upper_level_quantum_numbers", + "lower_level_quantum_numbers", + "lower_level_degeneracy"] + self.line_mixing = LineMixing(relaxation_matrix_path) + else: + self.line_mixing = None def update(self, spectral_lines, temperature, pressure, partial_pressure): """Calculate per-spectral-line pressure-broadened halfwidths. @@ -28,10 +37,18 @@ def update(self, spectral_lines, temperature, pressure, partial_pressure): pressure: Pressure [atm]. partial_pressure: Partial pressure [atm]. """ + if self.line_mixing is not None: + self.y = self.line_mixing.first_order_coefficients(spectral_lines, + temperature, + pressure, + partial_pressure/pressure) self.doppler_halfwidth = doppler_broadened_halfwidth(temperature, spectral_lines.mass, spectral_lines.v) - self.pressure_halfwidth = pressure_broadened_halfwidth(pressure, partial_pressure, temperature, - spectral_lines.n_air, spectral_lines.gamma_air, + self.pressure_halfwidth = pressure_broadened_halfwidth(pressure, partial_pressure, + temperature, + spectral_lines.n_air, + spectral_lines.n_self, + spectral_lines.gamma_air, spectral_lines.gamma_self) def profile(self, spectral_lines, v, index): @@ -45,8 +62,12 @@ def profile(self, spectral_lines, v, index): Returns: Line broadening [cm]. """ - return voigt_profile(v - spectral_lines.v[index], self.pressure_halfwidth[index], - self.doppler_halfwidth[index]) + x = voigt_profile(v - spectral_lines.v[index], self.pressure_halfwidth[index], + self.doppler_halfwidth[index]) + if self.line_mixing is None: + return x.real + else: + return ((1. - 1j*self.y[index])*x).real def voigt_profile(dv, pressure_halfwidth, doppler_halfwidth): @@ -61,4 +82,4 @@ def voigt_profile(dv, pressure_halfwidth, doppler_halfwidth): Voigt line profile broadening [cm]. """ sigma = doppler_halfwidth/sqrt(2*log(2)) - return wofz((dv + 1j*pressure_halfwidth)/sigma/sqrt(2)).real/sigma/sqrt(2*pi) + return wofz((dv + 1j*pressure_halfwidth)/sigma/sqrt(2))/sigma/sqrt(2*pi) diff --git a/pyrad/lbl/hitran/wigner3j.py b/pyrad/lbl/hitran/wigner3j.py new file mode 100644 index 0000000..6591760 --- /dev/null +++ b/pyrad/lbl/hitran/wigner3j.py @@ -0,0 +1,225 @@ +from numpy import power, sqrt +from scipy.special import factorial + + +def wigner_3j(j1, j2, j, m1, m2, m): + """Computes Wigner 3 j-symbols using + https://mathworld.wolfram.com/Wigner3j-Symbol.html + + Args: + j1: Top left 3 j-symbol value. + j2: Top center 3 j-symbol value. + j: Top right 3 j-symbol value. + m1: Bottom left 3 j-symbol value. + m2: Bottom center 3 j-symbol value. + m: Bottom right 3 j-symbol value. + + Returns: + Wigner 3 j-symbol value. + """ + if not _selection_rules(j1, j2, j, m1, m2, m): return 0. + try: + if j1 == j2 and m1 == -1*m2 and j == 0 and m == 0: + return _special_case1(j1, m1) + elif m1 == 0 and m2 == 0 and m == 0: + return _special_case2(j1, j2, j) + elif j == j1 + j2: + return _special_case3(j1, j2, m1, m2, m) + elif j1 == j2 + j: + if not _selection_rules(j2, j, j1, m2, m, m1): return 0. + return _special_case3(j2, j, m2, m, m1) + elif m1 == j1 and m2 == -j1: + return _special_case4(j1, j2, j, m) + elif m1 == j1 and m == -j1: + if not _selection_rules(j1, j, j2, m1, m, m2): return 0. + return power(-1., j1 + j2 + j)*_special_case4(j1, j, j2, m2) + else: + return _general_case(j1, j2, j, m1, m2, m) + except OverflowError: + print("{} {} {}\n{} {} {}".format(j1, j2, j, m1, m2, m)) + raise + + +def _selection_rules(j1, j2, j, m1, m2, m): + """Checks selection rules. + + Args: + j1: Top left 3 j-symbol value. + j2: Top center 3 j-symbol value. + j: Top right 3 j-symbol value. + m1: Bottom left 3 j-symbol value. + m2: Bottom center 3 j-symbol value. + m: Bottom right 3 j-symbol value. + + Returns: + True if values pass selection rules check, else False. + """ + if not (m1 + m2 == -m or (abs(j1 - j2) <= j and j <= j1 + j2)): + return False + for x, y in zip([j1, j2, j], [m1, m2, -m]): + if not (y >= -1*abs(x) and y <= abs(x)): + return False + return True + + +def _special_case1(l, m): + """Calculates wigner 3-j symbol of the form: + |l l 0| + |m -m 0| + + Args: + l: Top left 3 j-symbol value. + m: Bottom left 3 j-symbol value. + + Returns: + Wigner 3 j-symbol value. + """ + return power(-1., l - m)/sqrt(2*l + 1) + + +def _special_case2(j1, j2, j): + """Calculates wigner 3-j symbol of the form: + |j1 j2 j| + |0 0 0| + + Args: + j1: Top left 3 j-symbol value. + j2: Top center 3 j-symbol value. + j: Top right 3 j-symbol value. + + Returns: + Wigner 3 j-symbol value. + """ + J = j1 + j2 + j + if J % 2 == 0: + g = J/2 + term1 = power(-1., g) + term2 = sqrt((factorial(2*g - 2*j1, exact=True) * + factorial(2*g - 2*j2, exact=True) * + factorial(2*g - 2*j, exact=True)) / + factorial(2*g + 1, exact=True)) + term3 = factorial(g, exact=True) / \ + (factorial(g - j1, exact=True) * + factorial(g - j2, exact=True) * + factorial(g - j, exact=True)) + return term1*term2*term3 + else: + return 0. + + +def _special_case3(j1, j2, m1, m2, m): + """Calculates wigner 3-j symbol of the form: + |j1 j2 j1+j2| + |m1 m2 m | + + Args: + j1: Top left 3 j-symbol value. + j2: Top center 3 j-symbol value. + m1: Bottom left 3 j-symbol value. + m2: Bottom center 3 j-symbol value. + m: Bottom right 3 j-symbol value. + + Returns: + Wigner 3 j-symbol value. + """ + term1 = power(-1., j1 - j2 - m) + term2 = sqrt((factorial(2*j1, exact=True)*factorial(2*j2, exact=True)) / + factorial(2*j1 + 2*j2 + 1, exact=True)) + term3 = sqrt((factorial(j1 + j2 - m, exact=True) * + factorial(j1 + j2 + m, exact=True)) / + (factorial(j1 + m1, exact=True) * + factorial(j1 - m1, exact=True) * + factorial(j2 + m2, exact=True) * + factorial(j2 - m2, exact=True))) + return term1*term2*term3 + + +def _special_case4(j1, j2, j, m): + """Calculates wigner 3-j symbol of the form: + |j1 j2 j| + |j1 -j1 m| + + Args: + j1: Top left 3 j-symbol value. + j2: Top center 3 j-symbol value. + j: Top right 3 j-symbol value. + m: Bottom right 3 j-symbol value. + + Returns: + Wigner 3 j-symbol value. + """ + term1 = power(-1., -j1 + j2 - m) + term2 = sqrt((factorial(2*j1, exact=True) * + factorial(-j1 + j2 + j, exact=True)) / + (factorial(j1 + j2 + j + 1, exact=True) * + factorial(j1 - j2 + j, exact=True))) + term3 = sqrt((factorial(j1 + j2 - m, exact=True) * + factorial(j + m, exact=True)) / + (factorial(j1 + j2 - j, exact=True) * + factorial(-j1 + j2 + m, exact=True) * + factorial(j - m, exact=True))) + return term1*term2*term3 + + +def _general_case(a, b, c, alpha, beta, gamma): + """Calculates wigner 3-j symbol for the general case: + | a b c | + |alpha beta gamma| + + Args: + a: Top left 3 j-symbol value. + b: Top center 3 j-symbol value. + c: Top right 3 j-symbol value. + alpha: Bottom left 3 j-symbol value. + beta: Bottom center 3 j-symbol value. + gamma: Bottom right 3 j-symbol value. + + Returns: + Wigner 3 j-symbol value. + """ + term1 = power(-1., a - b - gamma) + numerator = sorted([x for x in [a + b - c, a - b + c, -a + b + c, a + alpha, + a - alpha, b + beta, b - beta, c + gamma, + c - gamma] if x > 1]) + denominator1 = [a + b + c + 1] + s = 0. + start = max(-c + a + beta, -c + b - alpha, 0) + stop = min(b + beta, a - alpha, a + b - c) + for t in range(start, stop + 1): + d = [x for x in [t, c - b + t + alpha, c - a + t - beta, + a + b - c - t, a - t - alpha, + b - t + beta] if x >= 0] + if len(d) < 6: continue + denominator = sorted([x for x in denominator1 + d + d if x > 1]) + top, bottom = [], [] + for x, y in zip(numerator, denominator): + if x == y: continue + if x > y: + top.append((y + 1, x)) + else: + bottom.append((x + 1, y)) + if len(numerator) > len(denominator): + top += [(1, x) for x in numerator[len(denominator):]] + elif len(numerator) < len(denominator): + bottom += [(1, x) for x in denominator[len(numerator):]] + f1 = 1 + for x in top: + f1 *= _fact(x[0], x[1]) + f2 = 1 + for x in bottom: + f2 *= _fact(x[0], x[1]) + s += power(-1., t)*sqrt(f1/f2) + return term1*s + + +def _fact(start, stop): + """Calculates the product of all integers in range start:stop (inclusive). + + Args: + start: Range lower bound. + stop: Range upper bound. + """ + x = start + for i in range(start + 1, stop + 1): + x *= i + return x diff --git a/pyrad/utils/database_utilities.py b/pyrad/utils/database_utilities.py index 2f3352a..06d8d72 100644 --- a/pyrad/utils/database_utilities.py +++ b/pyrad/utils/database_utilities.py @@ -3,8 +3,7 @@ from numpy import float64 -SQL_TYPES = {float64: "REAL", - int: "INTEGER"} +SQL_TYPES = {float64: "REAL", int: "INTEGER", str: "TEXT"} def ascii_table_records(response, block_size=512): diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..a86fb44 --- /dev/null +++ b/setup.py @@ -0,0 +1,23 @@ +from setuptools import find_packages, setup + + +with open("README.md", "r", encoding="utf-8") as readme: + long_description = readme.read() + +setup( + name="pyrad", + version="0.0.1", + author="R. Menzel", + author_email="author@example.com", + description="A is a simple (one-dimensional) all-sky atmospheric radiation package.", + long_description=long_description, + long_description_content_type="text/markdown", + url="https://github.com/menzel-gfdl/pylbl", + packages=find_packages(), + classifiers=[ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: LGPL-3.0 License", + "Operating System :: OS Independent", + ], + install_requires=["netCDF4", "numpy", "scipy"], + python_requires=">=3.5") diff --git a/tests/test_wigner3j.py b/tests/test_wigner3j.py new file mode 100644 index 0000000..6df9aa9 --- /dev/null +++ b/tests/test_wigner3j.py @@ -0,0 +1,36 @@ +from unittest import main, TestCase + +from numpy import seterr, sqrt + +from pyrad.lbl.hitran.wigner3j import wigner_3j + + +class TestWigner3j(TestCase): + + def test_wigner_3j(self): + inputs = [((24, 1, 23, 0, 1, -1), (1/7)*sqrt(23/94)), # Special case 3. + ((4, 1, 3, 1, 0, -1), -0.5*sqrt(5/21)), + ((18, 1, 19, 2, 0, -2), -sqrt(119/9139)), + ((23, 1, 22, 3, -1, -2), (1/3)*sqrt(65/1081)), + ((39, 1, 38, 3, 1, -4), sqrt(30/11297)), + ((1, 1, 1, 1, 0, -1), -sqrt(1/6)), # Special case 4. + ((1, 1, 1, 1, -1, 0), sqrt(1/6)), + ((2, 1, 2, 2, 0, -2), -sqrt(2/15)), + ((3, 1, 3, 3, 0, -3), -0.5*sqrt(3/7)), + ((4, 1, 4, 4, 0, -4), (-2/3)*sqrt(0.2)), + ((5, 1, 5, 5, 0, -5), -sqrt(5/66)), + ((2, 1, 2, 0, 1, -1), -sqrt(0.1)), # General case. + ((10, 1, 10, 0, 1, -1), -sqrt(1/42)), + ((24, 1, 24, 0, 1, -1), -(1/7)*sqrt(0.5)), + ((2, 1, 2, 1, 0, -1), sqrt(1/30)), + ((22, 1, 22, 1, 0, -1), (1/3)*sqrt(1/2530)), + ((12, 1, 12, 1, -1, 0), -0.2*sqrt(0.5)), + ((7, 1, 7, 3, -1, -2), 0.5*sqrt(5/42)), + ((29, 1, 29, 3, -1, -2), 6*sqrt(2/8555))] + for x in inputs: + self.assertAlmostEqual(wigner_3j(*x[0]), x[1], places=7) + + +if __name__ == "__main__": + seterr(divide="raise", over="raise", invalid="raise") + main()