diff --git a/stocproc/method_fft.py b/stocproc/method_fft.py
index 5ee3bdb..681d971 100644
--- a/stocproc/method_fft.py
+++ b/stocproc/method_fft.py
@@ -43,14 +43,43 @@ def find_integral_boundary(integrand, tol, ref_val, max_val, x0):
     """
     _max_num_iteration = 100
     _i = 0
+    x0 = float(x0)
     I_ref = integrand(ref_val)
     if I_ref < tol:
-        pass
+        # return find_integral_boundary(integrand = lambda x: 1/integrand(x),
+        #                               tol = 1/tol,
+        #                               ref_val=ref_val,
+        #                               max_val=max_val,
+        #                               x0=-x0)
+        x_old = ref_val
+        while True:
+            _i += 1
+            if _i > _max_num_iteration:
+                raise RuntimeError("max number of iteration reached")
+            x = ref_val - x0
+            try:
+                I_x = integrand(x)
+            except Exception as e:
+                raise RuntimeError("evaluation of integrand failed due to {}".format(e))
+
+            if I_x > tol:
+                break
+            x0 *= 2
+            x_old = x
+        a = brentq(lambda x: integrand(x) - tol, x_old, x)
+        return a
     elif I_ref > tol:
         x_old = ref_val
         while True:
+            _i += 1
+            if _i > _max_num_iteration:
+                raise RuntimeError("max number of iteration reached")
             x = ref_val + x0
-            I_x = integrand(x)
+            try:
+                I_x = integrand(x)
+            except Exception as e:
+                raise RuntimeError("evaluation of integrand failed due to {}".format(e))
+
             if I_x < tol:
                 break
             x0 *= 2
@@ -148,28 +177,33 @@ def _f_opt(x, integrand, a, b, N, t_max, ft_ref, diff_method, _f_opt_cache, b_on
         else:
             a_ = find_integral_boundary(integrand, tol=tol, ref_val=a, max_val=1e6, x0=-1)
             b_ = find_integral_boundary(integrand, tol=tol, ref_val=b, max_val=1e6, x0=1)
-    except ValueError:
-        a_ = -1
-        b_ = 1
+    except RuntimeError:
+        # in case 'find_integral_boundary' failes
+        d = 300
+        _f_opt_cache[key] = d, None, None
+        return d
+
 
     tau, ft_tau = fourier_integral_midpoint(integrand, a_, b_, N)
     idx = np.where(tau <= t_max)
     ft_ref_tau = ft_ref(tau[idx])
     d = diff_method(ft_ref_tau, ft_tau[idx])
     _f_opt_cache[key] = d, a_, b_
+    #log.debug("f_opt tol {} -> d {}".format(10**x, d))
     return np.log10(d)
 
 def _lower_contrs(x, integrand, a, b, N, t_max, ft_ref, diff_method, _f_opt_cache, b_only):
     _f_opt(x, integrand, a, b, N, t_max, ft_ref, diff_method, _f_opt_cache, b_only)
-    tol = 10**x
     d, a_, b_ = _f_opt_cache[float(x[0])]
+    if (a_ is None) or (b_ is None):
+        return -1
     v = N * np.pi / (b_ - a_) - t_max
-    log.debug("lower constr value {} for x {} (tol {})".format(v, x, tol))
+    #log.debug("lower constr value {} for x {} (tol {})".format(v, 10**x, tol))
     return v
 
 
 def _upper_contrs(x):
-    log.debug("upper constr value {}".format(-x))
+    #log.debug("upper constr value {}".format(-x))
     return -x
 
 
@@ -178,7 +212,8 @@ def opt_integral_boundaries(integrand, a, b, t_max, ft_ref, opt_b_only, N, diff_
 
     _f_opt_cache = dict()
     args = (integrand, a, b, N, t_max, ft_ref, diff_method, _f_opt_cache, opt_b_only)
-    r = minimize(_f_opt, x0 = [-0.1], args = args,
+    x0 = np.log10(0.1*integrand(b))
+    r = minimize(_f_opt, x0 = x0, args = args,
                  method='SLSQP',
                  constraints=[{"type": "ineq", "fun": _lower_contrs, "args": args},
                               {"type": "ineq", "fun": _upper_contrs}])
@@ -219,7 +254,7 @@ def get_N_a_b_for_accurate_fourier_integral(integrand, a, b, t_max, tol, ft_ref,
             raise RuntimeError("maximum number of points for Fourier Transform reached")
         i += 1    
 
-def get_dt_for_accurate_interpolation(t_max, tol, ft_ref):
+def get_dt_for_accurate_interpolation(t_max, tol, ft_ref, diff_method=_absDiff):
     N = 32
     sub_sampl = 8
     
@@ -231,43 +266,40 @@ def get_dt_for_accurate_interpolation(t_max, tol, ft_ref):
         ft_intp = ComplexInterpolatedUnivariateSpline(x = tau_sub, y = ft_ref_n[::sub_sampl], k=3)
         ft_intp_n = ft_intp(tau)
         
-        d = np.max(np.abs(ft_intp_n-ft_ref_n))
+        d = diff_method(ft_intp_n, ft_ref_n)
+        log.debug("acc interp N {} dt {:.3e} {:.3e} -> d {:.3e}".format(N, sub_sampl*tau[1], t_max/(N/sub_sampl), d))
         if d < tol:
             return t_max/(N/sub_sampl)
         N*=2
 
 
-def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, a, b, ft_ref, opt_b_only, N_max = 2**20):
+def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, a, b, ft_ref, opt_b_only, N_max = 2**20, diff_method=_absDiff):
     N, a, b = get_N_a_b_for_accurate_fourier_integral(integrand, a, b,
                                                       t_max  = t_max,
                                                       tol    = intgr_tol,
                                                       ft_ref = ft_ref,
                                                       opt_b_only=opt_b_only,
-                                                      N_max  = N_max)
+                                                      N_max  = N_max,
+                                                      diff_method=diff_method)
     dt_tol = get_dt_for_accurate_interpolation(t_max  = t_max,
                                                tol    = intpl_tol,
-                                               ft_ref = ft_ref)
-    
+                                               ft_ref = ft_ref,
+                                               diff_method=diff_method)
     dx = (b-a)/N
     dt = 2*np.pi/dx/N
     if dt <= dt_tol:
-        log.info("dt criterion fulfilled")
+        log.debug("dt criterion fulfilled")
         return a, b, N, dx, dt
     else:
-        log.info("down scale dx and dt to match new power of 2 N")
+        log.debug("down scale dx and dt to match new power of 2 N")
 
     N_min = 2*np.pi/dx/dt_tol
     N = 2**int(np.ceil(np.log2(N_min)))
-
-
-    #scale = np.sqrt(N_min/N)
-    #assert scale <= 1
-    scale = 1
-
+    scale = np.sqrt(N_min/N)
     dx_new = scale*dx
     b_minus_a = dx_new*N
-    
     dt_new = 2*np.pi/dx_new/N
+    assert dt_new < dt_tol
     if opt_b_only:
         b = a + b_minus_a
     else:
@@ -275,13 +307,4 @@ def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, a, b, ft_ref, opt_b_
         b += delta/2
         a -= delta/2
 
-    rd, a, b = opt_integral_boundaries(integrand=integrand, a=a, b=b, t_max=t_max, ft_ref=ft_ref,
-                                       opt_b_only=opt_b_only, N=N)
-
-    log.debug("rd after final optimization:{:.3e}".format(rd))
-
     return a, b, N, dx_new, dt_new
-        
-        
-    
-    
\ No newline at end of file
diff --git a/stocproc/method_kle.py b/stocproc/method_kle.py
index 28dfcb1..cea56fc 100644
--- a/stocproc/method_kle.py
+++ b/stocproc/method_kle.py
@@ -365,7 +365,7 @@ def auto_ng(corr, t_max, ngfac=2, meth=get_mid_point_weights_times, tol=1e-3, di
            the interpolation error is maximal when beeing in between the reference points.
         5) Now calculate the deviation :math:`\Delta(n)` for sequential n starting at n=0. Stop if
            :math:`\Delta(n) < tol`. If the deviation does not drop below tol for all :math:`0 \leq n < ng-1` increase
-           ng as follows :math:`ng = 2*ng-1` and start over at 1). (This update schema for ng asured that ng is odd
+           ng as follows :math:`ng = 2*ng-1` and start over at 1). (This update scheme for ng asured that ng is odd
            which is needed for the 'simpson' and 'fourpoint' integration weights)
 
     .. note::
diff --git a/stocproc/stocproc.py b/stocproc/stocproc.py
index 6e5781c..700c415 100644
--- a/stocproc/stocproc.py
+++ b/stocproc/stocproc.py
@@ -268,16 +268,16 @@ class StocProc_KLE(_absStocProc):
 class StocProc_FFT(_absStocProc):
     r"""Simulate Stochastic Process using FFT method
 
