From 9b58f01ba57d118fee9a7389ff6a2b7d51cffe23 Mon Sep 17 00:00:00 2001
From: Dario Panici <dpanici@princeton.edu>
Date: Mon, 2 Dec 2024 11:15:53 -0500
Subject: [PATCH 01/16] add NFP attribute to a few fields

---
 desc/magnetic_fields/_core.py      | 24 ++++++++++++++++++++++--
 desc/magnetic_fields/_dommaschk.py | 12 +++++++++++-
 2 files changed, 33 insertions(+), 3 deletions(-)

diff --git a/desc/magnetic_fields/_core.py b/desc/magnetic_fields/_core.py
index 29299d7a1..6312acb71 100644
--- a/desc/magnetic_fields/_core.py
+++ b/desc/magnetic_fields/_core.py
@@ -1948,12 +1948,22 @@ class ScalarPotentialField(_MagneticField):
         R,phi,Z are arrays of cylindrical coordinates.
     params : dict, optional
         default parameters to pass to potential function
+    NFP : int, optional
+        Whether the field has a discrete periodicity. This is only used when making
+        a ``SplineMagneticField`` from this field using its ``from_field`` method,
+        or when saving this field as an mgrid file using the ``save_mgrid`` method.
 
     """
 
-    def __init__(self, potential, params=None):
+    def __init__(self, potential, params=None, NFP=1):
         self._potential = potential
         self._params = params
+        self._NFP = NFP
+
+    @property
+    def NFP(self):
+        """int: Number of (toroidal) field periods."""
+        return self._NFP
 
     def compute_magnetic_field(
         self, coords, params=None, basis="rpz", source_grid=None, transforms=None
@@ -2042,12 +2052,22 @@ class VectorPotentialField(_MagneticField):
         R,phi,Z are arrays of cylindrical coordinates.
     params : dict, optional
         default parameters to pass to potential function
+    NFP : int, optional
+        Whether the field has a discrete periodicity. This is only used when making
+        a ``SplineMagneticField`` from this field using its ``from_field`` method,
+        or when saving this field as an mgrid file using the ``save_mgrid`` method.
 
     """
 
-    def __init__(self, potential, params=None):
+    def __init__(self, potential, params=None, NFP=1):
         self._potential = potential
         self._params = params
