diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml
index fa9829e..e31aa05 100644
--- a/.github/workflows/rustbca_compile_check.yml
+++ b/.github/workflows/rustbca_compile_check.yml
@@ -27,7 +27,7 @@ jobs:
         sudo apt-get install rustc
     - name: Install pip for Python-3
       run: |
-        sudo apt-get install python3-pip
+        sudo apt-get install python3-pip python3-dev
     - name: Install Python libraries
       run: |
         python3 -m pip install numpy shapely scipy
@@ -45,9 +45,16 @@ jobs:
         python3 setup.py install --root .
         cp -r usr/local/lib/python3.8/dist-packages/libRustBCA .
         python3 -c "from libRustBCA.pybca import *;"
-    - name: Test Fortran bindings
+    - name: Test Fortran and C bindings
       run : |
+        cargo build --release
+        cp examples/test_rustbca.f90 .
         gfortran -c rustbca.f90 target/release/liblibRustBCA.so
+        gfortran test_rustbca.f90 rustbca.f90 target/release/liblibRustBCA.so
+        ./a.out
+        cp examples/RustBCA.c .
+        g++ RustBCA.c RustBCA.h target/release/liblibRustBCA.so -Iexamples/ -I/usr/include/python3.8
+        ./a.out
     - name: Test RustBCA
       run: |
         sudo cargo test --features cpr_rootfinder_netlib,hdf5_input,distributions,parry3d
diff --git a/examples/RustBCA.c b/examples/RustBCA.c
new file mode 100644
index 0000000..a218531
--- /dev/null
+++ b/examples/RustBCA.c
@@ -0,0 +1,54 @@
+#include "RustBCA.h"
+#include <iostream>
+#include <vector>
+
+int main(int argc, char * argv[]) {
+  OutputTaggedBCA output;
+  double velocities[2][3] = {{500000.0, 0.1, 0.0}, {500000.0, 0.1, 0.0}};
+  double positions[2][3] = {{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}};
+  int tags[2] = {0, 1};
+  double weights[2] = {1.0, 1.0};
+  double Z[3] = {74.0, 74.0};
+  double m[2] = {184.0, 184.0};
+  double n[2] = {0.06306, 0.06306};
+  double Ec[2] = {1.0, 1.0};
+  double Es[2] = {8.79, 8.79};
+  double Eb[2] = {0.0, 0.0};
+
+  InputTaggedBCA input = {
+    2,
+    positions,
+    velocities,
+    1.0,
+    1.0,
+    1.0,
+    1.0,
+    2,
+    Z,
+    m,
+    n,
+    Ec,
+    Es,
+    Eb,
+    tags,
+    weights
+  };
+
+  //output = simple_bca_c(0., 0., 0., 0.5, 0.5, 0.00, 2000.0, 2.0, 4.0, 1.0, 0.0, 74.0, 184.0, 1.0, 8.79, 0.06306, 0.0);
+  //output = compound_bca_list_c(input);
+  output = compound_tagged_bca_list_c(input);
+
+  std::cout << "Particle 1 Z: ";
+  std::cout << output.particles[0][0];
+  std::cout << std::endl;
+  std::cout << "Particle 1 E [eV]: ";
+  std::cout << output.particles[0][2];
+  std::cout << std::endl;
+  std::cout << "Particle 2 Z: ";
+  std::cout << output.particles[1][0];
+  std::cout << std::endl;
+  std::cout << "Particle 2 E [eV]: ";
+  std::cout << output.particles[1][2];
+  std::cout << std::endl;
+  return 0;
+}
diff --git a/examples/test_rustbca.f90 b/examples/test_rustbca.f90
new file mode 100644
index 0000000..c0b8cae
--- /dev/null
+++ b/examples/test_rustbca.f90
@@ -0,0 +1,75 @@
+
+program test_rustbca
+
+    use rustbca
+    use, intrinsic :: iso_c_binding
+
+    integer :: N_ions
+    real(c_double), allocatable, dimension(:) :: ux, uy, uz, E, Z1, m1, Ec1, Es1
+    integer(c_int) :: num_species_target, num_emitted_particles
+    real(c_double), target :: Z2(2), m2(2), Ec2(2), Es2(2), Eb2(2), n2(2)
+    real(c_double) :: ux1, uy1, uz1, E1
+    type(c_ptr) :: bca_output_c
+    real(c_double), pointer, dimension(:,:) :: bca_output_f
+    real :: start, stop
+    logical(c_bool) :: track_recoils
+
+    !Initial ion conditions
+    N_ions = 100000
+    allocate(ux(N_ions), uy(N_ions), uz(N_ions), E(N_ions), Z1(N_ions), m1(N_ions), Ec1(N_ions), Es1(N_ions))
+    ux(:) = 0.999
+    uy(:) = sqrt(1.0 - 0.999*0.999)
+    uz(:) = 0.0
+    E(:) = 1000.0_8
+
+    !Hydrogen
+    Z1(:) = 1.0_8
+    m1(:) = 1.008_8
+    Ec1(:) = 1.0_8
+    Es1(:) = 1.5_8
+
+    !Titanium Hydride
+    num_species_target = 2
+    Z2(1) = 22.0_8
+    m2(1) = 47.867_8
+    Ec2(1) = 4.84_8
+    Es2(1) = 4.84_8
+    Eb2(1) = 3.0_8
+    n2(1) = 0.04527_8
+
+    Z2(2) = 1.0_8
+    m2(2) = 1.008_8
+    Ec2(2) = 1.5_8
+    Es2(2) = 1.5_8
+    Eb2(2) = 0.0_8
+    n2(2) = 0.09054_8
+
+    track_recoils = .false.
+
+    call cpu_time(start)
+    bca_output_c = compound_bca_list_fortran(N_ions, track_recoils, ux, uy, uz, E, &
+        Z1, m1, Ec1, Es1, &
+        num_species_target, Z2, m2, Ec2, Es2, Eb2, n2, &
+        num_emitted_particles)
+    call c_f_pointer(bca_output_c, bca_output_f, [num_emitted_particles, 6])
+    call cpu_time(stop)
+
+    write(*,*) "Elapsed time in seconds per ion per eV: ", (stop - start)/N_ions/1000.0
+
+    !write(*,*) bca_output_f
+
+    call cpu_time(start)
+    do i = 0, N_ions
+        !Test reflect_single_ion routine
+        ux1 = 0.999
+        uy1 = sqrt(1.0 - 0.999*0.999)
+        uz1 = 0.0
+        E1 = 1000.0
+        call reflect_single_ion_c(num_species_target, ux1, uy1, uz1, E1, Z1(1), m1(1), Ec1(1), Es1(1), Z2, m2, Ec2, Es2, Eb2, n2)
+    end do
+    call cpu_time(stop)
+    write(*,*) "Elapsed time in ions per eV per s: ", (stop - start)/N_ions/1000.0
+
+    !call exit(1)
+    
+end program test_rustbca