def METHOD_NAME(r=None, nr=100):
    """
    Returns the initial current density profile.
    :param r:      Minor radius (m). May be either scalar or numpy.ndarray.
    :param int nr: Radial grid resolution to use if ``r = None``.
    """
    global a, j0
    if r is None:
        r = np.linspace(0, a, nr)
    #return r, j0 * (1 - (r/a)**0.41)
    return r, j0 * (1 - (1-0.001**(1/0.41))*(r/a)**2)**0.41