Merge branch 'marcosertoli/read_spd' of github.com:ukaea/Indica into …

…marcosertoli/read_spd

indica/datatypes.py

@@ -13,29 +13,31 @@
 GeneralDataType = str
 SpecificDataType = str
 
-#: Information on the general datatypes, stored in a dictionary. Keys
-#  are the names of the datatypes, values are tuples where the first
-#  element is a description and the second is the units for the data
-#  (empty if unitless).
+"""
+Information on the Indica datatypes are stored in a dictionary.
+ - Keys:str = datatype name
+ - Values:tuple = (long_name, units) to be assigned to DataArray attrs.
+"""
+
 GENERAL_DATATYPES: Dict[GeneralDataType, Tuple[str, str]] = {
     "angular_freq": (
-        "Angular speed at which a species of ion completes a loop of the Tokamak",
-        "rad s^-1",
+        r"$\omega$",
+        r"$rad \; s^-1$",
     ),
     "asymmetry": (
-        "Parameter describing asymmetry between quantities on HFS and LFS",
+        "Asymmetry",
         "",
     ),
     "concentration": (
-        "Portion of nuclei which are the given type of ion",
+        "Concentration",
         "%",
     ),
     "effective_charge": (
-        "Zeff ratio of ion charge to electron charge in plasma",
+        "Zeff",
         "",
     ),
-    "emissivity": ("Radiation power produced per unit volume of space", "W m^-3"),
-    "line_emission": ("Line emission from excitation", "W m^3"),
+    "emissivity": ("Emissivity", r"$W \; m^{-3}$"),
+    "line_emission": ("Line", r"$W \; m^{-3}$"),
     "luminous_flux": (
         "Radiation power received per unit area at some point",
         "W m^-2",
@@ -51,16 +53,16 @@
         "m^3",
     ),
     "f_value": (
-        "Product of toroidal magnetic field strength and major radius",
+        "f",
         "Wb m",
     ),
-    "magnetic_flux": ("Poloidal component of magnetic flux", "Wb/2\\pi"),
+    "magnetic_flux": ("Poloidal flux", r"$Wb/2\pi$"),
     "norm_flux_pol": (
-        "Square root of normalised poloidal component of magnetic flux",
+        "Rho-poloidal",
         "",
     ),
     "norm_flux_tor": (
-        "Square root of normalised toroidal component of magnetic flux",
+        "Rho-toroidal",
         "",
     ),
     "toroidal_flux": ("Toroidal component of magnetic flux", "Wb"),

indica/models/bolometer_camera.py

@@ -57,6 +57,7 @@ def __call__(
         Lz: dict = None,
         t: LabeledArray = None,
         calc_rho=False,
+        **kwargs,
     ):
         """
         Calculate diagnostic measured values

indica/models/helike_spectroscopy.py

@@ -82,6 +82,8 @@ def __init__(
         self.measured_Ti: dict
         self.measured_Nimp: dict
         self.pos: dict
+        self.pos_err_in: dict
+        self.pos_err_out: dict
         self.spectra: DataArray
 
         self.Te: DataArray
@@ -237,6 +239,8 @@ def _moment_analysis(self):
         self.line_emission = line_emission
 
         rho_mean = {}
+        rho_err_in = {}
+        rho_err_out = {}
         measured_intensity = {}
         emission_los = {}
         measured_Te = {}
@@ -256,6 +260,23 @@ def _moment_analysis(self):
             rho_mean[line] = (emission_los * rho_los).sum(
                 "los_position", skipna=True
             ) / emission_sum
+            rho_err = rho_los  # rho_los.where(indx_err, np.nan,)
+            where_in = rho_err < rho_mean[line]
+            where_out = rho_err > rho_mean[line]
+            rho_in = xr.where(where_in, rho_los, np.nan)
+            rho_out = xr.where(where_out, rho_los, np.nan)
+            rho_err_in[line] = (
+                (emission_los * (rho_in - rho_mean[line]) ** 2).sum(
+                    "los_position", skipna=True
+                )
+                / emission_sum
+            ) ** 0.5
+            rho_err_out[line] = (
+                (emission_los * (rho_out - rho_mean[line]) ** 2).sum(
+                    "los_position", skipna=True
+                )
+                / emission_sum
+            ) ** 0.5
             measured_intensity[line] = los_integral
             emission_los[line] = emission_los
 
@@ -281,6 +302,8 @@ def _moment_analysis(self):
             ) / emission_sum
 
         self.pos = rho_mean
+        self.pos_err_in = rho_err_in
+        self.pos_err_out = rho_err_out
         self.measured_intensity = measured_intensity
         self.emission_los = emission_los
         self.measured_Te = measured_Te
@@ -333,6 +356,8 @@ def _build_bckc_dictionary(self):
 
                 if line in self.pos.keys():
                     self.bckc[quantity].attrs["pos"] = self.pos[line]
+                    self.bckc[quantity].attrs["pos_err_in"] = self.pos_err_in[line]
+                    self.bckc[quantity].attrs["pos_err_out"] = self.pos_err_out[line]
 
             if "int_k" in self.bckc.keys() and "int_w" in self.bckc.keys():
                 self.bckc["int_k/int_w"] = self.bckc["int_k"] / self.bckc["int_w"]
@@ -430,7 +455,9 @@ def __call__(
         return self.bckc
 
 
-def example_run(pulse: int = None, plasma=None, plot=False, **kwargs):
+def example_run(
+    pulse: int = None, plasma=None, plot=False, moment_analysis: bool = False, **kwargs
+):
     # TODO: LOS sometimes crossing bad EFIT reconstruction
     if plasma is None:
         plasma = example_plasma(
@@ -467,7 +494,7 @@ def example_run(pulse: int = None, plasma=None, plot=False, **kwargs):
     model.set_los_transform(los_transform)
     model.set_plasma(plasma)
 
-    bckc = model(**kwargs)
+    bckc = model(moment_analysis=moment_analysis, **kwargs)
 
     channels = model.los_transform.x1
     cols = cm.gnuplot2(np.linspace(0.1, 0.75, len(channels), dtype=float))

indica/models/sxr_camera.py

@@ -57,6 +57,7 @@ def __call__(
         Lz: dict = None,
         t: LabeledArray = None,
         calc_rho=False,
+        **kwargs,
     ):
         """
         Calculate diagnostic measured values

indica/models/thomson_scattering.py

@@ -63,6 +63,7 @@ def __call__(
         Te: DataArray = None,
         t: LabeledArray = None,
         calc_rho: bool = False,
+        **kwargs,
     ):
         """
         Calculate diagnostic measured values