diff --git a/indica/converters/time.py b/indica/converters/time.py index cd4ced5f..7c35bd72 100644 --- a/indica/converters/time.py +++ b/indica/converters/time.py @@ -329,7 +329,7 @@ def bin_in_time_dt(tstart: float, tend: float, dt: float, data: DataArray) -> Da ------- : Array like the input, but binned along the time axis. - + TODO: add possibility of doing 50% overlap of time bins! """ check_bounds_bin(tstart, tend, dt, data) tlabels = get_tlabels_dt(tstart, tend, dt) diff --git a/indica/equilibrium.py b/indica/equilibrium.py index 239505b2..04020285 100644 --- a/indica/equilibrium.py +++ b/indica/equilibrium.py @@ -65,12 +65,29 @@ def __init__( self.faxs = equilibrium_data["faxs"] self.fbnd = equilibrium_data["fbnd"] self.ftor = equilibrium_data["ftor"] + self.rmag = equilibrium_data["rmag"] + self.rbnd = equilibrium_data["rbnd"] + self.zmag = equilibrium_data["zmag"] + self.zbnd = equilibrium_data["zbnd"] + self.zx = self.zbnd.min("arbitrary_index") self.rhotor = np.sqrt( (self.ftor - self.ftor.sel(rho_poloidal=0.0)) / (self.ftor.sel(rho_poloidal=1.0) - self.ftor.sel(rho_poloidal=0.0)) ) self.psi = equilibrium_data["psi"] - self.rho = np.sqrt((self.psi - self.faxs) / (self.fbnd - self.faxs)) + + # Including workaround in case faxs or fbnd had messy data + rho: DataArray = np.sqrt((self.psi - self.faxs) / (self.fbnd - self.faxs)) + if np.any(np.isnan(rho.interp(R=self.rmag, z=self.zmag))): + self.faxs = self.psi.interp(R=self.rmag, z=self.zmag).drop(["R", "z"]) + self.fbnd = self.psi.interp(R=self.rbnd, z=self.zbnd).mean( + "arbitrary_index" + ) + rho = np.sqrt((self.psi - self.faxs) / (self.fbnd - self.faxs)) + if np.any(np.isnan(rho)): + rho = xr.where(rho > 0, rho, 0.0) + + self.rho = rho self.t = self.rho.t if "vjac" in equilibrium_data and "ajac" in equilibrium_data: self.psin = equilibrium_data["psin"] @@ -85,11 +102,6 @@ def __init__( if "rmji" and "rmjo" in equilibrium_data: self.rmji = equilibrium_data["rmji"] self.rmjo = equilibrium_data["rmjo"] - self.rmag = equilibrium_data["rmag"] - self.rbnd = equilibrium_data["rbnd"] - self.zmag = equilibrium_data["zmag"] - self.zbnd = equilibrium_data["zbnd"] - self.zx = self.zbnd.min("arbitrary_index") self.R_offset = R_shift self.z_offset = z_shift diff --git a/indica/models/diode_filters.py b/indica/models/diode_filters.py index 860a610a..f42f5253 100644 --- a/indica/models/diode_filters.py +++ b/indica/models/diode_filters.py @@ -31,7 +31,7 @@ def __init__( filter_fwhm: float = 1, # 1 filter_type: str = "boxcar", etendue: float = 1.0, - calibration: float = 2.0e-5, + calibration: float = 1, instrument_method="get_diode_filters", channel_mask: slice = None, # =slice(18, 28), ): diff --git a/indica/models/thomson_scattering.py b/indica/models/thomson_scattering.py index 94af48e8..030b8e52 100644 --- a/indica/models/thomson_scattering.py +++ b/indica/models/thomson_scattering.py @@ -166,7 +166,8 @@ def example_run( alpha=0.7, ) Ne = bckc["ne"].sel(t=t, method="nearest") - plt.scatter(Ne.rho_poloidal, Ne, color=cols_time[i], marker="o", alpha=0.7) + rho = Ne.transform.rho.sel(t=t, method="nearest") + plt.scatter(rho, Ne, color=cols_time[i], marker="o", alpha=0.7) plt.xlabel("Channel") plt.ylabel("Measured electron density (m^-3)") plt.legend() @@ -179,7 +180,8 @@ def example_run( alpha=0.7, ) Te = bckc["te"].sel(t=t, method="nearest") - plt.scatter(Te.rho_poloidal, Te, color=cols_time[i], marker="o", alpha=0.7) + rho = Te.transform.rho.sel(t=t, method="nearest") + plt.scatter(rho, Te, color=cols_time[i], marker="o", alpha=0.7) plt.xlabel("Channel") plt.ylabel("Measured electron temperature (eV)") plt.legend() diff --git a/indica/operators/spline_fit_easy.py b/indica/operators/spline_fit_easy.py index a012f5c0..d6549067 100644 --- a/indica/operators/spline_fit_easy.py +++ b/indica/operators/spline_fit_easy.py @@ -99,63 +99,86 @@ def residuals(yknots): def spline_fit_ts( - pulse: int, - tstart: float = 0.0, - tend: float = 0.2, - plot: bool = False, + pulse: int = 11314, + tstart: float = 0.03, + tend: float = 0.1, + dt: float = 0.01, + quantity: str = "te", + R_shift: float = 0.0, + knots: list = None, + plot: bool = True, ): + st40 = ReadST40(pulse, tstart=tstart, tend=tend, dt=dt) + st40(["ts"], R_shift=R_shift) + + if quantity == "te" and knots is None: + knots = [0, 0.3, 0.6, 0.8, 1.1] + if quantity == "ne" and knots is None: + knots = [0, 0.3, 0.6, 0.8, 0.95, 1.1] + data_all = st40.raw_data["ts"][quantity] + t = data_all.t + transform = data_all.transform + transform.convert_to_rho_theta(t=data_all.t) - ST40 = ReadST40(pulse, tstart=tstart, tend=tend) - ST40(["ts"]) - - Te_data_all = ST40.binned_data["ts"]["te"] - t = ST40.binned_data["ts"]["te"].t - transform = ST40.binned_data["ts"]["te"].transform R = transform.R Rmag = transform.equilibrium.rmag.interp(t=t) # Fit all available TS data - ind = np.full_like(Te_data_all, True) + ind = np.full_like(data_all, True) rho = xr.where(ind, transform.rho, np.nan) - Te_data = xr.where(ind, Te_data_all, np.nan) - Te_err = xr.where(ind, Te_data_all.error, np.nan) - Te_fit = fit_profile( - rho, Te_data, Te_err, knots=[0, 0.3, 0.5, 0.75, 0.95, 1.05], virtual_knots=False - ) + data = xr.where(ind, data_all, np.nan) + err = xr.where(ind, data_all.error, np.nan) + fit = fit_profile(rho, data, err, knots=knots, virtual_knots=False) # Use only HFS channels ind = R <= Rmag - rho = xr.where(ind, transform.rho, np.nan) - Te_data = xr.where(ind, Te_data_all, np.nan) - Te_err = xr.where(ind, Te_data_all.error, np.nan) - Te_fit_hfs = fit_profile(rho, Te_data, Te_err, virtual_knots=True) + rho_hfs = xr.where(ind, transform.rho, np.nan) + data_hfs = xr.where(ind, data_all, np.nan) + err_hfs = xr.where(ind, data_all.error, np.nan) + fit_hfs = fit_profile(rho_hfs, data_hfs, err_hfs, knots=knots, virtual_knots=True) + + # Use only LFS channels + ind = R >= Rmag + rho_lfs = xr.where(ind, transform.rho, np.nan) + data_lfs = xr.where(ind, data_all, np.nan) + err_lfs = xr.where(ind, data_all.error, np.nan) + fit_lfs = fit_profile(rho_lfs, data_lfs, err_lfs, knots=knots, virtual_knots=True) if plot: - for t in Te_data_all.t: + for t in data_all.t: plt.ioff() plt.errorbar( - Te_data_all.transform.rho.sel(t=t), - Te_data_all.sel(t=t), - Te_data_all.error.sel(t=t), + rho_hfs.sel(t=t), + data_hfs.sel(t=t), + err_hfs.sel(t=t), + marker="o", + label="data HFS", + color="blue", + ) + plt.errorbar( + rho_lfs.sel(t=t), + data_lfs.sel(t=t), + err_lfs.sel(t=t), marker="o", - label="data", + label="data LFS", + color="red", ) - Te_fit.sel(t=t).plot( - linewidth=5, alpha=0.5, color="orange", label="spline fit all" + fit.sel(t=t).plot( + linewidth=5, alpha=0.5, color="black", label="spline fit all" ) - Te_fit_hfs.sel(t=t).plot( - linewidth=5, alpha=0.5, color="red", label="spline fit HFS" + fit_lfs.sel(t=t).plot( + linewidth=5, alpha=0.5, color="red", label="spline fit LFS" + ) + fit_hfs.sel(t=t).plot( + linewidth=5, alpha=0.5, color="blue", label="spline fit HFS" ) plt.legend() plt.show() - return Te_data_all, Te_fit + return data_all, fit if __name__ == "__main__": - spline_fit_ts( - 10619, - tstart=0.0, - tend=0.2, - plot=True, - ) + plt.ioff() + spline_fit_ts(11089, quantity="ne") + plt.show() diff --git a/indica/plotters/plot_time_evolution.py b/indica/plotters/plot_time_evolution.py index bc4e897c..7bf950fa 100644 --- a/indica/plotters/plot_time_evolution.py +++ b/indica/plotters/plot_time_evolution.py @@ -1,11 +1,12 @@ +from copy import deepcopy import getpass import matplotlib.pyplot as plt import numpy as np import xarray as xr -from indica.numpy_typing import ArrayLike from indica.readers.read_st40 import ReadST40 +from indica.utilities import save_figure from indica.utilities import set_axis_sci from indica.utilities import set_plot_colors from indica.utilities import set_plot_rcparams @@ -14,99 +15,151 @@ CMAP, COLORS = set_plot_colors() LINESTYLES = ["solid", "dashed", "dotted"] set_plot_rcparams("profiles") +QUANTITIES: list = [ + "smmh:ne", + "xrcs:ti_w", + "xrcs:te_n3w", + "xrcs:te_kw", + "xrcs:spectra", + "sxrc_xy2:brightness", + "efit:ipla", + "efit:rmag", + "efit:zmag", + "sxr_spd:brightness", + # "cxff_pi:ti", + # "cxff_pi:vtor", + # "cxff_tws_c:ti", + # "cxff_tws_c:vtor", + "ts:te", + "ts:ne", + # "lines", +] +Y0 = {} +Y0["nirh1"] = True +Y0["smmh1"] = True +Y0["xrcs"] = True +Y0["sxr_diode_1"] = True +Y0["brems"] = True +Y0["efit"] = True +Y0["cxff_tws_c"] = True +Y0["cxff_pi"] = True -def plot_sxrc( - pulse: int = 10607, # 10605 - what_to_read: list = ["sxrc_xy2:brightness"], - tplot: ArrayLike = [0.067, 0.069], # [0.07, 0.09] - save_fig=False, - plot_raw: bool = False, - xvar: str = "t", - yvar: str = "channel", - tstart=0.02, - tend=0.1, - dt: float = 0.001, - data_key="binned_data", + +def plot_st40_data( + pulses: list, + tstart: float = 0.02, + tend: float = 0.1, + dt: float = 0.005, + quantities: list = [], + tplot: float = None, + save_fig: bool = False, + fig_path: str = None, + fig_style: str = "profiles", + plot_binned: bool = True, ): - instruments = [] - quantities = [] - linestyles = {} - for i, identifier in enumerate(what_to_read): - instr, quant = identifier.split(":") - instruments.append(instr) - quantities.append(quant) - linestyles[instr] = LINESTYLES[i] - - st40 = ReadST40(pulse, tstart, tend, dt=dt) - st40(instruments=instruments) - - data: dict = {} - for instr, quant in zip(instruments, quantities): - if instr not in data.keys(): - data[instr] = {} - data[instr][quant] = getattr(st40, data_key)[instr][quant].transpose(yvar, xvar) - _rho, _theta = data[instr][quant].transform.convert_to_rho_theta( - t=data[instr][quant].t - ) + if tplot is None: + tplot = np.mean([tstart, tend]) + if len(quantities) == 0: + quantities = QUANTITIES + instruments = list(np.unique([quant.split(":")[0] for quant in quantities])) + if fig_path is None: + fig_path = f"{FIG_PATH}" - for instr, quant in zip(instruments, quantities): - plt.figure() - plot = data[instr][quant].plot(label=f"{instr}:{quant}") - set_axis_sci(plot_object=plot) - plt.title(f"{instr.upper()} for pulse {pulse}") + set_plot_rcparams(fig_style) + xr.set_options(keep_attrs=True) + if len(pulses) > 1: + colors = CMAP(np.linspace(0.75, 0.1, len(pulses), dtype=float)) + else: + colors = ["blue"] - if tplot is not None: - cols = CMAP(np.linspace(0.75, 0.1, np.size(tplot), dtype=float)) + raw: dict = {quant: {} for quant in quantities} + binned = deepcopy(raw) + data = {"raw": raw, "binned": binned} + for i, pulse in enumerate(pulses): + fig_path += f"_{pulse}" + st40 = ReadST40(pulse, tstart=tstart, tend=tend, dt=dt) + st40(instruments) + for quantity in quantities: + instr, quant = quantity.split(":") + data["binned"][quantity][pulse] = st40.binned_data[instr][quant] + data["raw"][quantity][pulse] = st40.raw_data[instr][quant] + for quantity in quantities: + print(quantity) + instr, _ = quantity.split(":") plt.figure() - for icol, t in enumerate(tplot): - for instr, quant in zip(instruments, quantities): - _data = data[instr][quant].sel(t=t, method="nearest") - if "error" in _data.attrs: - _err = (_data.error + _data.stdev).sel(t=t, method="nearest") - else: - _err = xr.full_like(_data, 0.0) - _R_diff = ( - _data.transform.impact_parameter.value - - _data.transform.equilibrium.rmag - ).sel(t=t, method="nearest") - plt.fill_between( - _R_diff.values, - _data.values + _err.values, - _data.values - _err.values, - color=cols[icol], - alpha=0.5, - ) - plt.plot( - _R_diff.values, - _data.values, - label=f"t={t:.3f}", - linestyle=linestyles[instr], - color=cols[icol], + for i, pulse in enumerate(pulses): + color = colors[i] + if len(data["raw"][quantity].keys()) == 0: + continue + plot_data(data, quantity, pulse, tplot, key="raw", color=color) + if plot_binned: + plot_data(data, quantity, pulse, tplot, key="binned", color=color) + + set_axis_sci() + if instr in Y0.keys(): + plt.ylim( + 0, ) - plt.xlabel("Impact R - R$_{mag}$") - plt.ylabel(f"{_data.long_name} [{_data.units}]") - plt.title(f"{instr.upper()} for pulse {pulse}") - plt.legend() - set_axis_sci() - - return st40 - - -# def plot_time_surface( -# pulse: int = 10605, -# instruments: list = ["sxrc_xy1"], -# quantity: str = "brightness", -# tplot: ArrayLike = None, -# save_fig=False, -# plot_raw: bool = False, -# xvar: str = "time", -# yvar: str = "channel", -# ): + plt.legend() + plt.autoscale() + save_figure(fig_path, f"{quantity}", save_fig=save_fig) + + return data, st40 + + +def plot_data(data, quantity: str, pulse: int, tplot: float, key="raw", color=None): + str_to_add = "" + instr, quant = quantity.split(":") + if key == "raw": + marker = None + else: + marker = "o" + + _data = data[key][quantity][pulse] + tslice = slice(_data.t.min().values, _data.t.max().values) + if "error" not in _data.attrs: + _data.attrs["error"] = xr.full_like(_data, 0.0) + if "stdev" not in _data.attrs: + _data.attrs["stdev"] = xr.full_like(_data, 0.0) + _err = np.sqrt(_data.error**2 + _data.stdev**2) + _err = xr.where(_err / _data.values < 1.0, _err, 0.0) + if len(_data.dims) > 1: + str_to_add = f" @ {tplot:.3f} s" + tslice = _data.t.sel(t=tplot, method="nearest") + + _data = _data.sel(t=tslice) + _err = _err.sel(t=tslice) + if instr in "xrcs" and quant == "spectra": + bgnd = _data.sel(wavelength=slice(0.393, 0.388)).mean("wavelength") + _data -= bgnd + label = None + if key == "raw": + label = str(pulse) + alpha = 0.5 + plt.fill_between( + _data.coords[_data.dims[0]].values, + _data.values - _err.values, + _data.values + _err.values, + color=color, + alpha=alpha, + ) + if key == "binned": + alpha = 0.8 + plt.errorbar( + _data.coords[_data.dims[0]].values, + _data.values, + _err.values, + color=color, + alpha=alpha, + ) + _data.plot(label=label, color=color, alpha=alpha, marker=marker) + plt.title(f"{instr.upper()} {quant}" + str_to_add) + if __name__ == "__main__": plt.ioff() - plot_sxrc() + plot_st40_data([11226, 11225]) plt.show() diff --git a/indica/readers/abstractreader.py b/indica/readers/abstractreader.py index e34f253e..1bfeec52 100644 --- a/indica/readers/abstractreader.py +++ b/indica/readers/abstractreader.py @@ -862,6 +862,7 @@ def get_diode_filters( t = database_results["times"] t = DataArray(t, coords=[("t", t)], attrs={"long_name": "t", "units": "s"}) + label = database_results["labels"] coords = [("t", t)] if database_results["length"] > 1: coords.append(("channel", np.arange(database_results["length"]))) @@ -876,7 +877,7 @@ def get_diode_filters( coords, transform=transform, ) - data[quantity] = quant_data + data[quantity] = quant_data.assign_coords(label=("channel", label)) return data diff --git a/indica/readers/available_quantities.py b/indica/readers/available_quantities.py index bb5d9fb9..133d892e 100644 --- a/indica/readers/available_quantities.py +++ b/indica/readers/available_quantities.py @@ -1,5 +1,6 @@ """ Quantities that can be read with the current abstract reader implementation +TODO: change the tuple to DataArray (long_name, units) - see examples in abstractreader """ from typing import Dict diff --git a/indica/readers/read_st40.py b/indica/readers/read_st40.py index 8b2c305e..4536e4bb 100644 --- a/indica/readers/read_st40.py +++ b/indica/readers/read_st40.py @@ -12,24 +12,25 @@ from indica.readers import ST40Reader from indica.utilities import print_like -REVISIONS = { - "efit": 0, - "brems": -1, - "halpha": -1, - "nirh1": 0, - "smmh1": 0, - "cxff_pi": 0, - "cxff_tws_c": 0, - "cxqf_tws_c": 0, - "sxr_camera_4": 0, - "sxrc_xy1": 0, - "sxrc_xy2": 0, - "sxr_diode_1": 0, - "xrcs": 0, - "ts": 0, - "pi": 0, - "tws_c": 0, -} +INSTRUMENTS: list = [ + "efit", + "lines", + "nirh1", + "smmh1", + "smmh", + "cxff_pi", + "cxff_tws_c", + "cxqf_tws_c", + "sxr_spd", + "sxr_camera_4", + "sxrc_xy1", + "sxrc_xy2", + "sxr_diode_1", + "xrcs", + "ts", + "pi", + "tws_c", +] FILTER_LIMITS = { "cxff_pi": {"ti": (0, np.inf), "vtor": (0, np.inf)}, @@ -38,11 +39,13 @@ "xrcs": {"ti_w": (0, np.inf), "te_kw": (0, np.inf), "te_n3w": (0, np.inf)}, "brems": {"brightness": (0, np.inf)}, "halpha": {"brightness": (0, np.inf)}, + "sxr_spd": {"brightness": (0, np.inf)}, "sxr_diode_1": {"brightness": (0, np.inf)}, "sxr_camera_4": {"brightness": (0, np.inf)}, "sxrc_xy1": {"brightness": (0, np.inf)}, "sxrc_xy2": {"brightness": (0, np.inf)}, - "ts": {"te": (0, 10.0e3), "ne": (0, 1.0e21)}, + "blom_xy1": {"brightness": (0, np.inf)}, + "ts": {"te": (0, np.inf), "ne": (0, np.inf)}, "pi": {"spectra": (0, np.inf)}, "tws_c": {"spectra": (0, np.inf)}, } @@ -88,8 +91,12 @@ def __init__( self.tend = tend self.dt = dt - self.reader = ST40Reader(pulse, tstart - 0.05, tend + 0.05, tree=tree) - self.reader_equil = ST40Reader(pulse, tstart - 0.1, tend + 0.1, tree=tree) + _tend = tend + dt * 2 + _tstart = tstart - dt * 2 + if _tstart < 0: + _tstart = 0.0 + self.reader = ST40Reader(pulse, _tstart, _tend, tree=tree) + self.reader_equil = ST40Reader(pulse, _tstart, _tend, tree=tree) self.equilibrium: Equilibrium self.raw_data: dict = {} @@ -175,6 +182,7 @@ def map_diagnostics(self, instruments: list, map_raw: bool = False): break def filter_data(self, instruments: list): + if not hasattr(self, "binned_data"): raise ValueError( "Bin data before filtering. No action permitted on raw data structure!" @@ -296,7 +304,7 @@ def plot_profile( def __call__( self, - instruments: list = None, + instruments: list = [], revisions: dict = None, filters: dict = None, map_raw: bool = False, @@ -310,11 +318,15 @@ def __call__( debug: bool = False, ): self.debug = debug - if instruments is None: - instruments = list(REVISIONS.keys()) - if not revisions: - # TODO: fix default behaviour if missing key - revisions = REVISIONS + if len(instruments) == 0: + instruments = INSTRUMENTS + if revisions is None: + revisions = {instrument: 0 for instrument in instruments} + for instr in instruments: + if instr not in revisions.keys(): + revisions[instr] = 0 + if "efit" not in revisions: + revisions["efit"] = 0 if not filters: # TODO: fix default behaviour if missing key filters = FILTER_LIMITS @@ -326,7 +338,7 @@ def __call__( dt = self.dt self.reset_data() - self.get_equilibrium(R_shift=R_shift) + self.get_equilibrium(R_shift=R_shift, revision=revisions["efit"]) for instrument in instruments: print(f"Reading {instrument}") if debug: @@ -366,12 +378,17 @@ def astra_equilibrium(pulse: int, revision: RevisionLike): """Assign ASTRA to equilibrium class""" -def read_cxff_pi(): +def sxr_spd(): import indica.readers.read_st40 as read_st40 + from indica.utilities import set_axis_sci + + st40 = read_st40.ReadST40(11215) + st40(["sxr_spd"]) + + st40.raw_data["sxr_spd"]["brightness"].transform.plot() + + plt.figure() + st40.raw_data["sxr_spd"]["brightness"].sel(channel=0).plot() + set_axis_sci() - st40 = read_st40.ReadST40(10607) - st40(["cxff_pi"]) - st40.raw_data["cxff_pi"]["ti"].los_transform.set_equilibrium( - st40.raw_data["cxff_pi"]["ti"].transform.equilibrium - ) - st40.raw_data["cxff_pi"]["ti"].los_transform.plot() + return st40.raw_data["sxr_spd"]["brightness"] diff --git a/indica/readers/st40reader.py b/indica/readers/st40reader.py index 7855415e..92f567c2 100644 --- a/indica/readers/st40reader.py +++ b/indica/readers/st40reader.py @@ -75,12 +75,13 @@ class ST40Reader(DataReader): "efit": "get_equilibrium", "xrcs": "get_helike_spectroscopy", "princeton": "get_charge_exchange", - "brems": "get_diode_filters", - "halpha": "get_diode_filters", + "lines": "get_diode_filters", "nirh1": "get_interferometry", "nirh1_bin": "get_interferometry", "smmh1": "get_interferometry", + "smmh": "get_interferometry", "astra": "get_astra", + "sxr_spd": "get_diode_filters", "sxr_diode_1": "get_diode_filters", "sxr_diode_2": "get_diode_filters", "sxr_diode_3": "get_diode_filters", @@ -88,6 +89,7 @@ class ST40Reader(DataReader): "sxr_camera_4": "get_radiation", "sxrc_xy1": "get_radiation", "sxrc_xy2": "get_radiation", + "blom_xy1": "get_radiation", "cxff_pi": "get_charge_exchange", "cxff_tws_c": "get_charge_exchange", "cxqf_tws_c": "get_charge_exchange", @@ -95,15 +97,13 @@ class ST40Reader(DataReader): "tws_c": "get_spectrometer", "ts": "get_thomson_scattering", } + # TODO: this will not be necessary once the MDS+ standardisation is complete UIDS_MDS = { "xrcs": "sxr", "princeton": "spectrom", - "brems": "spectrom", - "halpha": "spectrom", "nirh1": "interferom", "nirh1_bin": "interferom", "smmh1": "interferom", - "astra": "", "sxr_diode_1": "sxr", "sxr_diode_2": "sxr", "sxr_diode_3": "sxr", @@ -158,11 +158,14 @@ class ST40Reader(DataReader): "smmh1": { "ne": ".line_int:ne", }, - "brems": { - "brightness": ".brem_mp1:intensity", + "smmh": { + "ne": ".global:ne_int", }, - "halpha": { - "brightness": ".h_alpha_mp1:intensity", + "lines": { + "brightness": ":emission", + }, + "sxr_spd": { + "brightness": ":emission", }, "sxr_camera_4": { "brightness": ".middle_head.filter_4:", @@ -180,6 +183,9 @@ class ST40Reader(DataReader): "sxrc_xy2": { "brightness": ".profiles:emission", }, + "blom_xy1": { + "brightness": ".profiles:emission", + }, "sxr_diode_1": { "brightness": ".filter_001:signal", }, @@ -292,7 +298,8 @@ class ST40Reader(DataReader): }, } - _IMPLEMENTATION_QUANTITIES = { # TODO: these will be different diagnostics!!!!!!! + # TODO: this can be deleted once MDS+ standardisation is complete + _IMPLEMENTATION_QUANTITIES = { "diode_arrays": { # GETTING THE DATA OF THE SXR CAMERA "sxr_camera_1": ("brightness", "total"), "sxr_camera_2": ("brightness", "50_Al_filtered"), @@ -301,6 +308,7 @@ class ST40Reader(DataReader): }, } + # TODO: this can be deleted once MDS+ standardisation is complete _RADIATION_RANGES = { "sxr_camera_1": (1, 20), "sxr_camera_2": (21, 40), @@ -313,7 +321,7 @@ def __init__( pulse: int, tstart: float, tend: float, - server: str = "192.168.1.7", # 192.168.1.7 10.0.40.13 + server: str = "smaug", tree: str = "ST40", default_error: float = 0.05, max_freq: float = 1e6, @@ -698,11 +706,14 @@ def _get_helike_spectroscopy( direction = np.array([direction]) results["times"], _ = self._get_signal(uid, instrument, ":time", revision) - results["wavelength"], _ = self._get_signal( - uid, instrument, ":wavelength", revision - ) + wavelength, _ = self._get_signal(uid, instrument, ":wavelength", revision) # TODO: change once wavelength in MDS+ has been fixed to nanometers! - results["wavelength"] /= 10.0 + wavelength /= 10.0 + if self.pulse >= 10307: + dlambda = float(np.abs(wavelength[1] - wavelength[0])) * 4 + wavelength += dlambda + results["wavelength"] = wavelength + for q in quantities: qval, q_path = self._get_signal( uid, instrument, self.QUANTITIES_MDS[instrument][q], revision @@ -761,6 +772,7 @@ def _get_charge_exchange( z, z_path = self._get_signal(uid, instrument, ":z", revision) R, R_path = self._get_signal(uid, instrument, ":R", revision) + # TODO: temporary fix until geometry sorted (especially pulse if statement..) try: location, location_path = self._get_signal( uid, instrument, ".geometry:location", revision @@ -772,7 +784,6 @@ def _get_charge_exchange( location = np.array([location]) direction = np.array([direction]) - # TODO: temporary fix until geometry sorted if location.shape[0] != x.shape[0]: if self.pulse > 10200: index = np.arange(18, 36) @@ -780,7 +791,6 @@ def _get_charge_exchange( index = np.arange(21, 36) location = location[index] direction = direction[index] - except TreeNNF: location = None direction = None @@ -805,10 +815,8 @@ def _get_charge_exchange( ) except TreeNNF: qval_err = np.full_like(qval, 0.0) - # q_path_err = "" dimensions, _ = self._get_signal_dims(q_path, len(qval.shape)) - results[q + "_records"] = q_path results[q] = qval results[f"{q}_error"] = qval_err @@ -861,13 +869,6 @@ def _get_spectrometer( location = np.array([location]) direction = np.array([direction]) - # if self.pulse > 10200: - # index = np.arange(18, 36) - # else: - # index = np.arange(21, 36) - # location = location[index] - # direction = direction[index] - for q in quantities: qval, q_path = self._get_signal( uid, @@ -885,7 +886,6 @@ def _get_spectrometer( ) except TreeNNF: qval_err = np.full_like(qval, 0.0) - # q_path_err = "" dimensions, _ = self._get_signal_dims(q_path, len(qval.shape)) @@ -910,47 +910,23 @@ def _get_diode_filters( revision: RevisionLike, quantities: Set[str], ) -> Dict[str, Any]: - + """ + TODO: labels are np.bytes_ type...is this correct? + """ if len(uid) == 0 and instrument in self.UIDS_MDS: uid = self.UIDS_MDS[instrument] - # TODO: change once new MDS+ standardisation has been completed - try: - location, location_path = self._get_signal( - uid, instrument, ".geometry:location", revision - ) - direction, position_path = self._get_signal( - uid, instrument, ".geometry:direction", revision - ) - except TreeNNF: - location = np.array( - [ - [1.0, 0, 0], - ] - ) - direction = ( - np.array( - [ - [0.17, 0, 0], - ] - ) - - location - ) + location, location_path = self._get_signal( + uid, instrument, ".geometry:location", revision + ) + direction, position_path = self._get_signal( + uid, instrument, ".geometry:direction", revision + ) if len(np.shape(location)) == 1: location = np.array([location]) direction = np.array([direction]) length = location[:, 0].size - if instrument == "brems": - _instrument = "lines" - revision = -1 - elif instrument == "halpha": - _instrument = "lines" - revision = -1 - elif "sxr_diode" in instrument: - _instrument = "diode_detr" - else: - raise ValueError(f"{instrument} not yet supported") results: Dict[str, Any] = { "length": length, @@ -958,23 +934,29 @@ def _get_diode_filters( } results["location"] = location results["direction"] = direction - results["revision"] = self._get_revision(uid, _instrument, revision) + results["revision"] = self._get_revision(uid, instrument, revision) results["revision"] = revision revision = results["revision"] quantity = "brightness" qval, q_path = self._get_signal( - uid, _instrument, self.QUANTITIES_MDS[instrument][quantity], revision + uid, instrument, self.QUANTITIES_MDS[instrument][quantity], revision ) - times, _ = self._get_signal_dims(q_path, len(qval.shape)) - times = times[0] + times, _ = self._get_signal(uid, instrument, ":time", revision) + _labels, _ = self._get_signal(uid, instrument, ":label", revision) + if type(_labels[0]) == np.bytes_: + labels = np.array([label.decode("UTF-8") for label in _labels]) + else: + labels = _labels + results["times"] = times + results["labels"] = labels results[quantity + "_records"] = q_path results[quantity] = qval try: qval_err, q_path_err = self._get_signal( uid, - _instrument, + instrument, self.QUANTITIES_MDS[instrument][quantity] + "_ERR", revision, ) @@ -993,6 +975,10 @@ def _get_interferometry( revision: RevisionLike, quantities: Set[str], ) -> Dict[str, Any]: + """ + TODO: SMMH 2023 launcher/receiver cross plasma on different poloidal paths! + Currently setting location and direction as average of the two!!!! + """ if len(uid) == 0 and instrument in self.UIDS_MDS: uid = self.UIDS_MDS[instrument] @@ -1011,6 +997,17 @@ def _get_interferometry( direction, direction_path = self._get_signal( uid, instrument, ".geometry:direction", revision ) + + if instrument == "smmh": + location_r, _ = self._get_signal( + uid, instrument, ".geometry:location_r", revision + ) + direction_r, _ = self._get_signal( + uid, instrument, ".geometry:direction_r", revision + ) + location = (location + location_r) / 2.0 + direction = (direction + direction_r) / 2.0 + if len(np.shape(location)) == 1: location = np.array([location]) direction = np.array([direction]) @@ -1083,6 +1080,8 @@ def _get_thomson_scattering( revision = results["revision"] times, times_path = self._get_signal(uid, instrument, ":time", revision) + # TODO: hardcoded correction to TS coordinates to be fixed in MDS+ + print("\n Hardcoded correction to TS coordinates to be fixed in MDS+ \n") # location, location_path = self._get_signal( # uid, instrument, ".geometry:location", revision # ) @@ -1171,6 +1170,7 @@ def get_revision_name(self, revision) -> str: """Return string defining RUN## or BEST if revision = 0""" if type(revision) == int: + rev_str = "" if revision < 0: rev_str = "" elif revision == 0: diff --git a/indica/workflows/load_modelling_plasma.py b/indica/workflows/load_modelling_plasma.py index 957577bc..2b530a42 100644 --- a/indica/workflows/load_modelling_plasma.py +++ b/indica/workflows/load_modelling_plasma.py @@ -11,6 +11,7 @@ from indica.converters.line_of_sight import LineOfSightTransform from indica.converters.transect import TransectCoordinates from indica.equilibrium import Equilibrium +from indica.models.bolometer_camera import Bolometer from indica.models.charge_exchange import ChargeExchange from indica.models.diode_filters import BremsstrahlungDiode from indica.models.helike_spectroscopy import HelikeSpectrometer @@ -29,6 +30,7 @@ CMAP, COLORS = set_plot_colors() DIAGNOSTIC_MODELS = { + "smmh": Interferometry, "smmh1": Interferometry, "nirh1": Interferometry, "xrcs": HelikeSpectrometer, @@ -37,10 +39,23 @@ "cxff_tws_c": ChargeExchange, "cxqf_tws_c": ChargeExchange, "brems": BremsstrahlungDiode, - "sxr_diode_1": SXRcamera, - "sxr_camera_4": SXRcamera, + "sxrc_xy1": Bolometer, + "sxrc_xy2": SXRcamera, + "sxr_spd": SXRcamera, + "blom_xy1": Bolometer, } -INSTRUMENTS = ["smmh1", "nirh1", "xrcs", "sxr_diode_1", "efit", "brems"] +INSTRUMENTS: list = [ + "smmh", + "smmh1", + "nirh1", + "xrcs", + "blom_xy1", + "sxrc_xy1", + "sxrc_xy2", + "efit", + "sxr_spd", +] +REVISIONS: dict = {instr: 0 for instr in INSTRUMENTS} FIG_PATH = f"/home/{getpass.getuser()}/figures/Indica/load_modelling_examples/" plt.ion() @@ -201,9 +216,9 @@ def initialize_diagnostic_models( def example_run( dt: float = 0.01, verbose: bool = True, - example: str = "predictive", + example: str = "alsu_11314", all_runs: bool = False, - plot: bool = False, + plot: bool = True, save_fig: bool = False, fig_style="profiles", alpha: float = 1.0, @@ -216,6 +231,7 @@ def example_run( pulse_code, pulse, equil, + equil_run, code, runs, comment, @@ -237,7 +253,8 @@ def example_run( if pulse is not None: print("Reading ST40 data") st40 = ReadST40(pulse, tstart, tend, dt=dt, tree="st40") - st40(instruments=INSTRUMENTS, map_diagnostics=False) + REVISIONS["efit"] = equil_run + st40(instruments=INSTRUMENTS, map_diagnostics=False, revisions=REVISIONS) if equil != code: equilibrium = {run: Equilibrium(st40.raw_data[equil]) for run in runs} @@ -283,7 +300,10 @@ def example_run( models[instrument].set_plasma(_plasma) - _bckc = models[instrument]() + if instrument == "xrcs": + _bckc = models[instrument](moment_analysis=True) + else: + _bckc = models[instrument]() bckc[run][instrument] = deepcopy(_bckc) if plot or save_fig: @@ -501,11 +521,13 @@ def plot_data_bckc_comparison( run_str += f"-{runs[1]}" fig_path = f"{FIG_PATH}{pulse_code}_{time}_{code}_{run_str}/" - norm = {} - norm["brems"] = True - norm["sxr_camera_4"] = True - norm["sxr_diode_1"] = True - norm["xrcs"] = True + norm: dict = {} + norm["xrcs"] = {} + norm["xrcs"]["spectra"] = True + # norm["brems"] = True + # norm["sxr_camera_4"] = True + # norm["sxrc_xy2"] = {} + # norm["sxrc_xy2"]["brightness"] = True y0 = {} y0["nirh1"] = True y0["smmh1"] = True @@ -518,39 +540,38 @@ def plot_data_bckc_comparison( for run in runs: bckc[run]["efit"] = {"wp": plasma[run].wp} - for instrument in st40.raw_data.keys(): - for quantity in bckc[run][instrument].keys(): - print(instrument) - print(f" {quantity}") - run = runs[0] - + instruments = st40.raw_data.keys() + for instrument in instruments: + quantities = bckc[runs[0]][instrument].keys() + for quantity in quantities: if ( quantity not in st40.binned_data[instrument].keys() or quantity not in st40.raw_data[instrument].keys() ): continue - plt.figure() + print(instrument) + print(f" {quantity}") _raw = st40.raw_data[instrument][quantity] _binned = st40.binned_data[instrument][quantity] - _bckc = bckc[run][instrument][quantity] + _bckc = bckc[runs[0]][instrument][quantity] + tslice = slice(_bckc.t.min().values, _bckc.t.max().values) + str_to_add = "" + tslice_raw = tslice + tslice_binned = tslice + if "error" not in _binned.attrs: _binned.attrs["error"] = xr.full_like(_binned, 0.0) if "stdev" not in _binned.attrs: _binned.attrs["stdev"] = xr.full_like(_binned, 0.0) - err = np.sqrt(_binned.error**2 + _binned.stdev**2) err = xr.where(err / _binned.values < 1.0, err, 0.0) - if len(_bckc.dims) > 1: + if len(_binned.dims) > 1: str_to_add = f" @ {time:.3f} s" tslice_binned = _binned.t.sel(t=time, method="nearest") tslice_raw = _raw.t.sel(t=time, method="nearest") - else: - str_to_add = "" - tslice_raw = tslice - tslice_binned = tslice _raw = _raw.sel(t=tslice_raw) _binned = _binned.sel(t=tslice_binned) @@ -562,6 +583,7 @@ def plot_data_bckc_comparison( _binned -= bgnd _raw -= bgnd + plt.figure() _raw.plot( label="Raw", color=COLORS["raw"], @@ -583,109 +605,178 @@ def plot_data_bckc_comparison( markersize=markersize, ) - label = None + label: str = "Model" for run in runs: - if run == runs[-1]: - label = "Model" - _bckc = bckc[run][instrument][quantity].sel(t=tslice_binned) if instrument in norm.keys(): - _bckc = _bckc * _binned.max() / _bckc.max() - - if label is not None: - label += " (scaled)" + if quantity in norm[instrument].keys(): + _bckc = _bckc / _bckc.max() * _binned.max() + if label is not None: + label += " (scaled)" (_bckc).plot( label=label, color=COLORS["bckc"], - linewidth=rcParams["lines.linewidth"] * 1.5, + linewidth=rcParams["lines.linewidth"], alpha=alpha, ) - set_axis_sci() - plt.title(f"{instrument.upper()} {quantity}" + str_to_add) - if instrument in y0.keys(): - plt.ylim( - 0, - ) + del label - if quantity == "spectra": - # TODO: wavelength axis is sorted from max to min... - plt.xlim(_bckc.wavelength.min(), _bckc.wavelength.max()) + set_axis_sci() + plt.title(f"{instrument.upper()} {quantity}" + str_to_add) + if instrument in y0.keys(): + plt.ylim( + 0, + ) + if quantity == "spectra": + plt.xlim(_bckc.wavelength.min(), _bckc.wavelength.max()) - plt.legend() - save_figure(fig_path, f"{instrument}_{quantity}", save_fig=save_fig) + plt.legend() + save_figure(fig_path, f"{instrument}_{quantity}", save_fig=save_fig) def example_params(example: str, all_runs: bool = False): - comment: str - pulse_code: int - pulse: int - equil: str - code: str - tstart: float - tend: float - tplot: float - runs = [f"RUN{run}" for run in (500 + np.arange(61, 77))] - - if example == "predictive": - comment = "Tests using fixed-boundary predictive ASTRA" - pulse_code = 13110009 - pulse = 10009 - equil = "astra" - code = "astra" - if not all_runs: - runs = ["RUN2621"] - tstart = 0.02 - tend = 0.08 - tplot = 0.06 - elif example == "interpretative_10009": - comment = "interpretative ASTRA using HDA/EFIT" - pulse_code = 13110009 - pulse = 10009 - equil = "efit" - code = "astra" - if not all_runs: - runs = ["RUN564"] - tstart = 0.03 - tend = 0.1 - tplot = 0.06 - elif example == "interpretative_9850": - comment = "ASTRA interpretative using HDA/EFIT" - pulse = 9850 - pulse_code = 13109850 - equil = "efit" - code = "astra" - if not all_runs: - runs = ["RUN564"] - tstart = 0.02 - tend = 0.1 - tplot = 0.08 - elif example == "interpretative_9229": - comment = "ASTRA interpretative using HDA/EFIT" - pulse = 9229 - pulse_code = 13109229 - equil = "efit" - code = "astra" - if not all_runs: - runs = ["RUN572"] - tstart = 0.03 - tend = 0.11 - tplot = 0.06 - elif example == "diverted": - comment = "predictive ASTRA using for diverted scenario" - pulse_code = 13000040 - pulse = 10009 - equil = "astra" - code = "astra" - if not all_runs: - runs = ["RUN292"] - tstart = 0.03 - tend = 0.11 - tplot = 0.1 - else: - raise ValueError(f"No parameters for example {example}") - - return pulse_code, pulse, equil, code, runs, comment, tstart, tend, tplot + runs_all: list = [f"RUN{run}" for run in (500 + np.arange(61, 77))] + + params = { + "predictive": dict( + comment="Tests using fixed-boundary predictive ASTRA", + pulse_code=13110009, + pulse=10009, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN2621"], + tstart=0.02, + tend=0.08, + tplot=0.06, + ), + "interpretative_10009": dict( + comment="interpretative ASTRA using HDA/EFIT", + pulse_code=13110009, + pulse=10009, + equil="efit", + equil_run=0, + code="astra", + runs=["RUN564"], + tstart=0.03, + tend=0.1, + tplot=0.06, + ), + "interpretative_9850": dict( + comment="ASTRA interpretative using HDA/EFIT", + pulse=9850, + pulse_code=13109850, + equil="efit", + equil_run=0, + code="astra", + runs=["RUN564"], + tstart=0.02, + tend=0.1, + tplot=0.08, + ), + "interpretative_9229": dict( + comment="ASTRA interpretative using HDA/EFIT", + pulse=9229, + pulse_code=13109229, + equil="efit", + equil_run=0, + code="astra", + runs=["RUN572"], + tstart=0.03, + tend=0.11, + tplot=0.06, + ), + "diverted": dict( + comment="predictive ASTRA using for diverted scenario", + pulse_code=13000040, + pulse=10009, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN292"], + tstart=0.03, + tend=0.11, + tplot=0.1, + ), + "michail_10014": dict( + comment="predictive ASTRA", + pulse_code=36010014, + pulse=10014, + equil="efit", + equil_run=0, + code="astra", + runs=["RUN24"], + tstart=0.03, + tend=0.1, + tplot=0.07, + ), + "aleksei_11228": dict( + comment="ASTRA using TS and invented Ti shapes", + pulse=11228, + pulse_code=13011228, + equil="efit", # "astra" + equil_run=1, + code="astra", + runs=["RUN610", "RUN611", "RUN612"], + tstart=0.03, + tend=0.11, + tplot=0.08, + ), + "alsu_11312": dict( + comment="ASTRA using TS and peaked Ti scaled to CXRS", + pulse=11312, + pulse_code=33011312, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN21"], + tstart=0.065, + tend=0.095, + tplot=0.075, + ), + "alsu_11314": dict( + comment="ASTRA using TS and peaked Ti scaled to CXRS", + pulse=11314, + pulse_code=33011314, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN12"], + tstart=0.065, + tend=0.095, + tplot=0.075, + ), + "alsu_11317": dict( + comment="ASTRA using TS and peaked Ti scaled to CXRS", + pulse=11317, + pulse_code=33011317, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN9"], + tstart=0.065, + tend=0.095, + tplot=0.075, + ), + } + + _params = params[example] + if all_runs: + _params["runs"] = runs_all + + return ( + _params["pulse_code"], + _params["pulse"], + _params["equil"], + _params["equil_run"], + _params["code"], + _params["runs"], + _params["comment"], + _params["tstart"], + _params["tend"], + _params["tplot"], + ) # def add_gacode_data( @@ -723,4 +814,6 @@ def example_params(example: str, all_runs: bool = False): if __name__ == "__main__": - example_run() + plt.ioff() + example_run(example="alsu_11314") + plt.show() diff --git a/indica/workflows/run_tomo_1d.py b/indica/workflows/run_tomo_1d.py index c4a73fa3..cc5f02b9 100644 --- a/indica/workflows/run_tomo_1d.py +++ b/indica/workflows/run_tomo_1d.py @@ -95,13 +95,15 @@ def phantom_examples( return inverted_emissivity, data_tomo, bckc_tomo -def experimental_examples( +def sxrc_example( + pulse: int = None, instrument: str = "sxrc_xy2", reg_level_guess: float = 0.5, phantom_data: bool = True, plot: bool = True, ): - pulse = PULSES[instrument] + if pulse is None: + pulse = PULSES[instrument] tstart = 0.02 tend = 0.1 @@ -112,17 +114,101 @@ def experimental_examples( quantity = list(st40.binned_data[instrument])[0] los_transform = st40.binned_data[instrument][quantity].transform los_transform.set_equilibrium(equilibrium, force=True) - if instrument == "pi": - BREMS_MODEL.set_los_transform(los_transform) - attrs = st40.binned_data[instrument]["spectra"].attrs - background, brightness = BREMS_MODEL.integrate_spectra( - st40.binned_data[instrument]["spectra"] + model = SXR_MODEL + model.set_los_transform(los_transform) + + if phantom_data: + plasma = example_plasma( + pulse, + tstart=tstart, + tend=tend, + dt=dt, ) - background.attrs = attrs - brightness.attrs = attrs - st40.binned_data[instrument]["background"] = background - st40.binned_data[instrument]["brightness"] = brightness - model = SXR_MODEL + plasma.build_atomic_data() + plasma.set_equilibrium(equilibrium) + model.set_plasma(plasma) + bckc = model() + emissivity = model.emissivity + brightness = bckc["brightness"] + else: + emissivity = None + brightness = st40.binned_data[instrument]["brightness"] + + z = los_transform.z + R = los_transform.R + dl = los_transform.dl + + data_t0 = brightness.isel(t=0).data + has_data = np.logical_not(np.isnan(brightness.isel(t=0).data)) & (data_t0 > 0) + rho_equil = equilibrium.rho.interp(t=brightness.t) + input_dict = dict( + brightness=brightness.data, + dl=dl, + t=brightness.t.data, + R=R, + z=z, + rho_equil=dict( + R=rho_equil.R.data, + z=rho_equil.z.data, + t=rho_equil.t.data, + rho=rho_equil.data, + ), + has_data=has_data, + debug=False, + ) + if emissivity is not None: + input_dict["emissivity"] = emissivity + + tomo = tomo_1D.SXR_tomography(input_dict, reg_level_guess=reg_level_guess) + tomo() + + if plot: + model.los_transform.plot() + tomo.show_reconstruction() + + inverted_emissivity = DataArray( + tomo.emiss, coords=[("t", tomo.tvec), ("rho_poloidal", tomo.rho_grid_centers)] + ) + inverted_error = DataArray( + tomo.emiss_err, + coords=[("t", tomo.tvec), ("rho_poloidal", tomo.rho_grid_centers)], + ) + inverted_emissivity.attrs["error"] = inverted_error + data_tomo = brightness + bckc_tomo = DataArray(tomo.backprojection, coords=data_tomo.coords) + + return inverted_emissivity, data_tomo, bckc_tomo + + +def pi_bremsstrahlung_example( + pulse: int = None, + instrument: str = "pi", + reg_level_guess: float = 0.5, + phantom_data: bool = True, + plot: bool = True, +): + if pulse is None: + pulse = PULSES[instrument] + + tstart = 0.02 + tend = 0.1 + dt = 0.01 + st40 = ReadST40(pulse, tstart, tend, dt=dt) + st40(instruments=[instrument, "efit"]) + equilibrium = Equilibrium(st40.raw_data["efit"]) + quantity = list(st40.binned_data[instrument])[0] + los_transform = st40.binned_data[instrument][quantity].transform + los_transform.set_equilibrium(equilibrium, force=True) + model = BREMS_MODEL + model.set_los_transform(los_transform) + attrs = st40.binned_data[instrument]["spectra"].attrs + background, brightness = BREMS_MODEL.integrate_spectra( + st40.binned_data[instrument]["spectra"] + ) + background.attrs = attrs + brightness.attrs = attrs + st40.binned_data[instrument]["background"] = background + st40.binned_data[instrument]["brightness"] = brightness if phantom_data: plasma = example_plasma( @@ -198,7 +284,7 @@ def sxrc_xy( reg_level_guess: float = 0.5, input_dict: dict = None, channels=slice(0, 15), - save_fig: bool = True, + save_fig: bool = False, instrument="sxrc_xy2", ): if input_dict is None: @@ -212,7 +298,7 @@ def sxrc_xy( equil = st40.equilibrium rho, theta = data.transform.convert_to_rho_theta(t=data.t) - rho = rho.sel(channel=channels) + # rho = rho.sel(channel=channels) brightness = data @@ -341,75 +427,3 @@ def old_camera( tomo.show_reconstruction() return input_dict - - -def fake_data( - pulse: int = 9229, - plasma=None, - model=None, - bckc=None, - debug=False, - exclude_bad_points=True, - plot=True, - reg_level_guess: float = 0.5, - nchannels=12, - input_dict: dict = None, -): - from indica.models.sxr_camera import example_run as sxr_example - - if input_dict is None: - if plasma is None or model is None or bckc is None: - plasma, model, bckc = sxr_example(pulse=pulse, nchannels=nchannels) - - tstart = plasma.t.min().values - tend = plasma.t.max().values - st40 = ReadST40(pulse, tstart, tend) - st40(instruments=["sxrc_xy2", "sxr_camera_4"]) - - try: - data = st40.binned_data["sxrc_xy2"]["brightness"] - except KeyError: - data = st40.binned_data["sxr_camera_4"]["brightness"] - model.set_los_transform(data.transform) - model.los_transform.set_equilibrium(plasma.equilibrium, force=True) - model() - - z = model.los_transform.z - R = model.los_transform.R - dl = model.los_transform.dl - - brightness = model.bckc["brightness"] - data_t0 = brightness.isel(t=0).data - if exclude_bad_points: - has_data = np.logical_not(np.isnan(data_t0)) & (data_t0 >= 1.0e3) - else: - has_data = np.logical_not(np.isnan(data_t0)) - - rho_equil = plasma.equilibrium.rho.interp(t=brightness.t) - input_dict = dict( - emissivity=model.emissivity, - brightness=brightness.data, - dl=dl, - t=brightness.t.data, - R=R, - z=z, - rho_equil=dict( - R=rho_equil.R.data, - z=rho_equil.z.data, - t=rho_equil.t.data, - rho=rho_equil.data, - ), - debug=debug, - has_data=has_data, - ) - - # return input_dict - - tomo = tomo_1D.SXR_tomography(input_dict, reg_level_guess=reg_level_guess) - - tomo() - - if plot: - tomo.show_reconstruction() - - return plasma diff --git a/indica/writers/bda_tree.py b/indica/writers/bda_tree.py index 25453c96..0a3bc894 100644 --- a/indica/writers/bda_tree.py +++ b/indica/writers/bda_tree.py @@ -238,7 +238,7 @@ def check_analysis_run( which_run, ): # Checker function to see if data already exists in a run - IP_address_smaug = "192.168.1.7:8000" + IP_address_smaug = "smaug" conn = Connection(IP_address_smaug) conn.openTree("BDA", pulseNo) diff --git a/tests/unit/readers/test_abstract_reader_pytest.py b/tests/unit/readers/test_abstract_reader_pytest.py index 894e2c92..cccd6aba 100644 --- a/tests/unit/readers/test_abstract_reader_pytest.py +++ b/tests/unit/readers/test_abstract_reader_pytest.py @@ -36,7 +36,7 @@ class Reader(DataReader): "radiation": "get_radiation", "helike_spectroscopy": "get_helike_spectroscopy", "interferometry": "get_interferometry", - "filters": "get_diode_filters", + "diode_filters": "get_diode_filters", } def __init__( @@ -368,6 +368,7 @@ def _get_diode_filters( results["location"] = np.array([[1.0, 2.0, 3.0]] * results["length"]) results["direction"] = np.array([[1.0, 2.0, 3.0]] * results["length"]) + results["labels"] = np.array(["label"] * results["length"]) for quantity in quantities: results[quantity] = np.random.uniform(0, 1.0e6, (nt, results["length"])) @@ -494,7 +495,7 @@ def test_get_helike_spectroscopy(): def test_get_diode_filters(): _test_get_methods( - instrument="filters", + instrument="diode_filters", nsamples=10, )