diff --git a/eoxmagmod/eoxmagmod/__init__.py b/eoxmagmod/eoxmagmod/__init__.py
index 78dbb0e..f2e09a3 100644
--- a/eoxmagmod/eoxmagmod/__init__.py
+++ b/eoxmagmod/eoxmagmod/__init__.py
@@ -65,7 +65,7 @@
     from .dipole_coords import (
         get_dipole_rotation_matrix, convert_to_dipole, vrot_from_dipole,
     )
-    from .sheval_dipole import sheval_dipole
+    from .dipole import sheval_dipole
     from .magnetic_time import mjd2000_to_magnetic_universal_time
     from .magnetic_model.loader_shc import load_model_shc, load_model_shc_combined
     from .magnetic_model.loader_igrf import load_model_igrf
diff --git a/eoxmagmod/eoxmagmod/sheval_dipole.py b/eoxmagmod/eoxmagmod/dipole.py
similarity index 98%
rename from eoxmagmod/eoxmagmod/sheval_dipole.py
rename to eoxmagmod/eoxmagmod/dipole.py
index f692e2f..ad1fa6a 100644
--- a/eoxmagmod/eoxmagmod/sheval_dipole.py
+++ b/eoxmagmod/eoxmagmod/dipole.py
@@ -33,6 +33,8 @@
 )
 from .dipole_coords import convert_to_dipole, vrot_from_dipole
 
+__all__ = ["sheval_dipole", "rotate_vectors_from_dipole"]
+
 
 def sheval_dipole(arr_in, coef, lat_ngp, lon_ngp,
                   coord_type_in=GEODETIC_ABOVE_WGS84,
diff --git a/eoxmagmod/eoxmagmod/dipole_coords.py b/eoxmagmod/eoxmagmod/dipole_coords.py
index 31e7f7a..ad22b73 100644
--- a/eoxmagmod/eoxmagmod/dipole_coords.py
+++ b/eoxmagmod/eoxmagmod/dipole_coords.py
@@ -33,6 +33,12 @@
     GEOCENTRIC_SPHERICAL, GEOCENTRIC_CARTESIAN,
 )
 
+__all__ = [
+    "get_dipole_rotation_matrix",
+    "convert_to_dipole",
+    "vrot_from_dipole",
+]
+
 DEG2RAD = pi / 180.0
 
 
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/coefficients.py b/eoxmagmod/eoxmagmod/magnetic_model/coefficients.py
index 3e5fae1..be96cda 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/coefficients.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/coefficients.py
@@ -37,6 +37,16 @@
     coeff_size, convert_value,
 )
 
+__all__ = [
+    "SHCoefficients",
+    "ComposedSHCoefficients",
+    "CombinedSHCoefficients",
+    "SparseSHCoefficients",
+    "SparseSHCoefficientsTimeDependent",
+    "SparseSHCoefficientsConstant",
+    "SparseSHCoefficientsTimeDependentDecimalYear",
+]
+
 
 class SHCoefficients:
     """ Abstract base class for all spherical harmonic coefficients. """
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/coefficients_mio.py b/eoxmagmod/eoxmagmod/magnetic_model/coefficients_mio.py
index 4bf448f..b32ed65 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/coefficients_mio.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/coefficients_mio.py
@@ -25,17 +25,20 @@
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 # THE SOFTWARE.
 #-------------------------------------------------------------------------------
-#pylint: disable=no-name-in-module
 
 from math import pi
 from collections import namedtuple
-from numpy import asarray, zeros, arange, broadcast_to, sin, cos
+from numpy import asarray, zeros
 from .coefficients import SparseSHCoefficients, coeff_size
 from ..time_util import (
     mjd2000_to_year_fraction as mjd2000_to_year_fraction_default,
 )
 from .._pymm import fourier2d
 
