Merge pull request #1226 from LourensVeen/issue_1221_protoplanetary_disk_broken

Fix protoplanetary disk parameter rename

Author: LourensVeen

src/amuse/ext/protodisk.py

@@ -1,5 +1,6 @@
 import numpy
 import warnings
+from amuse.support.helpers import rename_fn_par
 from amuse.ext.evrard_test import body_centered_grid_unit_cube
 from amuse.ext.evrard_test import regular_grid_unit_cube
 from amuse.ext.evrard_test import uniform_random_unit_cube
@@ -8,157 +9,148 @@
 
 from amuse.datamodel import Particles
 from amuse.datamodel import ParticlesWithUnitsConverted
+
+
 def approximate_inverse_error_function(x):
-  a=8*(numpy.pi-3)/3*numpy.pi*(4-numpy.pi)
-  return numpy.sign(x)*numpy.sqrt(
-    numpy.sqrt((2/numpy.pi/a+numpy.log(1-x**2)/2)**2-numpy.log(1-x**2)/a)-(2/numpy.pi/a+numpy.log(1-x**2)/2)
-   )
-
-class uniform_unit_cylinder(object):
-    def __init__(self,targetN, base_grid=None):
-        cube_cylinder_ratio=numpy.pi*0.5**2
-        self.targetN=targetN
-        self.estimatedN=targetN/cube_cylinder_ratio
+    a = 8 * (numpy.pi - 3) / 3 * numpy.pi * (4 - numpy.pi)
+
+    return numpy.sign(x) * numpy.sqrt(
+        numpy.sqrt(
+            (2 / numpy.pi / a + numpy.log(1 - x**2) / 2)**2 - numpy.log(1 - x**2) / a
+            ) - (2 / numpy.pi / a + numpy.log(1 - x**2) / 2)
+    )
+
+class uniform_unit_cylinder:
+    def __init__(self, targetN, base_grid=None):
+        cube_cylinder_ratio = numpy.pi*0.5**2
+        self.targetN = targetN
+        self.estimatedN = targetN / cube_cylinder_ratio
+
         if base_grid is None:
-            self.base_grid=uniform_random_unit_cube
+            self.base_grid = uniform_random_unit_cube
         else:
-            self.base_grid=base_grid
-   
-    def cutout_cylinder(self,x,y,z):
-        r=x**2+y**2
-        selection=r < numpy.ones_like(r)        
-        x=x.compress(selection)
-        y=y.compress(selection)
-        z=z.compress(selection)
-        return x,y,z
+            self.base_grid = base_grid
+
+    def cutout_cylinder(self, x, y, z):
+        r = x**2 + y**2
+        selection = r < numpy.ones_like(r)
+        x = x.compress(selection)
+        y = y.compress(selection)
+        z = z.compress(selection)
+        return x, y, z
 
     def make_xyz(self):
-        if(self.base_grid==uniform_random_unit_cube):
-            estimatedN=self.estimatedN
-            x=[]
+        if (self.base_grid == uniform_random_unit_cube):
+            estimatedN = self.estimatedN
+            x = []
             while len(x) < self.targetN:
-                estimadedN=estimatedN*1.1+1
-                x,y,z=self.cutout_cylinder(*(self.base_grid(estimatedN)).make_xyz())
-            return x[0:self.targetN],y[0:self.targetN],z[0:self.targetN]  
+                estimadedN = estimatedN * 1.1 + 1
+                x, y, z = self.cutout_cylinder(*(self.base_grid(estimatedN)).make_xyz())
+
+            return x[0:self.targetN], y[0:self.targetN], z[0:self.targetN]
         else:
             return self.cutout_cylinder(*(self.base_grid(self.estimatedN)).make_xyz())
 
 
 class ProtoPlanetaryDisk:
-
     def __init__(
         self, targetN, convert_nbody=None, discfraction=0.1,
-        densitypower=1., thermalpower=0.5, radius_min=1, radius_max=100,
-        gamma=1.,q_out=2.,base_grid=None, Rmin=None, Rmax=None,
+        densitypower=1., thermalpower=0.5, radius_min=None, radius_max=None,
+        gamma=1., q_out=2., base_grid=None, Rmin=None, Rmax=None,
     ):
