docs: def mete_upscale_iterative_alt(S, N, doublings):
      """
      This function is used to upscale from the anchor area.
  
      Parameters
      ----------
      S : int or float
          Number of species at anchor scale
      N : int or float
          Number of individuals at anchor scale
      doublings : int
          Number of doublings of A. Result vector will be length doublings + 1.
  
      Returns
      -------
      result : ndarray
          1D array of number of species at each doubling
  
      """
  
      # Arrays to store N and S at all doublings
      n_arr = np.empty(doublings+1)
      s_arr = np.empty(doublings+1)
  
      # Loop through all scales
      for i in xrange(doublings+1):
  
          # If this is first step (doubling 0), N and S are initial values
          if i == 0:
              n_arr[i] = N
              s_arr[i] = S
  
          # If not first step
          else:
  
              # Get previous S
              SA = s_arr[i-1]
  
              # N is double previous N
              n_arr[i] = 2 * n_arr[i-1]
              N2A = n_arr[i]
  
              # Eq 8 from Harte 2009, setup to return S2A given input of x
              # x is exp(-lam_phi, 2A)
              def S2A_calc(x, SA, N2A):
                  return ((SA +
                          N2A *
                          (1-x)/(x-x**(N2A+1)) *
                          (1 - (x**N2A)/(N2A+1))) /
                          x**-1)
  
              # Eq 9 from Harte 2009, setup to equal to zero, used to solve x
              # Note that two summations are replaced by known formulas for sum
              # of geometric and logarithmic series.
              # Note "offset" of 1e-23, which is needed because f(a) and f(b) do
              # not have the same sign in solver otherwise. This introduces no
              # more than a 1e-23 error in the calculation of x, which should not
              # cause a significant problem.
              def x_calc(x, SA, N2A):
                  return (S2A_calc(x,SA,N2A) /
                          N2A *
                          x*(x**N2A-1)/(x-1) -
                          (x**N2A * (-lerchphi(x,1,N2A+1))-np.log(1-x)) ) - 1e-23
  
              # Solve for x
              x = (optimize.brentq(x_calc, 1e-24, 1-1e-16, args=(SA,N2A),
                      xtol=1e-16, maxiter=1000, disp=True) + 1e-23)
  
              # Given x, calculate S2A
              s_arr[i] = S2A_calc(x,SA,N2A)
  
      return s_arr