+__all__ = [
+    "SparseSHCoefficientsMIO",
+]
+
 SCALE_SEASONAL = 2*pi
 SCALE_DIURNAL = 2*pi/24.
 
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/field_lines.py b/eoxmagmod/eoxmagmod/magnetic_model/field_lines.py
index 73ca11b..9aac5e2 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/field_lines.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/field_lines.py
@@ -30,6 +30,8 @@
 from .._pymm import GEOCENTRIC_CARTESIAN, GEOCENTRIC_SPHERICAL, convert
 from ..util import vnorm, vrotate
 
+__all__ = ["trace_field_line"]
+
 EARTH_RADIUS = 6371.2
 DEFAULT_MIN_RADIUS = 0.9977 * EARTH_RADIUS
 DEFAULT_MAX_RADIUS = 1e6 * EARTH_RADIUS
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/loader_emm.py b/eoxmagmod/eoxmagmod/magnetic_model/loader_emm.py
index 990d6da..b824cd6 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/loader_emm.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/loader_emm.py
@@ -34,6 +34,8 @@
 )
 from .parser_emm import combine_emm_coefficients, parse_emm_file
 
+__all__ = ["load_model_emm", "load_coeff_emm"]
+
 
 def load_model_emm(path_static, path_secvar):
     """ Load model from a EMM coefficient files. """
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/loader_igrf.py b/eoxmagmod/eoxmagmod/magnetic_model/loader_igrf.py
index c9a6d57..c705045 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/loader_igrf.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/loader_igrf.py
@@ -30,6 +30,8 @@
 from .coefficients import SparseSHCoefficientsTimeDependentDecimalYear
 from .parser_igrf import parse_igrf_file
 
+__all__ = ["load_model_igrf", "load_coeff_igrf"]
+
 
 def load_model_igrf(path):
     """ Load model from an IGRF coefficient file.
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/loader_mio.py b/eoxmagmod/eoxmagmod/magnetic_model/loader_mio.py
index d799226..0559fba 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/loader_mio.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/loader_mio.py
@@ -34,6 +34,13 @@
 from .coefficients_mio import SparseSHCoefficientsMIO
 from .parser_mio import parse_swarm_mio_file
 
+__all__ = [
+    "load_model_swarm_mio_internal",
+    "load_model_swarm_mio_external",
+    "load_coeff_swarm_mio_internal",
+    "load_coeff_swarm_mio_external",
+    "convert_external_mio_coeff",
+]
 
 def load_model_swarm_mio_internal(path):
     """ Load internal (secondary field) model from a Swarm MIO_SHA_2* product.
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/loader_mma.py b/eoxmagmod/eoxmagmod/magnetic_model/loader_mma.py
index ea7ef40..6ca6a73 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/loader_mma.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/loader_mma.py
@@ -42,6 +42,21 @@
     read_swarm_mma_2f_sm_internal, read_swarm_mma_2f_sm_external,
 )
 
+__all__ = [
+    "load_model_swarm_mma_2c_internal",
+    "load_model_swarm_mma_2c_external",
+    "load_model_swarm_mma_2f_geo_internal",
+    "load_model_swarm_mma_2f_geo_external",
+    "load_model_swarm_mma_2f_sm_internal",
+    "load_model_swarm_mma_2f_sm_external",
+    "load_coeff_swarm_mma_2c_internal",
+    "load_coeff_swarm_mma_2c_external",
+    "load_coeff_swarm_mma_2f_geo_internal",
+    "load_coeff_swarm_mma_2f_geo_external",
+    "load_coeff_swarm_mma_2f_sm_internal",
+    "load_coeff_swarm_mma_2f_sm_external",
+]
+
 # North geomagnetic coordinates used by the MMA products (IGRF-11, 2010.0)
 MMA2C_NGP_LATITUDE = 90 - 9.92   # deg.
 MMA2C_NGP_LONGITUDE = 287.78 - 360.0  # deg.
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/loader_shc.py b/eoxmagmod/eoxmagmod/magnetic_model/loader_shc.py
index 135b432..53a7f1c 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/loader_shc.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/loader_shc.py
@@ -39,16 +39,28 @@
 )
 from .parser_shc import parse_shc_file
 
