Skip to content

Commit

Permalink
Marcosertoli/read spd (#275)
Browse files Browse the repository at this point in the history 
* Fixed run_tomo_1d minor inconsistencies. Workflows should be refactored to make the code more modular

* Refactored readers for new LINES and SXR_SPD diagnostics

* Refactored abstractreader test for new diode_filters method

* Separated PI and SXRC examples

* Deleted commented-out line, generalised test call.

* Corrected bug for plotting TS modelled results vs. rho

* Simplified general call or read_st40.py wrapper & added example read and plot SXR_SPD

* Empty commit

---------

Co-authored-by: Marco Sertoli <[email protected]>
  • Loading branch information
@marcosertoli
marcosertoli and marcosertoli authored Sep 19, 2023
1 parent dd13acd commit 677dcd1
Show file tree
Hide file tree
Showing 6 changed files with 194 additions and 172 deletions.
6 changes: 4 additions & 2 deletions indica/models/thomson_scattering.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,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()
Expand All @@ -176,7 +177,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()
Expand Down
3 changes: 2 additions & 1 deletion indica/readers/abstractreader.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"])))
Expand All @@ -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

Expand Down
76 changes: 42 additions & 34 deletions indica/readers/read_st40.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,25 +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,
"blom_xy1": 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)},
Expand All @@ -39,6 +39,7 @@
"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)},
Expand Down Expand Up @@ -177,6 +178,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!"
Expand Down Expand Up @@ -298,8 +300,8 @@ def plot_profile(

def __call__(
self,
instruments: list = None,
revisions: dict = None,
instruments: list = [],
revisions: list = None,
map_raw: bool = False,
tstart: float = None,
tend: float = None,
Expand All @@ -311,10 +313,11 @@ def __call__(
debug: bool = False,
):
self.debug = debug
if instruments is None:
instruments = list(REVISIONS.keys())

if len(instruments) == 0:
instruments = INSTRUMENTS
if revisions is None:
revisions = REVISIONS
revisions = [0] * len(instruments)
if tstart is None:
tstart = self.tstart
if tend is None:
Expand All @@ -324,13 +327,13 @@ def __call__(

self.reset_data()
self.get_equilibrium(R_shift=R_shift)
for instrument in instruments:
for i, instrument in enumerate(instruments):
print(f"Reading {instrument}")
if debug:
self.get_raw_data("", instrument, revisions[instrument])
self.get_raw_data("", instrument, revisions[i])
else:
try:
self.get_raw_data("", instrument, revisions[instrument])
self.get_raw_data("", instrument, revisions[i])
except Exception as e:
print(f"Error reading {instrument}: {e}")

Expand Down Expand Up @@ -363,12 +366,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"]
90 changes: 43 additions & 47 deletions indica/readers/st40reader.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,12 +74,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",
Expand All @@ -98,12 +99,9 @@ class ST40Reader(DataReader):
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",
Expand Down Expand Up @@ -158,11 +156,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:",
Expand Down Expand Up @@ -913,71 +914,51 @@ 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,
"machine_dims": self.MACHINE_DIMS,
}
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])

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,
)
Expand All @@ -996,6 +977,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]
Expand All @@ -1014,6 +999,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])
Expand Down
Loading

0 comments on commit 677dcd1

Please sign in to comment.