-    This method works only for auto correlations functions of the form
+    This method uses the relation of the auto correlation to the spectral density
 
-    .. math:: \alpha(\tau) = \int_{\omega_\mathrm{min}}^{\omega_\mathrm{max}} \mathrm{d}\omega \, \frac{J(\omega)}{\pi} e^{-\mathrm{i}\omega \tau}
+    .. math:: \alpha(\tau) = \int \mathrm{d}\omega \, \frac{J(\omega)}{\pi} e^{-\mathrm{i}\omega \tau}
 
     where :math:`J(\omega)` is a real non negative function, usually called spectral density.
-    Then the integral can be approximated by a discrete integration schema:
+    Then the integral can be approximated by a discrete integration scheme:
 
     .. math:: \alpha(\tau) \approx \sum_{k=0}^{N-1} w_k \frac{J(\omega_k)}{\pi} e^{-\mathrm{i} k \omega_k \tau}
 
-    where the :math:`\omega_k` depend on the integration schema.
+    where the :math:`\omega_k` depend on the integration scheme.
 
     For a process defined by
 
@@ -294,14 +294,14 @@ class StocProc_FFT(_absStocProc):
                                            \approx & \alpha(t-s)
         \end{align}
 
-    To calculate :math:`Z(t)` the schema of the Discrete Fourier Transform (DFT) can be applied as follows:
+    To calculate :math:`Z(t)` the Discrete Fourier Transform (DFT) can be applied as follows:
 
     .. math:: Z(t_l) = e^{-\mathrm{i}\omega_\mathrm{min} t_l} \sum_{k=0}^{N-1} \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k  e^{-\mathrm{i} 2 \pi \frac{k l}{N} \frac{\Delta \omega \Delta t}{ 2 \pi} N}
 
     Here :math:`\omega_k` has to take the form :math:`\omega_k = \omega_\mathrm{min} + k \Delta \omega` and
     :math:`\Delta \omega = (\omega_\mathrm{max} - \omega_\mathrm{min}) / N-1` which limits