-        if Rmin is not None:
-            warnings.warn(
-                "Rmin is deprecated, use radius_min instead",
-                category=FutureWarning,
-            )
-            if radius_min is not None and radius_min != Rmin:
-                raise ValueError(
-                    "Rmin and radius_min have different values, "
-                    "this is only allowed if one of them is None"
-                )
-            radius_min = Rmin
-        if radius_min is None:
-            raise ValueError("radius_min must be set")
-        if Rmax is not None:
-            warnings.warn(
-                "Rmax is deprecated, use radius_max instead",
-                category=FutureWarning,
-            )
-            if radius_max is not None and radius_max != Rmax:
-                raise ValueError(
-                    "Rmax and radius_max have different values, "
-                    "this is only allowed if one of them is None"
-                )
-            radius_max = Rmax
-        if radius_max is None:
-            raise ValueError("radius_max must be set")
-
-        self.targetN=targetN
-        self.convert_nbody=convert_nbody
-        self.densitypower=densitypower
-        self.thermalpower=thermalpower
-        self.Rmin=radius_min
-        self.Rmax=radius_max
-        self.gamma=gamma
-        self.q_out=q_out
-        self.discfraction=discfraction
-
-        self.a=self.thermalpower
-        self.a2=self.thermalpower/2
-        self.g=densitypower
-        self.g2=2-densitypower
-        self.k_out=((1+discfraction)/Rmax**3)**0.5
-        self.sigma_out=self.g2*discfraction/(2*numpy.pi*Rmax**self.g*(Rmax**self.g2-Rmin**self.g2))
-        self.cs_out=self.q_out*numpy.pi*self.sigma_out/self.k_out
-        
-        self.base_cylinder=uniform_unit_cylinder(targetN,base_grid)   
-
-
-    def sigma(self,r):
-        return self.sigma_out*(self.Rmax/r)**self.g
-
-    def csound(self,r):
-        return self.cs_out*(self.Rmax/r)**self.a2
-
-    def cmass(self,r):
-        return self.discfraction*(r**self.g2-self.Rmin**self.g2)/(self.Rmax**self.g2-self.Rmin**self.g2)
-
-    def mass_encl(self,r):
-        return 1+self.cmass(r)
-
-    def kappa(self,r):
-        return (self.mass_encl(r)/r**3)**0.5
-  
-    def toomreQ(self,r):
-        return self.csound(r)*self.kappa(r)/numpy.pi/self.sigma(r)
-
-    def getradius(self,f):
-        return ((self.Rmax**self.g2-self.Rmin**self.g2)*f+self.Rmin**self.g2)**(1./self.g2)
-
-    def zscale(self,r):
-        return self.csound(r)/self.kappa(r) 
-
-    def u(self,r):
-        if self.gamma ==1.:
+        self.targetN = targetN
+        self.convert_nbody = convert_nbody
+        self.densitypower = densitypower
+        self.thermalpower = thermalpower
+        self.Rmin = rename_fn_par("radius_min", radius_min, "Rmin", Rmin, 1)
+        self.Rmax = rename_fn_par("radius_max", radius_max, "Rmax", Rmax, 100)
+        self.gamma = gamma
+        self.q_out = q_out
+        self.discfraction = discfraction
+
+        self.a = self.thermalpower
+        self.a2 = self.thermalpower/2
+        self.g = densitypower
+        self.g2 = 2 - densitypower
+        self.k_out = ((1 + discfraction) / self.Rmax**3)**0.5
+        self.sigma_out = self.g2 * discfraction / (
+                2 * numpy.pi * self.Rmax**self.g *
+                (self.Rmax**self.g2 - self.Rmin**self.g2))
+        self.cs_out = self.q_out * numpy.pi * self.sigma_out / self.k_out
+
+        self.base_cylinder = uniform_unit_cylinder(targetN, base_grid)
+
+    @property
+    def radius_min(self):
+        return self.Rmin
+
+    @property
+    def radius_max(self):
+        return self.Rmax
+
+    def sigma(self, r):
+        return self.sigma_out * (self.Rmax / r)**self.g
+
+    def csound(self, r):
+        return self.cs_out * (self.Rmax / r)**self.a2
+
+    def cmass(self, r):
+        return self.discfraction * (r**self.g2 - self.Rmin**self.g2) / (
+                self.Rmax**self.g2 - self.Rmin**self.g2)
+
+    def mass_encl(self, r):
+        return 1 + self.cmass(r)
+
+    def kappa(self, r):
+        return (self.mass_encl(r) / r**3)**0.5
+
+    def toomreQ(self, r):
+        return self.csound(r) * self.kappa(r) / numpy.pi / self.sigma(r)
+
+    def getradius(self, f):
+        return (
+                (self.Rmax**self.g2 - self.Rmin**self.g2) * f + self.Rmin**self.g2)**(
+                        1./self.g2)
+
+    def zscale(self, r):
+        return self.csound(r) / self.kappa(r)
+
+    def u(self, r):
+        if self.gamma == 1.:
             return self.csound(r)**2
         else:
-            return self.csound(r)**2/(self.gamma-1)
+            return self.csound(r)**2 / (self.gamma - 1)
+
+    def vcirc(self, r):
+        return (self.mass_encl(r) / r)**0.5
 
-    def vcirc(self,r):
-        return (self.mass_encl(r)/r)**0.5
-    
     def new_model(self):
