From ae21827d5a62dc6f3fd55d7014d295a428ef7fd2 Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@fb.com>
Date: Mon, 13 May 2024 15:33:26 -0700
Subject: [PATCH] update smoothing function

---
 Focusing2D/SubpixelSmoothing.ipynb | 27 +++++++++------------------
 1 file changed, 9 insertions(+), 18 deletions(-)

diff --git a/Focusing2D/SubpixelSmoothing.ipynb b/Focusing2D/SubpixelSmoothing.ipynb
index 71e0662..374549e 100644
--- a/Focusing2D/SubpixelSmoothing.ipynb
+++ b/Focusing2D/SubpixelSmoothing.ipynb
@@ -7,7 +7,7 @@
    "outputs": [],
    "source": [
     "using Pkg\n",
-    "Pkg.activate(\"/Users/ianhammond/GitHub/Emitter3DTopOpt\")\n",
+    "#Pkg.activate(\"/Users/ianhammond/GitHub/Emitter3DTopOpt\")\n",
     "using Gridap, Gridap.Geometry, Gridap.Fields, GridapGmsh\n",
     "using LinearAlgebra\n",
     "using NLopt\n",
@@ -285,24 +285,13 @@
     "    x_grad_norm = sqrt.(ifelse.(nonzero_norm, x_grad_helper, 1))\n",
     "    x_grad_norm_eff = ifelse.(nonzero_norm, x_grad_norm, 1)\n",
     "    d = (η .- x_smoothed) ./ x_grad_norm_eff\n",
+    "    \n",
     "    needs_smoothing = nonzero_norm .& (abs.(d) .<= pixel_radius)\n",
-    "    arccos_term = pixel_radius^2 * acos.(\n",
-    "        ifelse.(\n",
-    "            needs_smoothing,\n",
-    "            d / pixel_radius,\n",
-    "            0.0,\n",
-    "        )\n",
-    "    )\n",
-    "    sqrt_term = d .* sqrt.(\n",
-    "        ifelse.(\n",
-    "            needs_smoothing,\n",
-    "            pixel_radius^2 .- d.^2,\n",
-    "            1,\n",
-    "        )\n",
-    "    )\n",
+    "\n",
+    "    R = d ./ pixel_radius\n",
     "    fill_factor = ifelse.(\n",
     "        needs_smoothing,\n",
-    "        (1 / (pi * pixel_radius^2)) * (arccos_term .- sqrt_term),\n",
+    "        0.5 .- 15/16 * R .+ 5/8 * R .^ 3 .- 3 / 16 * R .^ 5,\n",
     "        1,\n",
     "    )\n",
     "    x_minus = x_smoothed .- x_grad_norm * pixel_radius\n",
@@ -319,7 +308,9 @@
     "        x_plus,\n",
     "        x_plus_eff_pert,\n",
     "    )\n",
+    "    \n",
     "    x_projected_smoothed = (1 .- fill_factor) .* x_minus_eff .+ (fill_factor) .* x_plus_eff\n",
+    "    \n",
     "    return flatten(ifelse.(\n",
     "        needs_smoothing,\n",
     "        x_projected_smoothed,\n",
@@ -1098,7 +1089,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Julia 1.10.1",
+   "display_name": "Julia 1.10.0",
    "language": "julia",
    "name": "julia-1.10"
   },
@@ -1106,7 +1097,7 @@
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
-   "version": "1.10.1"
+   "version": "1.10.0"
   }
  },
  "nbformat": 4,