Repository: PlasmaPy/PlasmaPy

Issue:
Drop support for Python 3.7
[NEP 29](https://numpy.org/neps/nep-0029-deprecation_policy.html) suggests a community schedule for dropping support old versions of Python & NumPy.  If we follow that schedule, we would drop support for Python 3.7 in PlasmaPy 0.8.0 which will likely be released in early 2022.  Python 3.10 is due to be released in ≲1 month, so we would be supporting the last three releases of Python.  Besides that, the next release of Astropy will not include support for Python 3.7, so we'd need to update our requirements fairly soon anyway.

The most significant changes in Python 3.8 are:
 - [f-string support for self-documenting expressions](https://docs.python.org/3.8/whatsnew/3.8.html#f-strings-support-for-self-documenting-expressions-and-debugging) (I'll probably start using these immediately)
 - [walrus operators](https://docs.python.org/3.8/whatsnew/3.8.html#assignment-expressions) (will likely improve readability here & there)
 - [positional-only parameters](https://docs.python.org/3.8/whatsnew/3.8.html#positional-only-parameters) (maybe appropriate but not necessary for mathematical stuff like the plasma dispersion function?)

Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1465?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) Report
> Merging [#1465](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1465?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) (327aa75) into [main](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/319bffe9bd91135cc6d42dfb108eae6077a7b54d?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) (319bffe) will **not change** coverage.
> The diff coverage is `50.00%`.

> :exclamation: Current head 327aa75 differs from pull request most recent head a3db08d. Consider uploading reports for the commit a3db08d to get more accurate results

[![Impacted file tree graph](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1465/graphs/tree.svg?width=650&height=150&src=pr&token=LzZ5mnXVad&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy)](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1465?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy)

```diff
@@           Coverage Diff           @@
##             main    #1465   +/-   ##
=======================================
  Coverage   96.49%   96.49%           
=======================================
  Files          82       82           
  Lines        7857     7857           
=======================================
  Hits         7582     7582           
  Misses        275      275           
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1465?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) | Coverage Δ | |
|---|---|---|
| [plasmapy/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1465/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy#diff-cGxhc21hcHkvX19pbml0X18ucHk=) | `34.48% <50.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1465?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1465?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy). Last update [319bffe...a3db08d](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1465?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy).

namurphy:  - [ ] Look into what is happening with the h5py in the Python 3.8 minimal imports behavior...
namurphy: I've been looking into the problem of why the tests with minimal dependencies were failing.  I tried requiring numpy ≥ 1.19.0, but the minimal dependencies tests didn't pass.  I removed the pinned version of h5py in the `tox.ini` env for minimal dependencies, and the tests passed (with h5py 3.6.x being used).  I'm thinking we'll need to specify a minimum version of h5py.  I'm trying h5py ≥ 3.5 now.
namurphy: Ready for review as long as this round of tests passes!


namurphy: One thing I'm noticed in this PR is that the tests with minimal dependencies seem to be taking a particularly long time to run.  I do not know why.  If that continues to be the case, it'd probably help to start them at the same time as the initial tests, but not require them to pass before the second round of tests is done.  
rocco8773: A note about `h5py` `v2.x` vs `v3.x`, right around the change over they made some kind of change that affected dataset indexing.  This affected me over a `bapsflib` since I was indexing something like `[i, j]` and in `v3.x` it got swapped to `[j, i]`.  I never discovered the root cause though.  With all that said, requiring `h5py>=3.0` is a good thing to do.
namurphy: Thanks!  I think I just addressed all of your suggestions.

Patch:
diff --git a/README.md b/README.md
index 89f57f46..4cef3119 100644
--- a/README.md
+++ b/README.md
@@ -25,7 +25,7 @@ # PlasmaPy
 [![Open Source Helpers](https://www.codetriage.com/plasmapy/plasmapy/badges/users.svg)](https://www.codetriage.com/plasmapy/plasmapy)
 
 [PlasmaPy](https://www.plasmapy.org/) is an open source, community-developed
-Python 3.7+ package for plasma science. PlasmaPy intends to be for plasma
+Python 3.8+ package for plasma science. PlasmaPy intends to be for plasma
 science what [Astropy](https://github.com/astropy/astropy) is for astronomy
 — a collection of functionality commonly used and shared between plasma
 scientists and researchers globally, running within and leveraging the
diff --git a/changelog/1427.trivial.rst b/changelog/1427.trivial.rst
new file mode 100644
index 00000000..d19c3b98
--- /dev/null
+++ b/changelog/1427.trivial.rst
@@ -0,0 +1,2 @@
+Added unit test cases for manual entry of vector values in order to improve code coverage
+in the null point finder.
diff --git a/changelog/1465.removal.rst b/changelog/1465.removal.rst
new file mode 100644
index 00000000..f7d59254
--- /dev/null
+++ b/changelog/1465.removal.rst
@@ -0,0 +1,3 @@
+Dropped support for Python 3.7 in accordance with the deprecation policy
+laid out in `NumPy Enhancement Proposal 29
+<https://numpy.org/neps/nep-0029-deprecation_policy.html>`__.
diff --git a/changelog/1465.trivial.1.rst b/changelog/1465.trivial.1.rst
new file mode 100644
index 00000000..ddeda7f9
--- /dev/null
+++ b/changelog/1465.trivial.1.rst
@@ -0,0 +1,3 @@
+Changed the raised exception to `ImportError` (from a general `Exception`)
+when attempting to import `plasmapy` from a Python version below the
+minimum supported version.
diff --git a/changelog/1465.trivial.2.rst b/changelog/1465.trivial.2.rst
new file mode 100644
index 00000000..0a4fc556
--- /dev/null
+++ b/changelog/1465.trivial.2.rst
@@ -0,0 +1,1 @@
+Bumped the minimum version of `h5py` to ``3.0.0``.
diff --git a/codemeta.json b/codemeta.json
index 4e8fca2c..422cce0e 100644
--- a/codemeta.json
+++ b/codemeta.json
@@ -39,15 +39,15 @@
         "science"
     ],
     "programmingLanguage": [
-        "Python 3.7",
         "Python 3.8",
-        "Python 3.9"
+        "Python 3.9",
+        "Python 3.10"
     ],
     "operatingSystem": [
         "OS Independent"
     ],
     "softwareRequirements": [
-        "python >= 3.7",
+        "python >= 3.8",
         "astropy >= 4.0",
         "h5py >= 2.8",
         "lmfit >= 1.0.1",
diff --git a/docs/conf.py b/docs/conf.py
index f4905350..f62ac877 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python3.7
+#!/usr/bin/env python3
 #
 # PlasmaPy documentation build configuration file, created by
 # sphinx-quickstart on Wed May 31 18:16:46 2017.
diff --git a/docs/development/code_guide.rst b/docs/development/code_guide.rst
index 3e23f1ed..e23e83d8 100644
--- a/docs/development/code_guide.rst
+++ b/docs/development/code_guide.rst
@@ -7,7 +7,7 @@ Code Development Guidelines
 This document describes the coding requirements and guidelines to be
 followed during the development of PlasmaPy and affiliated packages.
 
-Code written for PlasmaPy must be compatible with Python 3.7 and
+Code written for PlasmaPy must be compatible with Python 3.8 and
 later.
 
 Coding Style
diff --git a/docs/index.rst b/docs/index.rst
index 05b3cb9f..7b0b1e42 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -11,7 +11,7 @@
 PlasmaPy Documentation
 ######################
 
-PlasmaPy_ is an open source community-developed core Python_ 3.7+
+PlasmaPy_ is an open source community-developed core Python_ 3.8+
 package for plasma physics currently under development.
 
 Example highlights
diff --git a/docs/install.rst b/docs/install.rst
index 59e7f719..5ff0a237 100644
--- a/docs/install.rst
+++ b/docs/install.rst
@@ -15,7 +15,7 @@ Installing PlasmaPy
 Requirements
 ============
 
-PlasmaPy requires Python version 3.7 or newer.
+PlasmaPy requires Python version 3.8 or newer.
 PlasmaPy requires the following packages for installation:
 
 - `NumPy`_ — 1.18.1 or newer
diff --git a/plasmapy/__init__.py b/plasmapy/__init__.py
index 2d203b12..4f59583c 100644
--- a/plasmapy/__init__.py
+++ b/plasmapy/__init__.py
@@ -20,8 +20,8 @@
 # This is the same check as the one at the top of setup.py
 import sys
 
-if sys.version_info < (3, 7):
-    raise Exception("PlasmaPy does not support Python < 3.7")
+if sys.version_info < (3, 8):  # coverage: ignore
+    raise ImportError("PlasmaPy does not support Python < 3.8")
 
 # Packages may add whatever they like to this file, but
 # should keep this content at the top.
diff --git a/plasmapy/analysis/nullpoint.py b/plasmapy/analysis/nullpoint.py
index 07a7f330..6e0fc0e8 100644
--- a/plasmapy/analysis/nullpoint.py
+++ b/plasmapy/analysis/nullpoint.py
@@ -11,7 +11,7 @@
 import warnings
 
 # Declare Constants & global variables
-_ATOL = 1e-10
+_EQUALITY_ATOL = 1e-15
 _MAX_RECURSION_LEVEL = 10
 global _recursion_level
 _recursion_level = 0
@@ -54,11 +54,12 @@ def __eq__(self, point):
         Returns True if two null point objects have the same coordinates.
         False otherwise.
         """
-        return (
-            np.isclose(self.loc()[0], point.loc()[0], atol=_ATOL)
-            and np.isclose(self.loc()[1], point.loc()[1], atol=_ATOL)
-            and np.isclose(self.loc()[2], point.loc()[2], atol=_ATOL)
+        d = np.sqrt(
+            (self.loc[0] - point.loc[0]) ** 2
+            + (self.loc[1] - point.loc[1]) ** 2
+            + (self.loc[2] - point.loc[2]) ** 2
         )
+        return np.isclose(d, 0, atol=_EQUALITY_ATOL)
 
 
 def _vector_space(
@@ -456,9 +457,9 @@ def _reduction(vspace, cell):
     passY = False
     passZ = False
     # Check reduction criteria
-    sign_x = np.sign(u[cell[0], cell[1], cell[2]])
-    sign_y = np.sign(v[cell[0], cell[1], cell[2]])
-    sign_z = np.sign(w[cell[0], cell[1], cell[2]])
+    sign_x = np.sign(u[cell[0]][cell[1]][cell[2]])
+    sign_y = np.sign(v[cell[0]][cell[1]][cell[2]])
+    sign_z = np.sign(w[cell[0]][cell[1]][cell[2]])
     for point in corners:
         if (
             u[point[0]][point[1]][point[2]] == 0
@@ -514,41 +515,43 @@ def _bilinear_root(a1, b1, c1, d1, a2, b2, c2, d2):
     b = np.linalg.det(m2) + np.linalg.det(m3)
     c = np.linalg.det(m1)
 
-    if np.isclose(a, 0, atol=_ATOL):
-        if np.isclose(b, 0, atol=_ATOL):
-            return None, None
+    if np.isclose(a, 0, atol=_EQUALITY_ATOL):
+        if np.isclose(b, 0, atol=_EQUALITY_ATOL):
+            return np.array([])
         else:
             x1 = (-1.0 * c) / b
             x2 = (-1.0 * c) / b
+
     else:
         if (b ** 2 - 4.0 * a * c) < 0:
-            return None, None
+            return np.array([])
         else:
             x1 = (-1.0 * b + (b ** 2 - 4.0 * a * c) ** 0.5) / (2.0 * a)
             x2 = (-1.0 * b - (b ** 2 - 4.0 * a * c) ** 0.5) / (2.0 * a)
 
-    m1 = np.array([[a1, a2], [b1, b2]])
-    m2 = np.array([[a1, a2], [d1, d2]])
-    m3 = np.array([[b1, b2], [c1, c2]])
-    m4 = np.array([[c1, c2], [d1, d2]])
-    a = np.linalg.det(m4)
-    b = np.linalg.det(m2) - np.linalg.det(m3)
-    c = np.linalg.det(m1)
-
-    if np.isclose(a, 0, atol=_ATOL):
-        if np.isclose(b, 0, atol=_ATOL):
-            return None, None
-        else:
-            y1 = (-1.0 * c) / b
-            y2 = (-1.0 * c) / b
+    y1 = None
+    y2 = None
+    if not (np.isclose((c1 + d1 * x1), 0, atol=_EQUALITY_ATOL)):
+        y1 = (-a1 - b1 * x1) / (c1 + d1 * x1)
+    elif not (np.isclose((c2 + d2 * x1), 0, atol=_EQUALITY_ATOL)):
+        y1 = (-a2 - b2 * x1) / (c2 + d2 * x1)
+    if not (np.isclose((c1 + d1 * x2), 0, atol=_EQUALITY_ATOL)):
+        y2 = (-a1 - b1 * x2) / (c1 + d1 * x2)
+    elif not (np.isclose((c2 + d2 * x2), 0, atol=_EQUALITY_ATOL)):
+        y2 = (-a2 - b2 * x2) / (c2 + d2 * x2)
+
+    if y1 is None and y2 is None:
+        return np.array([])
+    elif y1 is None:
+        return np.array([(x2, y2)])
+    elif y2 is None:
+        return np.array([(x1, y1)])
     else:
-        if (b ** 2 - 4.0 * a * c) < 0:
-            return None, None
+        d = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
+        if np.isclose(d, 0, atol=_EQUALITY_ATOL):
+            return np.array([(x1, y1)])
         else:
-            y1 = (-1.0 * b - (b ** 2 - 4.0 * a * c) ** 0.5) / (2.0 * a)
-            y2 = (-1.0 * b + (b ** 2 - 4.0 * a * c) ** 0.5) / (2.0 * a)
-
-    return [(x1, y1), (x2, y2)]
+            return np.array([(x1, y1), (x2, y2)])
 
 
 def _trilinear_analysis(vspace, cell):
@@ -588,16 +591,6 @@ def _trilinear_analysis(vspace, cell):
         that a grid cell may contain more than one nullpoint.
     """
 
-    # Helper Function
-    def is_close(a, b):
-        arr = np.isclose(a, b, atol=_ATOL)
-        if type(arr) == np.bool_:
-            return arr
-        res = True
-        for b in arr:
-            res = res and b
-        return res
-
     # Critical Cell Corners
     f000 = cell
     f111 = [cell[0] + 1, cell[1] + 1, cell[2] + 1]
@@ -627,7 +620,7 @@ def is_close(a, b):
     ]  # y-coordinate of the front surface
     # Bx=By=0 Curve Endpoint
 
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ax + cx * yConst1,
         bx + ex * yConst1,
         dx + gx * yConst1,
@@ -637,16 +630,11 @@ def is_close(a, b):
         dy + gy * yConst1,
         fy + hy * yConst1,
     )
-
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxByEndpoints.append((root1[0], yConst1, root1[1]))
-        else:
-            BxByEndpoints.append((root1[0], yConst1, root1[1]))
-            BxByEndpoints.append((root2[0], yConst1, root2[1]))
+    for root in root_list:
+        BxByEndpoints.append((root[0], yConst1, root[1]))
 
     # Bx=BZ=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ax + cx * yConst1,
         bx + ex * yConst1,
         dx + gx * yConst1,
@@ -656,14 +644,12 @@ def is_close(a, b):
         dz + gz * yConst1,
         fz + hz * yConst1,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxBzEndpoints.append((root1[0], yConst1, root1[1]))
-        else:
-            BxBzEndpoints.append((root1[0], yConst1, root1[1]))
-            BxBzEndpoints.append((root2[0], yConst1, root2[1]))
+
+    for root in root_list:
+        BxBzEndpoints.append((root[0], yConst1, root[1]))
+
     # By=Bz=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ay + cy * yConst1,
         by + ey * yConst1,
         dy + gy * yConst1,
@@ -673,19 +659,16 @@ def is_close(a, b):
         dz + gz * yConst1,
         fz + hz * yConst1,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            ByBzEndpoints.append((root1[0], yConst1, root1[1]))
-        else:
-            ByBzEndpoints.append((root1[0], yConst1, root1[1]))
-            ByBzEndpoints.append((root2[0], yConst1, root2[1]))
+
+    for root in root_list:
+        ByBzEndpoints.append((root[0], yConst1, root[1]))
 
     # Back Surface
     yConst2 = vspace[0][1][f111[0]][f111[1]][
         f111[2]
     ]  # y-coordinate of the front surface
     # Bx=By=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ax + cx * yConst2,
         bx + ex * yConst2,
         dx + gx * yConst2,
@@ -695,14 +678,12 @@ def is_close(a, b):
         dy + gy * yConst2,
         fy + hy * yConst2,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxByEndpoints.append((root1[0], yConst2, root1[1]))
-        else:
-            BxByEndpoints.append((root1[0], yConst2, root1[1]))
-            BxByEndpoints.append((root2[0], yConst2, root2[1]))
+
+    for root in root_list:
+        BxByEndpoints.append((root[0], yConst2, root[1]))
+
     # Bx=Bz=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ax + cx * yConst2,
         bx + ex * yConst2,
         dx + gx * yConst2,
@@ -712,14 +693,12 @@ def is_close(a, b):
         dz + gz * yConst2,
         fz + hz * yConst2,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxBzEndpoints.append((root1[0], yConst2, root1[1]))
-        else:
-            BxBzEndpoints.append((root1[0], yConst2, root1[1]))
-            BxBzEndpoints.append((root2[0], yConst2, root2[1]))
+
+    for root in root_list:
+        BxBzEndpoints.append((root[0], yConst2, root[1]))
+
     # By=Bz=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ay + cy * yConst2,
         by + ey * yConst2,
         dy + gy * yConst2,
@@ -729,17 +708,14 @@ def is_close(a, b):
         dz + gz * yConst2,
         fz + hz * yConst2,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            ByBzEndpoints.append((root1[0], yConst2, root1[1]))
-        else:
-            ByBzEndpoints.append((root1[0], yConst2, root1[1]))
-            ByBzEndpoints.append((root2[0], yConst2, root2[1]))
+
+    for root in root_list:
+        ByBzEndpoints.append((root[0], yConst2, root[1]))
 
     # Right Surface
     xConst1 = vspace[0][0][f111[0]][f111[1]][f111[2]]
     # Bx=By=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ax + bx * xConst1,
         cx + ex * xConst1,
         dx + fx * xConst1,
@@ -750,15 +726,11 @@ def is_close(a, b):
         gy + hy * xConst1,
     )
 
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxByEndpoints.append((xConst1, root1[0], root1[1]))
-        else:
-            BxByEndpoints.append((xConst1, root1[0], root1[1]))
-            BxByEndpoints.append((xConst1, root2[0], root2[1]))
+    for root in root_list:
+        BxByEndpoints.append((xConst1, root[0], root[1]))
 
     # Bx=BZ=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ax + bx * xConst1,
         cx + ex * xConst1,
         dx + fx * xConst1,
@@ -768,15 +740,12 @@ def is_close(a, b):
         dz + fz * xConst1,
         gz + hz * xConst1,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxBzEndpoints.append((xConst1, root1[0], root1[1]))
-        else:
-            BxBzEndpoints.append((xConst1, root1[0], root1[1]))
-            BxBzEndpoints.append((xConst1, root2[0], root2[1]))
+
+    for root in root_list:
+        BxBzEndpoints.append((xConst1, root[0], root[1]))
 
     # By=Bz=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ay + by * xConst1,
         cy + ey * xConst1,
         dy + fy * xConst1,
@@ -786,17 +755,14 @@ def is_close(a, b):
         dz + fz * xConst1,
         gz + hz * xConst1,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            ByBzEndpoints.append((xConst1, root1[0], root1[1]))
-        else:
-            ByBzEndpoints.append((xConst1, root1[0], root1[1]))
-            ByBzEndpoints.append((xConst1, root2[0], root2[1]))
+
+    for root in root_list:
+        ByBzEndpoints.append((xConst1, root[0], root[1]))
 
     # Left Surface
     xConst2 = vspace[0][0][f000[0]][f000[1]][f000[2]]
     # Bx=By=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ax + bx * xConst2,
         cx + ex * xConst2,
         dx + fx * xConst2,
@@ -806,14 +772,12 @@ def is_close(a, b):
         dy + fy * xConst2,
         gy + hy * xConst2,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxByEndpoints.append((xConst2, root1[0], root1[1]))
-        else:
-            BxByEndpoints.append((xConst2, root1[0], root1[1]))
-            BxByEndpoints.append((xConst2, root2[0], root2[1]))
+
+    for root in root_list:
+        BxByEndpoints.append((xConst2, root[0], root[1]))
+
     # Bx=BZ=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ax + bx * xConst2,
         cx + ex * xConst2,
         dx + fx * xConst2,
@@ -823,15 +787,12 @@ def is_close(a, b):
         dz + fz * xConst2,
         gz + hz * xConst2,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxBzEndpoints.append((xConst2, root1[0], root1[1]))
-        else:
-            BxBzEndpoints.append((xConst2, root1[0], root1[1]))
-            BxBzEndpoints.append((xConst2, root2[0], root2[1]))
+
+    for root in root_list:
+        BxBzEndpoints.append((xConst2, root[0], root[1]))
 
     # By=Bz=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ay + by * xConst2,
         cy + ey * xConst2,
         dy + fy * xConst2,
@@ -841,17 +802,14 @@ def is_close(a, b):
         dz + fz * xConst2,
         gz + hz * xConst2,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            ByBzEndpoints.append((xConst2, root1[0], root1[1]))
-        else:
-            ByBzEndpoints.append((xConst2, root1[0], root1[1]))
-            ByBzEndpoints.append((xConst2, root2[0], root2[1]))
+
+    for root in root_list:
+        ByBzEndpoints.append((xConst2, root[0], root[1]))
 
     # Up Surface
     zConst1 = vspace[0][2][f111[0]][f111[1]][f111[2]]
     # Bx=By=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ax + dx * zConst1,
         bx + fx * zConst1,
         cx + gx * zConst1,
@@ -861,14 +819,12 @@ def is_close(a, b):
         cy + gy * zConst1,
         ey + hy * zConst1,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxByEndpoints.append((root1[0], root1[1], zConst1))
-        else:
-            BxByEndpoints.append((root1[0], root1[1], zConst1))
-            BxByEndpoints.append((root2[0], root2[1], zConst1))
-    # Bx=BZ=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+
+    for root in root_list:
+        BxByEndpoints.append((root[0], root[1], zConst1))
+
+    # Bx=Bz=0 Curve Endpoint
+    root_list = _bilinear_root(
         ax + dx * zConst1,
         bx + fx * zConst1,
         cx + gx * zConst1,
@@ -878,14 +834,12 @@ def is_close(a, b):
         cz + gz * zConst1,
         ez + hz * zConst1,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxBzEndpoints.append((root1[0], root1[1], zConst1))
-        else:
-            BxBzEndpoints.append((root1[0], root1[1], zConst1))
-            BxBzEndpoints.append((root2[0], root2[1], zConst1))
+
+    for root in root_list:
+        BxBzEndpoints.append((root[0], root[1], zConst1))
+
     # By=Bz=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ay + dy * zConst1,
         by + fy * zConst1,
         cy + gy * zConst1,
@@ -895,17 +849,14 @@ def is_close(a, b):
         cz + gz * zConst1,
         ez + hz * zConst1,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            ByBzEndpoints.append((root1[0], root1[1], zConst1))
-        else:
-            ByBzEndpoints.append((root1[0], root1[1], zConst1))
-            ByBzEndpoints.append((root2[0], root2[1], zConst1))
+
+    for root in root_list:
+        ByBzEndpoints.append((root[0], root[1], zConst1))
 
     # Down Surface
     zConst2 = vspace[0][2][f000[0]][f000[1]][f000[2]]
     # Bx=By=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ax + dx * zConst2,
         bx + fx * zConst2,
         cx + gx * zConst2,
@@ -915,14 +866,12 @@ def is_close(a, b):
         cy + gy * zConst2,
         ey + hy * zConst2,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxByEndpoints.append((root1[0], root1[1], zConst2))
-        else:
-            BxByEndpoints.append((root1[0], root1[1], zConst2))
-            BxByEndpoints.append((root2[0], root2[1], zConst2))
+
+    for root in root_list:
+        BxByEndpoints.append((root[0], root[1], zConst2))
+
     # Bx=Bz=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ax + dx * zConst2,
         bx + fx * zConst2,
         cx + gx * zConst2,
@@ -932,14 +881,12 @@ def is_close(a, b):
         cz + gz * zConst2,
         ez + hz * zConst2,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            BxBzEndpoints.append((root1[0], root1[1], zConst2))
-        else:
-            BxBzEndpoints.append((root1[0], root1[1], zConst2))
-            BxBzEndpoints.append((root2[0], root2[1], zConst2))
+
+    for root in root_list:
+        BxBzEndpoints.append((root[0], root[1], zConst2))
+
     # By=Bz=0 Curve Endpoint
-    root1, root2 = _bilinear_root(
+    root_list = _bilinear_root(
         ay + dy * zConst2,
         by + fy * zConst2,
         cy + gy * zConst2,
@@ -949,28 +896,31 @@ def is_close(a, b):
         cz + gz * zConst2,
         ez + hz * zConst2,
     )
-    if root1 is not None:
-        if is_close(root1, root2):
-            ByBzEndpoints.append((root1[0], root1[1], zConst2))
-        else:
-            ByBzEndpoints.append((root1[0], root1[1], zConst2))
-            ByBzEndpoints.append((root2[0], root2[1], zConst2))
+
+    for root in root_list:
+        ByBzEndpoints.append((root[0], root[1], zConst2))
 
     xbound = vspace[0][0][f111[0]][f111[1]][f111[2]]
     ybound = vspace[0][1][f111[0]][f111[1]][f111[2]]
     zbound = vspace[0][2][f111[0]][f111[1]][f111[2]]
 
     def bound(epoint):
-        a = (initial[0] < epoint[0] or is_close(initial[0], epoint[0])) and (
-            epoint[0] < xbound or is_close(epoint[0], xbound)
-        )
-        b = (initial[1] < epoint[1] or is_close(initial[1], epoint[1])) and (
-            epoint[1] < ybound or is_close(epoint[1], ybound)
-        )
-        c = (initial[2] < epoint[2] or is_close(initial[2], epoint[2])) and (
-            epoint[2] < zbound or is_close(epoint[2], zbound)
-        )
-        return a and b and c
+        """
+        Checks if the endpoints are located in or on the cube.
+        """
+        is_x_in_bound = (
+            initial[0] < epoint[0]
+            or np.isclose(initial[0], epoint[0], atol=_EQUALITY_ATOL)
+        ) and (epoint[0] < xbound or np.isclose(epoint[0], xbound, atol=_EQUALITY_ATOL))
+        is_y_in_bound = (
+            initial[1] < epoint[1]
+            or np.isclose(initial[1], epoint[1], atol=_EQUALITY_ATOL)
+        ) and (epoint[1] < ybound or np.isclose(epoint[1], ybound, atol=_EQUALITY_ATOL))
+        is_z_in_bound = (
+            initial[2] < epoint[2]
+            or np.isclose(initial[2], epoint[2], atol=_EQUALITY_ATOL)
+        ) and (epoint[2] < zbound or np.isclose(epoint[2], zbound, atol=_EQUALITY_ATOL))
+        return is_x_in_bound and is_y_in_bound and is_z_in_bound
 
     BxByEndpoints = list(filter(bound, BxByEndpoints))
     BxBzEndpoints = list(filter(bound, BxBzEndpoints))
@@ -978,16 +928,16 @@ def bound(epoint):
 
     tlApprox = trilinear_approx(vspace, cell)
 
-    # Check on the Surfaces
-    for p in BxByEndpoints:
-        if np.linalg.norm(tlApprox(p[0], p[1], p[2])) < _ATOL:
-            return True
-    for p in BxBzEndpoints:
-        if np.linalg.norm(tlApprox(p[0], p[1], p[2])) < _ATOL:
-            return True
-    for p in ByBzEndpoints:
-        if np.linalg.norm(tlApprox(p[0], p[1], p[2])) < _ATOL:
-            return True
+    # # Check on the Surfaces
+    # for p in BxByEndpoints:
+    #     if np.linalg.norm(tlApprox(p[0], p[1], p[2])) < _EQUALITY_ATOL:
+    #         return True
+    # for p in BxBzEndpoints:
+    #     if np.linalg.norm(tlApprox(p[0], p[1], p[2])) < _EQUALITY_ATOL:
+    #         return True
+    # for p in ByBzEndpoints:
+    #     if np.linalg.norm(tlApprox(p[0], p[1], p[2])) < _EQUALITY_ATOL:
+    #         return True
 
     # Check Grid Resolution
     if len(BxByEndpoints) == 0 and len(BxBzEndpoints) == 0 and len(ByBzEndpoints) == 0:
@@ -1120,18 +1070,27 @@ def _locate_null_point(vspace, cell, n, err):
         ]
     )
 
-    def inbound(pos):
+    def in_bound(pos):
+        """
+        Checks if the estimated position is located inside the cube.
+        """
         pos = pos.reshape(1, 3)[0]
-        A = (np.isclose(pos_000[0], pos[0], atol=_ATOL) or pos_000[0] < pos[0]) and (
-            np.isclose(pos[0], pos_111[0], atol=_ATOL) or pos[0] < pos_111[0]
+        is_x_in_bound = (
+            np.isclose(pos_000[0], pos[0], atol=_EQUALITY_ATOL) or pos_000[0] < pos[0]
+        ) and (
+            np.isclose(pos[0], pos_111[0], atol=_EQUALITY_ATOL) or pos[0] < pos_111[0]
         )
-        B = (np.isclose(pos_000[1], pos[1], atol=_ATOL) or pos_000[1] < pos[1]) and (
-            np.isclose(pos[1], pos_111[1], atol=_ATOL) or pos[1] < pos_111[1]
+        is_y_in_bound = (
+            np.isclose(pos_000[1], pos[1], atol=_EQUALITY_ATOL) or pos_000[1] < pos[1]
+        ) and (
+            np.isclose(pos[1], pos_111[1], atol=_EQUALITY_ATOL) or pos[1] < pos_111[1]
         )
-        C = (np.isclose(pos_000[2], pos[2], atol=_ATOL) or pos_000[2] < pos[2]) and (
-            np.isclose(pos[2], pos_111[2], atol=_ATOL) or pos[2] < pos_111[2]
+        is_z_in_bound = (
+            np.isclose(pos_000[2], pos[2], atol=_EQUALITY_ATOL) or pos_000[2] < pos[2]
+        ) and (
+            np.isclose(pos[2], pos_111[2], atol=_EQUALITY_ATOL) or pos[2] < pos_111[2]
         )
-        return A and B and C
+        return is_x_in_bound and is_y_in_bound and is_z_in_bound
 
     starting_pos = []
     # Adding the Corners
@@ -1162,13 +1121,20 @@ def inbound(pos):
             Bx0 = np.array([locx, locy, locz])
             Bx0 = Bx0.reshape(3, 1)
             prev_norm = np.linalg.norm(x0)
+            # Too many null points if the jacobian is zero
+            if np.isclose(
+                np.linalg.det(jcb(x0[0], x0[1], x0[2])), 0, atol=_EQUALITY_ATOL
+            ):
+                warnings.warn("Multiple null points")
+                return None
+            # Adjust position
             x0 = np.subtract(
                 x0, np.matmul(np.linalg.inv(jcb(x0[0], x0[1], x0[2])), Bx0)
             )
             norm = np.linalg.norm(x0)
-            if np.abs(norm - prev_norm) < err and inbound(x0):
+            if np.abs(norm - prev_norm) < err and in_bound(x0):
                 return x0
-        if inbound(x0):
+        if in_bound(x0):
             warnings.warn("Max Iterations Reached")
             return x0
 
diff --git a/requirements/environment.yml b/requirements/environment.yml
index eaa6072c..621e2eb8 100644
--- a/requirements/environment.yml
+++ b/requirements/environment.yml
@@ -5,11 +5,11 @@ channels:
 - conda-forge
 
 dependencies:
-- python=3.9
+- python=3.10
 - numpy
 - scipy
 - astropy >= 4.3.1
-- h5py
+- h5py >= 3.0.0
 - matplotlib >= 3.3.0
 - mpmath
 - pandas
diff --git a/requirements/extras.txt b/requirements/extras.txt
index c1c7ce01..4bb738ab 100644
--- a/requirements/extras.txt
+++ b/requirements/extras.txt
@@ -1,6 +1,6 @@
 # these are "extra" dependencies to run advanced package features
 # ought to mirror 'extras' under options.extras_require in setup.cfg
-h5py
+h5py >= 3.0.0
 lmfit
 mpmath
 numba
diff --git a/setup.cfg b/setup.cfg
index 21928f62..879dc59f 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -24,7 +24,6 @@ classifiers=
   License :: OSI Approved :: BSD License
   Operating System :: OS Independent
   Programming Language :: Python :: 3
-  Programming Language :: Python :: 3.7
   Programming Language :: Python :: 3.8
   Programming Language :: Python :: 3.8
   Programming Language :: Python :: 3.9
@@ -34,7 +33,7 @@ classifiers=
   Topic :: Scientific/Engineering :: Astronomy
 
 [options]
-python_requires = >=3.7
+python_requires = >=3.8
 packages = find:
 include_package_data = True
 setup_requires =
@@ -65,7 +64,7 @@ extras =
   setuptools_scm
   tox
   # for plasmapy.plasma
-  h5py
+  h5py >= 3.0.0
   # for plasmapy.formulary
   mpmath
   lmfit
diff --git a/tox.ini b/tox.ini
index 526b7bed..6a47df85 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,5 +1,5 @@
 [tox]
-envlist = clean,py37,build_docs
+envlist = clean,py38,build_docs
 isolated_build = True
 indexserver =
     NIGHTLY = https://pypi.anaconda.org/scipy-wheels-nightly/simple
@@ -63,8 +63,8 @@ setenv =
 commands = sphinx-build docs docs{/}_build{/}html -W -n --keep-going -b html {posargs}
 
 # This env requires tox-conda.
-[testenv:py37-conda]
-basepython = python3.7
+[testenv:py38-conda]
+basepython = python3.8
 extras =
 deps =
   lmfit
@@ -75,7 +75,7 @@ conda_deps =
   scipy>=1.2
   astropy>=4.3.1
   pytest>=5.1
-  h5py
+  h5py>=3.0.0
   matplotlib
   mpmath
   numpydoc
@@ -84,8 +84,8 @@ conda_deps =
   sphinx_rtd_theme
 
 # This env tests minimal versions of each dependency.
-[testenv:py37-all-minimal]
-basepython = python3.7
+[testenv:py38-all-minimal]
+basepython = python3.8
 extras =
 deps =
   lmfit==1.0.1
@@ -93,7 +93,7 @@ deps =
   numpy==1.18.1
   scipy==1.2
   astropy==4.3.1
-  h5py==2.8
+  h5py==3.0.0
   matplotlib==3.3.0
   pytest==5.1
   mpmath==1.0
@@ -119,8 +119,8 @@ commands =
   flake8 --bug-report
   flake8 {toxinidir}{/}plasmapy --count --show-source --statistics
 
-[testenv:py37-minimal-pypi-import]
-basepython = python3.7
+[testenv:py38-minimal-pypi-import]
+basepython = python3.8
 extras =
 deps =
 commands = python -c 'import plasmapy'


Test patch:
diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index f972224b..1205189e 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -52,15 +52,10 @@ jobs:
       matrix:
         include:
 
-        - name: Python 3.7 with minimal dependencies
+        - name: Python 3.8 with minimal dependencies
           os: ubuntu-latest
-          python: 3.7
-          toxenv: py37-all-minimal
-
-        - name: Python 3.7
-          os: ubuntu-latest
-          python: 3.7
-          toxenv: py37-all
+          python: 3.8
+          toxenv: py38-all-minimal
 
         - name: Python 3.8, all tests, code coverage
           os: ubuntu-latest
@@ -137,10 +132,10 @@ jobs:
       run: |
         git fetch --prune --unshallow
         git fetch --depth=1 origin +refs/tags/*:refs/tags/*
-    - name: Set up Python 3.7
+    - name: Set up Python 3.8
       uses: actions/setup-python@v3
       with:
-        python-version: 3.7
+        python-version: 3.8
     - name: Install requirements
       run: |
         pip install --progress-bar off --upgrade pip
diff --git a/plasmapy/analysis/tests/test_nullpoint.py b/plasmapy/analysis/tests/test_nullpoint.py
index 0c52e05e..59dd5d39 100644
--- a/plasmapy/analysis/tests/test_nullpoint.py
+++ b/plasmapy/analysis/tests/test_nullpoint.py
@@ -5,8 +5,8 @@
 import pytest
 
 from plasmapy.analysis.nullpoint import (
-    _ATOL,
     _bilinear_root,
+    _EQUALITY_ATOL,
     _locate_null_point,
     _reduction,
     _trilinear_analysis,
@@ -17,6 +17,10 @@
     trilinear_approx,
 )
 
+# Defining tolerance level for tests where the accuracy
+# level can not rise too high due to run time issues.
+_TESTING_ATOL = 1e-2
+
 
 def vspace_func_1(x, y, z):
     return [(y - 5.5), (z - 5.5), (x - 5.5)]
@@ -26,6 +30,26 @@ def vspace_func_2(x, y, z):
     return [2 * y - z - 5.5, 3 * x + z - 22, x ** 2 - 11 * x + y + 24.75]
 
 
+def vspace_func_3(x, y, z):
+    return [(y - 5.3) * (y - 5.5), (z - 5.5), (x - 5.5)]
+
+
+def vspace_func_4(x, y, z):
+    return [y - y, (z - 5.5), (x - 5.5)]
+
+
+def vspace_func_5(x, y, z):
+    return [(y - y), (z - z), (x - x)]
+
+
+def vspace_func_6(x, y, z):
+    return [(-1 - x - z), (-x - z - x * z), (y)]
+
+
+def vspace_func_7(x, y, z):
+    return [(y - 5.5) * (y + 5.5), (z - 5.5), (x - 5.5)]
+
+
 def test_trilinear_coeff_cal():
     r"""Test `~plasmapy.analysis.nullpoint.trilinear_coeff_cal`."""
     vspace1_args = {
@@ -75,19 +99,11 @@ def test_trilinear_jacobian():
     jcb = _trilinear_jacobian(vspace, [0, 0, 0])
     mtrx = jcb(0.5, 0.5, 0.5)
     exact_mtrx = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]])
-    assert np.allclose(mtrx, exact_mtrx, atol=_ATOL)
+    assert np.allclose(mtrx, exact_mtrx, atol=_EQUALITY_ATOL)
 
 
 def test_trilinear_approx():
     r"""Test `~plasmapy.analysis.nullpoint.trilinear_approx`."""
-    vspace1_args = {
-        "x_range": [0, 10],
-        "y_range": [0, 10],
-        "z_range": [0, 10],
-        "precision": [10 / 46, 10 / 46, 10 / 46],
-        "func": vspace_func_1,
-    }
-    vspace1 = _vector_space(**vspace1_args)
     vspace2_args = {
         "x_range": [0, 10],
         "y_range": [0, 10],
@@ -116,11 +132,13 @@ def test_trilinear_approx():
         approx = tlApprox(p[0], p[1], p[2])
         exact = vspace_func_2(p[0], p[1], p[2])
         approx = approx.reshape(1, 3)
-        assert np.allclose(approx, exact, atol=_ATOL)
+        assert np.allclose(approx, exact, atol=_EQUALITY_ATOL)
     # Testing Trilinear Approx function on a midpoint
     approx = tlApprox(mid[0], mid[1], mid[2])
     approx = approx.reshape(1, 3)
-    assert np.allclose(approx, [-5.39130435, -21.5652174, 23.68667299], atol=_ATOL)
+    assert np.allclose(
+        approx, [-5.39130435, -21.5652174, 23.68667299], atol=_EQUALITY_ATOL
+    )
 
 
 class Test_reduction:
@@ -188,10 +206,10 @@ def test_bilinear_root_vals(self, kwargs, expected):
         r"""Test expected values."""
         x1, y1 = _bilinear_root(**kwargs)[0]
         x2, y2 = _bilinear_root(**kwargs)[1]
-        assert np.isclose(x1, expected[0], atol=_ATOL)
-        assert np.isclose(y1, expected[1], atol=_ATOL)
-        assert np.isclose(x2, expected[2], atol=_ATOL)
-        assert np.isclose(y2, expected[3], atol=_ATOL)
+        assert np.isclose(x1, expected[0], atol=_EQUALITY_ATOL)
+        assert np.isclose(y1, expected[1], atol=_EQUALITY_ATOL)
+        assert np.isclose(x2, expected[2], atol=_EQUALITY_ATOL)
+        assert np.isclose(y2, expected[3], atol=_EQUALITY_ATOL)
 
 
 class Test_locate_null_point:
@@ -207,7 +225,7 @@ class Test_locate_null_point:
 
     test_locate_null_point_values = [
         (
-            {"vspace": vspace, "cell": [0, 0, 0], "n": 500, "err": _ATOL},
+            {"vspace": vspace, "cell": [0, 0, 0], "n": 500, "err": _EQUALITY_ATOL},
             np.array([5.5, 5.5, 5.5]),
         )
     ]
@@ -216,7 +234,7 @@ class Test_locate_null_point:
     def test_locate_null_point_vals(self, kwargs, expected):
         r"""Test expected values."""
         assert np.isclose(
-            _locate_null_point(**kwargs).reshape(1, 3), expected, atol=_ATOL
+            _locate_null_point(**kwargs).reshape(1, 3), expected, atol=_EQUALITY_ATOL
         ).all()
 
 
@@ -227,25 +245,122 @@ def test_null_point_find1():
         "x_range": [5, 6],
         "y_range": [5, 6],
         "z_range": [5, 6],
-        "precision": [1, 1, 1],
+        "precision": [0.1, 0.1, 0.1],
         "func": vspace_func_1,
     }
     npoints = null_point_find(**nullpoint_args)
     loc = npoints[0].loc.reshape(1, 3)
     assert len(npoints) == 1
-    assert np.isclose(loc, [5.5, 5.5, 5.5], atol=_ATOL).all()
+    assert np.isclose(loc, [5.5, 5.5, 5.5], atol=_EQUALITY_ATOL).all()
 
 
 def test_null_point_find2():
     r"""Test `~plasmapy.analysis.nullpoint.null_point_find`."""
     # Non-uniform grid
     nullpoint2_args = {
-        "x_arr": [0, 1, 2, 3, 4, 5, 6],
-        "y_arr": [0, 2, 4, 6, 8],
-        "z_arr": [0, 2, 4, 6, 8],
+        "x_arr": np.logspace(np.log10(5.48), np.log10(5.52), num=30),
+        "y_arr": np.logspace(np.log10(5.48), np.log10(5.52), num=30),
+        "z_arr": np.logspace(np.log10(5.48), np.log10(5.52), num=30),
         "func": vspace_func_1,
     }
     npoints2 = null_point_find(**nullpoint2_args)
     loc2 = npoints2[0].loc.reshape(1, 3)
     assert len(npoints2) == 1
-    assert np.isclose(loc2, [5.5, 5.5, 5.5], atol=_ATOL).all()
+    assert np.isclose(loc2, [5.5, 5.5, 5.5], atol=_EQUALITY_ATOL).all()
+
+
+def test_null_point_find3():
+    r"""Test `~plasmapy.analysis.nullpoint.null_point_find`."""
+    # Vector values passed by hand
+    nullpoint3_args = {
+        "x_range": [5, 6],
+        "y_range": [5, 6],
+        "z_range": [5, 6],
+        "precision": [1, 1, 1],
+        "u_arr": np.array([[[-0.5, -0.5], [0.5, 0.5]], [[-0.5, -0.5], [0.5, 0.5]]]),
+        "v_arr": np.array([[[-0.5, 0.5], [-0.5, 0.5]], [[-0.5, 0.5], [-0.5, 0.5]]]),
+        "w_arr": np.array([[[-0.5, -0.5], [-0.5, -0.5]], [[0.5, 0.5], [0.5, 0.5]]]),
+    }
+    npoints3 = null_point_find(**nullpoint3_args)
+    loc3 = npoints3[0].loc.reshape(1, 3)
+    assert len(npoints3) == 1
+    assert np.isclose(loc3, [5.5, 5.5, 5.5], atol=_EQUALITY_ATOL).all()
+    assert npoints3[0].get_type() == "N/A"
+
+
+def test_null_point_find4():
+    r"""Test `~plasmapy.analysis.nullpoint.null_point_find`."""
+    # Two null points
+    nullpoint4_args = {
+        "x_range": [5, 6],
+        "y_range": [5, 6],
+        "z_range": [5, 6],
+        "precision": [0.07, 0.003, 0.07],
+        "func": vspace_func_3,
+    }
+    npoints4 = null_point_find(**nullpoint4_args)
+    first_loc4 = npoints4[0].loc.reshape(1, 3)
+    second_loc4 = npoints4[1].loc.reshape(1, 3)
+    assert len(npoints4) == 2
+    assert np.isclose(first_loc4, [5.5, 5.3, 5.5], atol=_EQUALITY_ATOL).all()
+    assert np.isclose(second_loc4, [5.5, 5.5, 5.5], atol=_EQUALITY_ATOL).all()
+
+
+def test_null_point_find5():
+    r"""Test `~plasmapy.analysis.nullpoint.null_point_find`."""
+    # Many null points because a y vector dimension is zero
+    nullpoint5_args = {
+        "x_range": [5.4, 5.6],
+        "y_range": [5.4, 5.6],
+        "z_range": [5.4, 5.6],
+        "precision": [0.01, 0.01, 0.01],
+        "func": vspace_func_4,
+    }
+    npoints5 = null_point_find(**nullpoint5_args)
+    assert len(npoints5) == 0
+
+
+def test_null_point_find6():
+    r"""Test `~plasmapy.analysis.nullpoint.null_point_find`."""
+    # Many null points; All vector dimensions zero
+    nullpoint6_args = {
+        "x_range": [5, 6],
+        "y_range": [5, 6],
+        "z_range": [5, 6],
+        "precision": [0.08, 0.08, 0.08],
+        "func": vspace_func_5,
+    }
+    npoints6 = null_point_find(**nullpoint6_args)
+    assert len(npoints6) == 0
+
+
+def test_null_point_find7():
+    r"""Test `~plasmapy.analysis.nullpoint.null_point_find`."""
+    # No null points, discriminant less than zero
+    nullpoint7_args = {
+        "x_range": [-10, 10],
+        "y_range": [-10, 10],
+        "z_range": [-10, 10],
+        "precision": [1, 1, 1],
+        "func": vspace_func_6,
+    }
+    npoints7 = null_point_find(**nullpoint7_args)
+    assert len(npoints7) == 0
+
+
+def test_null_point_find8():
+    r"""Test `~plasmapy.analysis.nullpoint.null_point_find`."""
+    # Non-linear field
+    nullpoint8_args = {
+        "x_range": [-6, 6],
+        "y_range": [-6, 6],
+        "z_range": [-6, 6],
+        "precision": [0.3, 0.3, 0.3],
+        "func": vspace_func_7,
+    }
+    npoints8 = null_point_find(**nullpoint8_args)
+    assert len(npoints8) == 2
+    loc1 = npoints8[0].loc.reshape(1, 3)
+    loc2 = npoints8[1].loc.reshape(1, 3)
+    assert np.allclose(loc1, [5.5, -5.5, 5.5], atol=_TESTING_ATOL)
+    assert np.allclose(loc2, [5.5, 5.5, 5.5], atol=_TESTING_ATOL)