diff --git a/FDTD/operator_cylindermultigrid.cpp b/FDTD/operator_cylindermultigrid.cpp
index 8edd06b..b8ff83d 100644
--- a/FDTD/operator_cylindermultigrid.cpp
+++ b/FDTD/operator_cylindermultigrid.cpp
@@ -18,6 +18,7 @@
 #include "operator_cylindermultigrid.h"
 #include "engine_cylindermultigrid.h"
 #include "operator_ext_cylinder.h"
+#include "tools/useful.h"
 
 Operator_CylinderMultiGrid::Operator_CylinderMultiGrid(vector<double> Split_Radii) : Operator_Cylinder()
 {
@@ -127,22 +128,20 @@ void Operator_CylinderMultiGrid::Init()
 
 void Operator_CylinderMultiGrid::CalcStartStopLines(unsigned int &numThreads, vector<unsigned int> &start, vector<unsigned int> &stop) const
 {
-	if ((numLines[0] - m_Split_Pos + 1)<numThreads) //in case more threads requested as lines in r-direction, reduce number of worker threads
-		numThreads = numLines[0] - m_Split_Pos + 1;
+	vector<unsigned int> jpt = AssignJobs2Threads(numLines[0]- m_Split_Pos + 1, numThreads, true);
 
-	unsigned int linesPerThread = round((float)(numLines[0] - m_Split_Pos + 1) / (float)numThreads);
-	if ((numThreads-1) * linesPerThread >= (numLines[0] - m_Split_Pos + 1))
-		--numThreads;
+	numThreads = jpt.size();
 
 	start.resize(numThreads);
 	stop.resize(numThreads);
 
-	for (unsigned int n=0; n<numThreads; n++)
+	start.at(0)= m_Split_Pos-1;
+	stop.at(0)= jpt.at(0)-1  + m_Split_Pos-1;
+
+	for (unsigned int n=1; n<numThreads; n++)
 	{
-		start.at(n) = n * linesPerThread + m_Split_Pos - 1;
-		stop.at(n) = (n+1) * linesPerThread - 1 + m_Split_Pos - 1;
-		if (n == numThreads-1) // last thread
-			stop.at(n) = numLines[0]-1;
+		start.at(n) = stop.at(n-1)+1;
+		stop.at(n) = start.at(n) + jpt.at(n) - 1;
 	}
 }
 
diff --git a/FDTD/operator_multithread.cpp b/FDTD/operator_multithread.cpp
index e0f1a16..f058ea6 100644
--- a/FDTD/operator_multithread.cpp
+++ b/FDTD/operator_multithread.cpp
@@ -17,6 +17,7 @@
 
 #include "operator_multithread.h"
 #include "engine_multithread.h"
+#include "tools/useful.h"
 
 Operator_Multithread* Operator_Multithread::New(unsigned int numThreads)
 {
@@ -76,22 +77,20 @@ void Operator_Multithread::Reset()
 
 void Operator_Multithread::CalcStartStopLines(unsigned int &numThreads, vector<unsigned int> &start, vector<unsigned int> &stop) const
 {
-	if (numLines[0]<numThreads) //in case more threads requested as lines in x-direction, reduce number of worker threads
-		numThreads = numLines[0];
+	vector<unsigned int> jpt = AssignJobs2Threads(numLines[0], numThreads, true);
 
-	unsigned int linesPerThread = round((float)numLines[0] / (float)numThreads);
-	if ((numThreads-1) * linesPerThread >= numLines[0])
-		--numThreads;
+	numThreads = jpt.size();
 
 	start.resize(numThreads);
 	stop.resize(numThreads);
 
-	for (unsigned int n=0; n<numThreads; n++)
+	start.at(0)=0;
+	stop.at(0)=jpt.at(0)-1;
+
+	for (unsigned int n=1; n<numThreads; n++)
 	{
-		start.at(n) = n * linesPerThread;
-		stop.at(n) = (n+1) * linesPerThread - 1;
-		if (n == numThreads-1) // last thread
-			stop.at(n) = numLines[0]-1;
+		start.at(n) = stop.at(n-1)+1;
+		stop.at(n) = start.at(n) + jpt.at(n) - 1;
 	}
 }
 
diff --git a/tools/useful.cpp b/tools/useful.cpp
index 8baac9b..8124876 100644
--- a/tools/useful.cpp
+++ b/tools/useful.cpp
@@ -20,6 +20,9 @@
 #include <cstdlib>
 #include <cmath>
 #include <climits>
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
 
 unsigned int CalcNyquistNum(double fmax, double dT)
 {
@@ -29,3 +32,34 @@ unsigned int CalcNyquistNum(double fmax, double dT)
 	return floor(T0/2/dT);
 }
 
+std::vector<unsigned int> AssignJobs2Threads(unsigned int jobs, unsigned int nrThreads, bool RemoveEmpty)
+{
+	std::vector<unsigned int> jpt; //jobs per thread
+
+	unsigned int ui_jpt = jobs/nrThreads;
+	for (unsigned int n=0;n<nrThreads;++n)
+	{
+		jpt.push_back(ui_jpt);
+		jobs-=ui_jpt;
+	}
+
+	for (unsigned int n=0;n<nrThreads;++n)
+	{
+		if (jobs>0)
+		{
+			++jpt.at(n);
+			--jobs;
+		}
+	}
+
+	if (jobs>0)
+		std::cerr << "AssignJobs2Threads: Error, " << jobs << " remain to be assigned, this should not have happend..." << std::endl;
+
+	if (RemoveEmpty)
+	{
+		while (jpt.back()==0)
+			jpt.pop_back();
+	}
+
+	return jpt;
+}
diff --git a/tools/useful.h b/tools/useful.h
index a5a36a3..c6160cf 100644
--- a/tools/useful.h
+++ b/tools/useful.h
@@ -18,7 +18,12 @@
 #ifndef USEFUL_H
 #define USEFUL_H
 
+#include <vector>
+
 //! Calc the nyquist number of timesteps for a given frequency and timestep
 unsigned int CalcNyquistNum(double fmax, double dT);
 
+//! Calculate an optimal job distribution to a given number of threads. Will return a vector with the jobs for each thread.
+std::vector<unsigned int> AssignJobs2Threads(unsigned int jobs, unsigned int nrThreads, bool RemoveEmpty=false);
+
 #endif // USEFUL_H