Development

cherab/compass/equilibrium/compass_upgrade.py

@@ -68,7 +68,6 @@ def from_database(cls, shot_number):
         tmp.name = "Z_strike_point"
         tmp = tmp.to_dataset()
 
-
         equilibrium._data = xr.merge((equilibrium._data, tmp))
 
         equilibrium._data = equilibrium._data.merge(shot["f/fiesta_out"].to_dataset())
@@ -83,15 +82,11 @@ def from_database(cls, shot_number):
         equilibrium._data = equilibrium._data.merge(shot["Bt_vac_mag_axis/fiesta_out"].to_dataset())
         equilibrium._data = equilibrium._data.rename({"Bt_vac_mag_axis": "Btor_vacuum_radius"})
 
-
-        r_boundary_signal = cudb.get_signal("boundary_r/fiesta_out:{0}".format(shot_number))
-        z_boundary_signal = cudb.get_signal("boundary_z/fiesta_out:{0}".format(shot_number))
-
-        boundary = xr.Dataset({"R_lcfs": (("time", "lcfs_vertex"), r_boundary_signal.data.T),
-                               "Z_lcfs": (("time", "lcfs_vertex"), z_boundary_signal.data.T)},
-                              coords={"time": r_boundary_signal.time_axis.data})
-
-        equilibrium._data = xr.merge((equilibrium._data, boundary))
+        tmp = shot["boundary_closed_r/fiesta_out"].to_dataset().rename({"boundary_closed_r": "R_lcfs",
+                                                                 "boundary_closed_r_axis1": "lcfs_vertex"})
+        tmp = tmp.merge(shot["boundary_closed_z/fiesta_out"].to_dataset().rename({"boundary_closed_z": "Z_lcfs",
+                                                                 "boundary_closed_z_axis1": "lcfs_vertex"}))
+        equilibrium._data = equilibrium._data.merge(tmp)
 
         tmp =shot["Z_limiter/fiesta_out"]
         tmp = xr.Dataset({"Z_limiter": ("limiter_vertex", tmp.Z_limiter.values)})
@@ -105,16 +100,6 @@ def from_database(cls, shot_number):
 
         return equilibrium
 
-if __name__ == "__main__":
-    from cherab.tools.equilibrium.plot import plot_equilibrium
-    eq_shot = COMPASSUpgradeEquilibrium.from_database(7400)
-    time = 2.100
-
-    time_slice = eq_shot.time_slice(time)
-    plot_equilibrium(time_slice, detail=True)
-
-
-
 
 
 

cherab/compass/math/profile/electron.pxd

@@ -0,0 +1,15 @@
+from raysect.core.math.function cimport Function1D
+
+
+cdef class ModifiedTanhGaussian(Function1D):
+
+    cdef:
+        double _ped_height, _ped_sol, _ped_width, _ped_pos, _ped_slope, _core_height, _core_width, _core_exp
+
+    cdef double evaluate(self, double x) except? -1e999
+
+    cdef double pedestal(self, double x)
+
+    cdef double mtanh(self, double x)
+
+    cdef double core(self, double x)

cherab/compass/math/profile/electron.pyx

@@ -0,0 +1,157 @@
+from raysect.core.math.function cimport Function1D
+
+from libc.math cimport exp
+
+cdef class ModifiedTanhGaussian(Function1D):
+    """
+    Profile function used at the COMPASS tokamak for electron temperature and density profiles. Pedestal part
+    is described by a modified hyperbolic tanges. The core part is a Normal-like shape with power as a parameter.
+    Based on [1]_.
+
+
+    :param ped_height: Pedestal height
+    :param ped_sol: Scrapeoff layer position
+    :param ped_width: Width of the pedestal
+    :param ped_pos: Position of the pedestal
+    :param ped_slope: Slope of the pedestal
+    :param core_height: Value of the profile at the magnetic axis
+    :param core_width: Width or shape of the curvature of the core profile part
+    :param core_exp: Exponential parameter of the bell-shaped core profile
+
+
+    [1]Stefanikova, E., Peterka, M., Bohm, P., Bilkova, P., Aftanas, M., Sos, M., Urban, J., Hron, M. and Panek, R.,
+     2016. Fitting of the Thomson scattering density and temperature profiles on the COMPASS tokamak.
+    Review of Scientific Instruments, 87(11), p.11E536.
+    """
+
+    def __init__(self, ped_height, ped_sol, ped_width, ped_pos, ped_slope, core_height, core_width, core_exp):
+        super().__init__()
+
+        self.ped_height = ped_height
+        self.ped_sol = ped_sol
+        self.ped_width = ped_width
+        self.ped_pos = ped_pos
+        self.ped_slope = ped_slope
+        self.core_height = core_height
+        self.core_width = core_width
+        self.core_exp = core_exp
+
+    @property
+    def ped_height(self):
+        return self._ped_height
+
+    @ped_height.setter
+    def ped_height(self, value):
+        self._ped_height = value
+
+    @property
+    def ped_sol(self):
+        return self._ped_sol
+
+    @ped_sol.setter
+    def ped_sol(self, value):
+        self._ped_sol = value
+
+
+    @property
+    def ped_width(self):
+        return self._ped_width
+
+    @ped_width.setter
+    def ped_width(self, value):
+        self._ped_width = value
+
+    @property
+    def ped_pos(self):
+        return self._ped_pos
+
+    @ped_pos.setter
+    def ped_pos(self, value):
+        self._ped_pos = value
+
+    @property
+    def ped_slope(self):
+        return self._ped_slope
+
+    @ped_slope.setter
+    def ped_slope(self, value):
+        self._ped_slope = value
+
+    @property
+    def core_height(self):
+        return self._core_height
+
+    @core_height.setter
+    def core_height(self, value):
+        self._core_height = value
+
+
+    @property
+    def core_width(self):
+        return self._core_width
+
+    @core_width.setter
+    def core_width(self, value):
+        self._core_width = value
+
+    @property
+    def core_exp(self):
+        return self._core_exp
+
+    @core_exp.setter
+    def core_exp(self, value):
+        self._core_exp = value
+
+    cdef double evaluate(self, double x) except? -1e999:
+        """
+        Evaluate the function for x
+        :param x: Independent variable
+        :return: Profile value
+        """
+        cdef:
+            double core, ped
+
+        core = self.core(x)
+        ped = self.pedestal(x)
+
+        return ped + (self._core_height - ped) * core
+
+
+    cdef double pedestal(self, double x):
+        """
+        Pedestal part of the profile.
+        :param x: Independent variable
+        :return: Value of pedestal part of profile
+        """
+        cdef:
+            double mtanh
+
+        mtanh = self.mtanh(x)
+
+        return 0.5 * (self._ped_height - self._ped_sol) * (mtanh + 1)
+
+    cdef double mtanh(self, double x):
+        """
+        Modified hyperbolic tangens used in pedestal shape description
+        :param x: Independent variable
+        :return: 
+        """
+        cdef:
+            double xped
+
+        xped = (self._ped_pos - x ) / (2 * self._ped_width) # Independent variable for modified hyp. tangens
+        return ((1 + self._ped_slope * xped) * exp(xped) - exp(-1 * xped)) / (exp(xped) + exp(-1 * xped))
+
+
+    cdef double core(self, double x):
+        """
+        Core part of the profile described by a Normal-like distribution.
+        :param x: Independent variable
+        :return: Core part of the profile
+        """
+        cdef:
+            double exponent
+
+        exponent = -(x / self._core_width) ** self._core_exp
+        return exp(exponent)
+