+__all__ = [
+    "load_model_shc_combined",
+    "load_model_shc",
+    "load_coeff_shc_combined",
+    "load_coeff_shc_composed",
+    "load_coeff_shc",
+]
+
 
 def load_model_shc_combined(*paths, **kwargs):
-    """ Load model with coefficients combined from multiple SHC files. """
+    """ Load model with coefficients combined from multiple SHC files.
+    E.g., combining core and lithospheric models.
+    """
     return SphericalHarmomicGeomagneticModel(
         load_coeff_shc_combined(*paths, **kwargs)
     )
 
 
 def load_model_shc(*paths, **kwargs):
-    """ Load model from an SHC file. """
+    """ Load composed model from one or more SHC files.
+    E.g., composing multiple core models, each with a different time extent.
+    """
     return SphericalHarmomicGeomagneticModel(
         load_coeff_shc_composed(*paths, **kwargs)
     )
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/loader_wmm.py b/eoxmagmod/eoxmagmod/magnetic_model/loader_wmm.py
index 149a4e1..0339b95 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/loader_wmm.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/loader_wmm.py
@@ -30,6 +30,8 @@
 from .coefficients import SparseSHCoefficientsTimeDependentDecimalYear
 from .parser_wmm import parse_wmm_file
 
+__all__ = ["load_model_wmm", "load_coeff_wmm"]
+
 
 def load_model_wmm(path):
     """ Load model from a WMM COF file. """
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/model.py b/eoxmagmod/eoxmagmod/magnetic_model/model.py
index c1dfeda..ac40063 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/model.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/model.py
@@ -30,7 +30,7 @@
 
 from numpy import asarray, empty, full, nan, nditer
 from .._pymm import GRADIENT, GEOCENTRIC_SPHERICAL, sheval, shevaltemp
-from ..sheval_dipole import rotate_vectors_from_dipole
+from ..dipole import rotate_vectors_from_dipole
 from ..dipole_coords import convert_to_dipole
 from .util import reshape_times_and_coordinates
 from .coefficients import (
@@ -38,6 +38,12 @@
     CombinedSHCoefficients,
 )
 
+__all__ = [
+    "GeomagneticModel",
+    "SphericalHarmomicGeomagneticModel",
+    "DipoleSphericalHarmomicGeomagneticModel",
+]
+
 
 class GeomagneticModel:
     """ Abstract base class of the Earth magnetic field model. """
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/model_composed.py b/eoxmagmod/eoxmagmod/magnetic_model/model_composed.py
index 1dd3536..479e6a8 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/model_composed.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/model_composed.py
@@ -34,6 +34,7 @@
 )
 from .model import GeomagneticModel
 
+__all__ = ["ComposedGeomagneticModel"]
 
 Component = namedtuple("_Component", ["model", "scale", "parameters"])
 
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/model_mio.py b/eoxmagmod/eoxmagmod/magnetic_model/model_mio.py
index c234c6c..9c6df1f 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/model_mio.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/model_mio.py
@@ -25,7 +25,7 @@
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 # THE SOFTWARE.
 #-------------------------------------------------------------------------------
-# pylint: disable=too-many-arguments,too-many-locals
+# pylint: disable=too-many-arguments
 
 from numpy import nan, asarray, empty, isnan, full
 from .._pymm import GRADIENT, GEOCENTRIC_SPHERICAL, convert, sheval2dfs
@@ -34,6 +34,8 @@
 from .model import GeomagneticModel, DipoleSphericalHarmomicGeomagneticModel
 from .util import reshape_times_and_coordinates, reshape_array, mask_array
 
+__all__ = ["DipoleMIOGeomagneticModel", "MIOPrimaryGeomagneticModel"]
+
 MIO_HEIGHT = 110.0 # km
 MIO_EARTH_RADIUS = 6371.2 # km
 MIO_WOLF_RATIO = 0.014850
