diff --git a/python/MANIFEST.in b/python/MANIFEST.in
new file mode 100644
index 000000000..7c1688165
--- /dev/null
+++ b/python/MANIFEST.in
@@ -0,0 +1 @@
+include juliapkg.json
diff --git a/python/README.md b/python/README.md
new file mode 100644
index 000000000..41fd4b0e7
--- /dev/null
+++ b/python/README.md
@@ -0,0 +1,40 @@
+# testparticle-jl
+
+Python wrapper for [TestParticle.jl](https://github.com/JuliaSpacePhysics/TestParticle.jl) - trace charged particles in electromagnetic fields via Julia.
+
+## Installation
+
+```bash
+pip install testparticle-jl
+```
+
+The Julia package installs automatically on first use.
+
+## Usage
+
+```python
+import testparticle as tp
+from juliacall import Main as jl
+import numpy as np
+
+# Initialize a proton
+proton = tp.Proton
+
+# Define E and B fields (simple uniform fields for demonstration)
+def get_B(x):
+    return [0.0, 0.0, 1e-8] # 10 nT in z-direction
+
+def get_E(x):
+    return [0.0, 0.0, 0.0]
+
+# Prepare the trace problem
+# Note: Complex setup might require using Julia types directly via `jl` or helper functions
+# This is a placeholder for actual usage logic dependent on TestParticle.jl API
+```
+
+See the [TestParticle.jl documentation](https://juliaspacephysics.github.io/TestParticle.jl/dev/) for full details on the underlying Julia functions.
+
+## API
+
+Available functions mirror the Julia package.
+Functions ending with `!` in Julia are available with `_b` suffix (e.g., `trace!` -> `trace_b`).
diff --git a/python/juliapkg.json b/python/juliapkg.json
new file mode 100644
index 000000000..f68f71d5b
--- /dev/null
+++ b/python/juliapkg.json
@@ -0,0 +1,9 @@
+{
+    "julia": "1.10",
+    "packages": {
+        "TestParticle": {
+            "uuid": "953b605b-f162-4481-8f7f-a191c2bb40e3",
+            "path": ".."  
+        }
+    }
+}
diff --git a/python/pyproject.toml b/python/pyproject.toml
new file mode 100644
index 000000000..08c846989
--- /dev/null
+++ b/python/pyproject.toml
@@ -0,0 +1,29 @@
+[project]
+name = "testparticle-jl"
+version = "0.1.0"
+description = "Python wrapper for TestParticle.jl - trace charged particles in electromagnetic fields via Julia"
+readme = "README.md"
+license = "MIT"
+requires-python = ">=3.9"
+keywords = ["test-particle", "julia", "space-physics", "plasma-physics"]
+classifiers = [
+    "Development Status :: 4 - Beta",
+    "Intended Audience :: Science/Research",
+    "Programming Language :: Python :: 3",
+    "Topic :: Scientific/Engineering :: Physics",
+]
+dependencies = ["juliacall>=0.9.20"]
+
+[project.urls]
+Homepage = "https://henry2004y.github.io/TestParticle.jl/dev/"
+Repository = "https://github.com/henry2004y/TestParticle.jl"
+
+[build-system]
+requires = ["pdm-backend"]
+build-backend = "pdm.backend"
+
+[tool.pdm.build]
+includes = ["src/testparticle"]
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
diff --git a/python/src/testparticle/__init__.py b/python/src/testparticle/__init__.py
new file mode 100644
index 000000000..989263e88
--- /dev/null
+++ b/python/src/testparticle/__init__.py
@@ -0,0 +1,110 @@
+"""Python wrapper for TestParticle.jl."""
+
+from juliacall import Main as jl
+
+# Load the Julia package
+jl.seval("using TestParticle")
+
+# Re-export Julia functions and types
+prepare = jl.TestParticle.prepare
+prepare_gc = jl.TestParticle.prepare_gc
+get_gc = jl.TestParticle.get_gc
+get_gc_func = jl.TestParticle.get_gc_func
+
+# Functions with ! in Julia are accessed with _b suffix in Python via juliacall
+trace_b = jl.TestParticle.trace_b
+trace_relativistic_b = jl.TestParticle.trace_relativistic_b
+trace_normalized_b = jl.TestParticle.trace_normalized_b
+trace_relativistic_normalized_b = jl.TestParticle.trace_relativistic_normalized_b
+
+trace = jl.TestParticle.trace
+trace_relativistic = jl.TestParticle.trace_relativistic
+trace_normalized = jl.TestParticle.trace_normalized
+trace_relativistic_normalized = jl.TestParticle.trace_relativistic_normalized
+
+trace_gc_b = jl.TestParticle.trace_gc_b
+trace_gc_drifts_b = jl.TestParticle.trace_gc_drifts_b
+trace_gc_flr_b = jl.TestParticle.trace_gc_flr_b
+trace_gc_exb_b = jl.TestParticle.trace_gc_exb_b
+trace_fieldline_b = jl.TestParticle.trace_fieldline_b
+trace_fieldline = jl.TestParticle.trace_fieldline
+
+get_gc_velocity = jl.TestParticle.get_gc_velocity
+full_to_gc = jl.TestParticle.full_to_gc
+gc_to_full = jl.TestParticle.gc_to_full
+
+Proton = jl.TestParticle.Proton
+Electron = jl.TestParticle.Electron
+Ion = jl.TestParticle.Ion
+
+Maxwellian = jl.TestParticle.Maxwellian
+BiMaxwellian = jl.TestParticle.BiMaxwellian
+Kappa = jl.TestParticle.Kappa
+BiKappa = jl.TestParticle.BiKappa
+
+AdaptiveBoris = jl.TestParticle.AdaptiveBoris
+AdaptiveHybrid = jl.TestParticle.AdaptiveHybrid
+
+CurrentLoop = jl.TestParticle.CurrentLoop
+getB_loop = jl.TestParticle.getB_loop
+
+get_gyrofrequency = jl.TestParticle.get_gyrofrequency
+get_gyroperiod = jl.TestParticle.get_gyroperiod
+get_gyroradius = jl.TestParticle.get_gyroradius
+get_velocity = jl.TestParticle.get_velocity
+get_energy = jl.TestParticle.get_energy
+get_mean_magnitude = jl.TestParticle.get_mean_magnitude
+energy2velocity = jl.TestParticle.energy2velocity
+get_curvature_radius = jl.TestParticle.get_curvature_radius
+get_adiabaticity = jl.TestParticle.get_adiabaticity
+
+sample_unit_sphere = jl.TestParticle.sample_unit_sphere
+get_number_density_flux = jl.TestParticle.get_number_density_flux
+
+getB_zpinch = jl.TestParticle.getB_zpinch
+getB_bottle = jl.TestParticle.getB_bottle
+getB_mirror = jl.TestParticle.getB_mirror
+getB_tokamak_coil = jl.TestParticle.getB_tokamak_coil
+
+orbit = jl.TestParticle.orbit
+monitor = jl.TestParticle.monitor
+
+get_fields = jl.TestParticle.get_fields
+get_work = jl.TestParticle.get_work
+
+LazyTimeInterpolator = jl.TestParticle.LazyTimeInterpolator
+
+TraceProblem = jl.TestParticle.TraceProblem
+TraceGCProblem = jl.TestParticle.TraceGCProblem
+TraceHybridProblem = jl.TestParticle.TraceHybridProblem
+
+CartesianGrid = jl.TestParticle.CartesianGrid
+RectilinearGrid = jl.TestParticle.RectilinearGrid
+StructuredGrid = jl.TestParticle.StructuredGrid
+
+
+__version__ = "0.1.0"
+
+__all__ = [
+    "prepare", "prepare_gc", "get_gc", "get_gc_func",
+    "trace_b", "trace_relativistic_b", "trace_normalized_b", "trace_relativistic_normalized_b",
+    "trace", "trace_relativistic", "trace_normalized", "trace_relativistic_normalized",
+    "trace_gc_b",
+    "trace_gc_drifts_b", "trace_gc_flr_b", "trace_gc_exb_b", "trace_fieldline_b", "trace_fieldline",
+    "get_gc_velocity", "full_to_gc", "gc_to_full",
+    "Proton", "Electron", "Ion",
+    "Maxwellian", "BiMaxwellian", "Kappa", "BiKappa",
+    "AdaptiveBoris", "AdaptiveHybrid",
+    "CurrentLoop", "getB_loop",
+    "get_gyrofrequency",
+    "get_gyroperiod", "get_gyroradius", "get_velocity", "get_energy", "get_mean_magnitude",
+    "energy2velocity", "get_curvature_radius", "get_adiabaticity",
+    "sample_unit_sphere", "get_number_density_flux",
+    "getB_zpinch", "getB_bottle", "getB_mirror", "getB_tokamak_coil",
+    "orbit", "monitor",
+    "get_fields", "get_work",
+    "LazyTimeInterpolator",
+    "TraceProblem", "TraceGCProblem", "TraceHybridProblem",
+    "CartesianGrid", "RectilinearGrid", "StructuredGrid",
+    "__version__", "jl"
+]
diff --git a/python/tests/test_wrapper.py b/python/tests/test_wrapper.py
new file mode 100644
index 000000000..a270dc938
--- /dev/null
+++ b/python/tests/test_wrapper.py
@@ -0,0 +1,29 @@
+import pytest
+
+def test_import():
+    import testparticle as tp
+    assert hasattr(tp, "trace")
+    assert hasattr(tp, "trace_b")
+    assert hasattr(tp, "Proton")
+    assert hasattr(tp, "Electron")
+
+def test_basic_types():
+    import testparticle as tp
+    from juliacall import Main as jl
+    
+    # Check constants
+    m_p = tp.Proton.mass
+    q_p = tp.Proton.charge
+    assert m_p > 0
+    assert q_p > 0
+
+def test_gyrofrequency():
+    import testparticle as tp
+    
+    # B = 1.0 (scalar or magnitude), q = 1.0, m = 1.0
+    # The signature in Julia is usually get_gyrofrequency(B, species) or (B, q, m)
+    # Let's check (B, q, m)
+    
+    omega = tp.get_gyrofrequency(1.0, 1.0, 1.0)
+    assert abs(omega - 1.0) < 1e-6
+