deploy: 8d46f82

.buildinfo

@@ -0,0 +1,4 @@
+# Sphinx build info version 1
+# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
+config: 2af2662d93830d45f081373fcbe5308f
+tags: 645f666f9bcd5a90fca523b33c5a78b7

_downloads/041e21a5dcf4665a558ec24266dc2579/plot_reconstruct_velocity.py

@@ -0,0 +1,67 @@
+"""
+==================================
+Probes and velocity reconstruction
+==================================
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from pyudv.geometry import Probe, sketch_probes, reconstruct_velocity
+
+
+def U(z):
+    u = 5 * (5 - z) ** 2
+    v = u / 10
+    U = np.array([u, v])
+    return U
+
+
+# %%
+# Define probes and plot them
+# ===========================
+
+# define probes
+r = np.linspace(0, 5, 100)
+alpha1, alpha2 = -120, -70  # deg
+O1, O2 = np.array([1, 8]), np.array([-1, 7])
+probe1_pars = [r, alpha1, [0, O1]]
+probe2_pars = [r, alpha2, [0, O2]]
+#
+probe1 = Probe(*probe1_pars)
+probe2 = Probe(*probe2_pars)
+
+fig, ax = plt.subplots(1, 1, layout="constrained")
+sketch_probes(
+    [probe1, probe2],
+    combinations=[[0, 1]],
+    combination_colors=["k"],
+    ax=ax,
+)
+ax.set_xlabel("x")
+ax.set_ylabel("y")
+plt.show()
+
+# %%
+# Create fake signal
+# ==================
+
+u1 = U(probe1.z).T @ probe1.unit_vec
+u2 = U(probe2.z).T @ probe2.unit_vec
+
+# %%
+# Velocity reconstruction
+# =======================
+
+U_rec, z_interp, X, dx_1, dx_2 = reconstruct_velocity(u1, u2, probe1, probe2)
+U_th = U(z_interp)
+#
+fig, axarr = plt.subplots(1, 2, layout="constrained", sharey=True)
+for ax, u_th, u_rec in zip(axarr, U_th, U_rec):
+    ax.plot(u_th, z_interp, ".", label="base")
+    ax.plot(u_rec, z_interp, ".", label="reconstructed")
+    ax.legend()
+axarr[0].set_xlabel("u")
+axarr[1].set_xlabel("v")
+axarr[0].set_ylabel("z")
+plt.show()

_downloads/072ebe5500351177edd3cc97b2ba54e0/plot_direct_inversion.py

@@ -0,0 +1,182 @@
+"""
+==============================================
+Concentration inference from amplitude signals
+==============================================
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import pyudv.amplitude.direct_models as DM
+from pyudv.amplitude.inversion import explicit_inversion
+from pyudv.amplitude.sediment_acoustic_models import quartz_sand as quartz
+
+
+def C_to_phi(C, rho=2.65e3):
+    return C / rho
+
+
+def phi_to_C(phi, rho=2.65e3):
+    return phi * rho
+
+
+fig_width_small = 11.25
+
+# %%
+# Define the parameters for the direct model
+# ==========================================
+
+d = 100e-6  # grain mean diameter [m]
+rho = 2.65  # grain density [g/cm3]
+rho = rho * 1e-3 / 1e-6  # grain density [kg/m3]
+
+F = 2e6  # frequency [Hz]
+T = 20  # Temperature [Celsius degrees]
+a_w = DM.alpha_w(F, T)  # water absorption, [m-1]
+k = 2 * np.pi * F / DM.sound_velocity(T)  # wavenumber of the wave
+
+Ks, Kt = 1, 1  # sediment and transducer constants
+Xi = quartz.Xi(k, d / 2)  # sediment attenuation constant
+rn = 0.05  # near field distance [cm]
+r = np.linspace(0.001, 12, 330)  # radial coordinates [cm]
+
+# Volumic fraction of grain
+Phi = np.array([0.0001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1])
+
+
+# %%
+# Influence of the imposed point in the integration
+# =================================================
+#
+# without near field function
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+psi = DM.near_field_theoretical(r, rn) * 0 + 1  # near field function
+# defining a constant grain concentration profile
+C = Phi[:, None] * rho * (r[None, :] * 0 + 1)
+MSV = DM.create_MSvoltage(C, r, Xi, a_w, Ks, Kt, psi)
+#
+indexes = [2, 50, 150, 250, 320]
+
+color = []
+fig, ax = plt.subplots(
+    1, 1, figsize=(fig_width_small, fig_width_small), constrained_layout=True
+)
+for i, phi in enumerate(Phi):
+    if i == 0:
+        a = plt.plot(C[i, :], r, lw=4, alpha=0.3, label="imposed")
+    else:
+        a = plt.plot(C[i, :], r, lw=4, alpha=0.3)
+    color.append(a[0].get_color())
+for j, index in enumerate(indexes):
+    C0 = C[:, index]
+    r0 = r[index]
+    C_inferred = explicit_inversion(MSV, r, Xi, a_w, psi, C0, r0)
+    for i, phi in enumerate(Phi):
+        if i == 0:
+            plt.plot(
+                C_inferred[i, :], r, lw=2.5 - (j + 1) / 3, label=str(r0), color=color[i]
+            )
+        else:
+            plt.plot(C_inferred[i, :], r, lw=2.5 - (j + 1) / 3, color=color[i])
+plt.xlabel("Concentration~[kg/m3]")
+plt.ylabel("Distance from transducer~[m]")
+ax.set_xscale("log")
+plt.xlim([0.5 * C.min(), 2.5 * C.max()])
+plt.ylim([0, r.max()])
+ax.invert_yaxis()
+secax = ax.secondary_xaxis("top", functions=(C_to_phi, phi_to_C))
+secax.set_xlabel("Volumic fraction")
+plt.legend(title="Distance of imposed concentration [mm]")
+plt.title("Without near field")
+
+plt.show()
+
+# %%
+# with near field function
+# ^^^^^^^^^^^^^^^^^^^^^^^^
+
+psi = DM.near_field_theoretical(r, rn)  # near field function
+# defining a constant grain concentration profile
+C = Phi[:, None] * rho * (r[None, :] * 0 + 1)
+MSV = DM.create_MSvoltage(C, r[None, :], Xi, a_w, Ks, Kt, psi[None, :])
+#
+color = []
+fig, ax = plt.subplots(
+    1, 1, figsize=(fig_width_small, fig_width_small), constrained_layout=True
+)
+for i, phi in enumerate(Phi):
+    if i == 0:
+        a = plt.plot(C[i, :], r, lw=4, alpha=0.3, label="imposed")
+    else:
+        a = plt.plot(C[i, :], r, lw=4, alpha=0.3)
+    color.append(a[0].get_color())
+for j, index in enumerate(indexes):
+    C0 = C[:, index]
+    r0 = r[index]
+    C_inferred = explicit_inversion(MSV, r, Xi, a_w, psi * 0 + 1, C0, r0)
+    for i, phi in enumerate(Phi):
+        if i == 0:
+            plt.plot(
+                C_inferred[i, :], r, lw=2.5 - (j + 1) / 3, label=str(r0), color=color[i]
+            )
+        else:
+            plt.plot(C_inferred[i, :], r, lw=2.5 - (j + 1) / 3, color=color[i])
+plt.xlabel("Concentration~[kg/m3]")
+plt.ylabel("Distance from transducer~[m]")
+ax.set_xscale("log")
+plt.xlim([0.5 * C.min(), 2.5 * C.max()])
+plt.ylim([0, r.max()])
+ax.invert_yaxis()
+secax = ax.secondary_xaxis("top", functions=(C_to_phi, phi_to_C))
+secax.set_xlabel("Volumic fraction")
+plt.legend(title="Distance of imposed concentration [mm]")
+plt.title("With near field")
+
+plt.show()
+
+# %%
+# with real type concentration profiles
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+psi = DM.near_field_theoretical(r, rn)  # near field function
+# C = Phi[:, None]*rho*np.exp(-r[None, :]/1)  # defining an exponentially decreasing profile
+# C = Phi[:, None]*rho*np.exp(r[None, :]/10)  # defining an exponentially increasing profile
+C = Phi[:, None] * rho * (r[None, :] * 0 + 1)  # defining sedimentation-like profile
+C[..., :230] = 0
+#
+MSV = DM.create_MSvoltage(C, r[None, :], Xi, a_w, Ks, Kt, psi[None, :])
+#
+color = []
+fig, ax = plt.subplots(
+    1, 1, figsize=(fig_width_small, fig_width_small), constrained_layout=True
+)
+for i, phi in enumerate(Phi):
+    if i == 0:
+        a = plt.plot(C[i, :], r, lw=4, alpha=0.3, label="imposed")
+    else:
+        a = plt.plot(C[i, :], r, lw=4, alpha=0.3)
+    color.append(a[0].get_color())
+for j, index in enumerate(indexes):
+    C0 = C[:, index]
+    r0 = r[index]
+    C_inferred = explicit_inversion(MSV, r, Xi, a_w, psi * 0 + 1, C0, r0)
+    for i, phi in enumerate(Phi):
+        if i == 0:
+            plt.plot(
+                C_inferred[i, :], r, lw=2.5 - (j + 1) / 3, label=str(r0), color=color[i]
+            )
+        else:
+            plt.plot(C_inferred[i, :], r, lw=2.5 - (j + 1) / 3, color=color[i])
+plt.xlabel("Concentration~[kg/m3]")
+plt.ylabel("Distance from transducer~[m]")
+ax.set_xscale("log")
+plt.xlim([0.2 * C[:, -1].min(), 2.5 * C.max()])
+plt.ylim([0, r.max()])
+ax.invert_yaxis()
+secax = ax.secondary_xaxis("top", functions=(C_to_phi, phi_to_C))
+secax.set_xlabel("Volumic fraction")
+plt.legend(title="Distance of imposed concentration [mm]")
+plt.title("With exponential concentration profiles")
+
+plt.show()

