{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "296a21c0-86b3-498a-af11-23041d306a1f",
   "metadata": {},
   "source": [
    "# Dimensionally-inconsistent Calculations (Emperical formulae)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82ccbff3-3f88-403e-a9e5-c99db53c998d",
   "metadata": {},
   "source": [
    "Sometimes calculations and formulae are based on curve-fitting the results of various tests. The resulting formula from the curve-fitting may assume input units with a particular scaling and produce outputs assumed to be in a different scaling. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "6c82f362-cec3-4682-a6fd-cdd7e74e9dcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import forallpeople as si\n",
    "si.environment('default')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9a70767-f926-4e1e-9ff2-1cc6ac026957",
   "metadata": {},
   "source": [
    "## Dimensionally-inconsistent radicals"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "172026b3-8075-47e2-a952-c247d3a02198",
   "metadata": {},
   "source": [
    "Observed in many concrete design codes is a value expressed as $\\sqrt{f'_c}$ where $f'_c$ is in scaled units of stress (e.g. MPa) and $\\sqrt{f'_c}$ is also in units of stress (MPa in this example).\n",
    "\n",
    "The intention behind this particular expression is to compute the square-root of `35` and for that value to be in units of `MPa`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faa9f9f9-2501-455d-ae09-0e0a40c25732",
   "metadata": {},
   "source": [
    "### Strategy 1 - Manually add the units back on"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "f6d455f2-7509-4302-b994-c2e70919024c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "5.916 MPa"
      ],
      "text/latex": [
       "5.916\\ \\text{MPa}"
      ],
      "text/markdown": [
       "5.916 MPa"
      ],
      "text/plain": [
       "5.916 MPa"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import math\n",
    "MPa = 1e6 * si.Pa\n",
    "f_c = 35e6 * si.Pa\n",
    "\n",
    "math.sqrt(f_c) * MPa"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fa15a9f-a138-49f3-b700-99711a638a1f",
   "metadata": {},
   "source": [
    "### Strategy 2 - Define a custom `sqrt` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "acb6a9a7-d51e-4be8-bde2-d624ddbae09c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "35.000 MPa"
      ],
      "text/latex": [
       "35.000\\ \\text{MPa}"
      ],
      "text/markdown": [
       "35.000 MPa"
      ],
      "text/plain": [
       "35.000 MPa"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "5.916 MPa"
      ],
      "text/latex": [
       "5.916\\ \\text{MPa}"
      ],
      "text/markdown": [
       "5.916 MPa"
      ],
      "text/plain": [
       "5.916 MPa"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "35000.000 mm<sup>2</sup>"
      ],
      "text/latex": [
       "35000.000\\ \\text{mm}^{2}"
      ],
      "text/markdown": [
       "35000.000 mm<sup>2</sup>"
      ],
      "text/plain": [
       "35000.000 mm²"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "187.083 mm"
      ],
      "text/latex": [
       "187.083\\ \\text{mm}"
      ],
      "text/markdown": [
       "187.083 mm"
      ],
      "text/plain": [
       "187.083 mm"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "MPa = 1e6 * si.Pa\n",
    "\n",
    "def sqrt(x):\n",
    "    MPa = 1e6 * si.Pa\n",
    "    if (\n",
    "        isinstance(x, si.Physical) \n",
    "        and x.dimensions == (1, -1, -2, 0, 0, 0, 0)\n",
    "    ):\n",
    "        return math.sqrt(x) * MPa # Special case\n",
    "    else:\n",
    "        return x**(1/2) # General case\n",
    "    \n",
    "f_c = 35 * MPa # Example of special case\n",
    "area = 0.0350 * si.m**2 # Example of general case\n",
    "\n",
    "display(\n",
    "    f_c, \n",
    "    sqrt(f_c)\n",
    ")\n",
    "\n",
    "display(\n",
    "    area, \n",
    "    sqrt(area)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cffbcf8d-5b2d-40bb-8087-13a8a10de7e1",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e9ab264c-aa09-4181-9ca7-a093503f1d54",
   "metadata": {},
   "source": [
    "### Strategy 3\n",
    "\n",
    "Continuing with the example above, a calculation may intend to compute the square root of `35e6` (instead of `35`) and apply the same units to that resulting value.\n",
    "\n",
    "This strategy can be combined with **Strategy 2** to define a custom `sqrt` function for special cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "c643c3ab-9a6d-408e-823d-7324717efa70",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "5.916 kPa"
      ],
      "text/latex": [
       "5.916\\ \\text{kPa}"
      ],
      "text/markdown": [
       "5.916 kPa"
      ],
      "text/plain": [
       "5.916 kPa"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import math\n",
    "MPa = 1e6 * si.Pa\n",
    "f_c = 35 * MPa\n",
    "\n",
    "value, unit = f_c.split()\n",
    "math.sqrt(value) * unit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c2b11a4-7024-4e99-8ab9-1105c6af2786",
   "metadata": {},
   "source": [
    "## Hidden units"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a880206b-e4fa-420c-a9f7-316ca703d836",
   "metadata": {},
   "source": [
    "Another kind of _seemingly_ dimensionally-inconsistent formula may just be a case where units are assumed within a numerical value but not communicated.\n",
    "\n",
    "An example commonly used in Canadian flexural concrete design:\n",
    "\n",
    "$$A_s = 0.0015 (f'_c) (b) (d) - \\sqrt{d^2 - \\frac{3.85 M_f}{(f'_c) (b)}}$$\n",
    "\n",
    "This formula for calculating the required area of flexural tension steel should result in values of area (say, in mm<sup>2</sup>). In its current form, however, it is dimensionally-inconsistent because the value of `0.0015` contains hidden units of `1 / MPa`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "a8b14951-5ee1-4983-ae58-cbf08f3b0260",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "1.017 kN"
      ],
      "text/latex": [
       "1.017\\ \\text{kN}"
      ],
      "text/markdown": [
       "1.017 kN"
      ],
      "text/plain": [
       "1.017 kN"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "si.environment('default')\n",
    "MPa = 1e6 * si.Pa\n",
    "\n",
    "f_c = 35 * MPa\n",
    "b = 0.30 * si.m\n",
    "d = 0.60 * si.m\n",
    "M_f = 200e3 * si.N * si.m\n",
    "\n",
    "A_s = 0.0015 * f_c * b * (d - sqrt(d**2 - (3.85 * M_f)/(f_c * b)))\n",
    "A_s # Resulting units are in N, which are not expected"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24123956-69ca-469c-b10e-c4ed46e1702a",
   "metadata": {},
   "source": [
    "Here, we will manually add the units into the calculation as `(0.0015 / MPa)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "320e62c7-9dbd-495c-8cb2-8480ffcd1edd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "1017.251 mm<sup>2</sup>"
      ],
      "text/latex": [
       "1017.251\\ \\text{mm}^{2}"
      ],
      "text/markdown": [
       "1017.251 mm<sup>2</sup>"
      ],
      "text/plain": [
       "1017.251 mm²"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_s = (0.0015 / MPa) * f_c * b * (d - sqrt(d**2 - (3.85 * M_f)/(f_c * b)))\n",
    "A_s # These are the correct dimensions               "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa4417e0-8d47-41ab-9934-9d098eaa3029",
   "metadata": {},
   "source": [
    "## Completely inconsistent formulae"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bb56540-8b59-49c3-8589-af24cd66ff6f",
   "metadata": {},
   "source": [
    "Last, there are times where the formula essentially serves as a numeric transformation and the units simply need to be manually applied.\n",
    "\n",
    "An example from the Canadian wood design code:\n",
    "\n",
    "$$QS_i = 14t \\sqrt{\\frac{d_e}{1 - \\frac{d_e}{d}}}$$\n",
    "\n",
    "This formula for calculating tension resistance in the parallel-to-grain direction of a bolted timber connection is based on empirically-measured fracture test results.\n",
    "\n",
    "All symbols (`t`, `d_e`, and `d`) are numbers, scaled to units of `mm`, and the resulting value of `QS_i` is a number in units of force, `N`.\n",
    "\n",
    "Because this formula is a strict numeric computation with no dimensional consistency, it is best to treat it as such by manually removing units from the inputs and adding the units to the output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "4f6f4024-659f-4e00-afde-eca9a3f82416",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "228.499 kN"
      ],
      "text/latex": [
       "228.499\\ \\text{kN}"
      ],
      "text/markdown": [
       "228.499 kN"
      ],
      "text/plain": [
       "228.499 kN"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "si.environment('default')\n",
    "import math\n",
    "\n",
    "# Use .prefix to force the numerical value to scale in mm\n",
    "t = (0.354 * si.m).prefix('m')\n",
    "d = (0.620 * si.m).prefix('m')\n",
    "d_e = (0.480 * si.m).prefix('m')\n",
    "\n",
    "t = float(t)\n",
    "d = float(d)\n",
    "d_e = float(d_e)\n",
    "\n",
    "QS_i = 14 * t * math.sqrt(d_e / (1 - d_e/d)) * si.N\n",
    "QS_i"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "forallpeople",
   "language": "python",
   "name": "forallpeople"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}