-        x,y,z=self.base_cylinder.make_xyz()
-        self.actualN=len(x)
-        f=x**2+y**2
-        r=f**0.5
-        rtarget=self.getradius(f)
-        
-        mass=self.discfraction*numpy.ones_like(x)/self.actualN
-        internal_energy=self.u(rtarget)
-        zscale=self.zscale(rtarget)
-        r=r.clip(1.e-8,2.)
-        x=x/r
-        y=y/r
-
-        vx=-y*self.vcirc(rtarget)
-        vy=x*self.vcirc(rtarget)
-        vz=numpy.zeros_like(x)
-         
-        x=rtarget*x 
-        y=rtarget*y                  
-        z=approximate_inverse_error_function(z)*zscale*2.**0.5
+        x, y, z = self.base_cylinder.make_xyz()
+        self.actualN = len(x)
+        f = x**2 + y**2
+        r = f**0.5
+        rtarget = self.getradius(f)
+
+        mass = self.discfraction * numpy.ones_like(x) / self.actualN
+        internal_energy = self.u(rtarget)
+        zscale = self.zscale(rtarget)
+        r = r.clip(1.e-8, 2.)
+        x= x / r
+        y= y / r
+
+        vx = -y * self.vcirc(rtarget)
+        vy = x * self.vcirc(rtarget)
+        vz = numpy.zeros_like(x)
+
+        x = rtarget * x
+        y = rtarget * y
+        z = approximate_inverse_error_function(z) * zscale * 2.**0.5
 
         return (mass, x, y, z, vx, vy, vz, internal_energy)
-  
+
     @property
     def result(self):
-        masses, x,y,z, vx,vy,vz, internal_energies = self.new_model()
+        masses, x, y, z, vx, vy, vz, internal_energies = self.new_model()
         result = Particles(self.actualN)
         result.mass = nbody_system.mass.new_quantity(masses)
         result.x = nbody_system.length.new_quantity(x)
@@ -168,10 +160,10 @@ def result(self):
         result.vy = nbody_system.speed.new_quantity(vy)
         result.vz = nbody_system.speed.new_quantity(vz)
         result.u = nbody_system.specific_energy.new_quantity(internal_energies)
-                
+
         if not self.convert_nbody is None:
-            result = ParticlesWithUnitsConverted(result, self.convert_nbody.as_converter_from_si_to_generic())
+            result = ParticlesWithUnitsConverted(
+                    result, self.convert_nbody.as_converter_from_si_to_generic())
             result = result.copy()
-            
-        return result
 
+        return result

src/amuse/support/helpers.py

@@ -0,0 +1,55 @@
+import warnings
+
+
+def rename_fn_par(new_name, new_value, old_name, old_value, default_value):
+    """Get value of a renamed function parameter.
+
+    If you have a function f with parameter ``x`` and default value 1, like this:
+
+    .. code-block::
+
+        def f(x=1):
+            ...
+
+    and you want to rename ``x`` to ``y`` while staying backwards compatible, then you
+    can rename ``x`` to ``y``, add ``x`` back in at the end so that old keyword
+    arguments still work, give both a default value of ``None``, and use this function
+    like this:
+
+    .. code-block::
+
+        def f(y=None, x=None):
+            value = rename_fn_par("y", y, "x", x, 1)
+
+    Any callers using ``x`` explicitly will receive a warning to change their code to
+    use ``y`` in the future. If both ``x`` and ``y`` are set and the values are
+    different, then an exception is raised.
+
+    Args:
+        new_name (str): The new name of the variable
+        new_value: The value passed using the new name
+        old_name (str): The old name of the variable
+        old_value: The value passed using the old name
+        default_value: The default value if neither are set
+
+    Returns:
+        Either new_value or old_value if only one is set or they're set to the same
+        value, or default_value if neither is set.
+
+    Raises:
+        ValueError: If both new_value and old_value are set, and to different values.
+    """
+    if new_value is not None:
+        if old_value is not None and old_value != new_value:
+            raise ValueError(
+                    f"{old_name} and {new_name} have different values,"
+                    " which is not allowed because they represent the same thing.")
+        return new_value
+
+    if old_value is not None:
+        warnings.warn(
+                f"{old_name} is deprecated, please use {new_name} instead",
+                category=FutureWarning)
+        return old_value
+
+    return default_value

src/tests/core_tests/test_helpers.py

@@ -0,0 +1,40 @@
+from amuse.support.helpers import rename_fn_par
+
+import pytest
+
+
+@pytest.mark.filterwarnings("ignore: old is deprecated.*")
+def test_rename_fn_par():
+    assert rename_fn_par("new", None, "old", None, 1) == 1
+    assert rename_fn_par("new", None, "old", 1, 2) == 1
+    assert rename_fn_par("new", 1, "old", None, 2) == 1
+    assert rename_fn_par("new", 1, "old", 1, 2) == 1
+
+    with pytest.raises(ValueError):
+        rename_fn_par("new", 1, "old", 2, 3)
+
+
+class _RenamedMethod:
+    def __init__(self, new = None, a = 1, b = "test", old = None):
+        self.y = rename_fn_par("new", new, "old", old, 42)
+
+
+@pytest.mark.filterwarnings("ignore: old is deprecated.*")
+def test_rename_fn_par_usage():
+    rm = _RenamedMethod()
+    assert rm.y == 42
+
+    rm = _RenamedMethod(1)
+    assert rm.y == 1
+
+    rm = _RenamedMethod(new=1)
+    assert rm.y == 1
+
+    rm = _RenamedMethod(old=1)
+    assert rm.y == 1
+
+    rm = _RenamedMethod(1, 1, "test", 1)
+    assert rm.y == 1
+
+    with pytest.raises(ValueError):
+        _RenamedMethod(1, 2, "test", 3)