Marcosertoli/zeff profile (#272)

* New reader from spectroscopy pre-analysis (currently reading only emissivity and named "spectra") * Added gaunt factor approximation from Carson 1988 * first commit * first commit * second commit, updated background_fit * added **kwargs to call function * test * he_like * added pi diagnostic * new changes * commit * error with dimensions * now it works, but the data is not what is wanted to be * changed impurity from ar to c changed saving path * dict * updated dict, question about model_data in bayes * now bayes works * Temporarily stop caching (issues to be solved) * Copied Bremss bayes opt to new file * Marcosertoli/satakihino (#263) * Temporarily stop caching (issues to be solved) * Copied Bremss bayes opt to new file --------- Co-authored-by: Marco Sertoli <[email protected]> * CHanged .sel to .interp for plasma attributes as in other models. * Ran pre-commit * Fixed pre-commit errors * Fixed pre-commit errors * Reconverted interp to sel. Made transmission filter a DataArray. Started moving Aleksandra's spectra filtering to Diode model * Fixed pre-commit * Temporarily commented out t-interp case for diode filter (sel instead of interp...) * st40reader.py passing pre-commit * Moved background_fit to a method inside the diode filter model and generalised for spectral integration. * Fixed plotting axes ranges * Added generalised example that can use phantom data/equilibrium as well as experimental data from SXR cameras and diode filter models/data. * Added emissivity attribute = local Brems emission integrated over wavelength. Changed example geometry adding LOS to make inversion possible. * First stab at file with workflows to calculate Zeff using tomo_1d or Bayesian inference * Adapted tomo1d examples to run on "pi" bremsstrahlung. Working on phantom and real data. * Temporary fix to slash channels in spectrometer that are not acquired * Corrected factor multiplication for zeff calculation from bremsstrahlung * Initialize variables now easily callable after having instatiated the class * New example geometry and in standalone method so to be usable also from outside * Refactored methods to estimate Zeff from bremsstrahlung using inversion or bayesian inference. * Update bayesmodels.py adding zeff to saved quantites * Update bayes_workflow.py zeff added to plots * Update example_bayes_opt_brems.py zeff added to phantom profiles * Update bayesmodels.py Zeff -> zeff * Modified spectrometer reader to select good channels only * Zeff blob summed over elements * Zeff blob summed over elements * Set uncertainty to zero until correctly calculated * Changed filter limits to nested dictionary to be applicable to diagnostics with multiple quantities * Corrected minor typo * Formatting mod to pass pre-commit * Temporary fix to get abstract reader test to pass. * Modified workflows to accept both phantoms and experimental data, including fitting of TS. Four example methods to run at the end of the file. * Cleaned up so inference works better. Still some inconsistencies between Zeff and Impurity density... * Refactored to pass pre-commit * Marcosertoli/reorganize bremss calculation (#267) * Temporarily stop caching (issues to be solved) * Copied Bremss bayes opt to new file * CHanged .sel to .interp for plasma attributes as in other models. * Ran pre-commit * Fixed pre-commit errors * Fixed pre-commit errors * Reconverted interp to sel. Made transmission filter a DataArray. Started moving Aleksandra's spectra filtering to Diode model * Fixed pre-commit * Temporarily commented out t-interp case for diode filter (sel instead of interp...) * st40reader.py passing pre-commit * Moved background_fit to a method inside the diode filter model and generalised for spectral integration. * Fixed plotting axes ranges * Added generalised example that can use phantom data/equilibrium as well as experimental data from SXR cameras and diode filter models/data. * Added emissivity attribute = local Brems emission integrated over wavelength. Changed example geometry adding LOS to make inversion possible. * First stab at file with workflows to calculate Zeff using tomo_1d or Bayesian inference * Adapted tomo1d examples to run on "pi" bremsstrahlung. Working on phantom and real data. * Temporary fix to slash channels in spectrometer that are not acquired * Corrected factor multiplication for zeff calculation from bremsstrahlung * Initialize variables now easily callable after having instatiated the class * New example geometry and in standalone method so to be usable also from outside * Refactored methods to estimate Zeff from bremsstrahlung using inversion or bayesian inference. * Update bayesmodels.py adding zeff to saved quantites * Update bayes_workflow.py zeff added to plots * Update example_bayes_opt_brems.py zeff added to phantom profiles * Update bayesmodels.py Zeff -> zeff * Modified spectrometer reader to select good channels only * Zeff blob summed over elements * Zeff blob summed over elements * Set uncertainty to zero until correctly calculated * Changed filter limits to nested dictionary to be applicable to diagnostics with multiple quantities * Corrected minor typo * Formatting mod to pass pre-commit * Temporary fix to get abstract reader test to pass. * Modified workflows to accept both phantoms and experimental data, including fitting of TS. Four example methods to run at the end of the file. * Cleaned up so inference works better. Still some inconsistencies between Zeff and Impurity density... * Refactored to pass pre-commit --------- Co-authored-by: Marco Sertoli <[email protected]> Co-authored-by: michael-gemmell <[email protected]> * Fixed linting issues * Modified Bayesian phantoms and priors & changed plotting to include a rough uncertainty in Zeff * Limited Zeff plotting to (0, 10) * Moved start and end times to call kwargs * Refactored workflow to get all the steps in the right order. Both Bayes and Inv now working on phantoms. Inv also on data --> need to find other experimental cases with better EFIT. * Deleted unused example. * Fixed pre-commit error * Fixed issues tied to "c" vs "ar" * Profile assignment now element dependent * Deleted print statement in interferometry model * Deleted example_bayes_opt_brems.py as all zeff profile workflows are now contained in zeff_profile.py * Deleted background_fit.py that is now included in filter model * Reverted to previous figheader --------- Co-authored-by: Marco Sertoli <[email protected]> Co-authored-by: Aleksandra Alieva <[email protected]> Co-authored-by: satakihino <[email protected]> Co-authored-by: michael-gemmell <[email protected]>
indica-mcf · Jul 31, 2023 · ae49860 · ae49860
1 parent 60f81f1
commit ae49860
Show file tree

Hide file tree

Showing 18 changed files with 1,023 additions and 128 deletions.
diff --git a/Reconstruction.npz b/Reconstruction.npz
diff --git a/indica/bayesmodels.py b/indica/bayesmodels.py
@@ -76,11 +76,11 @@ def _ln_likelihood(self):
             # TODO: What to use as error?  Assume percentage error if none given...
             # Float128 is used since rounding of small numbers causes
             # problems when initial results are bad fits
-            model_data = self.bckc[key].values.astype("float128")
+            model_data = self.bckc[key].values.astype("float64")
             exp_data = (
                 self.data[key]
                 .sel(t=self.plasma.time_to_calculate)
-                .values.astype("float128")
+                .values.astype("float64")
             )
             _ln_likelihood = np.log(gaussian(model_data, exp_data, exp_data * 0.10))
             ln_likelihood += np.nanmean(_ln_likelihood)
@@ -185,6 +185,9 @@ def ln_posterior(self, parameters: dict, **kwargs):
             "impurity_density": self.plasma.impurity_density.sel(
                 t=self.plasma.time_to_calculate
             ),
+            "zeff": self.plasma.zeff.sum("element").sel(
+                t=self.plasma.time_to_calculate
+            ),
             # TODO: add Nh
         }
         blob = deepcopy({**self.bckc, **kin_profs})

diff --git a/indica/models/diode_filters.py b/indica/models/diode_filters.py
@@ -27,12 +27,13 @@ class BremsstrahlungDiode(DiagnosticModel):
     def __init__(
         self,
         name: str,
-        filter_wavelength: float = 530.0,
-        filter_fwhm: float = 10,
+        filter_wavelength: float = 531.5,  # 532
+        filter_fwhm: float = 1,  # 1
         filter_type: str = "boxcar",
         etendue: float = 1.0,
         calibration: float = 2.0e-5,
         instrument_method="get_diode_filters",
+        channel_mask: slice = None,  # =slice(18, 28),
     ):
         """
         Filtered diode diagnostic measuring Bremsstrahlung
@@ -63,6 +64,70 @@ def __init__(
         self.calibration = calibration
         self.instrument_method = instrument_method
         self.quantities = AVAILABLE_QUANTITIES[self.instrument_method]
+        self.channel_mask = channel_mask
+
+        wavelength = np.linspace(
+            self.filter_wavelength - self.filter_fwhm * 2,
+            self.filter_wavelength + self.filter_fwhm * 2,
+        )
+        self.wavelength = DataArray(wavelength, coords=[("wavelength", wavelength)])
+
+        # Transmission filter function
+        transmission = ph.make_window(
+            wavelength,
+            self.filter_wavelength,
+            self.filter_fwhm,
+            window=self.filter_type,
+        )
+        self.transmission = DataArray(transmission, coords=[("wavelength", wavelength)])
+
+    def integrate_spectra(self, spectra: DataArray, fit_background: bool = True):
+        """
+        Apply the diode transmission function to an input spectra
+
+        Parameters
+        ----------
+        spectra
+            Spectra with dimensions (channel, wavelength, time) in any order,
+            and with units of W/m**2
+        fit_background
+            If True - background fitted and then integrated
+            If False - spectra integrated using filter without any fitting
+
+        Returns
+        -------
+        Background emission & integral of spectra using the filter transmission
+        TODO: uncertainty on fit not calculated
+        """
+
+        # Interpolate transmission filter on spectral wavelength & restrict to > 0
+        _transmission = self.transmission.interp(wavelength=spectra.wavelength)
+        transmission = _transmission.where(_transmission > 1.0e-3, drop=True)
+        wavelength_slice = slice(
+            np.min(transmission.wavelength), np.max(transmission.wavelength)
+        )
+
+        # Take away neutron spikes in pixel intensity
+        _spectra = spectra.sortby("wavelength").sel(wavelength=wavelength_slice)
+        _spectra = _spectra.where(
+            xr.ufuncs.fabs(_spectra.diff("wavelength", n=2)) < 0.4e-3
+        )
+
+        # Fit spectra to calculate background emission, filter and integrate
+        if fit_background:
+            fit = _spectra.polyfit("wavelength", 0)
+            _spectra_to_integrate = fit.polyfit_coefficients.sel(degree=0)
+            spectra_to_integrate = _spectra_to_integrate.expand_dims(
+                dim={"wavelength": _spectra.wavelength}
+            )
+        else:
+            spectra_to_integrate = _spectra
+        integral = (spectra_to_integrate * transmission).sum("wavelength")
+
+        integral.attrs["error"] = integral * 0.0
+        spectra_to_integrate.attrs["error"] = spectra_to_integrate * 0.0
+
+        return spectra_to_integrate, integral
 
     def _build_bckc_dictionary(self):
         self.bckc = {}
@@ -81,7 +146,7 @@ def _build_bckc_dictionary(self):
                     "stdev": stdev,
                     "provenance": str(self),
                     "long_name": "Brightness",
-                    "units": "W $m^{-2}$",
+                    "units": "W m^{-2}",
                 }
             else:
                 print(f"{quant} not available in model for {self.instrument_method}")
@@ -94,6 +159,7 @@ def __call__(
         Zeff: DataArray = None,
         t: LabeledArray = None,
         calc_rho: bool = False,
+        **kwargs,
     ):
         """
         Calculate Bremsstrahlung emission and model measurement
@@ -110,14 +176,15 @@ def __call__(
             Total effective charge
         t
             time
+        TODO: emission needs a new name as it's in units [W m**-2 nm**-1]
         """
 
         if self.plasma is not None:
             if t is None:
-                t = self.plasma.t
-            Ne = self.plasma.electron_density.interp(t=t)
-            Te = self.plasma.electron_temperature.interp(t=t)
-            Zeff = self.plasma.zeff.interp(t=t).sum("element")
+                t = self.plasma.time_to_calculate
+            Ne = self.plasma.electron_density.sel(t=t)
+            Te = self.plasma.electron_temperature.sel(t=t)
+            Zeff = self.plasma.zeff.sel(t=t).sum("element")
         else:
             if Ne is None or Te is None or Zeff is None:
                 raise ValueError("Give inputs of assign plasma class!")
@@ -127,50 +194,51 @@ def __call__(
         self.Ne: DataArray = Ne
         self.Zeff: DataArray = Zeff
 
-        # Wavelength axis
-        wavelength = np.linspace(
-            self.filter_wavelength - self.filter_fwhm * 2,
-            self.filter_wavelength + self.filter_fwhm * 2,
-        )
-        self.wavelength = DataArray(wavelength, coords=[("wavelength", wavelength)])
-
-        # Transmission filter function
-        self.transmission = ph.make_window(
-            wavelength,
-            self.filter_wavelength,
-            self.filter_fwhm,
-            window=self.filter_type,
-        )
-
         # Bremsstrahlung emission for each time, radial position and wavelength
         wlength = deepcopy(self.wavelength)
         for dim in Ne.dims:
             wlength = wlength.expand_dims(dim={dim: self.Ne[dim]})
         self.emission = ph.zeff_bremsstrahlung(Te, Ne, wlength, zeff=Zeff)
-
+        self.emissivity = (self.emission * self.transmission).integrate("wavelength")
         los_integral = self.los_transform.integrate_on_los(
-            (self.emission * self.transmission).integrate("wavelength"),
+            self.emissivity,
             t=t,
             calc_rho=calc_rho,
         )
+        if self.channel_mask is not None:
+            los_integral = los_integral.where(
+                (los_integral.channel > self.channel_mask.start)
+                & (los_integral.channel < self.channel_mask.stop)
+            )
 
         self.los_integral = los_integral
 
         self._build_bckc_dictionary()
-
         return self.bckc
 
 
-def example_run(pulse: int = None, plasma=None, plot: bool = False):
+def example_geometry(nchannels: int = 12):
+
+    los_end = np.full((nchannels, 3), 0.0)
+    los_end[:, 0] = 0.0
+    los_end[:, 1] = np.linspace(-0.2, -1, nchannels)
+    los_end[:, 2] = 0.0
+    los_start = np.array([[1.5, 0, 0]] * los_end.shape[0])
+    origin = los_start
+    direction = los_end - los_start
+
+    return origin, direction
+
+
+def example_run(
+    pulse: int = None, nchannels: int = 12, plasma=None, plot: bool = False
+):
     if plasma is None:
         plasma = example_plasma(pulse=pulse)
 
     # Create new interferometers diagnostics
     diagnostic_name = "diode_brems"
-    los_start = np.array([[0.8, 0, 0], [0.8, 0, -0.1], [0.8, 0, -0.2]])
-    los_end = np.array([[0.17, 0, 0], [0.17, 0, -0.25], [0.17, 0, -0.2]])
-    origin = los_start
-    direction = los_end - los_start
+    origin, direction = example_geometry(nchannels=nchannels)
     los_transform = LineOfSightTransform(
         origin[:, 0],
         origin[:, 1],
@@ -214,8 +282,8 @@ def example_run(pulse: int = None, plasma=None, plot: bool = False):
         plt.figure()
         for i, t in enumerate(plasma.t.values):
             plt.plot(
-                model.emission.rho_poloidal,
-                model.emission.sel(t=t).integrate("wavelength"),
+                model.emissivity.rho_poloidal,
+                model.emissivity.sel(t=t),
                 color=cols_time[i],
                 label=f"t={t:1.2f} s",
             )

diff --git a/indica/models/helike_spectroscopy.py b/indica/models/helike_spectroscopy.py
@@ -351,7 +351,7 @@ def __call__(
         Nh: DataArray = None,
         t: LabeledArray = None,
         calc_rho: bool = False,
-        moment_analysis: bool = True,
+        moment_analysis: bool = False,
         **kwargs,
     ):
         """
@@ -427,7 +427,6 @@ def __call__(
             self._moment_analysis()
 
         self._build_bckc_dictionary()
-
         return self.bckc
 
 

diff --git a/indica/models/interferometry.py b/indica/models/interferometry.py
@@ -87,7 +87,6 @@ def __call__(
         self.los_integral_ne = los_integral_ne
 
         self._build_bckc_dictionary()
-
         return self.bckc