diff --git a/CMakeLists.txt b/CMakeLists.txt
index 97e21cde424a571fbda61d7851ca91e32f4a4a30..fe40da359def131c830b79d6f7779691b5c6ec5f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -23,12 +23,6 @@ add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/externals/autodiff")
 add_executable(main "${CMAKE_CURRENT_SOURCE_DIR}/src/main.cpp")
 target_link_libraries(main autodiff)
 
-add_executable(trace "${CMAKE_CURRENT_SOURCE_DIR}/src/trace_critical.cpp")
-if (MSVC)
-    target_sources(trace PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen/debug/msvc/eigen.natvis")
-endif()
-target_link_libraries(trace autodiff)
-
 add_executable(multifluid "${CMAKE_CURRENT_SOURCE_DIR}/src/multifluid.cpp")
 if (MSVC)
     target_sources(multifluid PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen/debug/msvc/eigen.natvis")
diff --git a/include/teqp/algorithms/critical_tracing.hpp b/include/teqp/algorithms/critical_tracing.hpp
index 7ec8ad98999fd92a930be38e23f4b7740e85817d..44baf516ff6ffed12169fe7cf82651fedbec3ecf 100644
--- a/include/teqp/algorithms/critical_tracing.hpp
+++ b/include/teqp/algorithms/critical_tracing.hpp
@@ -1,5 +1,7 @@
 #pragma once
 
+#include <fstream>
+
 #include <Eigen/Dense>
 #include "teqp/algorithms/rootfinding.hpp"
 
@@ -268,7 +270,7 @@ auto critical_polish_molefrac(const ModelType &model, const double T, const VecT
 }
 
 template<typename ModelType, typename VecType>
-void trace_critical_arclength_binary(const ModelType& model, VecType rhovec0, double T0, const std::string &filename) {
+void trace_critical_arclength_binary(const ModelType& model, double T0, const VecType &rhovec0, const std::string &filename) {
 
     double t = 0.0, dt = 100;
     std::valarray<double> last_drhodt;
@@ -292,7 +294,9 @@ void trace_critical_arclength_binary(const ModelType& model, VecType rhovec0, do
             std::cout << sout;
         };
         if (iter == 0) {
-            write_line();
+            if (!filename.empty()){
+                write_line();
+            }
         }
 
         auto drhodT = get_drhovec_dT_crit(model, T, rhovec);
@@ -332,6 +336,8 @@ void trace_critical_arclength_binary(const ModelType& model, VecType rhovec0, do
             break;
         }
         last_drhodt = c * drhodt;
-        write_line();
+        if (!filename.empty()) {
+            write_line();
+        }
     }
 }
\ No newline at end of file
diff --git a/src/multifluid.cpp b/src/multifluid.cpp
index 499229ec0ec2b1b7413568eb6b74e16df9ea68dc..6b90279e92b10bd8825f0c88c3bdaadc9235d368 100644
--- a/src/multifluid.cpp
+++ b/src/multifluid.cpp
@@ -57,7 +57,7 @@ void trace_critical_loci(const std::string &coolprop_root, const nlohmann::json
                 rhoc0 = rhovec.sum();
             }
             std::string filename = pp[0] + "_" + pp[1] + ".csv";
-            trace_critical_arclength_binary(model, rhovec, T0, filename);
+            trace_critical_arclength_binary(model, T0, rhovec, filename);
         }
     }
 }
diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx
index 240d88548855888170b58d2e8dbe2a743bf2f0a8..953bf365f2fbd6b7a6bfd3acfa236250835a9a8c 100644
--- a/src/tests/catch_tests.cxx
+++ b/src/tests/catch_tests.cxx
@@ -128,4 +128,21 @@ TEST_CASE("Check p three ways for vdW", "[virial][p]")
 
     CHECK(std::abs(pfromderiv - pexact)/pexact < 1e-15);
     CHECK(std::abs(pvir - pexact)/pexact < 1e-8);
+}
+
+TEST_CASE("Trace critical locus for vdW", "[vdW][crit]")
+{
+    // Argon + Xenon
+    std::valarray<double> Tc_K = { 150.687, 289.733 };
+    std::valarray<double> pc_Pa = { 4863000.0, 5842000.0 };
+    vdWEOS<double> vdW(Tc_K, pc_Pa);
+    auto Zc = 3.0/8.0;
+    auto rhoc0 = pc_Pa[0] / (vdW.R * Tc_K[0]) / Zc;
+    double T0 = Tc_K[0];
+    std::valarray<double> rhovec0 = { rhoc0, 0.0 };
+
+    auto tic0 = std::chrono::steady_clock::now();
+    std::string filename = "";
+    trace_critical_arclength_binary(vdW, T0, rhovec0, filename);
+    auto tic1 = std::chrono::steady_clock::now();
 }
\ No newline at end of file
diff --git a/src/trace_critical.cpp b/src/trace_critical.cpp
deleted file mode 100644
index dad53565a063d7bb96e29fbadb995f1ed1c027be..0000000000000000000000000000000000000000
--- a/src/trace_critical.cpp
+++ /dev/null
@@ -1,38 +0,0 @@
-#include <tuple>
-#include <valarray>
-#include <iostream>
-
-#include "teqp/core.hpp"
-
-void trace() { 
-    // Argon + Xenon
-    std::valarray<double> Tc_K = { 150.687, 289.733 };
-    std::valarray<double> pc_Pa = { 4863000.0, 5842000.0 };
-    vdWEOS<double> vdW(Tc_K, pc_Pa);
-    auto Zc = 3.0/8.0;
-    auto rhoc0 = pc_Pa[0]/(vdW.R*Tc_K[0]) / Zc; 
-    auto T = Tc_K[0];
-    const auto dT = 1;
-    std::valarray<double> rhovec = { rhoc0, 0.0 };
-    auto tic0 = std::chrono::steady_clock::now();
-    for (auto iter = 0; iter < 1000; ++iter){
-        auto drhovecdT = get_drhovec_dT_crit(vdW, T, rhovec);
-        rhovec += drhovecdT*dT;
-        T += dT;
-        int rr = 0;
-        auto z0 = rhovec[0] / rhovec.sum();
-        std::cout << z0 <<  "," << rhovec[0] << "," << T << "," << get_splus(vdW, T, rhovec) << std::endl;
-        if (z0 < 0) {
-            break;
-        }
-    }
-    auto tic1 = std::chrono::steady_clock::now();
-    std::cout << T << " " << rhovec[1] << std::endl;
-    std::cout << std::chrono::duration<double>(tic1 - tic0).count() << " s to trace" << std::endl;
-
-    int rr = 0;
-}
-int main() {   
-    trace();
-    return EXIT_SUCCESS;
-}
\ No newline at end of file