@@ -4854,6 +4854,74 @@ def test_NFW_virialsetup_wrtcrit():
48544854 return None
48554855
48564856
4857+ def test_EllipsoidalPotential_integration_scalar_and_array ():
4858+ # Test that the GL integration helpers for EllipsoidalPotential work
4859+ # correctly for both scalar and array inputs; this ensures that the
4860+ # scalar Python path is covered even when C extensions handle scalar
4861+ # calls in the main potential evaluation.
4862+ from galpy .potential .EllipsoidalPotential import (
4863+ _2ndDerivInt ,
4864+ _forceInt ,
4865+ _potInt ,
4866+ )
4867+
4868+ tp = potential .TriaxialNFWPotential (normalize = 1.0 , b = 0.7 , c = 0.5 )
4869+ glx , glw = tp ._glx , tp ._glw
4870+ b2 , c2 = tp ._b2 , tp ._c2
4871+ # Scalar inputs
4872+ x_s , y_s , z_s = 1.0 , 0.5 , 0.3
4873+ pot_s = _potInt (x_s , y_s , z_s , tp ._psi , b2 , c2 , glx = glx , glw = glw )
4874+ force_s = _forceInt (x_s , y_s , z_s , tp ._mdens , b2 , c2 , 0 , glx = glx , glw = glw )
4875+ deriv_s = _2ndDerivInt (
4876+ x_s , y_s , z_s , tp ._mdens , tp ._mdens_deriv , b2 , c2 , 0 , 0 , glx = glx , glw = glw
4877+ )
4878+ # Array inputs
4879+ x_a = numpy .array ([1.0 , 2.0 ])
4880+ y_a = numpy .array ([0.5 , 1.0 ])
4881+ z_a = numpy .array ([0.3 , 0.6 ])
4882+ pot_a = _potInt (x_a , y_a , z_a , tp ._psi , b2 , c2 , glx = glx , glw = glw )
4883+ force_a = _forceInt (x_a , y_a , z_a , tp ._mdens , b2 , c2 , 0 , glx = glx , glw = glw )
4884+ deriv_a = _2ndDerivInt (
4885+ x_a , y_a , z_a , tp ._mdens , tp ._mdens_deriv , b2 , c2 , 0 , 0 , glx = glx , glw = glw
4886+ )
4887+ # Scalar result should match first element of array result
4888+ assert numpy .fabs (pot_s - pot_a [0 ]) < 1e-10 , (
4889+ "Scalar and array _potInt results disagree"
4890+ )
4891+ assert numpy .fabs (force_s - force_a [0 ]) < 1e-10 , (
4892+ "Scalar and array _forceInt results disagree"
4893+ )
4894+ assert numpy .fabs (deriv_s - deriv_a [0 ]) < 1e-10 , (
4895+ "Scalar and array _2ndDerivInt results disagree"
4896+ )
4897+ return None
4898+
4899+
4900+ def test_EllipsoidalPotential_evaluate_array_inf ():
4901+ # Test that _evaluate handles array inputs with inf correctly
4902+ tp = potential .PerfectEllipsoidPotential (normalize = 1.0 , b = 0.7 , c = 0.5 )
4903+ # Scalar inf
4904+ val_inf = tp ._evaluate (numpy .inf , 0.0 )
4905+ # Array with inf
4906+ R_arr = numpy .array ([numpy .inf , 1.0 ])
4907+ z_arr = numpy .array ([0.0 , 0.1 ])
4908+ val_arr = tp ._evaluate (R_arr , z_arr )
4909+ assert numpy .fabs (val_inf - val_arr [0 ]) < 1e-10 , (
4910+ "Scalar and array inf _evaluate results disagree"
4911+ )
4912+ # Non-aligned with finite array
4913+ tp2 = potential .PerfectEllipsoidPotential (
4914+ normalize = 1.0 , b = 0.7 , c = 0.5 , zvec = [0.0 , numpy .sin (0.3 ), numpy .cos (0.3 )]
4915+ )
4916+ R_fin = numpy .array ([0.5 , 1.0 , 2.0 ])
4917+ z_fin = numpy .array ([0.1 , 0.2 , 0.3 ])
4918+ val_na = tp2 ._evaluate (R_fin , z_fin , phi = 0.5 )
4919+ assert not numpy .any (numpy .isnan (val_na )), (
4920+ "Non-aligned array _evaluate returned NaN"
4921+ )
4922+ return None
4923+
4924+
48574925def test_TriaxialNFW_virialsetup_wrtmeanmatter ():
48584926 H , Om , overdens , wrtcrit = 71.0 , 0.32 , 201.0 , False
48594927 ro , vo = 220.0 , 8.0
0 commit comments