diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 5421b8e..ee1b564 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -3,4 +3,9 @@ add_executable(basic basic.cpp)
 target_link_libraries(basic PUBLIC tensor)
 target_include_directories(basic PUBLIC "${CMAKE_SOURCE_DIR}")
 
-add_test(NAME basicTest COMMAND basic)
+add_test(NAME basic-test COMMAND basic)
+
+add_executable(two-flavour-osc two-flavour-osc.cpp)
+target_link_libraries(two-flavour-osc PUBLIC propagator)
+
+add_test(NAME two-flavour-osc-test COMMAND two-flavour-osc)
diff --git a/tests/two-flavour-osc.cpp b/tests/two-flavour-osc.cpp
new file mode 100644
index 0000000..f593853
--- /dev/null
+++ b/tests/two-flavour-osc.cpp
@@ -0,0 +1,78 @@
+#include <nuTens/propagator/vacuum-propagator.cpp>
+
+// Get absolute relative difference between two floats:
+//   | (f1 - f2) / f1 |
+float relativeDiff(float f1, float f2){
+    return std::abs((f1 - f2) / f1);
+}
+
+int main(){
+
+    float m1 = 0.1, m2 = 0.5;
+    float energy = 1.0;
+    float baseline = 0.5;
+
+    Tensor masses;
+    masses.ones({1, 2}, NTdtypes::kFloat).requiresGrad(false);
+    masses.setValue({0, 0}, m1);
+    masses.setValue({0, 1}, m2);
+    masses.requiresGrad(true);
+
+    Tensor energies;
+    energies.ones({1, 1}, NTdtypes::kFloat).requiresGrad(false);
+    energies.setValue({0, 0}, energy);
+    energies.requiresGrad(true);
+
+    VacuumPropagator propagator(2, baseline);
+    propagator.setMasses(masses);
+    propagator.setEnergies(energies);
+
+    
+    float theta = -M_PI;
+
+    // test that VacuumPropagator gives expected oscillation probabilites for a range of thetas
+    for( int i = 0; i < 20; i++){
+
+        theta += ( 2* M_PI / 20);
+
+        // construct the PMNS matrix for current theta value
+        Tensor PMNS;
+        PMNS.ones({1, 2, 2}, NTdtypes::kComplexFloat).requiresGrad(false);
+        PMNS.setValue({0, 0,0}, std::cos(theta));
+        PMNS.setValue({0, 0,1}, std::sin(theta));
+        PMNS.setValue({0, 1,0}, -std::sin(theta));
+        PMNS.setValue({0, 1,1}, std::cos(theta));
+        PMNS.requiresGrad(true);
+
+        propagator.setPMNS(PMNS);
+
+        Tensor probabilities = propagator.calculateProbs();
+
+        float sin2Theta = std::sin(2.0 * theta);
+        float sinPhi = std::sin((m1 *m1 - m2 * m2) * baseline / (4.0 * energy) );
+
+        float offAxis = sin2Theta * sin2Theta * sinPhi * sinPhi;
+        float onAxis = 1.0 - offAxis;
+
+        if ((relativeDiff(probabilities.getValue<float>({0, 0,0}), onAxis) > 0.00005)
+            || (relativeDiff(probabilities.getValue<float>({0, 1,1}), onAxis) > 0.00005)
+            || (relativeDiff(probabilities.getValue<float>({0, 0,1}), offAxis) > 0.00005)
+            || (relativeDiff(probabilities.getValue<float>({0, 1,0}), offAxis) > 0.00005)
+        ){
+            std::cerr << "ERROR: 2 flavour probabilities from VacuumPropagator do not match expected values for theta = " << theta << std::endl;
+            std::cerr << "Calculated probabilities: " << std::endl;
+            std::cerr << probabilities << std::endl;
+
+            std::cerr << std::endl;
+            std::cerr << "'True' probabilities: " << std::endl;
+            std::cerr << "  " << onAxis << "  " << offAxis << std::endl;
+            std::cerr << "  " << offAxis << "  " << onAxis << std::endl;
+            std::cerr << std::endl;
+            std::cerr << "sin(2 theta) = " << sin2Theta << std::endl;
+            std::cerr << "sin(Phi) = " << sinPhi << std::endl;
+
+            return 1;
+        }
+    }
+
+}
\ No newline at end of file