@@ -135,6 +137,7 @@ class DipoleMIOGeomagneticModel(DipoleSphericalHarmomicGeomagneticModel):
 
     def __init__(self, coefficients, north_pole, wolf_ratio=MIO_WOLF_RATIO,
                  height=MIO_HEIGHT, earth_radius=MIO_EARTH_RADIUS):
+        del height, earth_radius
 
         if not isinstance(coefficients, SparseSHCoefficientsMIO):
             raise TypeError(
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/parser_emm.py b/eoxmagmod/eoxmagmod/magnetic_model/parser_emm.py
index afadae3..a8c26d3 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/parser_emm.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/parser_emm.py
@@ -29,6 +29,10 @@
 import re
 from numpy import array, stack
 
+__all__ = [
+    "combine_emm_coefficients",
+    "parse_emm_file",
+]
 
 RE_HEADER_LINE = re.compile(r'^-+$')
 EMM_VALIDITY_PERIOD = 5.0 # years
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/parser_igrf.py b/eoxmagmod/eoxmagmod/magnetic_model/parser_igrf.py
index 2eaa2aa..f71bbd4 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/parser_igrf.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/parser_igrf.py
@@ -29,6 +29,8 @@
 
 from numpy import array
 
+__all__ = ["parse_igrf_file"]
+
 IGRF_EXTRAPOLATION_PERIOD = 5.0 # years
 
 
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/parser_mio.py b/eoxmagmod/eoxmagmod/magnetic_model/parser_mio.py
index 577f77a..ce3b310 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/parser_mio.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/parser_mio.py
@@ -29,6 +29,8 @@
 from numpy import array
 from .parser_shc import _strip_shc_comments
 
+__all__ = ["parse_swarm_mio_file"]
+
 
 def parse_swarm_mio_file(file_in):
     """ Parse Swarm MIO_SHA_2* product file format and return a dictionary
@@ -41,6 +43,7 @@ def parse_swarm_mio_file(file_in):
     data["degree_max"] = data["nm"][:, 0].max()
     return data
 
+
 def parse_swarm_mio_coefficients(lines, data):
     """ Parse the Swarm MIO_SHA_2* coefficients. """
     nm_idx = []
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/parser_mma.py b/eoxmagmod/eoxmagmod/magnetic_model/parser_mma.py
index 6a945eb..9eec631 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/parser_mma.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/parser_mma.py
@@ -29,6 +29,15 @@
 from numpy import array
 from spacepy import pycdf
 
+__all__ = [
+    "read_swarm_mma_2c_internal",
+    "read_swarm_mma_2c_external",
+    "read_swarm_mma_2f_geo_internal",
+    "read_swarm_mma_2f_geo_external",
+    "read_swarm_mma_2f_sm_internal",
+    "read_swarm_mma_2f_sm_external",
+]
+
 CDF_EPOCH_TYPE = pycdf.const.CDF_EPOCH.value
 CDF_EPOCH_2000 = 63113904000000.0
 CDF_EPOCH_TO_DAYS = 1.0/86400000.0
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/parser_shc.py b/eoxmagmod/eoxmagmod/magnetic_model/parser_shc.py
index 849edf5..dc89fd8 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/parser_shc.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/parser_shc.py
@@ -28,6 +28,8 @@
 
 from numpy import inf, array
 
+__all__ = ["parse_shc_file", "parse_shc_header"]
+
 
 def parse_shc_file(file_in):
     """ Parse SHC file and return a dictionary containing the parsed model data.
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/parser_wmm.py b/eoxmagmod/eoxmagmod/magnetic_model/parser_wmm.py
index de17369..b800e15 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/parser_wmm.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/parser_wmm.py
@@ -29,6 +29,8 @@
 import re
 from numpy import array
 
+__all__ = ["parse_wmm_file"]
+
 RE_TERMINATING_LINE = re.compile(r'^9+$')
 WMM_VALIDITY_PERIOD = 5.0 # years
 
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/tests/model_loaders.py b/eoxmagmod/eoxmagmod/magnetic_model/tests/model_loaders.py
index b33cf26..6adb065 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/tests/model_loaders.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/tests/model_loaders.py
@@ -80,7 +80,7 @@
     GEOCENTRIC_CARTESIAN, GEODETIC_ABOVE_WGS84, GEOCENTRIC_SPHERICAL, convert,
     GRADIENT, sheval,
 )
