{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "564d00b4",
   "metadata": {},
   "source": [
    "# Time Normalization Using the Integral Length Scale\n",
    "\n",
    "When comparing time series from a simulation against wind tunnel data or\n",
    "another simulation, raw seconds are not a useful axis. The natural time scale\n",
    "of an ABL flow is the **eddy turnover time** $T_u = L_u / U$, set by the\n",
    "streamwise integral length scale $L_u$ and a representative mean velocity $U$.\n",
    "Normalising time by $T_u$ removes the dependence on the specific reference\n",
    "velocity and length, so that one full unit of $t^* = t / T_u$ corresponds to\n",
    "the passage of one large-eddy structure past the probe.\n",
    "\n",
    "This notebook shows how to compute $L_u$ from a probe time series following\n",
    "Taylor's frozen-turbulence hypothesis, and how to use it to non-dimensionalise\n",
    "time for plotting, statistics windowing, and cross-case comparison.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "<strong>Note.</strong> The streamwise integral length scale $L_u$ is also a target quantity for\n",
    "the inlet itself. When validating an ABL case against a standard, it should be\n",
    "checked alongside $U(z)$ and $I_u(z)$. See the\n",
    "<a href=\"abl-velocity-turbulence-from-profile.ipynb\">ABL Velocity and Turbulence Intensity</a>\n",
    "notebook for those two profiles.\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bbbbac9",
   "metadata": {},
   "source": [
    "## Definition\n",
    "\n",
    "For a stationary, ergodic streamwise velocity signal $u(t)$ at a fixed height,\n",
    "decompose it as $u(t) = U + u'(t)$ with $U = \\overline{u(t)}$. The temporal\n",
    "autocorrelation of the fluctuating component is\n",
    "\n",
    "$$\n",
    "R_{uu}(\\tau) = \\frac{\\overline{u'(t)\\,u'(t+\\tau)}}{\\overline{u'(t)^2}}\n",
    "$$\n",
    "\n",
    "The integral time scale is\n",
    "\n",
    "$$\n",
    "T_{int} = \\int_{0}^{\\tau^*} R_{uu}(\\tau)\\, d\\tau\n",
    "$$\n",
    "\n",
    "where $\\tau^*$ is the first zero crossing of $R_{uu}$. This bounded integral\n",
    "avoids the slow-decay tail that contaminates the open-ended definition.\n",
    "\n",
    "Under Taylor's frozen-turbulence hypothesis (turbulent eddies advect past the\n",
    "probe faster than they evolve), the spatial integral length scale is\n",
    "\n",
    "$$\n",
    "L_u = U \\cdot T_{int}\n",
    "$$\n",
    "\n",
    "and the corresponding eddy turnover time is $T_u = L_u / U = T_{int}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28087999",
   "metadata": {},
   "source": [
    "## Step 1: Compute the Autocorrelation via FFT\n",
    "\n",
    "Direct evaluation of $R_{uu}(\\tau)$ scales as $O(N^2)$. The FFT-based form\n",
    "below is $O(N \\log N)$ and produces an unbiased estimator after dividing by\n",
    "the lag count.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a887b5ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "def autocorrelation_fft(signal: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"Normalised autocorrelation R_uu(tau) for tau >= 0 via FFT.\"\"\"\n",
    "    s = signal - signal.mean()\n",
    "    n = s.size\n",
    "    nfft = 1 << (2 * n - 1).bit_length()\n",
    "    F = np.fft.rfft(s, n=nfft)\n",
    "    acf_full = np.fft.irfft(F * np.conj(F), n=nfft)[:n]\n",
    "    norm = np.arange(n, 0, -1)        # unbiased lag count\n",
    "    acf = acf_full / norm\n",
    "    return acf / acf[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d6d176f",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "<strong>Tip.</strong> The zero-padding to <code>nfft = 1 &lt;&lt; (2 * n - 1).bit_length()</code>\n",
    "avoids circular-correlation wrap-around. The division by\n",
    "<code>np.arange(n, 0, -1)</code> corrects for the shrinking number of overlapping samples\n",
    "at large lag.\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1df1d90b",
   "metadata": {},
   "source": [
    "## Step 2: Integrate to the First Zero Crossing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4136cb9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def integral_length_scale(signal: np.ndarray, dt: float, u_mean: float) -> float:\n",
    "    \"\"\"Lu = u_mean * integral of R_uu(tau) up to the first zero crossing.\"\"\"\n",
    "    acf = autocorrelation_fft(signal)\n",
    "    sign_change = np.where(np.diff(np.sign(acf)))[0]\n",
    "    if sign_change.size == 0:\n",
    "        idx_zero = acf.size\n",
    "    else:\n",
    "        idx_zero = sign_change[0] + 1\n",
    "    T_int = np.trapz(acf[:idx_zero], dx=dt)\n",
    "    return u_mean * T_int"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14064b0f",
   "metadata": {},
   "source": [
    "If $R_{uu}$ does not cross zero within the signal length, the integration\n",
    "extends to the end of the available record; this typically means the signal is\n",
    "too short for a converged estimate, so flag it rather than trust the value.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccdda475",
   "metadata": {},
   "source": [
    "## Step 3: Apply to a Vertical Line Probe\n",
    "\n",
    "Given the probe CSV format produced by AeroSim (see the\n",
    "[ABL Velocity and Turbulence Intensity](abl-velocity-turbulence-from-profile.ipynb)\n",
    "notebook for the columns), compute $L_u$ at every probe point.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b43b68e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df_pts = pd.read_csv(\"line.line_profile line.points.csv\", index_col=\"idx\")\n",
    "df_ux  = pd.read_csv(\"line.line_profile line.ux.csv\")\n",
    "time_steps = df_ux[\"time_step\"].to_numpy(dtype=np.float64)\n",
    "df_ux.drop(columns=[\"time_step\"], inplace=True)\n",
    "df_ux.columns = df_ux.columns.astype(int)\n",
    "\n",
    "z  = df_pts[\"z\"].to_numpy()\n",
    "dt = time_steps[1] - time_steps[0]\n",
    "\n",
    "Lu = np.zeros(df_pts.shape[0])\n",
    "for i, point_idx in enumerate(df_pts.index):\n",
    "    u = df_ux[point_idx].to_numpy(dtype=float)\n",
    "    if u.std() < 1e-9:\n",
    "        Lu[i] = np.nan\n",
    "        continue\n",
    "    Lu[i] = integral_length_scale(u, dt=dt, u_mean=u.mean())\n",
    "\n",
    "for z_target in (10, 20, 50, 100):\n",
    "    j = int((np.abs(z - z_target)).argmin())\n",
    "    print(f\"z = {z[j]:6.2f} m  ->  Lu = {Lu[j]:7.2f} m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecc58af0",
   "metadata": {},
   "source": [
    "The EN 1991-1-4 Annex B.1 reference profile for $L_u(z)$ is\n",
    "\n",
    "$$\n",
    "L(z) = L_t \\left(\\frac{z}{z_t}\\right)^{\\alpha}, \\qquad \\alpha = 0.67 + 0.05 \\ln(z_0)\n",
    "$$\n",
    "\n",
    "with $L_t = 300$ m and $z_t = 200$ m, clamped to $L(z_{min})$ for $z < z_{min}$.\n",
    "Use it to overlay a target curve on the simulated $L_u(z)$, the same way\n",
    "$U(z)$ and $I_u(z)$ are validated.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8427f06",
   "metadata": {},
   "source": [
    "## Step 4: Pick a Reference Pair (Lu_ref, U_ref)\n",
    "\n",
    "For time normalisation, a single $T_u$ value must be chosen, not a profile.\n",
    "Common conventions:\n",
    "\n",
    "| Convention | $L_u$ | $U$ | When to use |\n",
    "| --- | --- | --- | --- |\n",
    "| Reference height | $L_u(z_{ref})$ | $U(z_{ref})$ | Building studies; $z_{ref}$ is the building reference height (e.g., $H$ or $2H/3$) |\n",
    "| Pedestrian level | $L_u(z = 10\\ \\text{m})$ | $U(z = 10\\ \\text{m})$ | Pedestrian comfort or terrain-only studies, matching standard 10 m reporting height |\n",
    "| Domain-mean | $\\overline{L_u}$ over the engineering height range | $\\overline{U}$ over the same range | Quick comparison; less reproducible across studies |\n",
    "\n",
    "State the convention explicitly when reporting results - the same dataset\n",
    "gives different non-dimensional times under different conventions.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "704973ac",
   "metadata": {},
   "source": [
    "## Step 5: Normalise Time\n",
    "\n",
    "Once a reference pair is chosen, the eddy turnover time and the\n",
    "non-dimensional time are\n",
    "\n",
    "$$\n",
    "T_u = \\frac{L_u}{U}, \\qquad t^* = \\frac{t}{T_u}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b620d6d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_ref = 30.0  # building reference height [m]\n",
    "j_ref = int(np.abs(z - z_ref).argmin())\n",
    "Lu_ref = Lu[j_ref]\n",
    "U_ref  = df_ux[df_pts.index[j_ref]].mean()\n",
    "T_u    = Lu_ref / U_ref\n",
    "\n",
    "t      = time_steps - time_steps[0]\n",
    "t_star = t / T_u\n",
    "\n",
    "print(f\"Lu_ref = {Lu_ref:.2f} m  |  U_ref = {U_ref:.2f} m/s  |  T_u = {T_u:.3f} s\")\n",
    "print(f\"Total non-dimensional time covered: t* = {t_star[-1]:.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "232e0c16",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "<strong>Important.</strong> For converged turbulence statistics on a high-rise study, the\n",
    "recommended sampling window is at least $t^* \\geq 100$. If the simulation\n",
    "duration falls short, autocorrelation tails will not flatten and quantities\n",
    "like $L_u$ itself will be underestimated.\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d12b6de8",
   "metadata": {},
   "source": [
    "## Practical Uses\n",
    "\n",
    "**Plotting time series.** Replace the seconds axis with $t^*$. Two simulations\n",
    "with different $U_{ref}$ or different reference heights become directly\n",
    "comparable.\n",
    "\n",
    "**Spectral analysis.** When plotting power spectral densities, use the reduced\n",
    "frequency $f^* = f \\cdot L_u / U$ rather than $f$ in Hz. This collapses the\n",
    "inertial range of independent runs onto a single $-5/3$ slope and makes\n",
    "Eurocode and von Karman target spectra applicable to any case.\n",
    "\n",
    "**Statistics windowing.** Discard the initial transient by trimming the first\n",
    "$\\sim 5$ to $10$ units of $t^*$ before computing means and standard deviations.\n",
    "\n",
    "**Comparing against wind tunnel data.** Tunnel records are reported at model\n",
    "scale with their own $L_u$ and $U_{ref}$. Plotting both datasets in $t^*$\n",
    "removes the scale factor and exposes whether spectra and integral statistics\n",
    "actually agree.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97eaf5b0",
   "metadata": {},
   "source": [
    "## Next Steps\n",
    "\n",
    "- [ABL Velocity and Turbulence Intensity](abl-velocity-turbulence-from-profile.ipynb) - mean profile and turbulence intensity from the same probe.\n",
    "- [ABL guided case post-processing](../guided-cases/ABL/post-processing.md) - full validation workflow against a target terrain category.\n",
    "- [Matching Custom Inlet](matching-custom-inlet.md) - reproducing a measured wind tunnel profile, where $L_u$ is also a calibration target.\n"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}