_downloads/25862d8e3051524fbbbd734856835c1b/plot_sketch_multiple_probes.ipynb

@@ -0,0 +1,43 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n# Sketch multiple probes\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import matplotlib.pyplot as plt\nimport numpy as np\n\nfrom pyudv.geometry import Probe, sketch_probes\n\n# ## define probes\nr = np.linspace(0, 5, 100)  # beam radial coordiantes, same for all probes\nalpha1, alpha2, alpha3, alpha4 = -120, -70, 70, -80  # probe inclinations in degree\n# probe origin points\nO1, O2, O3, O4 = (\n    np.array([1, 8]),\n    np.array([-1, 7]),\n    np.array([2, 2]),\n    np.array([3, 8.5]),\n)\n# probe parameters\nPROBE_PARS = [\n    [r, alpha1, [0, O1]],  # probe1\n    [r, alpha2, [0, O2]],  # probe2\n    [r, alpha3, [0, O3]],  # probe3\n    [r, alpha4, [0, O4]],  # probe4\n]\n#\nprobes = [Probe(*probe_pars) for probe_pars in PROBE_PARS]\nprobe_colors = [\"tab:blue\", \"tab:green\", \"tab:orange\", \"tab:red\"]\ncombinations = [(0, 1), (2, 3)]\ncombination_colors = [\"k\", \"grey\"]\n\nfig, ax = plt.subplots(1, 1, layout=\"constrained\")\nsketch_probes(\n    probes,\n    combinations=combinations,\n    probe_colors=probe_colors,\n    combination_colors=combination_colors,\n    ax=ax,\n)\n\n# adding another probe afterwards\nprobe_add = Probe(r, 70, [0, np.array([4, 2])])\nprobe_add.plot_probe(ax, color=\"tab:pink\")\n\nax.set_xlabel(\"x\")\nax.set_ylabel(\"y\")\nplt.show()"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.11.5"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}

