diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 6a18bae..7142d0a 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -17,7 +17,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest]
-        use_pymap: [true, false]
+        python_version: ['3.9', '3.12']
 
     runs-on: ${{ matrix.os }}
 
@@ -36,13 +36,10 @@ jobs:
 
     - uses: actions/setup-python@v5
       with:
-        python-version: '3.8'
+        python-version: ${{ matrix.python_version }}
 
     - run: pip install .[tests,lint]
 
-    - run: pip install pymap3d
-      if: ${{ matrix.use_pymap }}
-
     - run: flake8
     - run: mypy
 
diff --git a/pyproject.toml b/pyproject.toml
index 76112ad..1ff9661 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -15,7 +15,7 @@ classifiers = [
   "Topic :: Scientific/Engineering :: Atmospheric Science"
 ]
 dynamic = ["readme", "version"]
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 dependencies = ["python-dateutil", "numpy", "astropy", "xarray", "netcdf4"]
 
 [project.optional-dependencies]
diff --git a/scripts/PrintSourceRaDec.py b/scripts/PrintSourceRaDec.py
index 96eec0f..14b7596 100755
--- a/scripts/PrintSourceRaDec.py
+++ b/scripts/PrintSourceRaDec.py
@@ -3,19 +3,18 @@
 Utility program to print source Coordinates in RA,DEC that astrometry.net found
 use this with .rdls (rdls.fits) FITS files after solving an image
 """
-from pathlib import Path
-from astropy.io import fits
+
 import argparse
 
+from astrometry_azel.io import get_sources
+
 
 p = argparse.ArgumentParser(description="show RA DEC in .rdls files after solving an image")
 p.add_argument("fn", help=".rdls file from astrometry.net")
 p = p.parse_args()
 
-fn = Path(p.fn).expanduser().resolve(strict=True)
-with fits.open(fn, "readonly") as f:
-    radec = f[1].data
+radec_sources = get_sources(p.fn)
 
-print(radec.shape[0], "sources found in", fn, "with ra,dec coordinates:")
+print(radec_sources.shape[0], "sources found in", p.fn, "with ra,dec coordinates:")
 
-print(radec)
+print(radec_sources)
diff --git a/src/astrometry_azel/__init__.py b/src/astrometry_azel/__init__.py
index 4c8c83d..6c5a582 100644
--- a/src/astrometry_azel/__init__.py
+++ b/src/astrometry_azel/__init__.py
@@ -1,5 +1,204 @@
-from .base import fits2azel, fits2radec, radec2azel, doSolve
+from pathlib import Path
+from datetime import datetime
+from datetime import timezone as tz
+import functools
+import shutil
+import subprocess
+
+import numpy as np
+import xarray
+
+import logging
+from dateutil.parser import parse
+
+from astropy.time import Time
+from astropy.coordinates import AltAz, Angle, EarthLocation, SkyCoord
+from astropy import units as u
+from astropy.io import fits
+import astropy.wcs as awcs
 
 __all__ = ["fits2azel", "fits2radec", "radec2azel", "doSolve"]
 
-__version__ = "1.3.0"
+__version__ = "1.4.0"
+
+
+def fits2radec(fitsfn: Path, solve: bool = False, args: str = ""):
+    """
+    get RA, Decl from FITS file
+    """
+    fitsfn = Path(fitsfn).expanduser()
+
+    if solve:
+        doSolve(fitsfn, args)
+
+    with fits.open(fitsfn, mode="readonly") as f:
+        yPix, xPix = f[0].shape[-2:]
+
+    x, y = np.meshgrid(range(xPix), range(yPix))
+    # pixel indices to find RA/dec of
+    xy = np.column_stack((x.ravel(order="C"), y.ravel(order="C")))
+
+    if not (wcsfn := fitsfn.with_suffix(".wcs")).is_file():
+        if not (wcsfn := fitsfn.with_name("wcs.fits")).is_file():
+            raise FileNotFoundError(f"could not find WCS file for {fitsfn}")
+    with fits.open(wcsfn, mode="readonly") as f:
+        # %% use astropy.wcs to register pixels to RA/DEC
+        # http://docs.astropy.org/en/stable/api/astropy.wcs.WCS.html#astropy.wcs.WCS
+        # NOTE: it's normal to get this warning:
+        # WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
+        if f[0].header["WCSAXES"] == 2:
+            # greyscale image
+            radec = awcs.wcs.WCS(f[0].header).all_pix2world(xy, 0)
+        elif f[0].header["WCSAXES"] == 3:
+            # color image
+            radec = awcs.wcs.WCS(f[0].header, naxis=[0, 1]).all_pix2world(xy, 0)
+        else:
+            raise ValueError(f"{fitsfn} has {f[0].header['NAXIS']} axes -- expected 2 or 3")
+
+    ra = radec[:, 0].reshape((yPix, xPix), order="C")
+    dec = radec[:, 1].reshape((yPix, xPix), order="C")
+    # %% collect output
+    radec = xarray.Dataset(
+        {"ra": (("y", "x"), ra), "dec": (("y", "x"), dec)},
+        {"x": range(xPix), "y": range(yPix)},
+        attrs={"filename": str(fitsfn)},
+    )
+
+    radec["ra"].attrs["units"] = "Right Ascension degrees east"
+    radec["dec"].attrs["units"] = "Declination degrees north"
+    radec["x"].attrs["units"] = "pixel index"
+    radec["y"].attrs["units"] = "pixel index"
+
+    return radec
+
+
+def fits2azel(
+    fitsfn: Path,
+    *,
+    latlon: tuple[float, float],
+    time: datetime,
+    solve: bool = False,
+    args: str = "",
+):
+    fitsfn = Path(fitsfn).expanduser()
+
+    radec = fits2radec(fitsfn, solve, args)
+
+    return radec2azel(radec, latlon, time)
+
+
+def radec2azel(scale: xarray.Dataset, latlon: tuple[float, float], time: datetime):
+    """
+    right ascension/declination to azimuth/elevation
+    """
+
+    if isinstance(time, datetime):
+        pass
+    elif isinstance(time, (float, int)):  # assume UT1_Unix
+        time = datetime.fromtimestamp(time, tz=tz.UTC)
+    elif isinstance(time, str):
+        time = parse(time)
+    else:
+        raise TypeError(f"expected datetime, float, int, or str -- got {type(time)}")
+
+    print("image time:", time)
+    # %% knowing camera location, time, and sky coordinates observed, convert to az/el for each pixel
+    # .values is to avoid silently freezing AstroPy
+
+    az, el = pymap3d_radec2azel(scale["ra"].values, scale["dec"].values, *latlon, time)
+    if (el < 0).any():
+        Nbelow = (el < 0).nonzero()
+        logging.error(
+            f"{Nbelow} points were below the horizon."
+            "Currently this program assumed observer ~ ground level."
+        )
+
+    # %% collect output
+    scale["azimuth"] = (("y", "x"), az)
+    scale["azimuth"].attrs["units"] = "degrees clockwise from north"
+
+    scale["elevation"] = (("y", "x"), el)
+    scale["elevation"].attrs["units"] = "degrees above horizon"
+
+    scale["observer_latitude"] = latlon[0]
+    scale["observer_latitude"].attrs["units"] = "degrees north WGS84"
+
+    scale["observer_longitude"] = latlon[1]
+    scale["observer_longitude"].attrs["units"] = "degrees east WGS84"
+
+    # datetime64 can be saved to netCDF4, but Python datetime.datetime cannot
+    scale["time"] = np.datetime64(time)
+
+    return scale
+
+
+def pymap3d_radec2azel(
+    ra_deg,
+    dec_deg,
+    lat_deg,
+    lon_deg,
+    time: datetime,
+) -> tuple:
+    """
+    sky coordinates (ra, dec) to viewing angle (az, el)
+
+    Parameters
+    ----------
+    ra_deg : float
+         ecliptic right ascension (degress)
+    dec_deg : float
+         ecliptic declination (degrees)
+    lat_deg : float
+              observer latitude [-90, 90]
+    lon_deg : float
+              observer longitude [-180, 180] (degrees)
+    time : datetime.datetime
+           time of observation
+
+    Returns
+    -------
+    az_deg : float
+             azimuth [degrees clockwize from North]
+    el_deg : float
+             elevation [degrees above horizon (neglecting aberration)]
+    """
+
+    obs = EarthLocation(lat=lat_deg * u.deg, lon=lon_deg * u.deg)
+    points = SkyCoord(Angle(ra_deg, unit=u.deg), Angle(dec_deg, unit=u.deg), equinox="J2000.0")
+    altaz = points.transform_to(AltAz(location=obs, obstime=Time(time)))
+
+    return altaz.az.degree, altaz.alt.degree
+
+
+@functools.cache
+def get_solve_exe() -> str:
+    if not (solve := shutil.which("solve-field")):
+        raise FileNotFoundError("Astrometry.net solve-file exectuable not found")
+
+    return solve
+
+
+def doSolve(fitsfn: Path, args: str = "") -> None:
+    """
+    run Astrometry.net solve-field from Python
+    """
+
+    fitsfn = Path(fitsfn).expanduser()
+    if not fitsfn.is_file():
+        raise FileNotFoundError(f"{fitsfn} not found")
+
+    # %% build command
+    cmd = [get_solve_exe(), "--overwrite", str(fitsfn), "--verbose"]
+    if args:
+        # if args is a string, split it. Don't append an empty space or solve-field CLI fail
+        cmd += args.split(" ")
+    print("\n", " ".join(cmd), "\n")
+    # %% execute
+    # bufsize=1: line-buffered
+    with subprocess.Popen(cmd, stdout=subprocess.PIPE, bufsize=1, text=True) as p:
+        if p.stdout:
+            for line in p.stdout:
+                print(line, end="")
+
+    if "Did not solve" in line:
+        raise RuntimeError(f"could not solve {fitsfn}")
diff --git a/src/astrometry_azel/base.py b/src/astrometry_azel/base.py
deleted file mode 100644
index db1299d..0000000
--- a/src/astrometry_azel/base.py
+++ /dev/null
@@ -1,165 +0,0 @@
-from __future__ import annotations
-from pathlib import Path
-import subprocess
-import shutil
-import logging
-from dateutil.parser import parse
-from datetime import datetime
-from numpy import meshgrid, column_stack, datetime64
-import xarray
-from astropy.io import fits
-from astropy.wcs import wcs
-
-try:
-    import pymap3d
-except ImportError:
-    pymap3d = None
-
-__all__ = ["fits2radec", "fits2azel", "doSolve"]
-
-
-def fits2radec(fitsfn: Path, solve: bool = False, args: str = ""):
-    """
-    get RA, Decl from FITS file
-    """
-    fitsfn = Path(fitsfn).expanduser()
-
-    if solve:
-        doSolve(fitsfn, args)
-
-    with fits.open(fitsfn, mode="readonly") as f:
-        yPix, xPix = f[0].shape[-2:]
-
-    x, y = meshgrid(range(xPix), range(yPix))
-    # pixel indices to find RA/dec of
-    xy = column_stack((x.ravel(order="C"), y.ravel(order="C")))
-
-    if not (wcsfn := fitsfn.with_suffix(".wcs")).is_file():
-        if not (wcsfn := fitsfn.with_name("wcs.fits")).is_file():
-            raise FileNotFoundError(f"could not find WCS file for {fitsfn}")
-    with fits.open(wcsfn, mode="readonly") as f:
-        # %% use astropy.wcs to register pixels to RA/DEC
-        # http://docs.astropy.org/en/stable/api/astropy.wcs.WCS.html#astropy.wcs.WCS
-        # NOTE: it's normal to get this warning:
-        # WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
-        if f[0].header["WCSAXES"] == 2:
-            # greyscale image
-            radec = wcs.WCS(f[0].header).all_pix2world(xy, 0)
-        elif f[0].header["WCSAXES"] == 3:
-            # color image
-            radec = wcs.WCS(f[0].header, naxis=[0, 1]).all_pix2world(xy, 0)
-        else:
-            raise ValueError(f"{fitsfn} has {f[0].header['NAXIS']} axes -- expected 2 or 3")
-
-    ra = radec[:, 0].reshape((yPix, xPix), order="C")
-    dec = radec[:, 1].reshape((yPix, xPix), order="C")
-    # %% collect output
-    radec = xarray.Dataset(
-        {"ra": (("y", "x"), ra), "dec": (("y", "x"), dec)},
-        {"x": range(xPix), "y": range(yPix)},
-        attrs={"filename": str(fitsfn)},
-    )
-
-    radec["ra"].attrs["units"] = "Right Ascension degrees east"
-    radec["dec"].attrs["units"] = "Declination degrees north"
-    radec["x"].attrs["units"] = "pixel index"
-    radec["y"].attrs["units"] = "pixel index"
-
-    return radec
-
-
-def radec2azel(scale, latlon: tuple[float, float], time: datetime | None):
-    if pymap3d is None:
-        raise ImportError("azimuth, elevation computations require: pip install pymap3d")
-
-    assert isinstance(scale, xarray.Dataset)
-
-    if time is None:
-        with fits.open(scale.filename, mode="readonly") as f:
-            try:
-                t = f[0].header["FRAME"]
-                # NOTE: this only works from Solis FITS files
-            except KeyError:
-                raise ValueError("Could not determine time of image")
-        time = parse(t)
-        logging.info("using FITS header for time")
-    elif isinstance(time, datetime):
-        pass
-    elif isinstance(time, (float, int)):  # assume UT1_Unix
-        time = datetime.utcfromtimestamp(time)
-    else:  # user override of frame time
-        time = parse(time)
-
-    print("image time:", time)
-    # %% knowing camera location, time, and sky coordinates observed, convert to az/el for each pixel
-    # .values is to avoid silently freezing AstroPy
-
-    az, el = pymap3d.radec2azel(scale["ra"].values, scale["dec"].values, *latlon, time)
-    if (el < 0).any():
-        Nbelow = (el < 0).nonzero()
-        logging.error(
-            f"{Nbelow} points were below the horizon."
-            "Currently this program assumed observer ~ ground level."
-        )
-
-    # %% collect output
-    scale["azimuth"] = (("y", "x"), az)
-    scale["azimuth"].attrs["units"] = "degrees clockwise from north"
-
-    scale["elevation"] = (("y", "x"), el)
-    scale["elevation"].attrs["units"] = "degrees above horizon"
-
-    scale["observer_latitude"] = latlon[0]
-    scale["observer_latitude"].attrs["units"] = "degrees north WGS84"
-
-    scale["observer_longitude"] = latlon[1]
-    scale["observer_longitude"].attrs["units"] = "degrees east WGS84"
-
-    # datetime64 can be saved to netCDF4, but Python datetime.datetime cannot
-    scale["time"] = datetime64(time)
-
-    return scale
-
-
-def doSolve(fitsfn: Path, args: str = "") -> None:
-    """
-    run Astrometry.net solve-field from Python
-    """
-
-    fitsfn = Path(fitsfn).expanduser()
-    if not fitsfn.is_file():
-        raise FileNotFoundError(f"{fitsfn} not found")
-
-    if not (solve := shutil.which("solve-field")):
-        raise FileNotFoundError("Astrometry.net solve-file exectuable not found")
-
-    # %% build command
-    cmd = [solve, "--overwrite", str(fitsfn), "--verbose"]
-    if args:
-        # if args is a string, split it. Don't append an empty space or solve-field CLI fail
-        cmd += args.split(" ")
-    print("\n", " ".join(cmd), "\n")
-    # %% execute
-    # bufsize=1: line-buffered
-    with subprocess.Popen(cmd, stdout=subprocess.PIPE, bufsize=1, text=True) as p:
-        if p.stdout:
-            for line in p.stdout:
-                print(line, end="")
-
-    if "Did not solve" in line:
-        raise RuntimeError(f"could not solve {fitsfn}")
-
-
-def fits2azel(
-    fitsfn: Path,
-    *,
-    latlon: tuple[float, float],
-    time: datetime,
-    solve: bool = False,
-    args: str = "",
-):
-    fitsfn = Path(fitsfn).expanduser()
-
-    radec = fits2radec(fitsfn, solve, args)
-
-    return radec2azel(radec, latlon, time)
diff --git a/src/astrometry_azel/io.py b/src/astrometry_azel/io.py
index 0a36f29..226f480 100755
--- a/src/astrometry_azel/io.py
+++ b/src/astrometry_azel/io.py
@@ -29,6 +29,28 @@
     loadmat = None
 
 
+def get_sources(fn: Path) -> xarray.Dataset:
+    """
+    read source (star) coordinates from .rdls file
+
+    Parameters
+    ----------
+
+    fn: pathlib.Path
+        .rdls file computed by astrometry.net solve-field (astrometry_azel.doSolve())
+
+    Returns
+    -------
+
+    radec: xarray.Dataset
+        RA, Dec of sources
+    """
+
+    fn = Path(fn).expanduser().resolve(strict=True)
+    with fits.open(fn, "readonly") as f:
+        return f[1].data
+
+
 def rgb2grey(rgb_img):
     """
     rgb_img: ndarray
@@ -38,7 +60,7 @@ def rgb2grey(rgb_img):
 
     ndim = rgb_img.ndim
     if ndim == 2:
-        logging.info("assuming its already gray since ndim=2")
+        logging.info("assuming it's already gray since ndim=2")
         grey_img = rgb_img
     elif ndim == 3 and rgb_img.shape[-1] == 3:  # this is the normal case
         grey_img = np.around(rgb_img[..., :] @ [0.299, 0.587, 0.114]).astype(rgb_img.dtype)
@@ -174,7 +196,7 @@ def write_netcdf(ds: xarray.Dataset, out_file: Path) -> None:
             "zlib": True,
             "complevel": 3,
             "fletcher32": True,
-            "chunksizes": tuple(map(lambda x: x // 2, ds[k].shape))
+            "chunksizes": tuple(map(lambda x: x // 2, ds[k].shape)),
             # arbitrary, little impact on compression
         }
 
diff --git a/src/astrometry_azel/project.py b/src/astrometry_azel/project.py
index 98c7795..47655f6 100644
--- a/src/astrometry_azel/project.py
+++ b/src/astrometry_azel/project.py
@@ -6,7 +6,7 @@
 import numpy as np
 
 from .io import load_image, write_fits, write_netcdf
-from .base import fits2azel
+from . import fits2azel
 
 import pymap3d