From 4a5500b85bd732d477e356c304f8da58098c4722 Mon Sep 17 00:00:00 2001 From: Lourens Veen Date: Mon, 9 Mar 2026 13:46:35 +0100 Subject: [PATCH 1/2] Fix formatting I don't have my Fortran safety goggles with me, so this will avoid eye injury and OSHA violations. --- src/amuse/ext/protodisk.py | 222 +++++++++++++++++++------------------ 1 file changed, 117 insertions(+), 105 deletions(-) diff --git a/src/amuse/ext/protodisk.py b/src/amuse/ext/protodisk.py index ce7d00734b..4e2783dd55 100644 --- a/src/amuse/ext/protodisk.py +++ b/src/amuse/ext/protodisk.py @@ -8,48 +8,54 @@ 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, + gamma=1., q_out=2., base_grid=None, Rmin=None, Rmax=None, ): if Rmin is not None: warnings.warn( @@ -62,8 +68,10 @@ def __init__( "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", @@ -75,90 +83,94 @@ def __init__( "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 = 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.: 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 +180,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 From 09d8e90ec49e334653800a9cc6657e2bb1028328 Mon Sep 17 00:00:00 2001 From: Lourens Veen Date: Mon, 9 Mar 2026 18:20:02 +0100 Subject: [PATCH 2/2] Fix protodisk Rmin/radius_min rename --- src/amuse/ext/protodisk.py | 50 ++++++++----------------- src/amuse/support/helpers.py | 55 ++++++++++++++++++++++++++++ src/tests/core_tests/test_helpers.py | 40 ++++++++++++++++++++ 3 files changed, 110 insertions(+), 35 deletions(-) create mode 100644 src/amuse/support/helpers.py create mode 100644 src/tests/core_tests/test_helpers.py diff --git a/src/amuse/ext/protodisk.py b/src/amuse/ext/protodisk.py index 4e2783dd55..30da23a748 100644 --- a/src/amuse/ext/protodisk.py +++ b/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 @@ -54,45 +55,15 @@ def make_xyz(self): class ProtoPlanetaryDisk: def __init__( self, targetN, convert_nbody=None, discfraction=0.1, - densitypower=1., thermalpower=0.5, radius_min=1, radius_max=100, + 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.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 @@ -101,13 +72,22 @@ def __init__( self.a2 = self.thermalpower/2 self.g = densitypower self.g2 = 2 - densitypower - self.k_out = ((1 + discfraction) / Rmax**3)**0.5 + self.k_out = ((1 + discfraction) / self.Rmax**3)**0.5 self.sigma_out = self.g2 * discfraction / ( - 2 * numpy.pi * Rmax**self.g * (Rmax**self.g2 - Rmin**self.g2)) + 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 diff --git a/src/amuse/support/helpers.py b/src/amuse/support/helpers.py new file mode 100644 index 0000000000..0d1f67d2d4 --- /dev/null +++ b/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 diff --git a/src/tests/core_tests/test_helpers.py b/src/tests/core_tests/test_helpers.py new file mode 100644 index 0000000000..41d461ace9 --- /dev/null +++ b/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)