-    the itegration schemas to those with equidistant weights.
-    For the DFT schema to be applicable :math:`\Delta t` has to be chosen such that
+    the itegration schemes to those with equidistant weights.
+    For the DFT scheme to be applicable :math:`\Delta t` has to be chosen such that
 
     .. math:: 1 = \frac{\Delta \omega \Delta t}{2 \pi} N
 
diff --git a/tests/test_method_fft.py b/tests/test_method_fft.py
index 91d845f..3624ac1 100644
--- a/tests/test_method_fft.py
+++ b/tests/test_method_fft.py
@@ -1,22 +1,22 @@
+import os
+import pathlib
+import sys
+p = pathlib.PosixPath(os.path.abspath(__file__))
+sys.path.insert(0, str(p.parent.parent))
+
 import math
 import logging
 import numpy as np
-import os
-import pathlib
 import scipy.integrate as sp_int
 from scipy.special import gamma as gamma_func
 import stocproc as sp
-import sys
+from stocproc import method_fft
 
 try:
     import matplotlib.pyplot as plt
 except ImportError:
     print("matplotlib not found -> any plotting will crash")
 
-p = pathlib.PosixPath(os.path.abspath(__file__))
-sys.path.insert(0, str(p.parent.parent))
-
-
 def test_find_integral_boundary():
     def f(x):
         return np.exp(-(x)**2)