-from eoxmagmod.sheval_dipole import sheval_dipole
+from eoxmagmod.dipole import sheval_dipole
 
 
 class SHModelTestMixIn:
diff --git a/eoxmagmod/eoxmagmod/magnetic_model/util.py b/eoxmagmod/eoxmagmod/magnetic_model/util.py
index 5068ac2..2c93a9b 100644
--- a/eoxmagmod/eoxmagmod/magnetic_model/util.py
+++ b/eoxmagmod/eoxmagmod/magnetic_model/util.py
@@ -30,6 +30,17 @@
 from numpy import inf
 from numpy.lib.stride_tricks import as_strided
 
+__all__ = [
+    "coeff_size",
+    "convert_value",
+    "parse_file",
+    "reshape_times_and_coordinates",
+    "mask_array",
+    "reshape_array",
+    "get_nonoverlapping_intervals",
+    "aggregate_intersected_intervals",
+]
+
 
 def coeff_size(degree):
     """ Calculate size of the full coefficient array from the given degree. """
diff --git a/eoxmagmod/eoxmagmod/magnetic_time.py b/eoxmagmod/eoxmagmod/magnetic_time.py
index 46302de..967297f 100644
--- a/eoxmagmod/eoxmagmod/magnetic_time.py
+++ b/eoxmagmod/eoxmagmod/magnetic_time.py
@@ -30,6 +30,8 @@
 from numpy import sin, cos, arctan2
 from .solar_position import sunpos
 
+__all__ = ["mjd2000_to_magnetic_universal_time"]
+
 DEG2RAD = pi/180.0
 RAD2HOUR = 12.0/pi
 
diff --git a/eoxmagmod/eoxmagmod/quasi_dipole_coordinates.py b/eoxmagmod/eoxmagmod/quasi_dipole_coordinates.py
index 6278895..2e13b99 100644
--- a/eoxmagmod/eoxmagmod/quasi_dipole_coordinates.py
+++ b/eoxmagmod/eoxmagmod/quasi_dipole_coordinates.py
@@ -30,6 +30,13 @@
 from . import _pyqd
 from .data import APEX
 
+__all__ = [
+    "eval_qdlatlon",
+    "eval_qdlatlon_with_base_vectors",
+    "eval_mlt",
+    "eval_subsol",
+]
+
 
 def eval_qdlatlon(gclat, gclon, gcrad, time, fname=APEX):
     """
diff --git a/eoxmagmod/eoxmagmod/solar_position.py b/eoxmagmod/eoxmagmod/solar_position.py
index f5c0f43..e48b3cc 100644
--- a/eoxmagmod/eoxmagmod/solar_position.py
+++ b/eoxmagmod/eoxmagmod/solar_position.py
@@ -32,6 +32,8 @@
 
 from . import _pysunpos
 
+__all__ = ["sunpos"]
+
 
 def sunpos(time_mjd2k, lat, lon, rad=6371.2, dtt=0):
     """ Calculate solar equatorial and horizontal coordinates
diff --git a/eoxmagmod/eoxmagmod/tests/sheval_dipole.py b/eoxmagmod/eoxmagmod/tests/dipole.py
similarity index 99%
rename from eoxmagmod/eoxmagmod/tests/sheval_dipole.py
rename to eoxmagmod/eoxmagmod/tests/dipole.py
index 717bf2c..43c8403 100644
--- a/eoxmagmod/eoxmagmod/tests/sheval_dipole.py
+++ b/eoxmagmod/eoxmagmod/tests/dipole.py
@@ -38,7 +38,7 @@
     convert, relradpow, loncossin, legendre, spharpot, sphargrd,
 )
 from eoxmagmod.dipole_coords import convert_to_dipole, vrot_from_dipole