+        self._NFP = NFP
+
+    @property
+    def NFP(self):
+        """int: Number of (toroidal) field periods."""
+        return self._NFP
 
     def _compute_A_or_B(
         self,
diff --git a/desc/magnetic_fields/_dommaschk.py b/desc/magnetic_fields/_dommaschk.py
index 009c7fb72..a26dbff23 100644
--- a/desc/magnetic_fields/_dommaschk.py
+++ b/desc/magnetic_fields/_dommaschk.py
@@ -38,6 +38,10 @@ class DommaschkPotentialField(ScalarPotentialField):
         d_m_l coefficients of V_m_l terms, which multiply the sin(m*phi)*N_m_l-1 terms
     B0: float
         scale strength of the magnetic field's 1/R portion
+    NFP : int, optional
+        Whether the field has a discrete periodicity. This is only used when making
+        a ``SplineMagneticField`` from this field using its ``from_field`` method,
+        or when saving this field as an mgrid file using the ``save_mgrid`` method.
 
     """
 
@@ -50,6 +54,7 @@ def __init__(
         c_arr=jnp.array([0.0]),
         d_arr=jnp.array([0.0]),
         B0=1.0,
+        NFP=1,
     ):
         ms = jnp.atleast_1d(jnp.asarray(ms))
         ls = jnp.atleast_1d(jnp.asarray(ls))
@@ -68,6 +73,11 @@ def __init__(
             jnp.isscalar(B0) or jnp.atleast_1d(B0).size == 1
         ), "B0 should be a scalar value!"
 
+        ms_over_NFP = ms / NFP
+        assert jnp.allclose(
+            ms_over_NFP, ms_over_NFP.astype(int)
+        ), "To enforce desired NFP, `ms` should be all integer multiples of NFP"
+
         params = {}
         params["ms"] = ms
         params["ls"] = ls
@@ -77,7 +87,7 @@ def __init__(
         params["d_arr"] = d_arr
         params["B0"] = B0
 
-        super().__init__(dommaschk_potential, params)
+        super().__init__(dommaschk_potential, params, NFP)
 
     @classmethod
     def fit_magnetic_field(  # noqa: C901

From 65d562f3af3ea6e051a304282e829476b8360276 Mon Sep 17 00:00:00 2001
From: Dario Panici <dpanici@princeton.edu>
Date: Mon, 2 Dec 2024 11:36:37 -0500
Subject: [PATCH 02/16] add NFP check to from_field SplineMagneticField method

---
 desc/magnetic_fields/_core.py | 6 ++++--
 tests/test_magnetic_fields.py | 4 ++--
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/desc/magnetic_fields/_core.py b/desc/magnetic_fields/_core.py
index 6312acb71..59d40a9d4 100644
--- a/desc/magnetic_fields/_core.py
+++ b/desc/magnetic_fields/_core.py
@@ -1885,7 +1885,7 @@ def from_mgrid(cls, mgrid_file, extcur=None, method="cubic", extrap=False):
 
     @classmethod
     def from_field(
-        cls, field, R, phi, Z, params=None, method="cubic", extrap=False, NFP=1
+        cls, field, R, phi, Z, params=None, method="cubic", extrap=False, NFP=None
     ):
         """Create a splined magnetic field from another field for faster evaluation.
 
@@ -1903,7 +1903,8 @@ def from_field(
         extrap : bool
             whether to extrapolate splines beyond specified R,phi,Z
         NFP : int, optional
-            Number of toroidal field periods.
+            Number of toroidal field periods.  If not provided, will default to 1 or
+        the provided field's NFP, if it has that attribute.
 
         """
         R, phi, Z = map(np.asarray, (R, phi, Z))
@@ -1911,6 +1912,7 @@ def from_field(
         shp = rr.shape
         coords = np.array([rr.flatten(), pp.flatten(), zz.flatten()]).T
         BR, BP, BZ = field.compute_magnetic_field(coords, params, basis="rpz").T
+        NFP = setdefault(field.NFP, 1, hasattr(field, "_NFP"))
         try:
             AR, AP, AZ = field.compute_magnetic_vector_potential(
                 coords, params, basis="rpz"
diff --git a/tests/test_magnetic_fields.py b/tests/test_magnetic_fields.py
index 6484615bc..c6500f90c 100644
--- a/tests/test_magnetic_fields.py
+++ b/tests/test_magnetic_fields.py
@@ -1045,11 +1045,11 @@ def test_fourier_current_potential_field_coil_cut_warnings(self):
     @pytest.mark.unit
     def test_spline_field(self, tmpdir_factory):
         """Test accuracy of spline magnetic field."""
-        field1 = ScalarPotentialField(phi_lm, args)
+        field1 = ScalarPotentialField(phi_lm, args, NFP=5)
         R = np.linspace(0.5, 1.5, 20)
         Z = np.linspace(-1.5, 1.5, 20)
         p = np.linspace(0, 2 * np.pi / 5, 40)
-        field2 = SplineMagneticField.from_field(field1, R, p, Z, NFP=5)
+        field2 = SplineMagneticField.from_field(field1, R, p, Z)
 
         np.testing.assert_allclose(
             field1([1.0, 1.0, 1.0]), field2([1.0, 1.0, 1.0]), rtol=1e-2, atol=1e-2

From 612c74a170afad1fdc7ef457caff7e676992273d Mon Sep 17 00:00:00 2001
From: Dario Panici <dpanici@princeton.edu>
Date: Mon, 2 Dec 2024 16:07:14 -0500
Subject: [PATCH 03/16] fix attribute check

---
 desc/magnetic_fields/_core.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/desc/magnetic_fields/_core.py b/desc/magnetic_fields/_core.py
index 59d40a9d4..b8af799d8 100644
--- a/desc/magnetic_fields/_core.py
+++ b/desc/magnetic_fields/_core.py
@@ -1912,7 +1912,7 @@ def from_field(
         shp = rr.shape
         coords = np.array([rr.flatten(), pp.flatten(), zz.flatten()]).T
         BR, BP, BZ = field.compute_magnetic_field(coords, params, basis="rpz").T
-        NFP = setdefault(field.NFP, 1, hasattr(field, "_NFP"))
+        NFP = getattr(field, "_NFP", 1)
         try:
             AR, AP, AZ = field.compute_magnetic_vector_potential(
                 coords, params, basis="rpz"

From ed1a9c1276b2ea2b44a00094eabff4e21e5ec34c Mon Sep 17 00:00:00 2001
From: Dario Panici <dpanici@princeton.edu>
Date: Mon, 2 Dec 2024 16:11:30 -0500
Subject: [PATCH 04/16] update changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 58671a6bb..3545e2175 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,7 @@ New Features
     * use of both this and the ``QuadraticFlux`` objective allows for REGCOIL solutions to be obtained through the optimization framework, and combined with other objectives as well.
 - Changes local area weighting of Bn in QuadraticFlux objective to be the square root of the local area element (Note that any existing optimizations using this objective may need different weights to achieve the same result now.)
 - Adds a new tutorial showing how to use``REGCOIL`` features.
+- Adds an ``NFP`` attribute to ``ScalarPotentialField``, ``VectorPotentialField`` and ``DommaschkPotentialField``, to allow ``SplineMagneticField.from_field`` and ``MagneticField.save_mgrid`` to efficiently take advantage of the discrete toroidal symmetry of these fields, if present.
 
 
 Bug Fixes

From 6b378bbb24cc47817e778154e2254b422f8d6ab0 Mon Sep 17 00:00:00 2001
From: Dario Panici <dpanici@princeton.edu>
Date: Mon, 2 Dec 2024 16:14:27 -0500
Subject: [PATCH 05/16] add test

---
 tests/test_magnetic_fields.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/tests/test_magnetic_fields.py b/tests/test_magnetic_fields.py
index c6500f90c..134408bd4 100644
--- a/tests/test_magnetic_fields.py
+++ b/tests/test_magnetic_fields.py
@@ -1450,11 +1450,13 @@ def test_dommaschk_field_errors():
     b_arr = [1]
     c_arr = [1]
     d_arr = [1, 1]  # length is not equal to the rest
-    with pytest.raises(AssertionError):
+    with pytest.raises(AssertionError, match="size"):
         DommaschkPotentialField(ms, ls, a_arr, b_arr, c_arr, d_arr)
-    d_arr = [1]
+    d_arr = [1]  # test with incorrect NFP
+    with pytest.raises(AssertionError, match="desired NFP"):
+        DommaschkPotentialField(ms, ls, a_arr, b_arr, c_arr, d_arr, NFP=2)
     ms = [-1]  # negative mode number
-    with pytest.raises(AssertionError):
+    with pytest.raises(AssertionError, match=">= 0"):
         DommaschkPotentialField(ms, ls, a_arr, b_arr, c_arr, d_arr)
 
 

From 0d4c33143d67da5841ea2e15c8217e6cb83e6802 Mon Sep 17 00:00:00 2001
From: daniel-dudt <daniel.dudt@princetonstellarators.energy>
Date: Mon, 2 Dec 2024 17:13:24 -0500
Subject: [PATCH 06/16] bug fix for proximal scalar optimization

---
 desc/optimize/_constraint_wrappers.py | 19 ++++++++++++++
 tests/test_optimizer.py               | 38 +++++++++++++++++++++++++++
 2 files changed, 57 insertions(+)

diff --git a/desc/optimize/_constraint_wrappers.py b/desc/optimize/_constraint_wrappers.py
index dc2e3033c..7695a671b 100644
--- a/desc/optimize/_constraint_wrappers.py
+++ b/desc/optimize/_constraint_wrappers.py
@@ -836,6 +836,25 @@ def compute_scaled_error(self, x, constants=None):
         xopt, _ = self._update_equilibrium(x, store=False)
         return self._objective.compute_scaled_error(xopt, constants[0])
 
+    def compute_scalar(self, x, constants=None):
+        """Compute the sum of squares error.
+
+        Parameters
+        ----------
+        x : ndarray
+            State vector.
+        constants : list
+            Constant parameters passed to sub-objectives.
+
+        Returns
+        -------
+        f : float
+            Objective function scalar value.
+
+        """
+        f = jnp.sum(self.compute_scaled_error(x, constants=constants) ** 2) / 2
+        return f
+
     def compute_unscaled(self, x, constants=None):
         """Compute the raw value of the objective function.
 
diff --git a/tests/test_optimizer.py b/tests/test_optimizer.py
index 4f0cce0e0..9a4aac89b 100644
--- a/tests/test_optimizer.py
+++ b/tests/test_optimizer.py
@@ -324,6 +324,44 @@ def test_no_iterations():
     np.testing.assert_allclose(x0, out2["x"])
 
 
+@pytest.mark.regression
+@pytest.mark.optimize
+def test_proximal_scalar():
+    """Test that proximal scalar optimization works."""
+    # test fix for GH issue #1403
+
+    # optimize to reduce DSHAPE volume from 100 m^3 to 90 m^3
+    eq = desc.examples.get("DSHAPE")
+    optimizer = Optimizer("proximal-fmintr")  # proximal scalar optimizer
+    R_modes = np.vstack(
+        (
+            [0, 0, 0],
+            eq.surface.R_basis.modes[
+                np.max(np.abs(eq.surface.R_basis.modes), 1) > 1, :
+            ],
+        )
+    )
+    Z_modes = eq.surface.Z_basis.modes[
+        np.max(np.abs(eq.surface.Z_basis.modes), 1) > 1, :
+    ]
+    objective = ObjectiveFunction(Volume(eq=eq, target=90))  # scalar objective function
+    constraints = (
+        FixBoundaryR(eq=eq, modes=R_modes),
+        FixBoundaryZ(eq=eq, modes=Z_modes),
+        FixIota(eq=eq),
+        FixPressure(eq=eq),
+        FixPsi(eq=eq),
+        ForceBalance(eq=eq),  # force balance constraint for proximal projection
+    )
+    [eq], _ = optimizer.optimize(
+        things=eq,
+        objective=objective,
+        constraints=constraints,
+        verbose=3,
+    )
+    np.testing.assert_allclose(eq.compute("V")["V"], 90)
+
+
 @pytest.mark.regression
 @pytest.mark.slow
 @pytest.mark.optimize

From 6281c6e67396cbf51a8e39e3d2e2872f9d6d5294 Mon Sep 17 00:00:00 2001
From: daniel-dudt <daniel.dudt@princetonstellarators.energy>
Date: Tue, 3 Dec 2024 10:12:45 -0500
Subject: [PATCH 07/16] add new PoweredProfile class

---
 desc/profiles.py | 113 ++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 111 insertions(+), 2 deletions(-)

diff --git a/desc/profiles.py b/desc/profiles.py
index 064871e58..b4778eb03 100644
--- a/desc/profiles.py
+++ b/desc/profiles.py
@@ -252,10 +252,10 @@ class ScaledProfile(_Profile):
 
     Parameters
     ----------
-    profile : Profile
-        Base profile to scale.
     scale : float
         Scale factor.
+    profile : Profile
+        Base profile to scale.
 
     """
 
@@ -335,6 +335,115 @@ def __repr__(self):
         return s
 
 
+class PoweredProfile(_Profile):
+    """Profile raised to a power.
+
+    f_1(x) = f(x)**a
+
+    Parameters
+    ----------
+    power : float
+        Exponent of the new profile.
+    profile : Profile
+        Base profile to raise to a power.
+
+    """
+
+    _io_attrs_ = _Profile._io_attrs_ + ["_profile", "_power"]
+
+    def __init__(self, power, profile, **kwargs):
+        assert isinstance(
+            profile, _Profile
+        ), "profile in a PoweredProfile must be a Profile or subclass, got {}.".format(
+            str(profile)
+        )
+        assert np.isscalar(power), "power must be a scalar."
+
+        self._profile = profile.copy()
+        self._power = power
+
+        kwargs.setdefault("name", profile.name)
+        super().__init__(**kwargs)
+
+    @property
+    def params(self):
+        """ndarray: Parameters for computation [power, profile.params]."""
+        return jnp.concatenate([jnp.atleast_1d(self._power), self._profile.params])
+
+    @params.setter
+    def params(self, x):
+        self._power, self._profile.params = self._parse_params(x)
+
+    def _parse_params(self, x):
+        if x is None:
+            power = self._power
+            params = self._profile.params
+        elif isinstance(x, (tuple, list)) and len(x) == 2:
+            params = x[1]
+            power = x[0]
+        elif np.isscalar(x):
+            power = x
+            params = self._profile.params
+        elif len(x) == len(self._profile.params):
+            power = self._power
+            params = x
+        elif len(x) == len(self.params):
+            power = x[0]
+            params = x[1:]
+        else:
+            raise ValueError("Got wrong number of parameters for PoweredProfile")
+        return power, params
+
+    def compute(self, grid, params=None, dr=0, dt=0, dz=0):
+        """Compute values of profile at specified nodes.
+
+        Parameters
+        ----------
+        grid : Grid
+            locations to compute values at.
+        params : array-like
+            Parameters to use. If not given, uses the
+            values given by the self.params attribute.
+        dr, dt, dz : int
+            derivative order in rho, theta, zeta.
+
+        Returns
+        -------
+        values : ndarray
+            values of the profile or its derivative at the points specified.
+
+        """
+        if dt > 0 or dz > 0:
+            raise NotImplementedError(
+                "Poloidal and toroidal derivatives of PoweredProfile have not been "
+                + "implemented yet."
+            )
+        power, params = self._parse_params(params)
+        f0 = self._profile.compute(grid, params, 0, dt, dz)
+        if dr >= 1:
+            df1 = self._profile.compute(grid, params, 1, dt, dz)  # df/dr
+            fn1 = self.compute(grid, (power - 1, params), 0, dt, dz)  # f^(n-1)
+        if dr >= 2:
+            df2 = self._profile.compute(grid, params, 2, dt, dz)  # d^2f/dr^2
+            fn2 = self.compute(grid, (power - 2, params), 0, dt, dz)  # f^(n-2)
+        if dr == 0:
+            f = f0**power
+        elif dr == 1:
+            f = power * fn1 * df1
+        elif dr == 2:
+            f = power * ((power - 1) * fn2 * df1**2 + fn1 * df2)
+        else:
+            raise NotImplementedError("dr > 2 not implemented for PoweredProfile!")
+        return f
+
+    def __repr__(self):
+        """Get the string form of the object."""
+        s = super().__repr__()
+        s = s[:-1]
+        s += ", power={})".format(self._power)
+        return s
+
+
 class SumProfile(_Profile):
     """Sum of two or more Profiles.
 

From f106527526982058214781ef3ba4af5d235b0198 Mon Sep 17 00:00:00 2001
From: daniel-dudt <daniel.dudt@princetonstellarators.energy>
Date: Tue, 3 Dec 2024 10:37:30 -0500
Subject: [PATCH 08/16] update Changelog

---
 CHANGELOG.md | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 58671a6bb..1b784b5fb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@ Changelog
 =========
 
 New Features
+
 - Adds ``from_input_file`` method to ``Equilibrium`` class to generate an ``Equilibrium`` object with boundary, profiles, resolution and flux specified in a given DESC or VMEC input file
 - Adds function ``solve_regularized_surface_current`` to ``desc.magnetic_fields`` module that implements the REGCOIL algorithm (Landreman, (2017)) for surface current normal field optimization
     * Can specify the tuple ``current_helicity=(M_coil, N_coil)`` to determine if resulting contours correspond to helical topology (both ``(M_coil, N_coil)`` not equal to 0), modular (``N_coil`` equal to 0 and ``M_coil`` nonzero) or windowpane/saddle (``M_coil`` and ``N_coil`` both zero)
@@ -17,10 +18,10 @@ New Features
 
 Bug Fixes
 
-- Fixes bug that occurs when taking the gradient of ``root`` and ``root_scalar`` with newer versions of JAX (>=0.4.34) and unpins the JAX version
+- Fixes bug that occurs when taking the gradient of ``root`` and ``root_scalar`` with newer versions of JAX (>=0.4.34) and unpins the JAX version.
 - Changes ``FixLambdaGauge`` constraint to now enforce zero flux surface average for lambda, instead of enforcing lambda(rho,0,0)=0 as it was incorrectly doing before.
 - Fixes bug in ``softmin/softmax`` implementation.
-
+- Fixes bug that occured when using ``ProximalProjection`` with a scalar optimization algorithm.
 
 v0.12.3
 -------

From d00b4037d144e4f2def273ecdf13a429f3f6be9f Mon Sep 17 00:00:00 2001
From: daniel-dudt <daniel.dudt@princetonstellarators.energy>
Date: Tue, 3 Dec 2024 12:14:07 -0500
Subject: [PATCH 09/16] add tests

---
 tests/test_profiles.py | 42 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 42 insertions(+)

diff --git a/tests/test_profiles.py b/tests/test_profiles.py
index bf75489a6..a69f918ff 100644
--- a/tests/test_profiles.py
+++ b/tests/test_profiles.py
@@ -189,6 +189,7 @@ def test_repr(self):
         assert "SumProfile" in str(pp + zp)
         assert "ProductProfile" in str(pp * zp)
         assert "ScaledProfile" in str(2 * zp)
+        assert "PoweredProfile" in str(zp**2)
 
     @pytest.mark.unit
     def test_get_set(self):
@@ -347,6 +348,47 @@ def test_scaled_profiles(self):
         np.testing.assert_allclose(pp.params, [1, -2, 1])
         np.testing.assert_allclose(f(x), 8 * (pp(x)), atol=1e-3)
 
+    @pytest.mark.unit
+    def test_powered_profiles(self):
+        """Test raising profiles to a power."""
+        pp = PowerSeriesProfile(
+            modes=np.array([0, 1, 2, 4]), params=np.array([1, 0, -2, 1]), sym="auto"
+        )
+
+        f = pp**3
+        x = np.linspace(0, 1, 50)
+        np.testing.assert_allclose(f(x), (pp(x)) ** 3, atol=1e-3)
+
+        params = f.params
+        assert params[0] == 3
+        assert all(params[1:] == pp.params)
+
+        f.params = 2
+        np.testing.assert_allclose(f(x), (pp(x)) ** 2, atol=1e-3)
+
+        f.params = 0.5
+        np.testing.assert_allclose(f(x), np.sqrt(pp(x)), atol=1e-3)
+
+    @pytest.mark.unit
+    def test_powered_profiles_derivative(self):
+        """Test that powered profiles computes the derivative correctly."""
+        x = np.linspace(0, 1, 50)
+        p1 = PowerSeriesProfile(
+            modes=np.array([0, 1, 2, 4]), params=np.array([1, 3, -2, 4]), sym="auto"
+        )
+        p2 = p1 * p1
+        p3 = p1 * p2
+
+        f3 = p1**3
+        np.testing.assert_allclose(f3(x, dr=0), p3(x, dr=0))
+        np.testing.assert_allclose(f3(x, dr=1), p3(x, dr=1))
+        np.testing.assert_allclose(f3(x, dr=2), p3(x, dr=2))
+
+        f2 = f3 ** (2 / 3)
+        np.testing.assert_allclose(f2(x, dr=0), p2(x, dr=0))
+        np.testing.assert_allclose(f2(x, dr=1), p2(x, dr=1))
+        np.testing.assert_allclose(f2(x, dr=2), p2(x, dr=2))
+
     @pytest.mark.unit
     def test_profile_errors(self):
         """Test error checking when creating and working with profiles."""

From b792562a16c8b2e7b1a95e202ecefb00b6a7c2b2 Mon Sep 17 00:00:00 2001
From: daniel-dudt <daniel.dudt@princetonstellarators.energy>
Date: Tue, 3 Dec 2024 12:14:43 -0500
Subject: [PATCH 10/16] update profile docs

---
 desc/profiles.py         | 11 +++++++++++
 docs/api.rst             | 10 ++++++----
 docs/api_equilibrium.rst |  9 ++++++---
 3 files changed, 23 insertions(+), 7 deletions(-)

diff --git a/desc/profiles.py b/desc/profiles.py
index b4778eb03..7862a566d 100644
--- a/desc/profiles.py
+++ b/desc/profiles.py
@@ -244,6 +244,17 @@ def __sub__(self, x):
         """Subtract another profile from this one."""
         return self.__add__(-x)
 
+    def __pow__(self, x):
+        """Raise this profile to a power."""
+        if np.isscalar(x):
+            return PoweredProfile(x, self)
+        else:
+            raise NotImplementedError()
+
+    def __rpow__(self, x):
+        """Raise this profile to a power."""
+        return self.__pow__(x)
+
 
 class ScaledProfile(_Profile):
     """Profile times a constant value.
diff --git a/docs/api.rst b/docs/api.rst
index ca3484b5d..43a42adfd 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -290,14 +290,16 @@ Profiles
     :recursive:
     :template: class.rst
 
-    desc.profiles.PowerSeriesProfile
-    desc.profiles.TwoPowerProfile
-    desc.profiles.SplineProfile
+    desc.profiles.FourierZernikeProfile
     desc.profiles.HermiteSplineProfile
     desc.profiles.MTanhProfile
+    desc.profiles.PoweredProfile
+    desc.profiles.PowerSeriesProfile
+    desc.profiles.ProductProfile
     desc.profiles.ScaledProfile
+    desc.profiles.SplineProfile
     desc.profiles.SumProfile
-    desc.profiles.ProductProfile
+    desc.profiles.TwoPowerProfile
 
 Transform
 *********
diff --git a/docs/api_equilibrium.rst b/docs/api_equilibrium.rst
index ad2586262..39f89dbba 100644
--- a/docs/api_equilibrium.rst
+++ b/docs/api_equilibrium.rst
@@ -53,12 +53,15 @@ profiles together by addition, multiplication, or scaling.
     :recursive:
     :template: class.rst
 
-    desc.profiles.PowerSeriesProfile
-    desc.profiles.SplineProfile
+    desc.profiles.HermiteSplineProfile
     desc.profiles.MTanhProfile
+    desc.profiles.PoweredProfile
+    desc.profiles.PowerSeriesProfile
+    desc.profiles.ProductProfile
     desc.profiles.ScaledProfile
+    desc.profiles.SplineProfile
     desc.profiles.SumProfile
-    desc.profiles.ProductProfile
+    desc.profiles.TwoPowerProfile
 
 
 Utilities

From e4569136076fb884ef37689aec2fb0026a6c2aa4 Mon Sep 17 00:00:00 2001
From: daniel-dudt <daniel.dudt@princetonstellarators.energy>
Date: Tue, 3 Dec 2024 12:16:03 -0500
Subject: [PATCH 11/16] update Changelog

---
 CHANGELOG.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 58671a6bb..bf40b8624 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,7 +13,7 @@ New Features
     * use of both this and the ``QuadraticFlux`` objective allows for REGCOIL solutions to be obtained through the optimization framework, and combined with other objectives as well.
 - Changes local area weighting of Bn in QuadraticFlux objective to be the square root of the local area element (Note that any existing optimizations using this objective may need different weights to achieve the same result now.)
 - Adds a new tutorial showing how to use``REGCOIL`` features.
-
+- Adds a new profile class ``PoweredProfile`` for rasing profiles to a power.
 
 Bug Fixes
 

From a11567d68af8c78e1deb7043fad42a29376326c2 Mon Sep 17 00:00:00 2001
From: daniel-dudt <daniel.dudt@princetonstellarators.energy>
Date: Tue, 3 Dec 2024 15:22:46 -0500
Subject: [PATCH 12/16] reviewer suggested changes

---
 desc/profiles.py         | 32 +++++++++++++++++++++++++-------
 docs/api_equilibrium.rst | 12 ++++++------
 tests/test_profiles.py   |  4 ++++
 3 files changed, 35 insertions(+), 13 deletions(-)

diff --git a/desc/profiles.py b/desc/profiles.py
index 7862a566d..c278b1d16 100644
--- a/desc/profiles.py
+++ b/desc/profiles.py
@@ -373,9 +373,21 @@ def __init__(self, power, profile, **kwargs):
         self._profile = profile.copy()
         self._power = power
 
+        self._check_params()
+
         kwargs.setdefault("name", profile.name)
         super().__init__(**kwargs)
 
+    def _check_params(self, params=None):
+        """Check params and throw warnings or errors if necessary."""
+        params = self.params if params is None else params
+        power, params = self._parse_params(params)
+        warnif(
+            power < 0,
+            UserWarning,
+            "This profile may be undefined at some points because power < 0.",
+        )
+
     @property
     def params(self):
         """ndarray: Parameters for computation [power, profile.params]."""
@@ -383,6 +395,7 @@ def params(self):
 
     @params.setter
     def params(self, x):
+        self._check_params(x)
         self._power, self._profile.params = self._parse_params(x)
 
     def _parse_params(self, x):
@@ -844,17 +857,24 @@ def __init__(self, params=None, name=""):
             params = [0, 1, 1]
         self._params = np.atleast_1d(params)
 
+        self._check_params()
+
+    def _check_params(self, params=None):
+        """Check params and throw warnings or errors if necessary."""
+        params = self.params if params is None else params
         errorif(
-            self._params.size != 3, ValueError, "params must be an array of size 3."
+            params.size != 3,
+            ValueError,
+            f"params must be an array of size 3, got {len(params)}.",
         )
         warnif(
-            self._params[1] < 1,
+            params[1] < 1,
             UserWarning,
             "Derivatives of this profile will be infinite at rho=0 "
             + "because params[1] < 1.",
         )
         warnif(
-            self._params[2] < 1,
+            params[2] < 1,
             UserWarning,
             "Derivatives of this profile will be infinite at rho=1 "
             + "because params[2] < 1.",
@@ -868,10 +888,8 @@ def params(self):
     @params.setter
     def params(self, new):
         new = jnp.atleast_1d(jnp.asarray(new))
-        if new.size == 3:
-            self._params = jnp.asarray(new)
-        else:
-            raise ValueError(f"params should be an array of size 3, got {len(new)}.")
+        self._check_params(new)
+        self._params = new
 
     def compute(self, grid, params=None, dr=0, dt=0, dz=0):
         """Compute values of profile at specified nodes.
diff --git a/docs/api_equilibrium.rst b/docs/api_equilibrium.rst
index 39f89dbba..01bff7caf 100644
--- a/docs/api_equilibrium.rst
+++ b/docs/api_equilibrium.rst
@@ -53,15 +53,15 @@ profiles together by addition, multiplication, or scaling.
     :recursive:
     :template: class.rst
 
-    desc.profiles.HermiteSplineProfile
-    desc.profiles.MTanhProfile
-    desc.profiles.PoweredProfile
     desc.profiles.PowerSeriesProfile
-    desc.profiles.ProductProfile
-    desc.profiles.ScaledProfile
     desc.profiles.SplineProfile
-    desc.profiles.SumProfile
+    desc.profiles.MTanhProfile
     desc.profiles.TwoPowerProfile
+    desc.profiles.HermiteSplineProfile
+    desc.profiles.ScaledProfile
+    desc.profiles.ProductProfile
+    desc.profiles.SumProfile
+    desc.profiles.PoweredProfile
 
 
 Utilities
diff --git a/tests/test_profiles.py b/tests/test_profiles.py
index a69f918ff..548e63835 100644
--- a/tests/test_profiles.py
+++ b/tests/test_profiles.py
@@ -425,6 +425,10 @@ def test_profile_errors(self):
             tp.compute(grid, dr=3)
         with pytest.raises(NotImplementedError):
             mp.compute(grid, dr=3)
+        with pytest.raises(UserWarning):
+            tp.params = [1, 0.3, 0.7]
+        with pytest.raises(UserWarning):
+            a = sp**-1
 
     @pytest.mark.unit
     def test_default_profiles(self):

From 02ea2e1b9f428262acd84537c1592bb79379099b Mon Sep 17 00:00:00 2001
From: Daniel Dudt <33005725+ddudt@users.noreply.github.com>
Date: Tue, 3 Dec 2024 17:07:25 -0500
Subject: [PATCH 13/16] Update CHANGELOG.md

Co-authored-by: Dario Panici <37969854+dpanici@users.noreply.github.com>
---
 CHANGELOG.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 28448aa5b..2b2fd1353 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,7 +14,7 @@ New Features
     * use of both this and the ``QuadraticFlux`` objective allows for REGCOIL solutions to be obtained through the optimization framework, and combined with other objectives as well.
 - Changes local area weighting of Bn in QuadraticFlux objective to be the square root of the local area element (Note that any existing optimizations using this objective may need different weights to achieve the same result now.)
 - Adds a new tutorial showing how to use``REGCOIL`` features.
-- Adds a new profile class ``PoweredProfile`` for rasing profiles to a power.
+- Adds a new profile class ``PoweredProfile`` for raising profiles to a power.
 - Adds ``SurfaceQuadraticFlux`` objective which minimizes the quadratic magnetic flux through a ``FourierRZToroidalSurface`` object, allowing for optimizing for Quadratic flux minimizing (QFM) surfaces.
 - Allows ``ToroidalFlux`` objective to accept ``FourierRZToroidalSurface`` so it can be used to specify the toroidal flux through a QFM surface.
 - Adds ``eq_fixed`` flag to ``ToroidalFlux`` to allow for the equilibrium/QFM surface to vary during optimization, useful for single-stage optimizations.

From 95f856655b9e2e9c770c4b2fcebe24c0c5de3901 Mon Sep 17 00:00:00 2001
From: Dario Panici <dpanici@princeton.edu>
Date: Wed, 4 Dec 2024 01:59:32 -0500
Subject: [PATCH 14/16] add v0.13.0 marker to changelog

---
 CHANGELOG.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 14542a3e8..3aba729b7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,9 @@
 Changelog
 =========
 
+v0.13.0
+-------
+
 New Features
 
 - Adds ``from_input_file`` method to ``Equilibrium`` class to generate an ``Equilibrium`` object with boundary, profiles, resolution and flux specified in a given DESC or VMEC input file

From d6406de48c0af0a9b8e361425ebe792f097645b1 Mon Sep 17 00:00:00 2001
From: Dario Panici <dpanici@princeton.edu>
Date: Wed, 4 Dec 2024 12:22:49 -0500
Subject: [PATCH 15/16] add poincare_plot to docs

---
 docs/api_plotting.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/api_plotting.rst b/docs/api_plotting.rst
index f104735ad..4c15c760d 100644
--- a/docs/api_plotting.rst
+++ b/docs/api_plotting.rst
@@ -27,6 +27,7 @@ Plotting Flux Surfaces
     desc.plotting.plot_comparison
     desc.plotting.plot_boundary
     desc.plotting.plot_boundaries
+    desc.plotting.poincare_plot
 
 
 Contour Plots of 2-D data

From 0bc2cf1cf69df746172b194bd394d01e521d0ee3 Mon Sep 17 00:00:00 2001
From: daniel-dudt <daniel.dudt@princetonstellarators.energy>
Date: Wed, 4 Dec 2024 15:03:58 -0500
Subject: [PATCH 16/16] rename 'PoweredProfile' to 'PowerProfile'

---
 CHANGELOG.md             |  5 ++++-
 desc/profiles.py         | 12 ++++++------
 docs/api.rst             |  2 +-
 docs/api_equilibrium.rst |  2 +-
 tests/test_profiles.py   |  2 +-
 5 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7a2ced0af..c45684689 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,10 @@
 Changelog
 =========
 
+New Feature
+
+- Adds a new profile class ``PowerProfile`` for raising profiles to a power.
+
 v0.13.0
 -------
 
@@ -17,7 +21,6 @@ New Features
     * use of both this and the ``QuadraticFlux`` objective allows for REGCOIL solutions to be obtained through the optimization framework, and combined with other objectives as well.
 - Changes local area weighting of Bn in QuadraticFlux objective to be the square root of the local area element (Note that any existing optimizations using this objective may need different weights to achieve the same result now.)
 - Adds a new tutorial showing how to use``REGCOIL`` features.
-- Adds a new profile class ``PoweredProfile`` for raising profiles to a power.
 - Adds an ``NFP`` attribute to ``ScalarPotentialField``, ``VectorPotentialField`` and ``DommaschkPotentialField``, to allow ``SplineMagneticField.from_field`` and ``MagneticField.save_mgrid`` to efficiently take advantage of the discrete toroidal symmetry of these fields, if present.
 - Adds ``SurfaceQuadraticFlux`` objective which minimizes the quadratic magnetic flux through a ``FourierRZToroidalSurface`` object, allowing for optimizing for Quadratic flux minimizing (QFM) surfaces.
 - Allows ``ToroidalFlux`` objective to accept ``FourierRZToroidalSurface`` so it can be used to specify the toroidal flux through a QFM surface.
diff --git a/desc/profiles.py b/desc/profiles.py
index c278b1d16..ae1684abb 100644
--- a/desc/profiles.py
+++ b/desc/profiles.py
@@ -247,7 +247,7 @@ def __sub__(self, x):
     def __pow__(self, x):
         """Raise this profile to a power."""
         if np.isscalar(x):
-            return PoweredProfile(x, self)
+            return PowerProfile(x, self)
         else:
             raise NotImplementedError()
 
@@ -346,7 +346,7 @@ def __repr__(self):
         return s
 
 
-class PoweredProfile(_Profile):
+class PowerProfile(_Profile):
     """Profile raised to a power.
 
     f_1(x) = f(x)**a
@@ -365,7 +365,7 @@ class PoweredProfile(_Profile):
     def __init__(self, power, profile, **kwargs):
         assert isinstance(
             profile, _Profile
-        ), "profile in a PoweredProfile must be a Profile or subclass, got {}.".format(
+        ), "profile in a PowerProfile must be a Profile or subclass, got {}.".format(
             str(profile)
         )
         assert np.isscalar(power), "power must be a scalar."
@@ -415,7 +415,7 @@ def _parse_params(self, x):
             power = x[0]
             params = x[1:]
         else:
-            raise ValueError("Got wrong number of parameters for PoweredProfile")
+            raise ValueError("Got wrong number of parameters for PowerProfile")
         return power, params
 
     def compute(self, grid, params=None, dr=0, dt=0, dz=0):
@@ -439,7 +439,7 @@ def compute(self, grid, params=None, dr=0, dt=0, dz=0):
         """
         if dt > 0 or dz > 0:
             raise NotImplementedError(
-                "Poloidal and toroidal derivatives of PoweredProfile have not been "
+                "Poloidal and toroidal derivatives of PowerProfile have not been "
                 + "implemented yet."
             )
         power, params = self._parse_params(params)
@@ -457,7 +457,7 @@ def compute(self, grid, params=None, dr=0, dt=0, dz=0):
         elif dr == 2:
             f = power * ((power - 1) * fn2 * df1**2 + fn1 * df2)
         else:
-            raise NotImplementedError("dr > 2 not implemented for PoweredProfile!")
+            raise NotImplementedError("dr > 2 not implemented for PowerProfile!")
         return f
 
     def __repr__(self):
diff --git a/docs/api.rst b/docs/api.rst
index 8e675baa2..fa32edb86 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -295,7 +295,7 @@ Profiles
     desc.profiles.FourierZernikeProfile
     desc.profiles.HermiteSplineProfile
     desc.profiles.MTanhProfile
-    desc.profiles.PoweredProfile
+    desc.profiles.PowerProfile
     desc.profiles.PowerSeriesProfile
     desc.profiles.ProductProfile
     desc.profiles.ScaledProfile
diff --git a/docs/api_equilibrium.rst b/docs/api_equilibrium.rst
index b532edb04..30a04d730 100644
--- a/docs/api_equilibrium.rst
+++ b/docs/api_equilibrium.rst
@@ -61,7 +61,7 @@ profiles together by addition, multiplication, or scaling.
     desc.profiles.ScaledProfile
     desc.profiles.ProductProfile
     desc.profiles.SumProfile
-    desc.profiles.PoweredProfile
+    desc.profiles.PowerProfile
 
 
 Utilities
diff --git a/tests/test_profiles.py b/tests/test_profiles.py
index 548e63835..ef14c166f 100644
--- a/tests/test_profiles.py
+++ b/tests/test_profiles.py
@@ -189,7 +189,7 @@ def test_repr(self):
         assert "SumProfile" in str(pp + zp)
         assert "ProductProfile" in str(pp * zp)
         assert "ScaledProfile" in str(2 * zp)
-        assert "PoweredProfile" in str(zp**2)
+        assert "PowerProfile" in str(zp**2)
 
     @pytest.mark.unit
     def test_get_set(self):