{ "cells": [ { "cell_type": "markdown", "id": "61c6f09c-d285-42a9-b3d3-9c166ab5271f", "metadata": {}, "source": [ "We are going to learn hot to extract the main features (integral length scale, dissipation rate, Kolmogorov length scale) from hot wire and hot film (Eulerian) data. \n", "\n", "# Data from the largest wind tunnel Grid experiment in Modane S1MA.\n", "The details of the experimental setup can be found in the dedicated data description article [\"Investigation of the small scale statistics of turbulence in the Modane S1MA wind-tunnel\"](https://hal.science/hal-01578798/document).\n", "\n", "The Hot Wire data are [freely available](https://entrepot.recherche.data.gouv.fr/dataset.xhtml?persistentId=doi:10.57745/SL7ZCW), but in order to facilitate the class and avoid unnecessay network traffic, the subset of the data which we will be using is readily available on your virtual desktop." ] }, { "cell_type": "code", "execution_count": null, "id": "0bf0028f-c1e4-4b72-9c77-26057d43ad7a", "metadata": {}, "outputs": [], "source": [ "from ipyfilechooser import FileChooser\n", "\n", "base_path = ''\n", "fc = FileChooser(base_path)\n", "fc.filter_pattern = \"*.mat\"\n", "display(fc)" ] }, { "cell_type": "markdown", "id": "51955606-8fc6-430c-b360-257bfd6d983b", "metadata": {}, "source": [ "## Load the data\n", "The hot wire data have been sampled at 250 kHz and then translated to velocities. They are stored in Matlab 7 format and we will use the Scipy `loadmat` function to load them into the Python workspace." ] }, { "cell_type": "code", "execution_count": null, "id": "bcf8d3be-6d68-4214-83ca-99065ed6382a", "metadata": {}, "outputs": [], "source": [ "from scipy.io import loadmat\n", "data_path = fc.selected\n", "vv = loadmat (data_path)" ] }, { "cell_type": "code", "execution_count": null, "id": "6e424b78-e164-429d-a1cc-8fc1a2732b08", "metadata": {}, "outputs": [], "source": [ "v1 = vv['V4'].flatten()" ] }, { "cell_type": "markdown", "id": "f5291baa-3d3b-49c3-8d59-5b9f96972f4f", "metadata": {}, "source": [ "## Integral length scale\n", "Show a grid of sub-samples at various time scales (e.g. varying logarithmically from 1 ms to 10 s) and try to guess the order of magnitude of the large scale correlation time." ] }, { "cell_type": "code", "execution_count": null, "id": "9d6a6d6f-6a99-4df5-b0a0-7e83c96080e3", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "!jupyter nbconvert --to script TD1_HotWires.ipynb\n", "\n", "# Sample rate\n", "fps = 250e3\n", "\n", "# Time vector\n", "tt = np.arange(len(v1)) / fps\n", "tmax = tt[-1]\n", "\n", "# Create an array of 2-by-4 axes to show the signal at various magnification factors\n", "fig, axs = plt.subplots(2, 4,sharey=True)\n", "\n", "# Increase duration of samples logarithmically from 10**(-3) to 10**(1) seconds\n", "sample_time_s = np.logspace(-3,1,8)\n", "\n", "## Max number of samples shown onscreen\n", "n_samples_plot = 300\n", "\n", "for tau,ax in zip (sample_time_s, axs.flatten()):\n", " n_samples = int(np.round(tau*fps))\n", " stride = 1\n", " ## Don't show too much samples onscreen\n", " if n_samples_plot < n_samples:\n", " stride = int (np.round(n_samples/n_samples_plot))\n", " sample = v1[0:n_samples:stride]\n", " t = tt[0:n_samples:stride]\n", " ax.plot (t, sample)" ] }, { "cell_type": "markdown", "id": "e9a085c7-1f7f-4b90-a102-8aa502df9fee", "metadata": {}, "source": [ "The Integral length scale is defined as \n", "$$L = \\int_0^\\infty R(r) dr$$\n", "where $R$ is the longitudinal auto-correlation coefficient of the velocity, \n", "$$R(r) = \\frac{\\langle v'(x)v'(x+r)\\rangle}{\\langle v'(x)^2\\rangle}$$\n", "Here $v'$ is the fluctuating part of the velocity signal in the Reynolds decomposition $v = \\bar v + v'$.\n", "\n", "**Homework**: Compute and show the auto-correlation coefficient $R(x)$. Steps:\n", "1. First compute the average velocity, which will be used later to map time to space data through the Taylor Hypothesis.\n", "2. In this signal, like in many real life experimental data, there are parasitic long time trends that add up on turbulent velocity correlations. In short, simply subtracting the average velocity to obtain the fluctuating part of it, is not enough. More or less complicated methods have been developped to remove parasitic trends from experimental signals (see, e.g., [the empirical mode decomposition](https://www.sciencedirect.com/topics/engineering/empirical-mode-decomposition)), but here we will try highpass filtering the signal or cutting it into smaller peaces to remove the local mean velocity.\n", "3. Use `scipy.signal.correlate` to compute the correlation **function** and normalize it (with $\\sum v'^2$) to obtain the correlation coefficient $R$\n", "4. Use `scipy.signal.correlation_lags` to obtain time lags $\\tau$, and map to space lags $r$.\n", "5. Plot $R$ for positive values of $r$ (the function is even anyway)." ] }, { "cell_type": "code", "execution_count": null, "id": "33d90e86-c1d5-4842-b593-c0fbd8a68514", "metadata": {}, "outputs": [], "source": [ "from scipy.signal import correlate, correlation_lags,butter,filtfilt,detrend\n", "from scipy.interpolate import CubicSpline\n", "def autocorr (v):\n", " ## Compute the autocorrelation 'function' of the signal\n", " auto_corr = correlate (v, v, 'full')\n", " ## Normalize to obtain R\n", " norm = np.sum (v**2)\n", " R = (auto_corr/norm)\n", " lags = correlation_lags (v.size, v.size, 'full')\n", " ## Only retain positive times\n", " R = R[lags >= 0]\n", " lags = lags[lags >= 0]\n", " return R,lags\n", " \n", "## Several methods to remove the mean velocity and any eventual long term trends\n", "method = \"chunks\" # raw, raw-smoothchunks or highpass\n", "vm = np.mean (v1)\n", "\n", "if method == \"raw\":\n", " v = detrend(v1)\n", " sigmav = np.std (v)\n", " R,lags = autocorr(v1-vm)\n", "elif method == \"raw-smooth\":\n", " # Cut the signal into smaller peaces so that eventual large time trends\n", " # can be removed peace-wise\n", " tau_l = 1000 * (0.15/41);\n", " stride = int (tau_l * fps);\n", " nchunks = int(np.floor(v1.size / stride))\n", " tchunk = np.array(range(nchunks)).flatten() * stride/fps\n", " vchunk = np.zeros ((len(tchunk))).flatten()\n", " for ii in range(nchunks):\n", " print (f\"{ii} / {nchunks-1}\", end=\"\\r\")\n", " vchunk[ii] = np.mean (v1[(ii*stride):(ii*stride+stride)])\n", " spl = CubicSpline(tchunk, vchunk)\n", " vsmooth = spl(tt)\n", " v = v1-vsmooth\n", " sigmav = np.std (v)\n", " R,lags = autocorr(v1-vm)\n", " f, ax = plt.subplots()\n", " ax.plot (tchunk, vchunk)\n", "elif method == \"chunks\":\n", " # Cut the signal into smaller peaces so that eventual large time trends\n", " # can be removed peace-wise\n", " v = v1-vm\n", " tau_l = 1000 * (0.15/41);\n", " stride = int (tau_l * fps);\n", " nchunks = int(np.floor(v.size / stride))\n", " for ii in range(nchunks):\n", " print (f\"{ii} / {nchunks-1}\", end=\"\\r\")\n", " vchunk = v[(ii*stride):(ii*stride+stride)]\n", " if ii == 0:\n", " #R,lags = autocorr(vchunk - np.mean(vchunk))\n", " R,lags = autocorr(detrend (vchunk))\n", " sigmav = np.std(vchunk)\n", " else:\n", " R += autocorr(detrend (vchunk))[0]\n", " sigmav += np.std(vchunk)\n", " R /= nchunks\n", " sigmav /= nchunks\n", "elif method == \"highpass\":\n", " ## High pass filter at 2s (0.5Hz) ~ 80m\n", " fnyq = fps/2;\n", " tau_l = 560 * (0.18/41);\n", " fhp = 1/tau_l;#Hz\n", " b, a = butter(3, fhp/fnyq, 'highpass')\n", " v = filtfilt(b, a, v1)\n", " R, lags = autocorr(v)\n", " sigmav = np.std (v)\n", "\n", "## Taylor hypothesis to change time to space\n", "r = lags / fps * vm\n", "\n", "print (f'Umean = {vm:.2g}, Urms = {sigmav:.2g} m/s (intensity ~ {100*sigmav/vm:.2g} %)')" ] }, { "cell_type": "code", "execution_count": null, "id": "64c9902a-b974-4c3f-aacd-d41f8aa84712", "metadata": {}, "outputs": [], "source": [ "f, ax = plt.subplots()\n", "maxlen = 10\n", "\n", "ax.plot (r[r < maxlen], R[r < maxlen])\n", "ax.set_xlabel('$r$ [m]')\n", "ax.set_ylabel('$R(r)$')\n", "ax.grid()" ] }, { "cell_type": "markdown", "id": "04844771-1f11-4f08-b8da-d5d61f2c8c08", "metadata": {}, "source": [ "6. Now plot the cumulative sum of $R$ (try various max lengths) to try and guess a value for the integral length scale. Comments?" ] }, { "cell_type": "code", "execution_count": null, "id": "3bdaa092-5e02-4b38-87f0-be83adf5844e", "metadata": {}, "outputs": [], "source": [ "f, ax = plt.subplots()\n", "dx = r[1] - r[0]\n", "\n", "## Cumulative sum\n", "Rcum = np.cumsum(R.flatten() * dx)\n", "maxlen = 3\n", "\n", "ax.plot (r[r < maxlen], Rcum[r < maxlen])\n", "ax.set_xlabel('$r$[m]')\n", "ax.set_ylabel('$\\\\int_0^r R(r)dr$')\n", "ax.grid()" ] }, { "cell_type": "markdown", "id": "7a7b7f5b-94df-48c8-91d6-265ad84009f0", "metadata": {}, "source": [ "The order of magnitude of the L is $\\approx 15$ cm, a fraction of the mesh size (60 cm) as expected." ] }, { "cell_type": "markdown", "id": "bdfd8375-eef1-4d39-8f4b-6a20b9485def", "metadata": {}, "source": [ "## Taylor length scale from auto-correlation coeff.\n", "The Taylor length scale can be obtained by fitting $R$ at small $r$ with a parabola $p$ so that\n", "$$R(r) = 1-r^2/\\lambda^2$$\n", "\n", "and determining the value of $r$ where the curve crosses the x axis.\n", "\n", "**Homework:**\n", "Zoom-in $R(r)$, to see its shape at short distances (of order 10 cm). You should see two important features:\n", "1. between the first point, $R(0) = 1$, and the next one, there is a small drop (less than 1%) of $R$, while we would expect a null slope. This is due to the measurement noise adding up on velocity signals. Indeed, if we write the $v'_{meas}(t) = v'(t) + b(t)$ where $b(t)$ is a random noise (uncorrelated with the signal velocity), we obtain\n", " $$R_{meas}(t) = \\frac{\\langle v'(t)v'(t+\\tau)\\rangle + \\langle b(t)b(t+\\tau)\\rangle}{\\langle v'^2\\rangle + \\langle b^2\\rangle}$$ \n", "Here, the noise $b$ seems to have a correlation time smaller than the sampling frequency, because its contribution is only visible at $\\tau = 0 s$.\n", "2. Past the noise drop, the region which could eventually support a parabola is not so clear. We will adopt a *trial and error* method to define the best parabolic fit. " ] }, { "cell_type": "code", "execution_count": null, "id": "1dc5fa85-a633-489c-a1b1-15de3e972843", "metadata": {}, "outputs": [], "source": [ "f, ax = plt.subplots()\n", "maxlen = 0.01\n", "\n", "ax.plot (r[r < maxlen], R[r < maxlen])\n", "ax.set_xlabel('$r$ [m]')\n", "ax.set_ylabel('$R(r)$')\n", "ax.grid()" ] }, { "cell_type": "code", "execution_count": null, "id": "b794d42b-9356-48c1-96d8-4438182a104c", "metadata": {}, "outputs": [], "source": [ "f, ax = plt.subplots()\n", "maxlen = 0.01\n", "\n", "## maximum fit length\n", "l1 = 0.002 # m\n", "\n", "npts_fit = 3 # Number of samples to fit\n", "\n", "while npts_fit < (l1/vm*fps):\n", " ## Assume the second point is a good guess for R(0) intercept\n", " #a = R[1]\n", " #b = np.sum ((R[1:npts_fit] - a) * r[1:npts_fit]**2)/np.sum(r[1:npts_fit]**4)\n", " #lam = np.sqrt(-a/b)\n", " p = np.polyfit (r[1:npts_fit]**2, R[1:npts_fit], 1);\n", " a = p[1]\n", " lam = np.sqrt(-1/p[0])\n", " ax.plot (r[r < maxlen], a-r[r < maxlen]**2/lam**2,\n", " label = f'$\\\\lambda$ = {lam*1000:.2g} mm (n = {npts_fit})')\n", " npts_fit += 1\n", "\n", "ax.plot (r[r < maxlen], R[r < maxlen], linewidth=3, label = '$R(r)$')\n", "plt.xlabel('x[m]')\n", "plt.ylabel('$R$(x)')\n", "plt.grid()\n", "ax.legend()\n", "ax.set_ylim ([np.min(R[r < maxlen]), 1])" ] }, { "cell_type": "markdown", "id": "975eb120-5f5a-4e16-8f3f-ef99dc86d1be", "metadata": {}, "source": [ "The Taylor length scale is of order 10 mm here but we see that this method can hardly give us more than an order of magnitude. \n", "We will check later, using the dissipation rate $\\epsilon$, that this first estimate is coherent with the one obtained from classical HIT laws." ] }, { "cell_type": "markdown", "id": "8844d462-8281-4273-9034-b2b080432e29", "metadata": {}, "source": [ "## Dissipation rate from velocity increments structure function\n", "As explained in the original paper, *In the Kolmogorov 1941 framework, the scaling $S^{\\|}_n(r) \\propto (\\epsilon r)^{n/3}$ is expected for inertial range spacial increments* $r$, where $$S^{\\|}_n(r) = \\langle \\left[(\\mathbf{v}(x+r)-\\mathbf{v}(x)) \\cdot \\frac{\\mathbf{r}}{\\mathbf{r}}\\right]^n \\rangle $$ denotes the longitudinal structure function of order $n$ of the velocity increments across distances $r$ ($\\mathbf{r} = r\\mathbf{e_r}$).\n", "\n", "In this section we are going to use the 2nd and 3rd order structure functions to determine the dissipation rate $\\epsilon$ and then use classical HIT relation to derive the Taylor length scale and the Kolmogorov, dissipative, length scale.\n", "\n", "### Second order structure function\n", "For the second order structure function, the pre-factor of the scaling $S_2 = C_2 (\\epsilon r)^{2/3}$ has been determined empirically to be $C_2\\approx 2.1\\pm 0.1$.\n", "\n", "One can show that $S_2(r)$ can be obtained directly from the auto-correlation coefficient of the velocity:\n", "\n", "$$S_2(r) = 2\\langle v^2\\rangle\\left(1-R(r)\\right)$$\n" ] }, { "cell_type": "code", "execution_count": null, "id": "fda13276-3f4b-4b74-af5c-b6c1ac2d70d6", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "maxlen = 3\n", "idx = (r < maxlen) & (r > 1e-3)\n", "rs2 = r[idx]\n", "S2 = 2*sigmav**2*(1-R[idx])\n", "f,(ax1,ax2) = plt.subplots(1,2)\n", "ax1.loglog (r[idx], S2)\n", "r2_3 = np.logspace(-2,-1)\n", "two_third_law = 5*r2_3**(2/3)\n", "ax1.loglog (r2_3, two_third_law)\n", "\n", "ax1.set_xlabel ('$r$ [m]')\n", "ax1.set_ylabel ('$S_2$')\n", "ax1.grid(which='major', visible=True)\n", "ax1.grid(which='minor', visible=True)\n", "\n", "C2 = 2.1\n", "eps = (S2*rs2**(-2/3)/C2)**1.5\n", "epsilon = np.max(eps)\n", "ax2.loglog (rs2, eps)\n", "ax2.loglog (r2_3, epsilon*np.ones(len(r2_3)))\n", "ax2.text(r2_3[-1], epsilon, f'$\\\\epsilon$ = {epsilon:.2g} $m^2/s^3 (C_\\\\epsilon u^3/L_{1} \\\\approx {0.45*sigmav**3/0.15:.2g})$')\n", "ax2.set_xlabel ('$r$ [m]')\n", "ax2.set_ylabel ('$\\\\left(\\\\frac{S_2}{C_2r^{2/3}}\\\\right)^{3/2}$')\n", "ax2.yaxis.set_label_position(\"right\")\n", "ax2.yaxis.tick_right()\n", "ax2.grid(which='major', visible=True)\n", "ax2.grid(which='minor', visible=True)\n" ] }, { "cell_type": "markdown", "id": "42300311-c800-456c-9826-a79d5bcddb88", "metadata": {}, "source": [ "In the framework of HIT, one can show that the Taylor micro-scale can be obtained by\n", "$$\\lambda = \\sqrt{\\frac{15\\nu \\langle v^2 \\rangle}{\\epsilon}}$$\n", "With the above value of $\\epsilon$, we obtain" ] }, { "cell_type": "code", "execution_count": null, "id": "0b0ad759-b489-48d6-84dc-d046f968568b", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "nu = 1.5e-5\n", "lambd = sigmav*np.sqrt(15*nu/epsilon)\n", "print (f'λ ~ {lambd*1000:.2g} mm and the associated Reynolds number Rλ ~ {lambd*sigmav/nu:.3g}')" ] }, { "cell_type": "markdown", "id": "bc06596a-1ed1-48df-a608-bfcadc86e1a3", "metadata": {}, "source": [ "### Third order structure function\n", "The celebrated four-fith law is another way to devise the dissipation rate, with the advantage that it is an exact result from HIT theory that does not involve any empirical constant (like $C_2$ above).\n", "For sufficiently large Reynolds numbers this law writes\n", "\n", "$$S_3(r) = -\\frac{4}{5} \\epsilon r$$\n", "There is no simple method (other than loops) to obtain $S_3$ so we have to compute it manually:" ] }, { "cell_type": "code", "execution_count": null, "id": "1cb2d9c9-9072-48c7-a13c-0f623399b90c", "metadata": {}, "outputs": [], "source": [ "try:\n", " from numba import njit, prange\n", "except ModuleNotFoundError:\n", " from nonumba import njit\n", " prange = range\n", " \n", "@njit(parallel=True)\n", "def compute_moment_3 (v, n):\n", " S3 = np.zeros ((len(n)), dtype=np.float32).flatten()\n", " for i in prange(len(n)):\n", " nn = int (n[i])\n", " S3[i] = np.mean ((v[nn:]-v[:-nn])**3)\n", " return S3\n", "\n", "maxlen = 1 #m\n", "npts = round (maxlen/vm*fps);\n", "\n", "# Approximate number of distances r at which we compute S3\n", "ncompute = 100 \n", "n = np.unique (np.round (np.logspace (0, np.log10(npts), ncompute))).flatten()\n", "\n", "S3 = compute_moment_3 (v,n)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "195c3063-8829-405c-95d6-5c99b24489f2", "metadata": {}, "outputs": [], "source": [ "r = n / fps * vm\n", "\n", "f, ax = plt.subplots()\n", "eps = S3/(4/5*r)\n", "ax.semilogx (r, eps)\n", "\n", "r = [4e-3, 4e-2]\n", "epsilon = np.max(eps)\n", "ax.loglog (r, epsilon*np.ones(len(r)))\n", "ax.text(r[-1], epsilon, f'$\\\\epsilon$ = {epsilon:3.3g} $m^2/s^3$')\n", "ax.set_ylabel ('$\\\\frac{S_3}{-\\\\frac{4}{5}r}$')\n", "ax.grid(which='major')\n", "ax.grid(which='minor')" ] }, { "cell_type": "markdown", "id": "42c6999b-c549-451c-abd0-553fd01f188d", "metadata": {}, "source": [ "### Kolmogorov dissipation scale\n", "In the framework of HIT, the dissipation scale is defined as\n", "$$\\eta = \\left(\\nu^3/\\epsilon\\right)^{1/4}$$" ] }, { "cell_type": "code", "execution_count": null, "id": "7f11bcd9-88fa-44c0-aa00-37e46f900afa", "metadata": {}, "outputs": [], "source": [ "print (f'η ~ {(nu**3/epsilon)**0.25*1000:.2f} mm')" ] }, { "cell_type": "markdown", "id": "d83059bc-5f4c-4835-8e31-318b1d4de3c6", "metadata": {}, "source": [ "## Spectral representation\n", "For a stationnary process (what we hope is the case for our turbulent flow) the Wiener–Khintchine theorem states that the Power Spectral Density of a signal is equal to the Fourier transform of its auto-correlation function, $\\langle v'^2\\rangle R(t)$. Unfortunately simply applying FFT on the $R$ variable will lead to a very noisy PSD (see bellow). If one can afford to decrease the frequency resolution, it is possible to use more sophisticated algorithms, such as the Welch periodogram method, to devise a coarse grained PSD : conceptually, the signal is cut into smaller blocks (corresponding to time scales of the turbulence, in general smaller than the total duration of the signal) and the (squared module of the) Fourier transform on each block is averaged.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "77029aab-11e0-4fcb-a8c1-43e2b6033696", "metadata": {}, "outputs": [], "source": [ "## FFT method\n", "from scipy.fft import fft,fftfreq, ifft\n", "PSD = np.real(fft (sigmav**2*R))/fps * np.pi\n", "freq = fftfreq(R.size, d=1/fps)\n", "f, ax = plt.subplots()\n", "ax.loglog (freq[freq>0], PSD[freq>0])\n", "\n", "## Welch method\n", "from scipy.signal import welch\n", "freq, PSD = welch (v, fs=fps, nperseg=v1.size/50) ## Take blocks (segments) of 1/50th the total signal length\n", "ax.loglog (freq[freq>0], PSD[freq>0])\n", "ax.set_xlabel ('$f$ [Hz]')\n", "ax.set_ylabel ('PSD $ [m^2s^{-3}]$')\n", "ax.grid(True)\n", "print (f\"rms = {np.sqrt(np.sum(PSD*freq[1]))}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "1c2d3877-ff94-406e-a9bd-2f19c285e7c9", "metadata": {}, "outputs": [], "source": [ "Rpsd = sigmav**-2*np.real(ifft (PSD))*fps\n", "npsd = int(np.ceil(PSD.size/2))\n", "r_psd = np.arange(npsd)*1/(fps/2)*vm\n", "Rpsd = Rpsd[:npsd]/Rpsd[0]\n", "\n", "maxlen = 10\n", "\n", "_, (ax1,ax2) = plt.subplots(1,2)\n", "ax1.plot (r_psd[r_psd < maxlen],Rpsd[r_psd < maxlen])\n", "ax2.plot (r_psd[r_psd < maxlen], np.cumsum(Rpsd[r_psd < maxlen])*r_psd[1])" ] }, { "cell_type": "markdown", "id": "721e4b6e-b74d-47fe-a89c-359d71e3ea88", "metadata": {}, "source": [ "Again, using the Taylor hypothesis, one can transform (inverse) time scales into (inverse) spatial scales: $k[m^{-1}] = f /\\bar v$. Using this transformation, one can analyze the PSD in terms of spatial scales of the turbulence:" ] }, { "cell_type": "code", "execution_count": null, "id": "75b00ec7-2992-41c3-8835-70ce84f6d7b6", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "k = freq/vm\n", "f, ax = plt.subplots()\n", "ax.loglog (k[k>0], PSD[k>0])\n", "\n", "ax.set_xlabel ('$k$ [m$^{-1}$]')\n", "ax.set_ylabel ('$PSD(v)$')\n", "kl = 1/0.15\n", "ykL = PSD[np.argmax(k>kl)]\n", "ax.annotate ('1/$L$', (kl, ykL),\n", " arrowprops={'facecolor':'black'}) \n", "kl = 1/lambd\n", "ykL = PSD[np.argmax(k>kl)]\n", "ax.annotate ('1/$\\\\lambda$', (kl, ykL),\n", " arrowprops={'facecolor':'black'}) \n", "kl = 1/2E-4\n", "ykL = np.min(PSD)\n", "ax.annotate ('1/$\\\\eta$', (kl, ykL),\n", " arrowprops={'facecolor':'black'}) \n", "ax.grid(True)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }