diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
new file mode 100644
index 0000000..1bae4be
--- /dev/null
+++ b/.github/workflows/CI.yml
@@ -0,0 +1,120 @@
+# This file is autogenerated by maturin v1.4.0
+# To update, run
+#
+#    maturin generate-ci github
+#
+name: CI
+
+on:
+  push:
+    branches:
+      - main
+      - master
+    tags:
+      - '*'
+  pull_request:
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+jobs:
+  linux:
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        target: [x86_64, x86, aarch64, armv7, s390x, ppc64le]
+    steps:
+      - uses: actions/checkout@v3
+      - uses: actions/setup-python@v4
+        with:
+          python-version: '3.10'
+      - name: Build wheels
+        uses: PyO3/maturin-action@v1
+        with:
+          target: ${{ matrix.target }}
+          args: --release --out dist --find-interpreter
+          sccache: 'true'
+          manylinux: auto
+      - name: Upload wheels
+        uses: actions/upload-artifact@v3
+        with:
+          name: wheels
+          path: dist
+
+  windows:
+    runs-on: windows-latest
+    strategy:
+      matrix:
+        target: [x64, x86]
+    steps:
+      - uses: actions/checkout@v3
+      - uses: actions/setup-python@v4
+        with:
+          python-version: '3.10'
+          architecture: ${{ matrix.target }}
+      - name: Build wheels
+        uses: PyO3/maturin-action@v1
+        with:
+          target: ${{ matrix.target }}
+          args: --release --out dist --find-interpreter
+          sccache: 'true'
+      - name: Upload wheels
+        uses: actions/upload-artifact@v3
+        with:
+          name: wheels
+          path: dist
+
+  macos:
+    runs-on: macos-latest
+    strategy:
+      matrix:
+        target: [x86_64, aarch64]
+    steps:
+      - uses: actions/checkout@v3
+      - uses: actions/setup-python@v4
+        with:
+          python-version: '3.10'
+      - name: Build wheels
+        uses: PyO3/maturin-action@v1
+        with:
+          target: ${{ matrix.target }}
+          args: --release --out dist --find-interpreter
+          sccache: 'true'
+      - name: Upload wheels
+        uses: actions/upload-artifact@v3
+        with:
+          name: wheels
+          path: dist
+
+  sdist:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v3
+      - name: Build sdist
+        uses: PyO3/maturin-action@v1
+        with:
+          command: sdist
+          args: --out dist
+      - name: Upload sdist
+        uses: actions/upload-artifact@v3
+        with:
+          name: wheels
+          path: dist
+
+  release:
+    name: Release
+    runs-on: ubuntu-latest
+    if: "startsWith(github.ref, 'refs/tags/')"
+    needs: [linux, windows, macos, sdist]
+    steps:
+      - uses: actions/download-artifact@v3
+        with:
+          name: wheels
+      - name: Publish to PyPI
+        uses: PyO3/maturin-action@v1
+        env:
+          MATURIN_PYPI_TOKEN: ${{ secrets.PYPI_API_TOKEN }}
+        with:
+          command: upload
+          args: --non-interactive --skip-existing *
diff --git a/.gitignore b/.gitignore
index 571ce1c..035bd51 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,76 @@
-/target/
+/target
+
+# Byte-compiled / optimized / DLL files
+__pycache__/
+.pytest_cache/
+*.py[cod]
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+.venv/
+env/
+bin/
+build/
+develop-eggs/
+dist/
+eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+include/
+man/
+venv/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+pip-selfcheck.json
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.cache
+nosetests.xml
+coverage.xml
+
+# Translations
+*.mo
+
+# Mr Developer
+.mr.developer.cfg
+.project
+.pydevproject
+
+# Rope
+.ropeproject
+
+# Django stuff:
+*.log
+*.pot
+
+.DS_Store
+
+# Sphinx documentation
+docs/_build/
+
+# PyCharm
+.idea/
+
+# VSCode
+.vscode/
+
+# Pyenv
+.python-version
+
 **/*.rs.bk
 Cargo.lock
 .criterion
-.vscode/
-.DS_Store
diff --git a/Cargo.toml b/Cargo.toml
index 25c7684..9d48faf 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,17 +1,22 @@
 [package]
 edition = "2021"
-name = "score-rs"
-version = "0.0.0"
+name = "score_rs"
+version = "0.0.1"
 authors = ["Moritz Althaus <moritzalthaus@gmx.de>"]
 license = "MIT"
 repository = "https://github.com/moldhouse/score-rs.git"
 homepage = "https://github.com/moldhouse/score-rs.git"
 readme = "README.md"
 
+[lib]
+name = "score_rs"
+crate-type = ["cdylib"]
+
 [features]
 default = ["rayon"]
 
 [dependencies]
+pyo3 = "0.20.0"
 cfg-if = "1.0.0"
 failure = "^0.1.1"
 flat_projection = "0.4.0"
@@ -20,13 +25,9 @@ ordered-float = "2.0.1"
 ord_subset = "^3.1.0"
 rayon = { version = "^1.0", optional = true }
 itertools = "0.10.0"
+numpy = "0.20.0"
 
 [dev-dependencies]
 assert_approx_eq = "^1.0.0"
-criterion = "^0.3.0"
 igc = "0.2.2"
 env_logger = "0.8.2"
-
-[[bench]]
-name = "free"
-harness = false
\ No newline at end of file
diff --git a/README.md b/README.md
index 9b282b4..adc7129 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # score-rs
 
-Find the seven GPS points of a single flight whose straight connection gives the maximum total distance.
+Find `n` points out of a set of possible tens of thousands of GPS points, such that the straight distance between them is maximized.
 There is one constraint: The finish altitude must not be more than 1000 m less than the start altitude.
 
 The algorithm does the same optimization that [WeGlide](https://www.weglide.org) does to assign a distance to every flight:
@@ -17,14 +17,38 @@ The code is based on the excellent [aeroscore-rs](https://github.com/glide-rs/ae
 2. If the 1000 m altitude is satisfied by the best result, the optimization is similar. If not, this library uses a caching system to quickly determine if start candidates can give a better solution than the current best without traversing the whole graph.
 3. Also look for potential solutions by adjusting the start- and end points of a given solution and keeping the middle points constant. This is not used to find the actual solution (as it does not guarantee optimality), but it speeds up the optimization by helping to find better intermediate results and discard candidates that do not offer a better solution
 
-## Test
+## Develop
+
+Python bindings are generated with [maturin](https://github.com/PyO3/maturin/). Create a virtual env first with
 
 ```bash
-cargo test
+python -m venv ./env && source .env/bin/activate
 ```
 
-## Bench
+and install numpy in your virtual env:
 
 ```bash
-cargo bench
+pip install numpy
+```
+
+To develop, run
+
+```bash
+maturin develop
+``` 
+
+or for a (faster) release version
+
+```bash
+maturin build --release
+pip install .
+```
+
+## Test
+
+You can run the tests with
+
+```bash
+cargo test
+python -m pytest
 ```
\ No newline at end of file
diff --git a/benches/free.rs b/benches/free.rs
deleted file mode 100644
index 682847f..0000000
--- a/benches/free.rs
+++ /dev/null
@@ -1,51 +0,0 @@
-#[macro_use]
-extern crate criterion;
-
-extern crate igc;
-
-use criterion::Criterion;
-use igc::util::Time;
-use score_rs::free;
-use score_rs::point::PointImpl;
-
-const LEGS: usize = 6;
-
-fn criterion_benchmark(c: &mut Criterion) {
-    c.bench_function("free", |b| {
-        b.iter(|| {
-            let release = Time::from_hms(8, 12, 29);
-            let fixes = include_str!("../tests/fixtures/2023-06-17_288167.igc")
-                .lines()
-                .filter(|l| l.starts_with('B'))
-                .filter_map(|line| {
-                    igc::records::BRecord::parse(&line)
-                        .ok()
-                        .map_or(None, |record| {
-                            if record.timestamp.seconds_since_midnight()
-                                >= release.seconds_since_midnight()
-                            {
-                                Some(PointImpl {
-                                    latitude: record.pos.lat.into(),
-                                    longitude: record.pos.lon.into(),
-                                    altitude: record.pressure_alt,
-                                })
-                            } else {
-                                None
-                            }
-                        })
-                })
-                .collect::<Vec<_>>();
-
-            free::optimize(&fixes, 0.0, LEGS).unwrap()
-        })
-    });
-}
-
-criterion_group! {
-    name = benches;
-    config = Criterion::default()
-        .sample_size(10);
-
-    targets = criterion_benchmark
-}
-criterion_main!(benches);
diff --git a/tests/fixtures/2023-06-17_288167.igc b/fixtures/2023-06-17_288167.igc
similarity index 100%
rename from tests/fixtures/2023-06-17_288167.igc
rename to fixtures/2023-06-17_288167.igc
diff --git a/tests/fixtures/schunk_1000m.igc b/fixtures/schunk_1000m.igc
similarity index 100%
rename from tests/fixtures/schunk_1000m.igc
rename to fixtures/schunk_1000m.igc
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..6e5f8cd
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,16 @@
+[build-system]
+requires = ["maturin>=1.4,<2.0"]
+build-backend = "maturin"
+
+[project]
+name = "score-rs"
+requires-python = ">=3.8"
+classifiers = [
+    "Programming Language :: Rust",
+    "Programming Language :: Python :: Implementation :: CPython",
+    "Programming Language :: Python :: Implementation :: PyPy",
+]
+dynamic = ["version"]
+
+[tool.maturin]
+features = ["pyo3/extension-module"]
diff --git a/pyscore-rs/free_test.py b/pyscore-rs/free_test.py
new file mode 100644
index 0000000..6fe218f
--- /dev/null
+++ b/pyscore-rs/free_test.py
@@ -0,0 +1,59 @@
+from typing import NamedTuple
+import datetime as dt
+
+import score_rs
+import numpy as np
+from numpy.testing import assert_almost_equal
+
+
+class Fix(NamedTuple):
+    lon: float
+    lat: float
+    pressure_alt: int
+    time: int
+
+
+def parse_line(line: str) -> Fix:
+    time = ((int(line[1:3]) * 60) + int(line[3:5])) * 60 + int(line[5:7])
+    lat = float(line[7:9]) + (float(line[9:11]) + float(line[11:14]) / 1000.0) / 60.0
+    if line[14] == "S":
+        lat = -lat
+
+    lon = float(line[15:18]) + (float(line[18:20]) + float(line[20:23]) / 1000.0) / 60.0
+    if line[23] == "W":
+        lon = -lon
+
+    pressure_alt = int(float(line[25:30]))
+    return Fix(lon, lat, pressure_alt, time)
+
+
+def read_igc(
+    file_path: str, release: dt.time
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    lon, lat, pressure_alt = [], [], []
+    release_seconds_since_midnight = (
+        release.hour * 60 + release.minute
+    ) * 60 + release.second
+    with open(file_path, "r") as file:
+        for line in file.read().split("\n"):
+            if line.startswith("B"):
+                fix = parse_line(line)
+                if fix.time < release_seconds_since_midnight:
+                    continue
+                lon.append(fix.lon)
+                lat.append(fix.lat)
+                pressure_alt.append(fix.pressure_alt)
+
+    return np.array(lon), np.array(lat), np.array(pressure_alt)
+
+
+def test_free():
+    release = dt.time(8, 12, 29)
+    data = read_igc("fixtures/2023-06-17_288167.igc", release)
+    res = score_rs.optimize(data[0], data[1], data[2], 6)
+    assert_almost_equal(res[1], 1018.54, 2)
+    assert res[0] == [0, 936, 2847, 3879, 5048, 7050, 8128]
+
+    res = score_rs.optimize(data[0], data[1], data[2], 2)
+    assert_almost_equal(res[1], 804.95, 2)
+    assert res[0] == [887, 3886, 7801]
diff --git a/src/free.rs b/src/free.rs
index a3dfb1d..fc10f6d 100644
--- a/src/free.rs
+++ b/src/free.rs
@@ -613,3 +613,56 @@ fn calculate_distance<T: Point>(points: &[T], path: &Path) -> f32 {
         .map(|(fix1, fix2)| vincenty_distance(fix1, fix2))
         .sum()
 }
+
+#[cfg(test)]
+mod tests {
+    use crate::free;
+    use crate::free::OptimizationResult;
+    use crate::point::PointImpl;
+    use assert_approx_eq::assert_approx_eq;
+    use igc::records::BRecord;
+    use igc::util::Time;
+
+    const LEGS: usize = 6;
+
+    #[test]
+    fn free_distance() {
+        let release = Time::from_hms(8, 12, 29);
+        let result = run_free_test(include_str!("../fixtures/2023-06-17_288167.igc"), release);
+        assert_approx_eq!(result.distance, 1018.5, 0.1);
+        assert_eq!(result.path, [0, 936, 2847, 3879, 5048, 7050, 8128]);
+    }
+
+    #[test]
+    fn free_distance_with_1000m() {
+        let release = Time::from_hms(8, 16, 30);
+        let result = run_free_test(include_str!("../fixtures/schunk_1000m.igc"), release);
+        assert_approx_eq!(result.distance, 1158.61, 0.1);
+        assert_eq!(result.path, [335, 10099, 14740, 15482, 24198, 34160, 35798]);
+    }
+
+    fn run_free_test(file: &str, release: Time) -> OptimizationResult {
+        env_logger::try_init().ok();
+
+        let fixes = file
+            .lines()
+            .filter(|l| l.starts_with('B'))
+            .filter_map(|line| {
+                BRecord::parse(&line).ok().map_or(None, |record| {
+                    if record.timestamp.seconds_since_midnight() >= release.seconds_since_midnight()
+                    {
+                        Some(PointImpl {
+                            latitude: record.pos.lat.into(),
+                            longitude: record.pos.lon.into(),
+                            altitude: record.pressure_alt,
+                        })
+                    } else {
+                        None
+                    }
+                })
+            })
+            .collect::<Vec<_>>();
+
+        free::optimize(&fixes, 0.0, LEGS).unwrap()
+    }
+}
diff --git a/src/lib.rs b/src/lib.rs
index d035e16..9d0f2ad 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,6 +1,36 @@
+use numpy::PyReadonlyArray1;
+use pyo3::prelude::*;
+
 pub mod flat;
 pub mod free;
 pub mod parallel;
 pub mod point;
 pub mod utils;
-pub mod vincenty;
\ No newline at end of file
+pub mod vincenty;
+
+#[pymodule]
+fn score_rs(_py: Python, m: &PyModule) -> PyResult<()> {
+    #[pyfn(m)]
+    #[pyo3(name = "optimize")]
+    fn optimize_py<'py>(
+        longitude: PyReadonlyArray1<'py, f64>,
+        latitude: PyReadonlyArray1<'py, f64>,
+        alt: PyReadonlyArray1<'py, i64>,
+        legs: usize,
+    ) -> PyResult<(Vec<usize>, f32)> {
+        let mut points = Vec::new();
+        let longitude = longitude.as_slice().unwrap();
+        let latitude = latitude.as_slice().unwrap();
+        let alt = alt.as_slice().unwrap();
+        for i in 0..longitude.len() {
+            points.push(point::PointImpl {
+                longitude: longitude[i] as f32,
+                latitude: latitude[i] as f32,
+                altitude: alt[i] as i16,
+            });
+        }
+        let result = free::optimize(&points, 0.0, legs).unwrap();
+        Ok((result.path, result.distance))
+    }
+    Ok(())
+}
diff --git a/tests/free.rs b/tests/free.rs
deleted file mode 100644
index abad4a7..0000000
--- a/tests/free.rs
+++ /dev/null
@@ -1,53 +0,0 @@
-#[macro_use]
-extern crate assert_approx_eq;
-extern crate igc;
-
-use igc::util::Time;
-use score_rs::free;
-use score_rs::free::OptimizationResult;
-use score_rs::point::PointImpl;
-
-const LEGS: usize = 6;
-
-#[test]
-fn free_distance() {
-    let release = Time::from_hms(8, 12, 29);
-    let result = run_free_test(include_str!("fixtures/2023-06-17_288167.igc"), release);
-    assert_approx_eq!(result.distance, 1018.5, 0.1);
-    assert_eq!(result.path, [0, 936, 2847, 3879, 5048, 7050, 8128]);
-}
-
-#[test]
-fn free_distance_with_1000m() {
-    let release = Time::from_hms(8, 16, 30);
-    let result = run_free_test(include_str!("fixtures/schunk_1000m.igc"), release);
-    assert_approx_eq!(result.distance, 1158.61, 0.1);
-    assert_eq!(result.path, [335, 10099, 14740, 15482, 24198, 34160, 35798]);
-}
-
-fn run_free_test(file: &str, release: Time) -> OptimizationResult {
-    env_logger::try_init().ok();
-
-    let fixes = file
-        .lines()
-        .filter(|l| l.starts_with('B'))
-        .filter_map(|line| {
-            igc::records::BRecord::parse(&line)
-                .ok()
-                .map_or(None, |record| {
-                    if record.timestamp.seconds_since_midnight() >= release.seconds_since_midnight()
-                    {
-                        Some(PointImpl {
-                            latitude: record.pos.lat.into(),
-                            longitude: record.pos.lon.into(),
-                            altitude: record.pressure_alt,
-                        })
-                    } else {
-                        None
-                    }
-                })
-        })
-        .collect::<Vec<_>>();
-
-    free::optimize(&fixes, 0.0, LEGS).unwrap()
-}