Adding support for the ADF11 Ionisation and Recombination rates.

cherab/openadas/install.py

@@ -28,6 +28,12 @@
 def install_files(configuration, download=False, repository_path=None, adas_path=None):
 
     for adf in configuration:
+        if adf.lower() == 'adf11scd':
+            for args in configuration[adf]:
+                install_adf11scd(*args, download=download, repository_path=repository_path, adas_path=adas_path)
+        if adf.lower() == 'adf11acd':
+            for args in configuration[adf]:
+                install_adf11acd(*args, download=download, repository_path=repository_path, adas_path=adas_path)
         if adf.lower() == 'adf12':
             for args in configuration[adf]:
                 install_adf12(*args, download=download, repository_path=repository_path, adas_path=adas_path)
@@ -45,6 +51,48 @@ def install_files(configuration, download=False, repository_path=None, adas_path
                 install_adf22bme(*args, download=download, repository_path=repository_path, adas_path=adas_path)
 
 
+def install_adf11scd(element, file_path, download=False, repository_path=None, adas_path=None):
+    """
+    Adds the ionisation rate defined in an ADF11 file to the repository.
+
+    :param element: The element described by the rate file.
+    :param file_path: Path relative to ADAS root.
+    :param download: Attempt to download file if not present (Default=True).
+    :param repository_path: Path to the repository in which to install the rates (optional).
+    :param adas_path: Path to ADAS files repository (optional).
+    """
+
+    print('Installing {}...'.format(file_path))
+    path = _locate_adas_file(file_path, download, adas_path)
+    if not path:
+        raise ValueError('Could not locate the specified ADAS file.')
+
+    # decode file and write out rates
+    rate = parse_adf11(element, path)
+    repository.update_ionisation_rates(rate, repository_path)
+
+
+def install_adf11acd(element, file_path, download=False, repository_path=None, adas_path=None):
+    """
+    Adds the recombination rate defined in an ADF11 file to the repository.
+
+    :param element: The element described by the rate file.
+    :param file_path: Path relative to ADAS root.
+    :param download: Attempt to download file if not present (Default=True).
+    :param repository_path: Path to the repository in which to install the rates (optional).
+    :param adas_path: Path to ADAS files repository (optional).
+    """
+
+    print('Installing {}...'.format(file_path))
+    path = _locate_adas_file(file_path, download, adas_path)
+    if not path:
+        raise ValueError('Could not locate the specified ADAS file.')
+
+    # decode file and write out rates
+    rate = parse_adf11(element, path)
+    repository.update_recombination_rates(rate, repository_path)
+
+
 # todo: move print calls to logging
 def install_adf12(donor_ion, donor_metastable, receiver_ion, receiver_ionisation, file_path, download=False, repository_path=None, adas_path=None):
     """

cherab/openadas/openadas.py

@@ -56,7 +56,37 @@ def wavelength(self, ion, ionisation, transition):
             ion = ion.element
         return repository.get_wavelength(ion, ionisation, transition)
 
-    def beam_cx_rate(self, donor_ion, receiver_ion, receiver_ionisation, transition):
+    def ionisation_rate(self, ion, ionisation):
+
+        if isinstance(ion, Isotope):
+            ion = ion.element
+
+        try:
+            data = repository.get_ionisation_rate(ion, ionisation)
+
+        except RuntimeError:
+            if self._missing_rates_return_null:
+                return NullIonisationRate()
+            raise
+
+        return IonisationRate(data, extrapolate=self._permit_extrapolation)
+
+    def recombination_rate(self, ion, ionisation):
+
+        if isinstance(ion, Isotope):
+            ion = ion.element
+
+        try:
+            data = repository.get_recombination_rate(ion, ionisation)
+
+        except RuntimeError:
+            if self._missing_rates_return_null:
+                return NullRecombinationRate()
+            raise
+
+        return RecombinationRate(data, extrapolate=self._permit_extrapolation)
+
+    def beam_cx_pec(self, donor_ion, receiver_ion, receiver_ionisation, transition):
         """
 
         :param donor_ion:
@@ -80,13 +110,13 @@ def beam_cx_rate(self, donor_ion, receiver_ion, receiver_ionisation, transition)
 
         except RuntimeError:
             if self._missing_rates_return_null:
-                return [NullBeamCXRate()]
+                return [NullBeamCXPEC()]
             raise
 
         # load and interpolate the relevant transition data from each file
         rates = []
         for donor_metastable, rate_data in data:
-            rates.append(BeamCXRate(donor_metastable, wavelength, rate_data, extrapolate=self._permit_extrapolation))
+            rates.append(BeamCXPEC(donor_metastable, wavelength, rate_data, extrapolate=self._permit_extrapolation))
         return rates
 
     def beam_stopping_rate(self, beam_ion, plasma_ion, ionisation):
@@ -146,7 +176,7 @@ def beam_population_rate(self, beam_ion, metastable, plasma_ion, ionisation):
         # load and interpolate data
         return BeamPopulationRate(data, extrapolate=self._permit_extrapolation)
 
-    def beam_emission_rate(self, beam_ion, plasma_ion, ionisation, transition):
+    def beam_emission_pec(self, beam_ion, plasma_ion, ionisation, transition):
         """
 
         :param beam_ion:
@@ -170,13 +200,13 @@ def beam_emission_rate(self, beam_ion, plasma_ion, ionisation, transition):
 
         except RuntimeError:
             if self._missing_rates_return_null:
-                return NullBeamEmissionRate()
+                return NullBeamEmissionPEC()
             raise
 
         # load and interpolate data
-        return BeamEmissionRate(data, wavelength, extrapolate=self._permit_extrapolation)
+        return BeamEmissionPEC(data, wavelength, extrapolate=self._permit_extrapolation)
 
-    def impact_excitation_rate(self, ion, ionisation, transition):
+    def impact_excitation_pec(self, ion, ionisation, transition):
         """
 
         :param ion:
@@ -194,12 +224,12 @@ def impact_excitation_rate(self, ion, ionisation, transition):
 
         except RuntimeError:
             if self._missing_rates_return_null:
-                return NullImpactExcitationRate()
+                return NullImpactExcitationPEC()
             raise
 
-        return ImpactExcitationRate(wavelength, data, extrapolate=self._permit_extrapolation)
+        return ImpactExcitationPEC(wavelength, data, extrapolate=self._permit_extrapolation)
 
-    def recombination_rate(self, ion, ionisation, transition):
+    def recombination_pec(self, ion, ionisation, transition):
         """
 
         :param ion:
@@ -217,10 +247,10 @@ def recombination_rate(self, ion, ionisation, transition):
 
         except (FileNotFoundError, KeyError):
             if self._missing_rates_return_null:
-                return NullRecombinationRate()
+                return NullRecombinationPEC()
             raise
 
-        return RecombinationRate(wavelength, data, extrapolate=self._permit_extrapolation)
+        return RecombinationPEC(wavelength, data, extrapolate=self._permit_extrapolation)
 
     # def stage_resolved_line_ radiation_rate(self, ion, ionisation):
     #

cherab/openadas/parse/__init__.py

@@ -1,3 +1,5 @@
+
+from .adf11 import parse_adf11
 from .adf12 import parse_adf12
 from .adf15 import parse_adf15
 from .adf21 import parse_adf21

cherab/openadas/parse/adf11.py

@@ -16,10 +16,11 @@
 # See the Licence for the specific language governing permissions and limitations
 # under the Licence.
 
-import numpy as np
 import re
+import numpy as np
 
-# todo: to be implemented in a future release
+from cherab.core.atomic import Element
+from cherab.core.utility import RecursiveDict, Cm3ToM3, PerCm3ToPerM3
 
 
 # def parse_adf11acd(ion, ionisation, adf_file_path):
@@ -58,6 +59,110 @@
 #     pass
 
 
+def parse_adf11(element, adf_file_path):
+    """
+    Reads contents of open adas adf11 files
+
+    :param element: Element described by ADF file.
+    :param adf_file_path: Path to ADF11 file from ADAS root.
+    :return: temperature, density, rates as numpy array
+    """
+
+    if not isinstance(element, Element):
+        raise TypeError('The element must be an Element object.')
+
+    with open(adf_file_path, "r") as source_file:
+
+        lines = source_file.readlines()  # read file contents by lines
+        tmp = re.split("\s{2,}", lines[0].strip())  # split into relevant variables
+        # exctract variables
+        z_nuclear = int(tmp[0])
+        n_densities = int(tmp[1])
+        n_temperatures = int(tmp[2])
+        z_min = int(tmp[3])
+        z_max = int(tmp[4])
+        element_name = tmp[5].strip('/').lower()
+        projectname = tmp[6]
+
+        if element.atomic_number != z_nuclear or element.name != element_name:
+            raise ValueError("The requested element '{}' does not match the element description in the"
+                             "specified ADF11 file, '{}'.".format(element.name, element_name))
+
+        # check if it is a resolved file
+        if re.match("\s*[0-9]+", lines[3]):  # is it unresolved?
+            startsearch = 2
+        else:
+            startsearch = 4  # skip vectors with info about resolved states
+
+        # get temperature and density vectors
+        for i in range(startsearch, len(lines)):
+            if re.match("^\s*C{0}-{2,}", lines[i]):
+                tmp = re.sub("\n*\s+", "\t",
+                             "".join(lines[startsearch:i]).strip())  # replace unwanted chars
+                tmp = np.fromstring(tmp, sep="\t", dtype=float)  # put into nunpy array
+                densities = tmp[:n_densities]  # read density values
+                densities = 10**densities  # convert from log10(cm^-3) to cm^-3
+                densities = PerCm3ToPerM3.to(densities)  # convert units from cm^-3 to m^-3
+                temperatures = tmp[n_densities:]  # read temperature values
+                temperatures = 10**temperatures  # convert from log10(eV) to eV
+                startsearch = i
+                break
+
+        # process rate data
+        rates = RecursiveDict()
+
+        # get beginnig and end of requested rates data block and add it to xarray
+        blockrates_start = None
+        blockrates_stop = None
+        for i in range(startsearch, len(lines)):
+
+            if re.match("^\s*C*-{2,}", lines[i]):  # is it a rates block header?
+
+                # is it a first data block found?
+                if not blockrates_start is None:
+                    blockrates_stop = i  # end of the requested block
+
+                    rates_table = re.sub("\n*\s+", "\t",
+                                         "".join(lines[
+                                                 blockrates_start:blockrates_stop]).strip())  # replace unwanted chars
+                    rates_table = np.fromstring(rates_table, sep="\t",
+                                                dtype=float).reshape((n_temperatures,
+                                                                      n_densities))  # transform into an array
+
+                    # convert units from cm^-3 to m^-3
+                    rates_table = 10**rates_table
+                    rates_table = Cm3ToM3.to(rates_table)
+
+                    rates[element][ion_charge]['ne'] = densities
+                    rates[element][ion_charge]['te'] = temperatures
+                    rates[element][ion_charge]['rates'] = np.swapaxes(rates_table, 0, 1)
+
+                    print()
+                    print("density", densities.shape)
+                    print("temperatures", temperatures.shape)
+                    print("rates", rates_table.shape)
+                    print()
+
+                    # if end of data block beak the loop or reassign start of data block for next stage
+                    if re.match("^\s*C{1}-{2,}", lines[i]) or re.match("^\s*C{0,1}-{2,}", lines[i]) and \
+                            re.match("^\s*C\n", lines[i + 1]):
+                        break
+
+                z1_pos = re.search("Z1\s*=*\s*[0-9]+\s*", lines[i]).group()  # get Z1 part
+                ion_charge = int(re.sub("Z1[\s*=]", "", z1_pos))  # remove Z1 to avoid getting 1  later
+                if not re.search("IGRD\s*=*\s*[0-9]+\s*", lines[i]) is None:  # get the IGRD part
+                    igrd_pos = re.search("IGRD\s*=*\s*[0-9]+\s*", lines[i]).group()  # get the IGRD part
+                else:
+                    igrd_pos = "No spec"
+                if not re.search("IPRT\s*=*\s*[0-9]+\s*", lines[i]) is None:
+                    iptr_pos = re.search("IPRT\s*=*\s*[0-9]+\s*", lines[i]).group()  # get the IPRT part
+                else:
+                    iptr_pos = "No spec"
+                blockrates_start = i + 1  # if block start not known, check if we are at the right position
+
+        return rates
+
+
 # def parse_adf11xxx(file_path, ion, ionisation):
 #
 #     adf11_fh = open(file_path, 'r')

cherab/openadas/rates/__init__.pxd

@@ -19,3 +19,4 @@
 from cherab.openadas.rates.beam cimport *
 from cherab.openadas.rates.cx cimport *
 from cherab.openadas.rates.pec cimport *
+from cherab.openadas.rates.atomic cimport *

cherab/openadas/rates/__init__.py

@@ -18,4 +18,5 @@
 
 from .beam import *
 from .cx import *
-from .pec import *
+from .pec import *
+from .atomic import *

cherab/openadas/rates/atomic.pxd

@@ -0,0 +1,47 @@
+# Copyright 2016-2018 Euratom
+# Copyright 2016-2018 United Kingdom Atomic Energy Authority
+# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas
+#
+# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the
+# European Commission - subsequent versions of the EUPL (the "Licence");
+# You may not use this work except in compliance with the Licence.
+# You may obtain a copy of the Licence at:
+#
+# https://joinup.ec.europa.eu/software/page/eupl5
+#
+# Unless required by applicable law or agreed to in writing, software distributed
+# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR
+# CONDITIONS OF ANY KIND, either express or implied.
+#
+# See the Licence for the specific language governing permissions and limitations
+# under the Licence.
+
+from cherab.core cimport IonisationRate as CoreIonisationRate
+from cherab.core cimport RecombinationRate as CoreRecombinationRate
+from cherab.core.math cimport Interpolate2DCubic
+
+
+cdef class IonisationRate(CoreIonisationRate):
+
+    cdef:
+        readonly dict raw_data
+        readonly double wavelength
+        readonly tuple density_range, temperature_range
+        Interpolate2DCubic _rate
+
+
+cdef class NullImpactExcitationRate(CoreIonisationRate):
+    pass
+
+
+cdef class RecombinationRate(CoreRecombinationRate):
+
+    cdef:
+        readonly dict raw_data
+        readonly double wavelength
+        readonly tuple density_range, temperature_range
+        Interpolate2DCubic _rate
+
+
+cdef class NullRecombinationRate(CoreRecombinationRate):
+    pass