-from eoxmagmod.sheval_dipole import sheval_dipole
+from eoxmagmod.dipole import sheval_dipole
 from eoxmagmod.tests.data import mma_external, mma_internal
 
 
diff --git a/eoxmagmod/eoxmagmod/tests/time_util.py b/eoxmagmod/eoxmagmod/tests/time_util.py
index ff89793..50b7e9e 100644
--- a/eoxmagmod/eoxmagmod/tests/time_util.py
+++ b/eoxmagmod/eoxmagmod/tests/time_util.py
@@ -28,16 +28,39 @@
 # pylint: disable=missing-docstring, invalid-name, too-few-public-methods
 
 from unittest import TestCase, main
+from datetime import date, datetime
 from numpy import vectorize, inf, nan
 from numpy.random import uniform
 from numpy.testing import assert_allclose
+from eoxmagmod.tests.util import FunctionTestMixIn
 from eoxmagmod.time_util import (
     decimal_year_to_mjd2000_simple,
     mjd2000_to_decimal_year_simple,
     mjd2000_to_year_fraction_simple,
+    datetime_to_decimal_year,
 )
 
 
+class TestDatetimeToDecimalYear(FunctionTestMixIn, TestCase):
+    NAME = "datetime_to_decimal_year"
+
+    @staticmethod
+    def eval(input_):
+        return datetime_to_decimal_year(input_)
+
+    ACCEPTED = [
+        (datetime(2001, 1, 12), 2001.0301369863014),
+        (datetime(2012, 8, 31), 2012.6639344262296),
+        (datetime(2014, 8, 31), 2014.66301369863),
+        (datetime(2024, 12, 31, 23, 59, 59, 999), 2024.9999999684085),
+    ]
+
+    REJECTED = [
+        (None, TypeError),
+        (date(2001, 1, 12), TypeError),
+    ]
+
+
 class TestMjd2000ToYearFractionSimple(TestCase):
 
     @staticmethod
diff --git a/eoxmagmod/eoxmagmod/tests/util.py b/eoxmagmod/eoxmagmod/tests/util.py
index 5f403a0..e8d37d5 100644
--- a/eoxmagmod/eoxmagmod/tests/util.py
+++ b/eoxmagmod/eoxmagmod/tests/util.py
@@ -28,10 +28,9 @@
 # pylint: disable=missing-docstring
 
 from unittest import main, TestCase
-from datetime import date, datetime
 from numpy import array
 from numpy.testing import assert_allclose
-from eoxmagmod.util import vnorm, vincdecnorm, datetime_to_decimal_year
+from eoxmagmod.util import vnorm, vincdecnorm
 
 
 class FunctionTestMixIn:
@@ -62,26 +61,6 @@ def test_rejected(self):
                 raise AssertionError(f"\nInput: {input_}\n{exc}") from None
 
 
-class TestDatetimeToDecimalYear(FunctionTestMixIn, TestCase):
-    NAME = "datetime_to_decimal_year"
-
-    @staticmethod
-    def eval(input_):
-        return datetime_to_decimal_year(input_)
-
-    ACCEPTED = [
-        (datetime(2001, 1, 12), 2001.0301369863014),
-        (datetime(2012, 8, 31), 2012.6639344262296),
-        (datetime(2014, 8, 31), 2014.66301369863),
-        (datetime(2024, 12, 31, 23, 59, 59, 999), 2024.9999999684085),
-    ]
-
-    REJECTED = [
-        (None, TypeError),
-        (date(2001, 1, 12), TypeError),
-    ]
-
-
 class TestVnorm(FunctionTestMixIn, TestCase):
     NAME = "vnorm"
 
