diff --git a/data/x2s5823/_info.yaml b/data/x2s5823/_info.yaml index 3a30f42..f17a315 100644 --- a/data/x2s5823/_info.yaml +++ b/data/x2s5823/_info.yaml @@ -49,6 +49,7 @@ probe-info: c2c: 5.6 # mm gauge-diam: 3.05 # mm gauge-c2c: 4 #mm + dist-uncert: 2 #mm gauges: - 1 @@ -60,6 +61,7 @@ probe-info: - type: "scope" config: "eProbe-Scope.txt" data: "eProbe-Scope.csv" + time-uncert: 0 # s trigger: # Redundant? type: "channel" diff --git a/data/x2s5824/_info.yaml b/data/x2s5824/_info.yaml index bc3cb44..c4da523 100644 --- a/data/x2s5824/_info.yaml +++ b/data/x2s5824/_info.yaml @@ -49,6 +49,7 @@ probe-info: c2c: 5.6 # mm gauge-diam: 0.8 # mm gauge-c2c: 1.8 #mm + dist-uncert: 0.025 #mm gauges: - 1 @@ -61,6 +62,7 @@ probe-info: - type: "scope" config: "eProbe-Scope.txt" data: "eProbe-Scope.csv" + time-uncert: 0 # s trigger: # Redundant? type: "channel" diff --git a/data/x2s5827/_info.yaml b/data/x2s5827/_info.yaml index 5e59d9d..8895a08 100644 --- a/data/x2s5827/_info.yaml +++ b/data/x2s5827/_info.yaml @@ -52,6 +52,7 @@ probe-info: c2c: 5.6 # mm gauge-diam: 0.8 # mm gauge-c2c: 1.8 #mm + dist-uncert: 0.025 #mm gauges: - 1 @@ -64,6 +65,7 @@ probe-info: - type: "scope" config: "eProbe-Scope.txt" data: "eProbe-Scope.csv" + time-uncert: 0 # s trigger: # Redundant? type: "channel" diff --git a/data/x2s5829/_info.yaml b/data/x2s5829/_info.yaml index 9f678fa..7a530a2 100644 --- a/data/x2s5829/_info.yaml +++ b/data/x2s5829/_info.yaml @@ -54,6 +54,7 @@ probe-info: c2c: 5.6 # mm gauge-diam: 0.8 # mm gauge-c2c: 1.8 #mm + dist-uncert: 0.025 #mm gauges: - 1 @@ -66,6 +67,7 @@ probe-info: - type: "scope" config: "eProbe-Scope.txt" data: "eProbe-Scope.csv" + time-uncert: 0 # s trigger: # Redundant? type: "channel" diff --git a/data/x2s5830/_info.yaml b/data/x2s5830/_info.yaml index 5f37ad5..6d33d80 100644 --- a/data/x2s5830/_info.yaml +++ b/data/x2s5830/_info.yaml @@ -52,6 +52,7 @@ probe-info: c2c: 5.6 # mm gauge-diam: 0.8 # mm gauge-c2c: 1.8 #mm + dist-uncert: 0.025 #mm gauges: - 1 @@ -64,6 +65,7 @@ probe-info: - type: "scope" config: "eProbe-Scope.txt" data: "eProbe-Scope.csv" + time-uncert: 0 # s trigger: # Redundant? type: "channel" diff --git a/data/x2s5831/_info.yaml b/data/x2s5831/_info.yaml index f13b577..9536d70 100644 --- a/data/x2s5831/_info.yaml +++ b/data/x2s5831/_info.yaml @@ -55,11 +55,13 @@ probe-info: c2c: 5.6 # mm gauge-diam: 0.8 # mm gauge-c2c: 1.8 #mm + dist-uncert: 0.025 #mm data-records: - type: "scope" config: "eProbe-Scope.txt" data: "eProbe-Scope.csv" + time-uncert: 0 # s trigger: # Redundant? type: "channel" diff --git a/data/x2s5832/_info.yaml b/data/x2s5832/_info.yaml index 0df85c8..0b5bbed 100644 --- a/data/x2s5832/_info.yaml +++ b/data/x2s5832/_info.yaml @@ -58,11 +58,13 @@ probe-info: c2c: 5.6 # mm gauge-diam: 0.8 # mm gauge-c2c: 1.8 #mm + dist-uncert: 0.025 #mm data-records: - type: "scope" config: "eProbe-Scope.txt" data: "eProbe-Scope.csv" + time-uncert: 0 # s trigger: # Redundant? type: "channel" diff --git a/main.py b/main.py index e04636e..61dcbd6 100644 --- a/main.py +++ b/main.py @@ -39,6 +39,21 @@ data_to_load = [ "x2s5832" ] +# ==== Uncerts ==== +# Taken from DOI: 10.1007/s00193-017-0763-3 (Implementation of a state-to-state analytical framework for the calculation of expansion tube flow properties) + +UNCERTS = TUNNEL_INFO["uncertainties"] + +def deltaX(delta_x_1: float, delta_x_2: float): + return np.sqrt(np.pow(delta_x_1, 2) + np.pow(delta_x_2, 2)) + +def deltaT(delta_t_1: float, delta_t_2: float, delta_t_sr: float): + return np.sqrt(np.pow(delta_t_1, 2) + np.pow(delta_t_2, 2) + np.pow(delta_t_sr, 2)) + +def deltaVs(V: float, dx: float, dt: float, delta_x: tuple[float, float], delta_t: tuple[float, float, float]): + return V * np.sqrt(np.pow(deltaX(*delta_x) / dx, 2) + np.pow(deltaT(*delta_t) / dt, 2)) + + # ==== Data Loading & Processing ==== def load_data(data_path: str, data={}) -> dict: data_info_path = data_path + DATA_INFO @@ -114,6 +129,7 @@ def load_data(data_path: str, data={}) -> dict: "x2": None, "probes": None, # This may be x2 but may not - ie a scope was used "trigger_index": None, + "probe_uncert": None, #s }, "data": { "x2": {}, # Only pop channels with a voltage scale in ./tunnel-info.yaml @@ -196,6 +212,7 @@ def load_data(data_path: str, data={}) -> dict: if "scope" in data_locs: data[x2_shot]["data"]["probes"] = [data[x2_shot]["data"]["scope"][1], data[x2_shot]["data"]["scope"][2]] data[x2_shot]["time"]["probes"] = data[x2_shot]["time"]["scope"] + data[x2_shot]["time"]["probe_uncert"] = scope_data_info["time-uncert"] # Find Shock Times @@ -215,7 +232,7 @@ def load_data(data_path: str, data={}) -> dict: first_value, first_value_uncertainty, _, _ = canny_shock_finder(x2_time_us, refData, sigma=sigma, post_suppression_threshold=post_sup_thresh, plot=False, print_func=None) shock_point = np.where(x2_time_us >= first_value)[0][0] # [BUG] Seems to give n+1 - data[x2_shot]["shock-point"][ref] = shock_point, first_value + data[x2_shot]["shock-point"][ref] = shock_point, first_value, first_value_uncertainty @@ -250,7 +267,7 @@ def load_data(data_path: str, data={}) -> dict: shock_point = np.where(scope_time >= first_value)[0][0] # [BUG] Seems to give n+1 - data[x2_shot]["shock-point"][f"{probe}-g1"] = shock_point, first_value + data[x2_shot]["shock-point"][f"{probe}-g1"] = shock_point, first_value, first_value_uncertainty if 2 in dataInfo["probe-info"]["gauges"]: # Do the same for G2 @@ -264,7 +281,7 @@ def load_data(data_path: str, data={}) -> dict: raise ValueError(f"{probe}-g2 not detected") shock_point = np.where(scope_time >= first_value)[0][0] # [BUG] Seems to give n+1 - data[x2_shot]["shock-point"][f"{probe}-g2"] = shock_point, first_value + data[x2_shot]["shock-point"][f"{probe}-g2"] = shock_point, first_value, first_value_uncertainty @@ -276,11 +293,19 @@ def load_data(data_path: str, data={}) -> dict: if i == 0: continue p1_time = data[x2_shot]["shock-point"][refProbe][1] / 1e6 # Convert to seconds p2_time = data[x2_shot]["shock-point"][dataInfo["pcb-refs"][i-1]][1] / 1e6 # Convert to seconds - p2p_dist = (TUNNEL_INFO["distance"][refProbe] - TUNNEL_INFO["distance"][dataInfo["pcb-refs"][i-1]]) / 1000 # convert to m - probe_velocity = p2p_dist / abs(p2_time - p1_time) # m/s - print(f"{refProbe}-{dataInfo["pcb-refs"][i-1]} Measured a shock speed of {probe_velocity:.2f} m/s ({probe_velocity/1000:.2f} km/s)") - data[x2_shot]["shock-speed"][f"{refProbe}-{dataInfo["pcb-refs"][i-1]}"] = probe_velocity, True # Speed, Ref + p2p_dist = (TUNNEL_INFO["distance"][refProbe] - TUNNEL_INFO["distance"][dataInfo["pcb-refs"][i-1]]) / 1000 # convert to m + p2p_time = abs(p2_time - p1_time) + + probe_velocity = p2p_dist / p2p_time # m/s + + p1_time_uncert = data[x2_shot]["shock-point"][refProbe][2] / 1e6 # Convert to seconds + p2_time_uncert = data[x2_shot]["shock-point"][dataInfo["pcb-refs"][i-1]][2] / 1e6 # Convert to seconds + + uncert = deltaVs(probe_velocity, p2p_dist, p2p_time, (UNCERTS["probe-dist"][refProbe], UNCERTS["probe-dist"][dataInfo["pcb-refs"][i-1]]), (p1_time_uncert, p2_time_uncert, UNCERTS["time"]["x2-daq"])) + + print(f"{refProbe}-{dataInfo["pcb-refs"][i-1]} Measured a shock speed of {probe_velocity:.2f} +/- {uncert:.2f} m/s ({probe_velocity/1000:.2f} +/- {uncert/1000:.2f} km/s [{uncert/probe_velocity * 100 :.2f}%])") + data[x2_shot]["shock-speed"][f"{refProbe}-{dataInfo["pcb-refs"][i-1]}"] = probe_velocity, uncert, True # Speed, Ref print() for probe in dataInfo["probe-info"]["locations"]: @@ -288,11 +313,17 @@ def load_data(data_path: str, data={}) -> dict: g1_time = data[x2_shot]["shock-point"][f"{probe}-g1"][1] / 1e6 # Convert to seconds g2_time = data[x2_shot]["shock-point"][f"{probe}-g2"][1] / 1e6 # Convert to seconds c2c_dist = dataInfo["probe-info"]["c2c"] / 1000 # convert to m + c2c_time = abs(g2_time - g1_time) - probe_velocity = c2c_dist / abs(g2_time - g1_time) # m/s + probe_velocity = c2c_dist / c2c_time # m/s - print(f"{probe} Measured a shock speed of {probe_velocity:.2f} m/s ({probe_velocity/1000:.2f} km/s)") - data[x2_shot]["shock-speed"][probe] = probe_velocity, False # Speed, Ref # m/s + g1_time_uncert = data[x2_shot]["shock-point"][f"{probe}-g1"][2] / 1e6 # Convert to seconds + g2_time_uncert = data[x2_shot]["shock-point"][f"{probe}-g2"][2] / 1e6 # Convert to seconds + + uncert = deltaVs(probe_velocity, p2p_dist, p2p_time, (0.05/1000, 0.05/1000), (g1_time_uncert, g2_time_uncert, data[x2_shot]["time"]["probe_uncert"])) + + print(f"{probe} Measured a shock speed of {probe_velocity:.2f} +/- {uncert:.2f} m/s ({probe_velocity/1000:.2f} +/- {uncert/1000:.2f} km/s)") + data[x2_shot]["shock-speed"][probe] = probe_velocity, uncert, False # Speed, Ref # m/s else: print(f"Unable to calculate probe velocity, only have one probe: {f"{probe}-g2" if f"{probe}-g2" in data[x2_shot]["shock-point"] else f"{probe}-g1"}") @@ -300,23 +331,37 @@ def load_data(data_path: str, data={}) -> dict: for i in range(len(dataInfo["probe-info"]["locations"]) - 1): probe_locs = dataInfo["probe-info"]["locations"] - p2p = (TUNNEL_INFO["distance"][probe_locs[1]] - TUNNEL_INFO["distance"][probe_locs[0]]) / 1000 # convert to m + p2p_dist = (TUNNEL_INFO["distance"][probe_locs[1]] - TUNNEL_INFO["distance"][probe_locs[0]]) / 1000 # convert to m if f"{probe_locs[i]}-g1" in data[x2_shot]["shock-point"] and f"{probe_locs[i+1]}-g1" in data[x2_shot]["shock-point"]: p1_g1_time = data[x2_shot]["shock-point"][f"{probe_locs[i]}-g1"][1] / 1e6 # Convert to seconds p2_g1_time = data[x2_shot]["shock-point"][f"{probe_locs[i+1]}-g1"][1] / 1e6 # Convert to seconds - p2p_1 = p2p / abs(p2_g1_time - p1_g1_time) # m/s + + p2p_time = abs(p2_g1_time - p1_g1_time) + p2p_1 = p2p_dist / p2p_time # m/s - print(f"{probe_locs[i]}-{probe_locs[i + 1]} - G1 - Measured a shock speed of {p2p_1:.2f} m/s ({p2p_1/1000:.2f} km/s)") - data[x2_shot]["shock-speed"][f"{probe_locs[i]}-{probe_locs[i + 1]}-g1"] = p2p_1, False # Speed, Ref + p1_time_uncert = data[x2_shot]["shock-point"][f"{probe_locs[i]}-g1"][2] / 1e6 # Convert to seconds + p2_time_uncert = data[x2_shot]["shock-point"][f"{probe_locs[i+1]}-g1"][2] / 1e6 # Convert to seconds + + uncert = deltaVs(p2p_1, p2p_dist, p2p_time, (UNCERTS["probe-dist"][probe_locs[i]], UNCERTS["probe-dist"][probe_locs[i+1]]), (p1_time_uncert, p2_time_uncert, data[x2_shot]["time"]["probe_uncert"])) + + print(f"{probe_locs[i]}-{probe_locs[i + 1]} - G1 - Measured a shock speed of {p2p_1:.2f} +/- {uncert:.2f} m/s ({p2p_1/1000:.2f} +/- {uncert/1000:.2f} [{uncert/p2p_1 * 100:.2f}%] km/s)") + data[x2_shot]["shock-speed"][f"{probe_locs[i]}-{probe_locs[i + 1]}-g1"] = p2p_1, uncert, False # Speed, Ref if f"{probe_locs[i]}-g2" in data[x2_shot]["shock-point"] and f"{probe_locs[i+1]}-g2" in data[x2_shot]["shock-point"]: p1_g2_time = data[x2_shot]["shock-point"][f"{probe_locs[i]}-g2"][1] / 1e6 # Convert to seconds p2_g2_time = data[x2_shot]["shock-point"][f"{probe_locs[i+1]}-g2"][1] / 1e6 # Convert to seconds - p2p_2 = p2p / abs(p2_g2_time - p1_g2_time) # m/s + + p2p_time = abs(p2_g2_time - p1_g2_time) + p2p_2 = p2p_dist / p2p_time # m/s - print(f"{probe_locs[i]}-{probe_locs[i + 1]} - G2 - Measured a shock speed of {p2p_2:.2f} m/s ({p2p_2/1000:.2f} km/s)") - data[x2_shot]["shock-speed"][f"{probe_locs[i]}-{probe_locs[i + 1]}-g2"] = p2p_2, False # Speed, Ref + p1_time_uncert = data[x2_shot]["shock-point"][f"{probe_locs[i]}-g2"][2] / 1e6 # Convert to seconds + p2_time_uncert = data[x2_shot]["shock-point"][f"{probe_locs[i+1]}-g2"][2] / 1e6 # Convert to seconds + + uncert = deltaVs(p2p_2, p2p_dist, p2p_time, (UNCERTS["probe-dist"][probe_locs[i]], UNCERTS["probe-dist"][probe_locs[i+1]]), (p1_time_uncert, p2_time_uncert, data[x2_shot]["time"]["probe_uncert"])) + + print(f"{probe_locs[i]}-{probe_locs[i + 1]} - G2 - Measured a shock speed of {p2p_2:.2f} +/- {uncert:.2f} m/s ({p2p_2/1000:.2f} +/- {uncert/1000:.2f} [{uncert/p2p_2 * 100:.2f}%] km/s)") + data[x2_shot]["shock-speed"][f"{probe_locs[i]}-{probe_locs[i + 1]}-g2"] = p2p_2, uncert, False # Speed, Ref print() # Return the data & the successfully loaded data keys @@ -344,6 +389,11 @@ def genGraph(gData: dict, showPlot: bool = True, doLimits: bool = True, forcePlo "plots": [] } + #if forcePlots or not doLimits: graphData["title"] += "\n" + #if forcePlots: graphData["title"] += "(All Data Shown)" + #if not doLimits: graphData["title"] += () + "Full Re" + + lims = [] for label in gData["info"]["pcb-refs"]: # + ["trigbox"]: @@ -398,13 +448,14 @@ def genGraph(gData: dict, showPlot: bool = True, doLimits: bool = True, forcePlo for shock_speed_loc in gData["shock-speed"]: probeText += "\n" #probeText += "(Reference) " if gData["shock-speed"][shock_speed_loc][1] else "" - probeText += f"{shock_speed_loc} - {gData["shock-speed"][shock_speed_loc][0]/1000:.2f} km/s" + probeText += f"{shock_speed_loc} - {gData["shock-speed"][shock_speed_loc][0]/1000:.2f} $\\pm${gData["shock-speed"][shock_speed_loc][1]/1000:.2f} [{gData["shock-speed"][shock_speed_loc][1]/gData["shock-speed"][shock_speed_loc][0]*100:.2f}%] km/s" graphData["plots"].append({ "type": "text", "text": f"Measured Shock Speeds{probeText}", "align": ("top", "right"), - "x": 0.94, "y": 0.94 + "x": 0.94 if len(gData["info"]["probe-info"]["locations"]) < 3 else 0.90, + "y": 0.94 }) if doLimits and len(lims) > 1: @@ -420,7 +471,7 @@ print("Graphing Data") for shot in loaded_data: #if shot != loaded_data[-2]: continue genGraph(data[shot], showPlot=False) - genGraph(data[shot], showPlot=False, doLimits=False, forcePlots=True) + #genGraph(data[shot], showPlot=False, doLimits=False, forcePlots=True) # This forces matplotlib to hang until I tell it to close all windows pltKeyClose() diff --git a/tunnel-info.yaml b/tunnel-info.yaml index 0e01073..2994ad7 100644 --- a/tunnel-info.yaml +++ b/tunnel-info.yaml @@ -35,3 +35,37 @@ volt-scale: trigbox: 0.001 #V / mV trigbox_delay: 0.001 #V / mV + + +uncertainties: + probe-dist: + sd1: 0.002 # +/- m + sd2: 0.002 # +/- m + sd3: 0.002 # +/- m + + st1: 0.002 # +/- m + st2: 0.002 # +/- m + st3: 0.002 # +/- m + + at1: 0.002 # +/- m + at2: 0.002 # +/- m + at3: 0.002 # +/- m + + at4: 0.002 # +/- m + at5: 0.002 # +/- m + at6: 0.002 # +/- m + + # [TODO] This could be better + st1-g1: 0.002 # +/- m + st1-g2: 0.002 # +/- m + + st2-g1: 0.002 # +/- m + st2-g2: 0.002 # +/- m + + st3-g1: 0.002 # +/- m + st3-g2: 0.002 # +/- m + + + time: + x2-daq: 0 +