@@ -41,14 +41,14 @@ def test_find_integral_boundary():
     b = sp.method_fft.find_integral_boundary(integrand=f2, tol=tol, ref_val=1, x0=+1, max_val=1e6)
     a = sp.method_fft.find_integral_boundary(integrand=f2, tol=tol, ref_val=-1, x0=-1, max_val=1e6)
     assert a != b
-    assert abs(f2(a)/f2(1) -tol) < 1e-14
-    assert abs(f2(b)/f2(-1)-tol) < 1e-14
+    assert abs(f2(a) -tol) < 1e-14
+    assert abs(f2(b)-tol) < 1e-14
     
     b = sp.method_fft.find_integral_boundary(integrand=f2, tol=tol, ref_val=1, x0=b+5, max_val=1e6)
     a = sp.method_fft.find_integral_boundary(integrand=f2, tol=tol, ref_val=-1, x0=a-5, max_val=1e6)
     assert a != b
-    assert abs(f2(a)/f2(1)-tol) < 1e-14, "diff {}".format(abs(f2(a)/f2(1)-tol))
-    assert abs(f2(b)/f2(-1)-tol) < 1e-14, "diff {}".format(abs(f2(b)/f2(-1)-tol))
+    assert abs(f2(a)-tol) < 1e-14, "diff {}".format(abs(f2(a)/f2(1)-tol))
+    assert abs(f2(b)-tol) < 1e-14, "diff {}".format(abs(f2(b)/f2(-1)-tol))
     
     def f3(x):
         return np.exp(-(x-5)**2)*x**2
@@ -57,16 +57,30 @@ def test_find_integral_boundary():
     b = sp.method_fft.find_integral_boundary(integrand=f3, tol=tol, ref_val=5, x0=+1, max_val=1e6)
     a = sp.method_fft.find_integral_boundary(integrand=f3, tol=tol, ref_val=5, x0=-1, max_val=1e6)
     assert a != b
-    f3_refval = f3(5)
-    assert abs(f3(a)/f3_refval-tol) < 1e-14
-    assert abs(f3(b)/f3_refval-tol) < 1e-14
+    assert abs(f3(a)-tol) < 1e-14
+    assert abs(f3(b)-tol) < 1e-14
     
     b = sp.method_fft.find_integral_boundary(integrand=f3, tol=tol, ref_val=5, x0=b+5, max_val=1e6)
     a = sp.method_fft.find_integral_boundary(integrand=f3, tol=tol, ref_val=5, x0=a-5, max_val=1e6)
     assert a != b