cherab/compass/metis/cudb.py

@@ -1,11 +1,14 @@
 from cherab.metis.data_source import MetisDataSource_base
-import cdb_extras.xarray_support as cdbxr
+
 from pyCDB import client
+import cdb_extras.xarray_support as cdbxr
+
 import numpy as np
 
+
 class METISFromCUDB(MetisDataSource_base):
 
-    def __init__(self, shot_number:int = None):
+    def __init__(self, shot_number: int = None):
         super().__init__()
 
         self._cudb = client.CDBClient(host="cudb.tok.ipp.cas.cz", data_root="/compass/CC19_COMPASS-U_data/")
@@ -15,44 +18,68 @@ def __init__(self, shot_number:int = None):
         if shot_number is not None:
             self.shot_number = shot_number
 
-
     @property
     def shot_number(self):
         return self._shot_number
 
     @shot_number.setter
-    def shot_number(self, value:int):
+    def shot_number(self, value: int):
         self._shot_number = value
         self._shot = cdbxr.Shot(self._shot_number)
 
         self.get_data()
 
-    def _get_data(self):
+    @property
+    def profil0d_signals(self):
+        return self._profil0d_signals
+
+    @property
+    def zerod_signals(self):
+        return self._zerod_signals
+
+    def update_signal_names(self, shot_number=None):
+        if shot_number is None:
+            shot_number = self._shot_number
+
+        # get all generic signals belonging to the metis output data source and available for the given shot number
         metis_signals = self._cudb.query(
             "SELECT gs.generic_signal_name FROM generic_signals gs INNER JOIN data_signals ds USING(generic_signal_id) WHERE gs.data_source_id = {1} AND ds.record_number = {0}".format(
                 int(self._shot_number), int(self._metis_data_source_id)))
+        # get the metis output signal name
         metis_output = self._cudb.query(
             "SELECT name FROM data_sources WHERE data_source_id = {0}".format(self._metis_data_source_id))
 
-        zerod = {}
-        profil0d = {}
-
+        # get dictionariers with profile names as keys and signal names as items
+        self._profil0d_signals = {}
+        self._zerod_signals = {}
         for i in metis_signals:
             if "profil0d" in i["generic_signal_name"]:
                 key = i["generic_signal_name"].replace("profil0d_", "")
                 signal_name = i["generic_signal_name"] + "/" + metis_output[0]["name"]
-                try:
-                    profil0d[key] = self._shot[signal_name].values.T
-                except Exception:
-                    pass
-
+                self._profil0d_signals[key] = signal_name
             elif "zerod" in i["generic_signal_name"]:
                 key = i["generic_signal_name"].replace("zerod_", "")
                 signal_name = i["generic_signal_name"] + "/" + metis_output[0]["name"]
-                try:
-                    zerod[key] = self._shot[signal_name].values
-                except Exception:
-                    pass
+                self._zerod_signals[key] = signal_name
+
+    def _get_data(self):
+        # update signal names just in case
+        self.update_signal_names()
+
+        zerod = {}
+        profil0d = {}
+
+        for key, item in self._profil0d_signals.items():
+            try:
+                profil0d[key] = self._shot[item].values.T
+            except Exception:
+                pass
+
+        for key, item in self._zerod_signals.items():
+            try:
+                zerod[key] = self._shot[item].values
+            except Exception:
+                pass
 
         profil0d["psin"] = np.divide((profil0d["psi"] - profil0d["psi"][0, :]),
                                      profil0d["psi"][-1, :] - profil0d["psi"][0, :])