diff --git a/eoxmagmod/eoxmagmod/time_util.py b/eoxmagmod/eoxmagmod/time_util.py
index dfb7f57..721d01b 100644
--- a/eoxmagmod/eoxmagmod/time_util.py
+++ b/eoxmagmod/eoxmagmod/time_util.py
@@ -28,6 +28,7 @@
 # THE SOFTWARE.
 #-------------------------------------------------------------------------------
 
+from datetime import datetime
 from numpy import asarray, floor
 from ._pytimeconv import (
     decimal_year_to_mjd2000,
@@ -35,6 +36,16 @@
     mjd2000_to_year_fraction,
 )
 
+__all__ = [
+    "decimal_year_to_mjd2000",
+    "mjd2000_to_decimal_year",
+    "mjd2000_to_year_fraction",
+    "mjd2000_to_decimal_year_simple",
+    "mjd2000_to_year_fraction_simple",
+    "decimal_year_to_mjd2000_simple",
+    "datetime_to_decimal_year",
+]
+
 
 def mjd2000_to_decimal_year_simple(mjd2000):
     """ Convert Modified Julian Date since 2000 to (Julian) decimal year
@@ -59,3 +70,19 @@ def decimal_year_to_mjd2000_simple(decimal_year):
         mjd2000 = (decimal_year - 2000.0) * 365.25
     """
     return (asarray(decimal_year) - 2000.0) * 365.25
+
+
+def datetime_to_decimal_year(time):
+    """ Convert time given by a `datetime.datetime` object to a decimal year
+    value.
+    """
+    if not isinstance(time, datetime):
+        raise TypeError("The input must be a datetime object.")
+
+    year_start = datetime(year=time.year, month=1, day=1)
+    next_year_start = datetime(year=time.year+1, month=1, day=1)
+
+    year_elapsed = (time - year_start).total_seconds()
+    year_total = (next_year_start - year_start).total_seconds()
+
+    return time.year + year_elapsed / year_total
diff --git a/eoxmagmod/eoxmagmod/util.py b/eoxmagmod/eoxmagmod/util.py
index c9b9a55..aee6e8f 100644
--- a/eoxmagmod/eoxmagmod/util.py
+++ b/eoxmagmod/eoxmagmod/util.py
@@ -27,7 +27,6 @@
 #-------------------------------------------------------------------------------
 # pylint: disable=no-name-in-module
 
-from datetime import datetime
 from numpy import sqrt, asarray
 from ._pymm import (
     GEODETIC_ABOVE_WGS84, GEOCENTRIC_SPHERICAL, GEOCENTRIC_CARTESIAN,
@@ -36,6 +35,13 @@
 
 SPHERICAL_COORD_TYPES = (GEODETIC_ABOVE_WGS84, GEOCENTRIC_SPHERICAL)
 
+__all__ = [
+    "vrotate",
+    "vnorm",
+    "vincdecnorm",
+]
+
+
 def vrotate(arr, coord_in, coord_out, coord_type_in, coord_type_out):
     """ Rotate vectors from one coordinate system to another.
         Input:
@@ -96,19 +102,3 @@ def vincdecnorm(arr):
     """
     tmp = convert(arr, GEOCENTRIC_CARTESIAN, GEOCENTRIC_SPHERICAL)
     return -tmp[..., 0], tmp[..., 1], tmp[..., 2]
-
-
-def datetime_to_decimal_year(time):
-    """ Convert time given by a `datetime.datetime` object to a decimal year
-    value.
-    """
-    if not isinstance(time, datetime):
-        raise TypeError("The input must be a datetime object.")
-
-    year_start = datetime(year=time.year, month=1, day=1)
-    next_year_start = datetime(year=time.year+1, month=1, day=1)
-
-    year_elapsed = (time - year_start).total_seconds()
-    year_total = (next_year_start - year_start).total_seconds()
-
-    return time.year + year_elapsed / year_total