diff --git a/CMakeLists.txt b/CMakeLists.txt
index 82d74b1d..2403993c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -50,6 +50,10 @@ option (TEQP_MULTICOMPLEX_ENABLED
         "Enable the use of multi-complex arithmetic for taking derivatives"
         OFF)
         
+option (TEQP_COMPLEXSTEP_ENABLED
+        "Enable the use of complex-step derivatives for taking first derivatives"
+        ON)
+
 
 ####  SETUP
 set(CMAKE_CXX_STANDARD 17)
@@ -192,6 +196,7 @@ if (NOT TEQP_NO_TESTS)
 
   target_compile_definitions(catch_tests PRIVATE -DTEQPC_CATCH)
   target_compile_definitions(catch_tests PRIVATE -DTEQP_MULTICOMPLEX_ENABLED)
+  target_compile_definitions(catch_tests PRIVATE -DTEQP_COMPLEXSTEP_ENABLED)
   target_compile_definitions(catch_tests PRIVATE -DTEQP_MULTIPRECISION_ENABLED)
   target_link_libraries(catch_tests PUBLIC teqpcpp PRIVATE autodiff PRIVATE teqpinterface PRIVATE Catch2WithMain)
   add_test(normal_tests catch_tests)
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 878f00ce..512dd0bb 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -1,6 +1,8 @@
 #pragma once
 
+#if defined(TEQP_COMPLEXSTEP_ENABLED)
 #include <complex>
+#endif
 #include <map>
 #include <tuple>
 #include <numeric>
@@ -119,27 +121,29 @@ enum class ADBackends { autodiff
 #if defined(TEQP_MULTICOMPLEX_ENABLED)
     ,multicomplex
 #endif
+#if defined(TEQP_COMPLEXSTEP_ENABLED)
     ,complex_step
+#endif
 };
 
 template<typename Model, typename Scalar = double, typename VectorType = Eigen::ArrayXd>
 struct TDXDerivatives {
-
+    
     static auto get_Ar00(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
         return model.alphar(T, rho, molefrac);
     }
-
+    
     /**
-    * Calculate the derivative \f$\Lambda_{xy}\f$, where 
-    * \f[
-    * \Lambda^{\rm r}_{ij} = (1/T)^i\rho^j\left(\frac{\partial^{i+j}(\alpha^*)}{\partial(1/T)^i\partial\rho^j}\right)
-    * \f]
-    * 
-    * Note: none of the intermediate derivatives are returned, although they are calculated
-    */
+     * Calculate the derivative \f$\Lambda_{xy}\f$, where
+     * \f[
+     * \Lambda^{\rm r}_{ij} = (1/T)^i\rho^j\left(\frac{\partial^{i+j}(\alpha^*)}{\partial(1/T)^i\partial\rho^j}\right)
+     * \f]
+     *
+     * Note: none of the intermediate derivatives are returned, although they are calculated
+     */
     template<int iT, int iD, ADBackends be = ADBackends::autodiff, class AlphaWrapper>
     static auto get_Agenxy(const AlphaWrapper& w, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
-
+        
         static_assert(iT > 0 || iD > 0);
         if constexpr (iT == 0 && iD > 0) {
             if constexpr (be == ADBackends::autodiff) {
@@ -148,11 +152,13 @@ struct TDXDerivatives {
                 auto f = [&w, &T, &molefrac](const auto& rho__) { return w.alpha(T, rho__, molefrac); };
                 return powi(rho, iD)*derivatives(f, along(1), at(rho_))[iD];
             }
+#if defined(TEQP_COMPLEXSTEP_ENABLED)
             else if constexpr (iD == 1 && be == ADBackends::complex_step) {
                 double h = 1e-100;
                 auto rho_ = std::complex<Scalar>(rho, h);
                 return powi(rho, iD) * w.alpha(T, rho_, molefrac).imag() / h;
             }
+#endif
 #if defined(TEQP_MULTICOMPLEX_ENABLED)
             else if constexpr (be == ADBackends::multicomplex) {
                 using fcn_t = std::function<mcx::MultiComplex<Scalar>(const mcx::MultiComplex<Scalar>&)>;
@@ -173,11 +179,13 @@ struct TDXDerivatives {
                 auto f = [&w, &rho, &molefrac](const auto& Trecip__) {return w.alpha(forceeval(1.0/Trecip__), rho, molefrac); };
                 return powi(Trecip, iT)*derivatives(f, along(1), at(Trecipad))[iT];
             }
+#if defined(TEQP_COMPLEXSTEP_ENABLED)
             else if constexpr (iT == 1 && be == ADBackends::complex_step) {
                 double h = 1e-100;
                 auto Trecipcsd = std::complex<Scalar>(Trecip, h);
                 return powi(Trecip, iT)* w.alpha(1/Trecipcsd, rho, molefrac).imag()/h;
             }
+#endif
 #if defined(TEQP_MULTICOMPLEX_ENABLED)
             else if constexpr (be == ADBackends::multicomplex) {
                 using fcn_t = std::function<mcx::MultiComplex<Scalar>(const mcx::MultiComplex<Scalar>&)>;
@@ -1054,9 +1062,11 @@ struct IsochoricDerivatives{
             return build_Psir_gradient_multicomplex(model, T, rho);
         }
 #endif
+#if defined(TEQP_COMPLEXSTEP_ENABLED)
         else if constexpr (be == ADBackends::complex_step) {
             return build_Psir_gradient_complex_step(model, T, rho);
         }
+#endif
     }
 
     /***