2 files changed, 40 insertions(+), 6 deletions(-)

M clusterfunc.py
M test_clusterfunc.py

@@ 55,6 55,9 @@ def unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
         4) Not enough points to reach unit length at full stretch
         """
 
+    if N < 2:
+        raise ClusteringException(f'Need at least two points to cluster, N={N}')
+
     if Dmax <= Dmin:
         raise ClusteringException(f'Dmax={Dmax} should be > Dmin={Dmin}')
 


@@ 76,25 79,29 @@ def unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
         f'Not enough points N={N} to reach unit length for ERmax={ERmax}, Dmin={Dmin}, regardless of Dmax., L={Lfull}, Nfull={Nfull}'
         )
 
-
     # Assume we do not need to cap the spacing
     # What expansion ratio is needed, and what is the max spacing?
     Nm1, Nm2, Nm3 = N-1, N-2, N-3
+    Np1 = N+1
     Nm1Nm2 = Nm1*Nm2
     def _iter_ER(ERi):
-        y = ERi - 1. + Dmin*(1.-ERi**Nm1)
-        dy = 1. - Dmin*Nm1*ERi**Nm2
-        d2y = -Dmin*Nm1Nm2*ERi**Nm3
-        return y, dy, d2y
+        # y = ERi - 1. + Dmin*(1.-ERi**Nm1)
+        # dy = 1. - Dmin*Nm1*ERi**Nm2
+        # d2y = -Dmin*Nm1Nm2*ERi**Nm3
+        y = Dmin*(1.-ERi**Nm1)/(1.-ERi) - 1.
+        dy = Dmin*(Nm2*ERi**Np1 - Nm1*ERi**N + ERi**2)/(1.-ERi)**2/ERi**2
+        return y, dy#, d2y
 
     # If we cannot reach irregardless of capping, give up
+    # print(N,_iter_ER(ERmax)[0])
+    # quit()
     if _iter_ER(ERmax)[0] > 0.:
         raise ClusteringException(
         f'Insufficient points N={N} to reach unit length with Dmin={Dmin}, ERmax={ERmax}, no capping to Dmax needed.'
         )
 
     # Find root 
-    soln = root_scalar(_iter_ER, x0=ERmax, fprime2=True, maxiter=200)
+    soln = root_scalar(_iter_ER, x0=ERmax, fprime=True, maxiter=200)
     assert soln.converged
     ERnocap = soln.root
     Dnocap = Dmin*ERnocap**Nm2


@@ 194,7 201,25 @@ def unit_free_N(Dmin, Dmax, ERmax, mult=8, check=True):
     return x
 
 
+def fixed_N(x0, x1, Dmin, Dmax, ERmax, N, check=True):
+    """Single-sided clustering between two values with with fixed number of points."""
+    return x0 + (x1-x0) * unit_fixed_N(Dmin, Dmax, ERmax, N, check)
+
+def free_N(x0, x1, Dmin, Dmax, ERmax, mult=8, check=True):
+    """Single-sided clustering between two values with with free number of points."""
+    dx = x1-x0
+    return x0 + (x1-x0) * unit_free_N(Dmin/dx, Dmax/dx, ERmax, mult, check)
+
+def symmetric_unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
+    Na = N//2
+    Nb = N - Na + 1
+    xa = fixed_N(0., 0.5, Dmin, Dmax, ERmax, Na, check)
+    xb = np.flip(fixed_N(1., 0.5, Dmin, Dmax, ERmax, Na, check))
+    x = np.concatenate((xa, xb[1:]))
+    assert len(x) == N
+
 
 def cumsum0(x, axis=None):
     return np.insert(np.cumsum(x, axis=axis), 0, 0.0, axis=axis)
 
+# symmetric_unit_fixed_N(1e-3,1e-2, 1.2, 1000)

@@ 21,6 21,15 @@ def test_cluster():
                 for N in range(Nmin,Nmax+1):
                     clusterfunc.unit_fixed_N(Dmin, Dmax, ER, N)
 
+def test_unif():
+
+    N = np.arange(2, 101).astype(int)
+    Dmax = 2.
+    ER = 1.2
+    for n in N:
+        Dmin = 1./(n-1)
+        clusterfunc.unit_fixed_N(Dmin, Dmax, ER, n)
 
 if __name__=='__main__':
     test_cluster()
+    # test_unif()