_downloads/39a7d881e7c4cd005f75283a91b4a81d/plot_read_plot_mfprof_data.py

@@ -0,0 +1,69 @@
+"""
+==========================
+Read and plot .mfprof data
+==========================
+"""
+
+import matplotlib.colors as colors
+import matplotlib.pyplot as plt
+import numpy as np
+
+from pyudv.read_mfprof import (
+    amplitude_from_mfprof_reading,
+    read_mfprof,
+    velocity_from_mfprof_reading,
+)
+
+#
+path_data = "src/data_sample.mfprof"
+
+# #### Loading data
+Data, Parameters, Info, Units = read_mfprof(path_data)
+Amplitude_data = amplitude_from_mfprof_reading(Data, Parameters)
+Velocity_data = velocity_from_mfprof_reading(Data, Parameters)
+time = Data["profileTime"]
+z_coordinates = Data["DistanceAlongBeam"] * 1e2
+indmax_time = np.argwhere(Data["transducer"] == 0)[1][0] - 1
+indmax_z = -1
+ind_bottom = 572
+
+# #### plotting velocity
+divnorm = colors.TwoSlopeNorm(vcenter=0, vmin=-0.2, vmax=0.2)
+fig, ax = plt.subplots(1, 1, constrained_layout=True)
+c = ax.pcolormesh(
+    time[:indmax_time],
+    z_coordinates[:indmax_z],
+    Velocity_data[:indmax_z, :indmax_time],
+    cmap="PuOr",
+    norm=divnorm,
+    rasterized=True,
+    shading="auto",
+)
+ax.axhline(y=z_coordinates[ind_bottom], color="k", lw="0.5", ls="--")
+ax.invert_yaxis()
+fig.colorbar(c, label="Velocity [m/s]")
+ax.set_xlabel("Time [s]")
+ax.set_ylabel("DistanceAlongBeam [cm]")
+# fig.draw_without_rendering()
+plt.savefig("plots/Spatio_temporal_velocity.pdf", dpi=600)
+plt.show()
+
+# #### plotting amplitude
+divnorm = colors.TwoSlopeNorm(vcenter=0, vmin=-0.1, vmax=0.1)
+fig, ax = plt.subplots(1, 1, constrained_layout=True)
+c = ax.pcolormesh(
+    time[:indmax_time],
+    z_coordinates[:indmax_z],
+    Amplitude_data[:indmax_z, :indmax_time],
+    rasterized=True,
+    shading="auto",
+    cmap="PuOr",
+    norm=divnorm,
+)
+ax.axhline(y=z_coordinates[ind_bottom], color="k", lw="0.5", ls="--")
+ax.invert_yaxis()
+fig.colorbar(c, label="Amplitude [V]")
+ax.set_xlabel("Time [s]")
+ax.set_ylabel("DistanceAlongBeam [cm]")
+plt.savefig("plots/Spatio_temporal_amplitude.pdf", dpi=600)
+plt.show()