{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3"
  }
 },
 "cells": [
  {
   "cell_type": "code",
   "id": "cell-0",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# --- house style (the Press palette: accent/navy/gold/parchment) ---\nimport matplotlib as mpl\nACCENT, NAVY, GOLD, PARCH = '#7a1f1f', '#1f3a5f', '#b8860b', '#f7f2e7'\nmpl.rcParams.update({'figure.facecolor': PARCH, 'axes.facecolor': '#fffdf6',\n                     'axes.edgecolor': '#c9bfa3', 'font.family': 'serif'})\n\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom matplotlib.collections import PolyCollection\n\n# --- The computational opus [0,1]^2: uniform in xi, logarithmic in eta -----\nNTH, NR = 64, 16\ndelta, L = 0.15, 3.0                               # wall scale, domain radius ratio\nth = np.linspace(0, 2*np.pi, NTH+1)\ny  = delta*((1 + L/delta)**(np.arange(NR+1)/NR) - 1)     # Thm 13.3 closed form\nR  = 1.0 + y                                       # radial stations off the unit circle\n\n# --- Joukowski map about an offset circle ----------------------------------\nc, eps = 1.0, 0.09                                 # map constant, offset (camber/thickness)\nZc = -eps + 1j*eps                                 # circle center -> Joukowski airfoil\nZT = np.exp(1j*th)                                 # cylinder boundary\n\nmesh = np.empty((NR+1, NTH+1), dtype=complex)\nfor i in range(NR+1):\n    zeta = Zc + R[i]*ZT\n    mesh[i] = zeta + c*c/zeta                      # z = zeta + c^2/zeta\n\n# --- Cell tesserae: area and the Jacobian tiling rule |J| ------------------\ndef quad_area(z0, z1, z2, z3):                     # shoelace on a quad\n    zs = np.array([z0, z1, z2, z3])\n    return 0.5*abs(np.sum(zs*np.roll(zs.conj(), -1)).imag)\n\nareas, jacs, polys = [], [], []\nfor i in range(NR):\n    for j in range(NTH):\n        z0, z1 = mesh[i, j],   mesh[i, j+1]\n        z2, z3 = mesh[i+1, j+1], mesh[i+1, j]\n        areas.append(quad_area(z0, z1, z2, z3))\n        zeta_c = Zc + 0.5*(R[i]+R[i+1])*np.exp(1j*0.5*(th[j]+th[j+1]))\n        jacs.append(abs(1 - (c/zeta_c)**2))        # |J| from the analytic derivative\n        polys.append(np.column_stack([[z0.real, z1.real, z2.real, z3.real],\n                                      [z0.imag, z1.imag, z2.imag, z3.imag]]))\nareas, jacs = np.array(areas), np.array(jacs)\n\n# --- Thm 13.2: |J|-share equals volume share, ring by ring -----------------\ncomp_area = (2*np.pi/NTH) * (y[1:] - y[:-1])       # d(xi) x d(eta) per ring\nnum = (jacs.reshape(NR, NTH).mean(axis=1) * comp_area)\njac_share = np.repeat(num/num.sum(), NTH)\nvol_share = areas/areas.sum()\nprint(f\"Thm 13.2 ring check - max |J-share - vol-share|: \"\n      f\"{np.abs(jac_share - vol_share).max():.2e}\")\n\nnear_wall = areas[:NTH].sum()/areas.sum()\nprint(f\"near-wall ring holds {near_wall:.4f} of the area with 1/{NR} of the cells\")\nprint(f\"tessera areas: min {areas.min():.2e}, max {areas.max():.2e}, \"\n      f\"ratio {areas.max()/areas.min():.0f}:1 - unequal pieces, allotted shares\")\n\n# --- Draw: cells colored by |J| (navy -> gold -> accent) -------------------\nfig, ax = plt.subplots(figsize=(8.5, 4.6))\nnorm = mpl.colors.Normalize(jacs.min(), jacs.max())\ncmap = mpl.colors.LinearSegmentedColormap.from_list('press', [NAVY, GOLD, ACCENT])\nax.add_collection(PolyCollection(polys, array=jacs, cmap=cmap, norm=norm,\n                                 edgecolor='k', lw=0.15))\nax.autoscale_view(); ax.set_aspect('equal')\nax.set_title('Joukowski airfoil mesh - a tesseraction in two directions; color = |J| (Thm 13.2)')\nplt.tight_layout(); plt.show()\n\nfig, ax = plt.subplots(figsize=(7, 3))\nax.hist(np.log10(areas), bins=40, color=NAVY, edgecolor='k', lw=0.4)\nax.set(xlabel='log10(tessera area)', ylabel='count',\n       title='the tesserae of an airfoil mesh are wildly unequal - by design')\nplt.tight_layout(); plt.show()\n\n# --- The morph: cylinder -> airfoil, slider-driven (fixed blend) -----------\n# Blend: z = zeta + m*c^2/zeta about the offset center m*(-eps + i*eps).\n# The map strength itself scales with m, so every intermediate frame is a\n# smooth, singularity-free deformation -- no figure-8 pinch.\nfrom ipywidgets import interact, FloatSlider\n\ndef show_morph(m=1.0):\n    \"\"\"m = 0: cylinder mesh.  m = 1: Joukowski airfoil mesh.\"\"\"\n    Zc_m = m*(-eps + 1j*eps)\n    fig, ax = plt.subplots(figsize=(8.5, 4.6))\n    polys_m, jacs_m = [], []\n    def P(rr, tt):\n        zeta = Zc_m + rr*np.exp(1j*tt)\n        z = zeta + m*c*c/zeta                    # map strength scales with m\n        return np.array([z.real, z.imag])\n    for i in range(NR):\n        for j in range(NTH):\n            z0, z1 = P(R[i],   th[j]),   P(R[i],   th[j+1])\n            z2, z3 = P(R[i+1], th[j+1]), P(R[i+1], th[j])\n            polys_m.append(np.array([z0, z1, z2, z3]))\n            zc = Zc_m + 0.5*(R[i]+R[i+1])*np.exp(1j*0.5*(th[j]+th[j+1]))\n            jacs_m.append(abs(1 - m*(c/zc)**2))\n    norm_m = mpl.colors.Normalize(min(jacs_m), max(jacs_m))\n    ax.add_collection(PolyCollection(polys_m, array=np.array(jacs_m),\n                                     cmap=cmap, norm=norm_m, edgecolor='k', lw=0.15))\n    ax.autoscale_view(); ax.set_aspect('equal')\n    ax.set_title(f\"morph = {m:.2f}  {'(cylinder)' if m < 0.5 else '(airfoil)'}  -  color = |J|\")\n    plt.show()\n\ninteract(show_morph,\n         m=FloatSlider(value=1.0, min=0.0, max=1.0, step=0.02,\n                       description='morph', continuous_update=False));\n"
  }
 ]
}