diff --git a/C++/nonLin/libNonLin.cc b/C++/nonLin/libNonLin.cc
index 5c5d63075c0d5df6926a12d2802cf49b23fc325f..3e93a2634ede2863661d2f4aaa0709494cf1f57f 100644
--- a/C++/nonLin/libNonLin.cc
+++ b/C++/nonLin/libNonLin.cc
@@ -79,19 +79,20 @@ void declare_nonLinCF_harmonicBalance(py::module &m){
 	.def(py::init(
 			[](shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, 
                         const vector<shared_ptr<CoefficientFunction>> ptrInput_i, 
+                        const vector<double> omega_i, 
                         vector<double> H_KL, 
                         vector<double> B_KL,
                         vector<double> ti,  
                         const std::string cFieldIn, const std::string cFieldOut, const shared_ptr<CoefficientFunction>& ptrMask, unsigned int order, 
                         bool usePreisBiro)
 	{
-		return new nonLinCF_harmonicBalance<T, DIM>(mesh, intrule, ptrInput_i, H_KL, B_KL, ti, cFieldIn, cFieldOut, ptrMask, order, usePreisBiro);
-	}), py::arg("mesh"), py::arg("IntegrationRule"), py::arg("InputFunctions"), py::arg("H_KL"), py::arg("B_KL"), py::arg("ti"), py::arg("field_in") = "H", py::arg("field_out") = "B",py::arg("mask")=nullptr, py::arg("order")=3, py::arg("usePreisBiro")=false)
+		return new nonLinCF_harmonicBalance<T, DIM>(mesh, intrule, ptrInput_i, omega_i, H_KL, B_KL, ti, cFieldIn, cFieldOut, ptrMask, order, usePreisBiro);
+	}), py::arg("mesh"), py::arg("IntegrationRule"), py::arg("InputFunctions"), py::arg("omegas"), py::arg("H_KL"), py::arg("B_KL"), py::arg("ti"), py::arg("field_in") = "H", py::arg("field_out") = "B",py::arg("mask")=nullptr, py::arg("order")=3, py::arg("usePreisBiro")=false)
     .def("Update",  overload_cast_<>()(&nonLinCF_harmonicBalance<T, DIM>::Update), "Updates the output CF according to the current input CF")
     .def("Update",  overload_cast_<const vector<shared_ptr<CoefficientFunction>> >()(&nonLinCF_harmonicBalance<T, DIM>::Update), "Updates the output CF according to the given input CF")
     .def("getParts", &nonLinCF_harmonicBalance<T, DIM>::getParts)
     .def("getTi", &nonLinCF_harmonicBalance<T, DIM>::getTi)
-    .def_property("omega0",  &nonLinCF_harmonicBalance<T, DIM>::getOmega, &nonLinCF_harmonicBalance<T, DIM>::setOmega)
+    .def_property("omega",  &nonLinCF_harmonicBalance<T, DIM>::getOmegas, &nonLinCF_harmonicBalance<T, DIM>::setOmegas)
 
     ;
 }
diff --git a/C++/nonLin/nonLinCF.cc b/C++/nonLin/nonLinCF.cc
index bf1cdee84f59821f31994a5b446a048534576bf2..e75e42c314a0262acc4c2fb96e27c5b630ab8e80 100644
--- a/C++/nonLin/nonLinCF.cc
+++ b/C++/nonLin/nonLinCF.cc
@@ -411,13 +411,14 @@ unsigned int nonLinCF<T, DIM>::getOrder() const{
 template <typename T, unsigned int DIM>
 nonLinCF_harmonicBalance<T, DIM>::nonLinCF_harmonicBalance(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, 
             const vector<shared_ptr<CoefficientFunction>> ptrInput_i, 
+            const vector<double> omega_i, 
             vector<double> H_KL, 
             vector<double> B_KL, 
             vector<double> ti, 
             const std::string cFieldIn, const std::string cFieldOut, 
             const shared_ptr<CoefficientFunction>& ptrMask, unsigned int order, bool usePreisBiro):
             nonLinCF<T, DIM>{mesh, intrule, nullptr, H_KL, B_KL, cFieldIn, cFieldOut, ptrMask, order, usePreisBiro}, 
-            ptrInput_i{ptrInput_i}, ti{ti}, omega_0{2.0*M_PI/(ti.back() - ti.front())}{
+            ptrInput_i{ptrInput_i}, omega_i{omega_i}, ti{ti}{
                 unsigned int i;
 
                 for(i = 0; i!= ti.size(); ++i){
@@ -464,9 +465,10 @@ void nonLinCF_harmonicBalance<T, DIM>::Update(const vector<shared_ptr<Coefficien
                 for(j = 0; j!= ptrInput_i.size(); ++j){
                     ptrInput_i[j]->Evaluate(this->mips[el][ip], fV_tmp);
 
-                    val_in_c += fV_tmp(0) * std::exp( j * 1i * omega_0 * ti[i]);
+                    val_in_c += fV_tmp(0) * std::exp( 1i * omega_i[j] * ti[i]);
                 }
 
+                // TODO: let the user decide.... 
                 val_in = std::abs(val_in_c.real());
 
                 if(!std::strcmp(this->sFieldIn.c_str(), "H")){       
diff --git a/C++/nonLin/nonLinCF.h b/C++/nonLin/nonLinCF.h
index 46ffcd0af2b3614c0ca6f6453dec0400dcba2fb1..8e6668cec5ebfedd5e37a313f7911732c2c1a230 100644
--- a/C++/nonLin/nonLinCF.h
+++ b/C++/nonLin/nonLinCF.h
@@ -144,6 +144,7 @@ class nonLinCF_harmonicBalance : public nonLinCF<T, DIM>{
     public:
         nonLinCF_harmonicBalance(shared_ptr<MeshAccess>& mesh, IntegrationRule& intrule, 
                     const vector<shared_ptr<CoefficientFunction>> ptrInput_i, 
+                    const vector<double> omega_i, 
                     vector<double> H_KL, 
                     vector<double> B_KL, 
                     vector<double> ti, 
@@ -155,12 +156,12 @@ class nonLinCF_harmonicBalance : public nonLinCF<T, DIM>{
         virtual void Update(const vector<shared_ptr<CoefficientFunction>> ptrInput);
 
         
-        double getOmega()const{
-            return omega_0;
+        vector<double> getOmegas()const{
+            return omega_i;
         }
         
-        void setOmega(const double val){
-            omega_0 = val;
+        void setOmegas(const vector<double> val){
+            omega_i = val;
         }
 
         vector<double> getTi() const{
@@ -174,8 +175,8 @@ class nonLinCF_harmonicBalance : public nonLinCF<T, DIM>{
 
     private:
         vector<shared_ptr<CoefficientFunction>> ptrInput_i;
+        vector<double> omega_i;
         vector<double> ti;
-        double omega_0 = 0;
         vector<std::shared_ptr<nonLinCF<T, DIM>>> nonLin_parts;
         
 };