-    assert abs(f3(a)/f3_refval-tol) < 1e-14, "diff {}".format(abs(f3(a)-tol))
-    assert abs(f3(b)/f3_refval-tol) < 1e-14, "diff {}".format(abs(f3(b)-tol))
-    
+    assert abs(f3(a)-tol) < 1e-14, "diff {}".format(abs(f3(a)-tol))
+    assert abs(f3(b)-tol) < 1e-14, "diff {}".format(abs(f3(b)-tol))
+
+
+    ##################################
+    ## the case where f(xref) < tol ##
+    ##################################
+
+    def f(x):
+        return np.exp(-x ** 2)
+
+    tol = 1e-3
+    b = sp.method_fft.find_integral_boundary(integrand=f, tol=tol, ref_val=10, x0=+1, max_val=1e6)
+    assert abs(f(b) - tol) < 1e-14
+    a = sp.method_fft.find_integral_boundary(integrand=f, tol=tol, ref_val=-10, x0=-1., max_val=1e6)
+    assert abs(f(a) - tol) < 1e-14
+
+
 def fourier_integral_trapz(integrand, a, b, N):
     """
         approximates int_a^b dx integrand(x) by the riemann sum with N terms
@@ -324,44 +338,44 @@ def test_calc_abN():
     
     tol = 1e-3
     tmax=40
+    diff_method = method_fft._absDiff
 
-    a,b = sp.method_fft.find_integral_boundary_auto(integrand=intg, tol=1e-2, ref_val=1)
-    a, b, N, dx, dt = sp.method_fft.calc_ab_N_dx_dt(integrand = intg,
-                                                  intgr_tol = tol, 
-                                                  intpl_tol = tol,
-                                                  t_max = tmax,
-                                                  a    = 0, 
-                                                  b    = b, 
-                                                  ft_ref = bcf_ref,
-                                                  opt_b_only = True)
+    a, b = sp.method_fft.find_integral_boundary_auto(integrand=intg, tol=1e-2, ref_val=1)
+    a, b, N, dx, dt = sp.method_fft.calc_ab_N_dx_dt(integrand=intg,
+                                                    intgr_tol=tol,
+                                                    intpl_tol=tol,
+                                                    t_max=tmax,
+                                                    a=0,
+                                                    b=b,
+                                                    ft_ref=bcf_ref,
+                                                    opt_b_only=True,
+                                                    diff_method=diff_method)
 
     tau, ft_tau = sp.method_fft.fourier_integral_midpoint(intg, a, b, N)
-    print(dt, tau[1])
     idx = np.where(tau <= tmax)
     ft_ref_tau = bcf_ref(tau[idx])
-    rd = np.max(np.abs(ft_tau[idx] - ft_ref_tau) / np.abs(ft_ref_tau))
-    print("rd {:.3e} <= tol {:.3e}".format(rd, tol))
-    assert rd < tol
+    rd = diff_method(ft_tau[idx], ft_ref_tau)
 
     tau_fine = np.linspace(0, tmax, 1500)
     ft_ref_n = bcf_ref(tau_fine)
     ft_intp = sp.tools.ComplexInterpolatedUnivariateSpline(x=tau[idx], y=ft_tau[idx], k=3)
     ft_intp_n = ft_intp(tau_fine)
-    d = np.max(np.abs(ft_intp_n - ft_ref_n))
-    print("d {:.3e} <= tol {:.3e}".format(d, tol))
-    assert d < tol
+    d = diff_method(ft_intp_n, ft_ref_n)
 
+    assert rd < tol
+    assert d < tol
     assert (np.abs(dx*dt*N - np.pi*2)) < 1e-15
 
     
     
 if __name__ == "__main__":
     logging.basicConfig(level=logging.INFO)
+    #logging.basicConfig(level=logging.DEBUG)
     test_find_integral_boundary()
-    # test_fourier_integral_finite_boundary()
-    # test_fourier_integral_infinite_boundary(plot=True)
-    # test_get_N_a_b_for_accurate_fourier_integral()
-    # test_get_N_a_b_for_accurate_fourier_integral_b_only()
-    # test_get_dt_for_accurate_interpolation()
-    # test_sclicing()
-    # test_calc_abN()
\ No newline at end of file
+    test_fourier_integral_finite_boundary()
+    test_fourier_integral_infinite_boundary(plot=False)
+    test_get_N_a_b_for_accurate_fourier_integral()
+    test_get_N_a_b_for_accurate_fourier_integral_b_only()
+    test_get_dt_for_accurate_interpolation()
+    test_sclicing()
+    test_calc_abN()
\ No newline at end of file