diff --git a/setup.py b/setup.py
index acf5e4fe..9a0d3489 100644
--- a/setup.py
+++ b/setup.py
@@ -68,6 +68,7 @@
             "dcore_bse = dcore.dcore_bse:run",
             "dcore_gk = dcore.dcore_gk:run",
             "dcore_mpicheck = dcore.dcore_mpicheck:run",
+            "dcore_anacont = dcore.dcore_anacont:run",
             "dcore_anacont_pade = dcore.dcore_anacont_pade:run",
             "dcore_anacont_spm = dcore.dcore_anacont_spm:run",
             "dcore_anacont_spm_interactive = dcore.dcore_anacont_spm_interactive:run",
diff --git a/src/dcore/anacont/__init__.py b/src/dcore/anacont/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/src/dcore/anacont/pade.py b/src/dcore/anacont/pade.py
new file mode 100644
index 00000000..7f30c99c
--- /dev/null
+++ b/src/dcore/anacont/pade.py
@@ -0,0 +1,62 @@
+#
+# DCore -- Integrated DMFT software for correlated electrons
+# Copyright (C) 2017 The University of Tokyo
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+#
+
+import numpy
+
+from dcore._dispatcher import MeshImFreq, GfReFreq, GfImFreq
+from dcore.program_options import create_parser, parse_parameters
+
+
+def _set_n_pade(omega_cutoff, beta, n_min, n_max):
+    """
+    Return (int)n_pade: the number of Matsubara frequencies below the cutoff frequency.
+    n_pade is bounded between n_min and n_max
+    """
+    n_pade = int((beta * omega_cutoff + numpy.pi) / (2.0 * numpy.pi))
+    print("n_pade = {} (evaluated from omega_pade)".format(n_pade))
+    n_pade = max(n_pade, n_min)
+    n_pade = min(n_pade, n_max)
+    print("n_pade = {}".format(n_pade))
+    return n_pade
+
+
+def anacont(sigma_iw_npz, beta, mesh_w, params_pade):
+    iw_max = params_pade["iomega_max"]
+    n_min = params_pade["n_min"]
+    n_max = params_pade["n_max"]
+    eta = params_pade["eta"]
+    n_pade = _set_n_pade(iw_max, beta, n_min=n_min, n_max=n_max)
+
+    data_w = {}
+    num_data = numpy.sum([key.startswith("data") for key in sigma_iw_npz.keys()])
+    for idata in range(num_data):
+        key = f"data{idata}"
+        data = sigma_iw_npz[key]
+        mesh_iw = MeshImFreq(beta, "Fermion", data.shape[0] // 2)
+        sigma_iw = GfImFreq(data=data, beta=beta, mesh=mesh_iw)
+        sigma_w = GfReFreq(mesh=mesh_w, target_shape=data.shape[1:])
+        sigma_w.set_from_pade(sigma_iw, n_points=n_pade, freq_offset=eta)
+        data_w[key] = sigma_w.data
+    return data_w
+
+
+def parameters_from_ini(inifile):
+    parser = create_parser(["post.anacont.pade"])
+    parser.read(inifile)
+    params = parser.as_dict()
+    return params["post.anacont.pade"]
diff --git a/src/dcore/anacont/spm.py b/src/dcore/anacont/spm.py
new file mode 100644
index 00000000..775f9492
--- /dev/null
+++ b/src/dcore/anacont/spm.py
@@ -0,0 +1,351 @@
+#
+# DCore -- Integrated DMFT software for correlated electrons
+# Copyright (C) 2017 The University of Tokyo
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/li/show_fitcenses/>.
+#
+
+import sys
+
+import numpy as np
+import cvxpy as cp
+from scipy.interpolate import interp1d
+from scipy.sparse.linalg import svds
+
+from dcore._dispatcher import MeshImFreq, GfImFreq, GfReFreq
+from dcore.program_options import create_parser
+
+def set_default_values(dictionary, default_values_dict):
+    for key, value in default_values_dict.items():
+        if key not in dictionary:
+            dictionary[key] = value
+    return dictionary
+
+
+def _find_sum_rule_const(matsubara_frequencies, gf_wn, ntail, show_fit=False):
+    wntail = matsubara_frequencies[-ntail:]
+    gftail = gf_wn[-ntail:]
+    const_imagpart = -np.mean(np.imag(gftail) * wntail)
+    const_realpart = np.mean(np.real(gftail))
+    if show_fit:
+        import matplotlib.pyplot as plt
+
+        z = np.linspace(wntail[0], wntail[-1], num=1000)
+        plt.scatter(wntail, np.imag(gftail), zorder=5, color="C0", label="data")
+        plt.plot(z, -const_imagpart / z, zorder=10, color="C1", label="fit")
+        plt.xlabel(r"$\omega_n$")
+        plt.ylabel(r"Im $G( \omega_n )$")
+        plt.legend()
+        plt.tight_layout()
+        plt.show()
+    return const_realpart, const_imagpart
+
+
+def _calc_gf_tau(
+    matsubara_frequencies, gf_wn, beta, const_real_tail, const_imag_tail, ntau
+):
+    tau_grid = np.linspace(0, beta, num=ntau)
+    gf_tau = -0.5 * const_imag_tail * np.ones(tau_grid.shape, dtype=np.float64)
+    kernel = lambda tau: np.sum(
+        (np.real(gf_wn) - const_real_tail) * np.cos(matsubara_frequencies * tau)
+        + (np.imag(gf_wn) + const_imag_tail / matsubara_frequencies)
+        * np.sin(matsubara_frequencies * tau)
+    )
+    gf_tau += 2 / beta * np.vectorize(kernel)(tau_grid)
+    gf_tau *= -1
+    return tau_grid, gf_tau
+
+
+def calc_gf_tau_from_gf_matsubara(
+    matsubara_frequencies, gf_wn, ntau, ntail, beta, show_fit=False
+):
+    const_real_tail, const_imag_tail = _find_sum_rule_const(
+        matsubara_frequencies, gf_wn, ntail, show_fit=show_fit
+    )
+    print("Determined tail constants: {}, {}".format(const_real_tail, const_imag_tail))
+    tau_grid, gf_tau = _calc_gf_tau(
+        matsubara_frequencies, gf_wn, beta, const_real_tail, const_imag_tail, ntau
+    )
+    return tau_grid, gf_tau, const_real_tail, const_imag_tail
+
+
+def get_single_continuation(
+    tau_grid,
+    gf_tau,
+    nsv,
+    beta,
+    emin,
+    emax,
+    num_energies,
+    sum_rule_const,
+    lambd,
+    verbose=True,
+    max_iters=100,
+    solver="ECOS",
+):
+    U, S, Vt, delta_energy, energies_extract = _get_svd_for_continuation(
+        tau_grid, nsv, beta, emin, emax, num_energies
+    )
+    rho_prime, gf_tau_fit, chi2 = _solveProblem(
+        delta_energy,
+        U,
+        S,
+        Vt,
+        gf_tau,
+        sum_rule_const,
+        lambd,
+        verbose=verbose,
+        max_iters=max_iters,
+        solver=solver,
+    )
+    rho = np.dot(Vt.T, rho_prime)
+    rho_integrated = np.trapz(y=rho, x=energies_extract)
+    return rho, gf_tau_fit, energies_extract, rho_integrated, chi2
+
+
+def get_multiple_continuations(
+    tau_grid,
+    gf_tau,
+    nsv,
+    beta,
+    emin,
+    emax,
+    num_energies,
+    sum_rule_const,
+    lambdas,
+    verbose=True,
+    max_iters=100,
+    solver="ECOS",
+):
+    U, S, Vt, delta_energy, energies_extract = _get_svd_for_continuation(
+        tau_grid, nsv, beta, emin, emax, num_energies
+    )
+    continued_densities = []
+    chi2_values = []
+    rho_values = []
+    for lambd in lambdas:
+        rho_prime, _, chi2 = _solveProblem(
+            delta_energy,
+            U,
+            S,
+            Vt,
+            gf_tau,
+            sum_rule_const,
+            lambd,
+            verbose=verbose,
+            max_iters=max_iters,
+            solver=solver,
+        )
+        rho = np.dot(Vt.T, rho_prime)
+        rho_integrated = np.trapz(y=rho, x=energies_extract)
+        continued_densities.append(rho)
+        chi2_values.append(chi2)
+        rho_values.append(rho_integrated)
+    return energies_extract, continued_densities, chi2_values, rho_values
+
+
+def _get_kernel_matrix(
+    energies, tau_grid, beta, delta_energy, z_critical=20
+):  # with z_critical=20 error of asymptote is below 1e-10
+    assert tau_grid[0] == 0
+    assert tau_grid[-1] == beta
+    i_lc = np.searchsorted(energies, -z_critical / beta)
+    i_rc = np.searchsorted(energies, +z_critical / beta)
+    kernel = np.zeros((len(tau_grid), len(energies)))
+    kernel[:, i_lc : i_rc + 1] = (
+        0.5
+        * np.exp(np.multiply.outer(0.5 * beta - tau_grid, energies[i_lc : i_rc + 1]))
+        / np.cosh(0.5 * beta * energies[i_lc : i_rc + 1])
+    )
+    if i_lc > 0:  # asymptote for kernel for w -> -inf
+        kernel[:, 0:i_lc] = np.exp(np.multiply.outer(beta - tau_grid, energies[0:i_lc]))
+    if i_rc < len(energies) - 1:  # asymptote for kernel for w -> +inf
+        kernel[:, i_rc + 1 : len(energies)] = np.exp(
+            np.multiply.outer(-tau_grid, energies[i_rc + 1 : len(energies)])
+        )
+    kernel[:, 0] *= 0.5
+    kernel[:, -1] *= 0.5
+    kernel *= delta_energy
+    return kernel
+
+
+def _getSVD(matrix, nsv=None):
+    if nsv == None or nsv >= min(matrix.shape) or nsv <= 0:
+        raise Exception(
+            "nsv must be 0 < nsv < {} (min(matrix.shape))".format(min(matrix.shape))
+        )
+    U, S, V = svds(matrix, k=nsv, which="LM")
+    S = np.flip(S, axis=0)
+    U = np.flip(U, axis=1)
+    Vt = np.flip(V, axis=0)
+    return U, S, Vt
+
+
+def _solveProblem(
+    delta_energy,
+    U,
+    S,
+    Vt,
+    gf_tau,
+    sum_rule_const,
+    lambd,
+    verbose=True,
+    max_iters=100,
+    solver="ECOS",
+):
+    Q = len(S)
+    Smat = np.diag(S)
+    rho_prime = cp.Variable(Q)
+    errTerm = gf_tau - U @ Smat @ rho_prime
+    objective = cp.Minimize(0.5 * cp.norm2(errTerm) ** 2 + lambd * cp.norm1(rho_prime))
+    V_mod = np.copy(Vt.T)
+    V_mod[0, :] *= 0.5
+    V_mod[-1, :] *= 0.5
+    constraints = [
+        Vt.T @ rho_prime >= 0,
+        cp.sum(delta_energy * V_mod @ rho_prime) == sum_rule_const,
+    ]  # uniform real energy grid is assumed here
+    prob = cp.Problem(objective, constraints)
+    _ = prob.solve(verbose=verbose, solver=solver, max_iters=max_iters)
+    gf_tau_fit = np.dot(U, np.dot(Smat, rho_prime.value))
+    chi2 = 0.5 * np.linalg.norm(gf_tau - gf_tau_fit, ord=2) ** 2
+    return rho_prime.value, gf_tau_fit, chi2
+
+
+def _get_svd_for_continuation(tau_grid, nsv, beta, emin, emax, num_energies):
+    energies_extract = np.linspace(emin, emax, num=num_energies)
+    delta_energy = energies_extract[1] - energies_extract[0]
+    kernel = _get_kernel_matrix(energies_extract, tau_grid, beta, delta_energy)
+    U, S, Vt = _getSVD(kernel, nsv=nsv)
+    return U, S, Vt, delta_energy, energies_extract
+
+
+def _integral_kramers_kronig(
+    energies_imagpart, energy_realpart, gf_imag_interp, energy_threshold
+):
+    enum1 = gf_imag_interp(energies_imagpart) - gf_imag_interp(energy_realpart)
+    enum2 = np.gradient(gf_imag_interp(energies_imagpart), energies_imagpart)
+    den = energies_imagpart - energy_realpart
+    mask = np.where(np.abs(den) < energy_threshold, 1, 0)  # 1 if den is below threshold
+    # avoid divide by zero with mask * energy_threshold, in this case (1 - mask == 0)
+    term1 = enum1 / (den + mask * energy_threshold) * (1 - mask)
+    term2 = enum2 * mask
+    kernel = term1 + term2
+    integral = np.trapz(y=kernel, x=energies_imagpart)
+    return integral
+
+
+def get_kramers_kronig_realpart(
+    energies, gf_imag, energy_threshold=1e-10, dos_threshold=1e-5
+):
+    if np.abs(gf_imag[0]) > dos_threshold:
+        print("Warning! DOS at left interval end exceeds {}.".format(dos_threshold))
+    if np.abs(gf_imag[-1]) > dos_threshold:
+        print("Warning! DOS at right interval end exceeds {}.".format(dos_threshold))
+    gf_real = np.zeros(energies.shape, dtype=np.float64)
+    gf_imag_interp = interp1d(
+        energies, gf_imag, kind="linear", bounds_error=False, fill_value=0.0
+    )
+    energies_noend = energies[1:-1]
+    a, b = energies[0], energies[-1]
+    gf_real[1:-1] = (
+        gf_imag_interp(energies_noend)
+        / np.pi
+        * np.log((b - energies_noend) / (energies_noend - a))
+    )  # assume zero DOS at endpoints
+    integral_func = lambda y: _integral_kramers_kronig(
+        energies, y, gf_imag_interp, energy_threshold
+    )
+    gf_real += np.vectorize(integral_func)(energies) / np.pi
+    return energies, gf_real, gf_imag
+
+
+def dos_to_gf_imag(dos):
+    return -np.pi * dos
+
+
+def check_parameters(params):
+    n_tau = params["n_tau"]
+    if n_tau <= 0:
+        sys.exit("n_tau must be positive.")
+    pass
+
+def parameters_from_ini(inifile):
+    parser = create_parser(["post.anacont.spm"])
+    parser.read(inifile)
+    params = parser.as_dict()
+    return params["post.anacont.spm"]
+
+def anacont(sigma_iw_npz, beta, mesh_w, params_spm):
+
+    n_tau = params_spm["n_tau"]
+    n_tail = params_spm["n_tail"]
+    show_fit = params_spm["show_fit"]
+    
+    n_sv = params_spm["n_sv"]
+    L1_coeff = params_spm["lambda"]
+    max_iters = params_spm["max_iters_opt"]
+    solver_opt = params_spm["solver_opt"]
+    verbose_opt = params_spm["verbose_opt"]
+    
+    data_w = {}
+    num_data = np.sum([key.startswith("data") for key in sigma_iw_npz.keys()])
+    for idata in range(num_data):
+        key = f"data{idata}"
+        data = sigma_iw_npz[key]
+        n_matsubara = data.shape[0] // 2
+        n_matsubara_retain = min(n_matsubara, params_spm["n_matsubara"])
+        mesh_iw = MeshImFreq(beta, "Fermion", n_matsubara_retain)
+        gf_iw = GfImFreq(data=data, beta=beta, mesh=mesh_iw)
+        n_orbitals = gf_iw.data.shape[1]
+        matsubara_frequencies = np.imag(gf_iw.mesh.values()[n_matsubara_retain:])
+        sigma_w_data = np.zeros(
+            (mesh_w.size, n_orbitals, n_orbitals), dtype=np.complex128
+        )
+        for i_orb in range(n_orbitals):
+            print(
+                f"Performing analytic continuation for data index {idata} and orbital index {i_orb}..."
+            )
+            gf_imag_matsubara = gf_iw.data[
+                n_matsubara : n_matsubara + n_matsubara_retain, i_orb, i_orb
+            ]
+            
+            tau_grid, gf_tau, const_real_tail, const_imag_tail = calc_gf_tau_from_gf_matsubara(
+                matsubara_frequencies, gf_imag_matsubara, n_tau, n_tail, beta, show_fit=show_fit
+            )
+
+            density, gf_tau_fit, energies_extract, sum_rule_const, chi2 = get_single_continuation(
+                tau_grid, gf_tau, n_sv, beta, mesh_w.values()[0], mesh_w.values()[-1], mesh_w.size, sum_rule_const=const_imag_tail, lambd=L1_coeff, verbose=verbose_opt, max_iters=max_iters, solver=solver_opt
+            )
+            energies, gf_real, gf_imag = get_kramers_kronig_realpart(
+                energies_extract, dos_to_gf_imag(density)
+            )
+            gf_real += const_real_tail
+
+            sigma_w_data[:, i_orb, i_orb] = gf_real + 1j * gf_imag
+            if params_spm["show_result"]:
+                import matplotlib.pyplot as plt
+
+                plt.axhline(y=0, xmin=energies[0], xmax=energies[-1], color="lightgrey")
+                plt.plot(energies, gf_real, label=r"Re $G( \omega )$")
+                plt.plot(energies, gf_imag, label=r"Im $G( \omega )$")
+                plt.xlim(energies[0], energies[-1])
+                plt.xlabel(r"$\omega$")
+                plt.legend()
+                plt.tight_layout()
+                plt.show()
+                plt.close()
+        sigma_w = GfReFreq(mesh=mesh_w, data=sigma_w_data)
+        data_w[key] = sigma_w.data
+    return data_w
\ No newline at end of file
diff --git a/src/dcore/anacont_spm.py b/src/dcore/anacont_spm.py
deleted file mode 100644
index 54070cf1..00000000
--- a/src/dcore/anacont_spm.py
+++ /dev/null
@@ -1,157 +0,0 @@
-#
-# DCore -- Integrated DMFT software for correlated electrons
-# Copyright (C) 2017 The University of Tokyo
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <https://www.gnu.org/licenses/>.
-#
-
-import numpy as np
-import cvxpy as cp
-from scipy.interpolate import interp1d
-from scipy.sparse.linalg import svds
-
-def set_default_values(dictionary, default_values_dict):
-    for key, value in default_values_dict.items():
-        if key not in dictionary:
-            dictionary[key] = value
-    return dictionary
-
-def _find_sum_rule_const(matsubara_frequencies, gf_wn, ntail, show_fit=False):
-    wntail = matsubara_frequencies[-ntail:]
-    gftail = gf_wn[-ntail:]
-    const_imagpart = -np.mean(np.imag(gftail) * wntail)
-    const_realpart = np.mean(np.real(gftail))
-    if show_fit:
-        import matplotlib.pyplot as plt
-        z = np.linspace(wntail[0], wntail[-1], num=1000)
-        plt.scatter(wntail, np.imag(gftail), zorder=5, color='C0', label='data')
-        plt.plot(z, -const_imagpart / z, zorder=10, color='C1', label='fit')
-        plt.xlabel(r'$\omega_n$')
-        plt.ylabel(r'Im $G( \omega_n )$')
-        plt.legend()
-        plt.tight_layout()
-        plt.show()
-    return const_realpart, const_imagpart
-
-def _calc_gf_tau(matsubara_frequencies, gf_wn, beta, const_real_tail, const_imag_tail, ntau):
-    tau_grid = np.linspace(0, beta, num=ntau)
-    gf_tau = -0.5 * const_imag_tail * np.ones(tau_grid.shape, dtype=np.float64)
-    kernel = lambda tau : np.sum((np.real(gf_wn) - const_real_tail) * np.cos(matsubara_frequencies * tau) + (np.imag(gf_wn) + const_imag_tail / matsubara_frequencies) * np.sin(matsubara_frequencies * tau))
-    gf_tau += 2 / beta * np.vectorize(kernel)(tau_grid)
-    gf_tau *= -1
-    return tau_grid, gf_tau
-
-def calc_gf_tau_from_gf_matsubara(matsubara_frequencies, gf_wn, ntau, ntail, beta, show_fit=False):
-    const_real_tail, const_imag_tail = _find_sum_rule_const(matsubara_frequencies, gf_wn, ntail, show_fit=show_fit)
-    print('Determined tail constants: {}, {}'.format(const_real_tail, const_imag_tail))
-    tau_grid, gf_tau = _calc_gf_tau(matsubara_frequencies, gf_wn, beta, const_real_tail, const_imag_tail, ntau)
-    return tau_grid, gf_tau, const_real_tail, const_imag_tail
-
-def get_single_continuation(tau_grid, gf_tau, nsv, beta, emin, emax, num_energies, sum_rule_const, lambd, verbose=True, max_iters=100, solver='ECOS'):
-    U, S, Vt, delta_energy, energies_extract = _get_svd_for_continuation(tau_grid, nsv, beta, emin, emax, num_energies)
-    rho_prime, gf_tau_fit, chi2 = _solveProblem(delta_energy, U, S, Vt, gf_tau, sum_rule_const, lambd, verbose=verbose, max_iters=max_iters, solver=solver)
-    rho = np.dot(Vt.T, rho_prime)
-    rho_integrated = np.trapz(y=rho, x=energies_extract)
-    return rho, gf_tau_fit, energies_extract, rho_integrated, chi2
-
-def get_multiple_continuations(tau_grid, gf_tau, nsv, beta, emin, emax, num_energies, sum_rule_const, lambdas, verbose=True, max_iters=100, solver='ECOS'):
-    U, S, Vt, delta_energy, energies_extract = _get_svd_for_continuation(tau_grid, nsv, beta, emin, emax, num_energies)
-    continued_densities = []
-    chi2_values = []
-    rho_values = []
-    for lambd in lambdas:
-        rho_prime, _, chi2 = _solveProblem(delta_energy, U, S, Vt, gf_tau, sum_rule_const, lambd, verbose=verbose, max_iters=max_iters, solver=solver)
-        rho = np.dot(Vt.T, rho_prime)
-        rho_integrated = np.trapz(y=rho, x=energies_extract)
-        continued_densities.append(rho)
-        chi2_values.append(chi2)
-        rho_values.append(rho_integrated)
-    return energies_extract, continued_densities, chi2_values, rho_values
-
-def _get_kernel_matrix(energies, tau_grid, beta, delta_energy, z_critical=20): #with z_critical=20 error of asymptote is below 1e-10
-    assert tau_grid[0] == 0
-    assert tau_grid[-1] == beta
-    i_lc = np.searchsorted(energies, -z_critical / beta)
-    i_rc = np.searchsorted(energies, +z_critical / beta)
-    kernel = np.zeros((len(tau_grid), len(energies)))
-    kernel[:, i_lc:i_rc + 1] = 0.5 * np.exp(np.multiply.outer(0.5 * beta - tau_grid, energies[i_lc:i_rc + 1])) / np.cosh(0.5 * beta * energies[i_lc:i_rc + 1])
-    if i_lc > 0: #asymptote for kernel for w -> -inf
-        kernel[:, 0:i_lc] = np.exp(np.multiply.outer(beta - tau_grid, energies[0:i_lc]))
-    if i_rc < len(energies) - 1: #asymptote for kernel for w -> +inf
-        kernel[:, i_rc+1:len(energies)] = np.exp(np.multiply.outer(-tau_grid, energies[i_rc+1:len(energies)]))
-    kernel[:, 0] *= 0.5
-    kernel[:, -1] *= 0.5
-    kernel *= delta_energy
-    return kernel
-
-def _getSVD(matrix, nsv=None):
-    if nsv == None or nsv >= min(matrix.shape) or nsv <= 0:
-        raise Exception("nsv must be 0 < nsv < {} (min(matrix.shape))".format(min(matrix.shape)))
-    U, S, V = svds(matrix, k=nsv, which='LM')
-    S = np.flip(S, axis=0)
-    U = np.flip(U, axis=1)
-    Vt = np.flip(V, axis=0)
-    return U, S, Vt
-
-def _solveProblem(delta_energy, U, S, Vt, gf_tau, sum_rule_const, lambd, verbose=True, max_iters=100, solver='ECOS'):
-    Q = len(S)
-    Smat = np.diag(S)
-    rho_prime = cp.Variable(Q)
-    errTerm = gf_tau - U @ Smat @ rho_prime
-    objective = cp.Minimize(0.5 * cp.norm2(errTerm) ** 2 + lambd * cp.norm1(rho_prime))
-    V_mod = np.copy(Vt.T)
-    V_mod[0, :] *= 0.5
-    V_mod[-1, :] *= 0.5
-    constraints = [Vt.T @ rho_prime >= 0, cp.sum(delta_energy * V_mod @ rho_prime) == sum_rule_const] #uniform real energy grid is assumed here
-    prob = cp.Problem(objective, constraints)
-    _ = prob.solve(verbose=verbose, solver=solver, max_iters=max_iters)
-    gf_tau_fit = np.dot(U, np.dot(Smat, rho_prime.value))
-    chi2 = 0.5 * np.linalg.norm(gf_tau - gf_tau_fit, ord=2) ** 2
-    return rho_prime.value, gf_tau_fit, chi2
-
-def _get_svd_for_continuation(tau_grid, nsv, beta, emin, emax, num_energies):
-    energies_extract = np.linspace(emin, emax, num=num_energies)
-    delta_energy = energies_extract[1] - energies_extract[0]
-    kernel = _get_kernel_matrix(energies_extract, tau_grid, beta, delta_energy)
-    U, S, Vt = _getSVD(kernel, nsv=nsv)
-    return U, S, Vt, delta_energy, energies_extract
-
-def _integral_kramers_kronig(energies_imagpart, energy_realpart, gf_imag_interp, energy_threshold):
-    enum1 = gf_imag_interp(energies_imagpart) - gf_imag_interp(energy_realpart)
-    enum2 = np.gradient(gf_imag_interp(energies_imagpart), energies_imagpart)
-    den = energies_imagpart - energy_realpart
-    mask = np.where(np.abs(den) < energy_threshold, 1, 0) #1 if den is below threshold
-    #avoid divide by zero with mask * energy_threshold, in this case (1 - mask == 0)
-    term1 = enum1 / (den + mask * energy_threshold) * (1 - mask)
-    term2 = enum2 * mask 
-    kernel = term1 + term2
-    integral = np.trapz(y=kernel, x=energies_imagpart)
-    return integral
-
-def get_kramers_kronig_realpart(energies, gf_imag, energy_threshold=1e-10, dos_threshold=1e-5):
-    if np.abs(gf_imag[0]) > dos_threshold:
-        print('Warning! DOS at left interval end exceeds {}.'.format(dos_threshold))
-    if np.abs(gf_imag[-1]) > dos_threshold:
-        print('Warning! DOS at right interval end exceeds {}.'.format(dos_threshold))
-    gf_real = np.zeros(energies.shape, dtype=np.float64)
-    gf_imag_interp = interp1d(energies, gf_imag, kind='linear', bounds_error=False, fill_value=0.0)
-    energies_noend = energies[1:-1]
-    a, b = energies[0], energies[-1]
-    gf_real[1:-1] = gf_imag_interp(energies_noend) / np.pi * np.log((b - energies_noend) / (energies_noend - a)) #assume zero DOS at endpoints
-    integral_func = lambda y : _integral_kramers_kronig(energies, y, gf_imag_interp, energy_threshold)
-    gf_real += np.vectorize(integral_func)(energies) / np.pi
-    return energies, gf_real, gf_imag
-
-def dos_to_gf_imag(dos):
-    return -np.pi * dos
diff --git a/src/dcore/dcore_anacont.py b/src/dcore/dcore_anacont.py
new file mode 100644
index 00000000..5134e53f
--- /dev/null
+++ b/src/dcore/dcore_anacont.py
@@ -0,0 +1,92 @@
+#
+# DCore -- Integrated DMFT software for correlated electrons
+# Copyright (C) 2017 The University of Tokyo
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+#
+
+import argparse
+import os.path
+import numpy
+from dcore._dispatcher import MeshReFreq
+from dcore.version import version, print_header
+from dcore.program_options import create_parser, parse_parameters
+from dcore.anacont import pade, spm
+
+def dcore_anacont(inifile):
+    parser = create_parser(["system", "model", "post", "post.anacont"])
+    parser.read(inifile)
+    params = parser.as_dict()
+    parse_parameters(params)
+
+    beta = params["system"]["beta"]
+
+    omega_min = params["post"]["omega_min"]
+    omega_max = params["post"]["omega_max"]
+    if omega_min >= omega_max:
+        # ToDo: stop the program properly
+        assert omega_min < omega_max
+
+    Nomega = params["post"]["Nomega"]
+    mesh_w = MeshReFreq(omega_min, omega_max, Nomega)
+
+    seedname = params["model"]["seedname"]
+    file_sigma_iw = seedname + "_sigma_iw.npz"
+    # file_sigma_iw = params["post"]["file_sigma_iw"]
+    if not os.path.exists(file_sigma_iw):
+        assert False, "File not found: " + file_sigma_iw
+    sigma_iw_npz = numpy.load(file_sigma_iw)
+
+    dir_work = params["post"]["dir_work"]
+    if not os.path.exists(dir_work):
+        os.makedirs(dir_work)
+
+    solver = params["post.anacont"]["solver"]
+    if solver == "pade":
+        params_ac = pade.parameters_from_ini(inifile)
+        data_w = pade.anacont(sigma_iw_npz, beta, mesh_w, params_ac)
+    elif solver == "spm":
+        params_ac = spm.parameters_from_ini(inifile)
+        data_w = spm.anacont(sigma_iw_npz, beta, mesh_w, params_ac)
+    else:
+        assert False, "Unknown solver: " + solver
+
+    file_sigma_w = os.path.join(dir_work, "sigma_w.npz")
+    print("Writing to", file_sigma_w + "...")
+    numpy.savez(file_sigma_w, **data_w)
+
+def run():
+
+    print_header()
+
+    parser = argparse.ArgumentParser(
+        prog='dcore_anacont.py',
+        description='DCore post script -- analytic continuation.',
+        usage='$ dcore_anacont input.ini',
+        add_help=True,
+        formatter_class=argparse.RawTextHelpFormatter,
+        #epilog=generate_all_description()
+    )
+    parser.add_argument('path_input_files',
+                        action='store',
+                        default=None,
+                        type=str,
+                        nargs='*',
+                        help="Input filename(s)",
+                        )
+    parser.add_argument('--version', action='version', version='DCore {}'.format(version))
+
+    args = parser.parse_args()
+
+    dcore_anacont(args.path_input_files)
\ No newline at end of file
diff --git a/src/dcore/dcore_anacont_pade.py b/src/dcore/dcore_anacont_pade.py
index ce2149c8..16c2bfe9 100644
--- a/src/dcore/dcore_anacont_pade.py
+++ b/src/dcore/dcore_anacont_pade.py
@@ -17,70 +17,59 @@
 #
 
 import argparse
-from dcore._dispatcher import MeshReFreq, MeshImFreq, GfReFreq, GfImFreq
-from dcore.version import version, print_header
 import numpy
 import toml
 
-def _set_n_pade(omega_cutoff, beta, n_min, n_max):
-    """
-    Return (int)n_pade: the number of Matsubara frequencies below the cutoff frequency.
-    n_pade is bounded between n_min and n_max
-    """
-    n_pade = int((beta * omega_cutoff + numpy.pi) / (2.0 * numpy.pi))
-    print("n_pade = {} (evaluated from omega_pade)".format(n_pade))
-    n_pade = max(n_pade, n_min)
-    n_pade = min(n_pade, n_max)
-    print("n_pade = {}".format(n_pade))
-    return n_pade
+from dcore._dispatcher import MeshReFreq
+from dcore.version import version, print_header
+from dcore.anacont.pade import anacont
+
 
-def dcore_anacont_pade(seedname):
+def dcore_anacont_pade_from_seedname(seedname):
     print("Reading ", seedname + "_anacont.toml...")
     with open(seedname + "_anacont.toml", "r") as f:
         params = toml.load(f)
- 
+
     print("Reading ", seedname + "_sigma_iw.npz...")
-    npz = numpy.load(seedname + "_sigma_iw.npz")
+    sigma_iw_npz = numpy.load(seedname + "_sigma_iw.npz")
 
     assert params["omega_min"] < params["omega_max"]
     mesh_w = MeshReFreq(params["omega_min"], params["omega_max"], params["Nomega"])
 
-    n_pade = _set_n_pade(params["pade"]["omega_max"], params['beta'], n_min=params["pade"]["n_min"], n_max=params["pade"]["n_max"])
+    beta = params["beta"]
+    params_pade = params.get("pade", {})
+    params_pade["iomega_max"] = params_pade.get("omega_max", -1.0)
+    params_pade["eta"] = params_pade.get("eta", 0.01)
+    params_pade["n_min"] = params_pade.get("n_min", 0)
+    params_pade["n_max"] = params_pade.get("n_max", 100000000)
 
-    data_w = {}
+    data_w = anacont(
+        sigma_iw_npz, beta=beta, mesh_w=mesh_w, params_pade=params_pade
+    )
 
-    num_data = numpy.sum([key.startswith("data") for key in npz.keys()])
-    for idata in range(num_data):
-        key = f"data{idata}"
-        data = npz[key]
-        mesh_iw = MeshImFreq(params["beta"], "Fermion", data.shape[0]//2)
-        sigma_iw = GfImFreq(data=data, beta=params["beta"], mesh=mesh_iw)
-        sigma_w = GfReFreq(mesh=mesh_w, target_shape=data.shape[1:])
-        sigma_w.set_from_pade(sigma_iw, n_points=n_pade, freq_offset=params["pade"]["eta"])
-        data_w[key] = sigma_w.data
     print("Writing to", seedname + "_sigma_w.npz...")
     numpy.savez(seedname + "_sigma_w.npz", **data_w)
 
+
 def run():
 
     print_header()
 
     parser = argparse.ArgumentParser(
-        prog='dcore_anacont_pade.py',
-        description='pre script for dcore.',
-        usage='$ dcore_anacont_pade input',
+        prog="dcore_anacont_pade.py",
+        description="pre script for dcore.",
+        usage="$ dcore_anacont_pade input",
         add_help=True,
         formatter_class=argparse.RawTextHelpFormatter,
-        #epilog=generate_all_description()
+        # epilog=generate_all_description()
+    )
+    parser.add_argument(
+        "seedname", action="store", default=None, type=str, help="seedname"
+    )
+    parser.add_argument(
+        "--version", action="version", version="DCore {}".format(version)
     )
-    parser.add_argument('seedname',
-                        action='store',
-                        default=None,
-                        type=str,
-                        help="seedname"
-                        )
-    parser.add_argument('--version', action='version', version='DCore {}'.format(version))
 
     args = parser.parse_args()
 
-    dcore_anacont_pade(args.seedname)
+    dcore_anacont_pade_from_seedname(args.seedname)
diff --git a/src/dcore/dcore_anacont_spm.py b/src/dcore/dcore_anacont_spm.py
index 476a955e..f45752d0 100644
--- a/src/dcore/dcore_anacont_spm.py
+++ b/src/dcore/dcore_anacont_spm.py
@@ -19,77 +19,48 @@
 import argparse
 import toml
 import numpy as np
-from dcore._dispatcher import MeshReFreq, MeshImFreq, GfReFreq, GfImFreq
+from dcore._dispatcher import MeshReFreq
 from dcore.version import version, print_header
-from dcore.anacont_spm import set_default_values, calc_gf_tau_from_gf_matsubara, get_single_continuation, get_kramers_kronig_realpart, dos_to_gf_imag
+from dcore.anacont.spm import set_default_values, anacont
 
-def _anacont_spm_per_gf(params, matsubara_frequencies, gf_matsubara):
-    tau_grid, gf_tau, const_real_tail, const_imag_tail = calc_gf_tau_from_gf_matsubara(matsubara_frequencies, gf_matsubara, params['spm']['n_tau'], params['spm']['n_tail'], params['beta'], show_fit=params['spm']['show_fit'])
-    density, gf_tau_fit, energies_extract, density_integrated, chi2 = get_single_continuation(tau_grid, gf_tau, params['spm']['n_sv'], params['beta'], params['omega_min'], params['omega_max'], params['Nomega'], const_imag_tail, params['spm']['lambda'], verbose=params['spm']['verbose_opt'], max_iters=params['spm']['max_iters_opt'], solver=params['spm']['solver_opt'])
-    energies, gf_real, gf_imag = get_kramers_kronig_realpart(energies_extract, dos_to_gf_imag(density))
-    gf_real += const_real_tail
-    return energies, gf_real, gf_imag
 
 def set_default_config(params):
-    default_values = {'show_fit' : False,
-                      'show_result' : False,
-                      'show_fit': False,
-                      'verbose_opt' : False,
-                      'max_iters_opt' : 100,
-                      'solver_opt' : 'ECOS'}
-    params['spm'] = set_default_values(params['spm'], default_values)
+    default_values = {
+        "show_fit": False,
+        "show_result": False,
+        "show_fit": False,
+        "verbose_opt": False,
+        "max_iters_opt": 100,
+        "solver_opt": "ECOS",
+    }
+    params["spm"] = set_default_values(params["spm"], default_values)
     return params
 
-def dcore_anacont_spm(seedname):
-    print('Reading ', seedname + '_anacont.toml...')
-    with open(seedname + '_anacont.toml', 'r') as f:
+
+def dcore_anacont_spm_from_seedname(seedname):
+    print("Reading ", seedname + "_anacont.toml...")
+    with open(seedname + "_anacont.toml", "r") as f:
         params = toml.load(f)
     params = set_default_config(params)
-    print('Using configuration: ', params)
- 
-    print('Reading ', seedname + '_sigma_iw.npz...')
-    npz = np.load(seedname + '_sigma_iw.npz')
-
-    assert params['omega_min'] < params['omega_max']
-    mesh_w = MeshReFreq(params['omega_min'], params['omega_max'], params['Nomega'])
-
-    assert params['spm']['n_matsubara'] > 0
-
-    data_w = {}
-    num_data = np.sum([key.startswith('data') for key in npz.keys()])
-    for idata in range(num_data):
-        key = f'data{idata}'
-        data = npz[key]
-        n_matsubara = data.shape[0]//2
-        n_matsubara_retain = min(n_matsubara, params['spm']['n_matsubara'])
-        mesh_iw = MeshImFreq(params['beta'], 'Fermion', n_matsubara_retain)
-        gf_iw = GfImFreq(data=data, beta=params['beta'], mesh=mesh_iw)
-        n_orbitals = gf_iw.data.shape[1]
-        matsubara_frequencies = np.imag(gf_iw.mesh.values()[n_matsubara_retain:])
-        sigma_w_data = np.zeros((params['Nomega'], n_orbitals, n_orbitals), dtype=np.complex128)
-        for i_orb in range(n_orbitals):
-            print(f'Performing analytic continuation for data index {idata} and orbital index {i_orb}...')
-            gf_imag_matsubara = gf_iw.data[n_matsubara:n_matsubara + n_matsubara_retain, i_orb, i_orb]
-            energies, gf_real, gf_imag = _anacont_spm_per_gf(params, matsubara_frequencies, gf_imag_matsubara)
-            sigma_w_data[:, i_orb, i_orb] = gf_real + 1j * gf_imag
-            if params['spm']['show_result']:
-                import matplotlib.pyplot as plt
-                plt.axhline(y=0, xmin=energies[0], xmax=energies[-1], color='lightgrey')
-                plt.plot(energies, gf_real, label=r'Re $G( \omega )$')
-                plt.plot(energies, gf_imag, label=r'Im $G( \omega )$')
-                plt.xlim(energies[0], energies[-1])
-                plt.xlabel(r'$\omega$')
-                plt.legend()
-                plt.tight_layout()
-                plt.show()
-                plt.close()
-        sigma_w = GfReFreq(mesh=mesh_w, data=sigma_w_data)
-        data_w[key] = sigma_w.data
-
-    print('Writing to', seedname + '_sigma_w.npz...')
-    np.savez(seedname + '_sigma_w.npz', **data_w)
-
-#example file for 'seedname_anacont.toml'
+    print("Using configuration: ", params)
+
+    print("Reading ", seedname + "_sigma_iw.npz...")
+    npz = np.load(seedname + "_sigma_iw.npz")
+
+    assert params["omega_min"] < params["omega_max"]
+    mesh_w = MeshReFreq(params["omega_min"], params["omega_max"], params["Nomega"])
+
+    assert params["spm"]["n_matsubara"] > 0
+
+    beta = params["beta"]
+    params_spm = params["spm"]
+    data_w = anacont(npz, beta, mesh_w, params_spm)
+
+    print("Writing to", seedname + "_sigma_w.npz...")
+    np.savez(seedname + "_sigma_w.npz", **data_w)
+
+
+# example file for 'seedname_anacont.toml'
 # beta = 40.0
 # Nomega = 4000
 # omega_min = -6.0
@@ -107,25 +78,25 @@ def dcore_anacont_spm(seedname):
 # max_iters_opt = 100
 # solver_opt = 'ECOS'
 
+
 def run():
     print_header()
 
     parser = argparse.ArgumentParser(
-        prog='dcore_anacont_spm.py',
-        description='pre script for dcore.',
-        usage='$ dcore_anacont_spm input',
+        prog="dcore_anacont_spm.py",
+        description="pre script for dcore.",
+        usage="$ dcore_anacont_spm input",
         add_help=True,
         formatter_class=argparse.RawTextHelpFormatter,
-        #epilog=generate_all_description()
+        # epilog=generate_all_description()
+    )
+    parser.add_argument(
+        "seedname", action="store", default=None, type=str, help="seedname"
+    )
+    parser.add_argument(
+        "--version", action="version", version="DCore {}".format(version)
     )
-    parser.add_argument('seedname',
-                        action='store',
-                        default=None,
-                        type=str,
-                        help='seedname'
-                        )
-    parser.add_argument('--version', action='version', version='DCore {}'.format(version))
 
     args = parser.parse_args()
 
-    dcore_anacont_spm(args.seedname)
+    dcore_anacont_spm_from_seedname(args.seedname)
diff --git a/src/dcore/dcore_anacont_spm_interactive.py b/src/dcore/dcore_anacont_spm_interactive.py
index d7bd0d97..133fc3e5 100644
--- a/src/dcore/dcore_anacont_spm_interactive.py
+++ b/src/dcore/dcore_anacont_spm_interactive.py
@@ -23,7 +23,7 @@
 from matplotlib.ticker import MultipleLocator
 from dcore._dispatcher import MeshReFreq, MeshImFreq, GfReFreq, GfImFreq
 from dcore.version import version, print_header
-from dcore.anacont_spm import set_default_values, calc_gf_tau_from_gf_matsubara, get_single_continuation, get_multiple_continuations, get_kramers_kronig_realpart, dos_to_gf_imag
+from dcore.anacont.spm import set_default_values, calc_gf_tau_from_gf_matsubara, get_single_continuation, get_multiple_continuations, get_kramers_kronig_realpart, dos_to_gf_imag
 
 def _plot_overview(lambdas, chi2_values, energies, densities, nrows=3, ncols=5):
     fig, axes = plt.subplots(nrows=nrows, ncols=ncols, sharex=False, figsize=(15, 10))
diff --git a/src/dcore/program_options.py b/src/dcore/program_options.py
index f6fba8b4..53312425 100644
--- a/src/dcore/program_options.py
+++ b/src/dcore/program_options.py
@@ -46,7 +46,7 @@ def create_parser(target_sections=None):
     Create a parser for all program options of DCore
     """
     if target_sections is None:
-        parser = TypedParser(['mpi', 'model', 'system', 'impurity_solver', 'control', 'tool', 'bse', 'vertex', 'sparse_bse'])
+        parser = TypedParser(['mpi', 'model', 'system', 'impurity_solver', 'control', 'tool', 'post', 'bse', 'vertex', 'sparse_bse'])
     else:
         parser = TypedParser(target_sections)
 
@@ -143,6 +143,36 @@ def create_parser(target_sections=None):
         "Each line after the first line must contain two integer numbers: The first one is an sequential index of a sampling frequency (start from 0), and the second one is a fermionic sampling frequency (i.e, -1, 0, 1)."
     )
 
+    # [post]
+    parser.add_option("post", "omega_min", float, -1, "Minimum value of real frequency")
+    parser.add_option("post", "omega_max", float, 1, "Max value of real frequency")
+    parser.add_option("post", "Nomega", int, 100, "Number of real frequencies")
+    # parser.add_option("post", "file_sigma_iw", str, "sigma_iw.npz", "Filename of the self-energy in Matsubara frequency")
+    parser.add_option("post", "dir_post", str, "post", "Directory to which results of dcore_post are stored.")
+    parser.add_option("post", "dir_work", str, "work/post", "Directory to which temporary files of dcore_post are stored.")
+
+    # [post.anacont]
+    parser.add_option("post.anacont", "solver", str, "algorithm", "Algorithm for analytic continuation")
+
+    # [post.anacont.pade]
+    parser.add_option( "post.anacont.pade", "iomega_max", float, -1.0, "Cut-off frequency of the Matsubara frequency",)
+    parser.add_option( "post.anacont.pade", "n_min", int, 0, "lower bound of the order of Pade approximation",)
+    parser.add_option( "post.anacont.pade", "n_max", int, 100000000, "upper bound of the order of Pade approximation",)
+    parser.add_option( "post.anacont.pade", "eta", float, 0.01, "Imaginary Frequency shift to avoid divergence",)
+
+    # [post.anacont.spm]
+    parser.add_option("post.anacont.spm", "n_matsubara", int, 100000, "number of tau points")
+    parser.add_option("post.anacont.spm", "n_tau", int, -1, "number of tau points")
+    parser.add_option("post.anacont.spm", "n_tail", int, 10, "number of matsubara points for tail-fitting")
+    parser.add_option("post.anacont.spm", "n_sv", int, 50, "number of singular values to be used")
+    parser.add_option("post.anacont.spm", "lambda", float, 1e-5, "coefficient of L1 regularization")
+    parser.add_option("post.anacont.spm", "max_iters_opt", int, 100, "maximum number of iterations")
+    parser.add_option("post.anacont.spm", "solver_opt", str, "ECOS", "solver to be used")
+
+    parser.add_option("post.anacont.spm", "verbose_opt", bool, False, "show optimization progress")
+    parser.add_option("post.anacont.spm", "show_fit", bool, False, "plot result of tail-fitting")
+    parser.add_option("post.anacont.spm", "show_result", bool, False, "plot result of analytic continuation")
+
     # [bse]
     parser.add_option("bse", "num_wb", int, 1, "Number of bosonic frequencies (>0)")
     parser.add_option("bse", "num_wf", int, 10, "Number of fermionic frequencies (>0)")
diff --git a/tests/non-mpi/anacont_spm/anacont_spm.py b/tests/non-mpi/anacont_spm/anacont_spm.py
index ab94c406..93c4d530 100644
--- a/tests/non-mpi/anacont_spm/anacont_spm.py
+++ b/tests/non-mpi/anacont_spm/anacont_spm.py
@@ -32,7 +32,7 @@ def _get_dos_semicircular(energies, r):
     return 2.0 / np.pi / r * np.sqrt(np.maximum(1 - np.square(energies / r), 0))
 
 def test_find_sum_rule_const():
-    from dcore.anacont_spm import _find_sum_rule_const
+    from dcore.anacont.spm import _find_sum_rule_const
     a = 1.123
     b = 2.234
     beta = 40
@@ -45,7 +45,7 @@ def test_find_sum_rule_const():
     assert np.allclose(b_test, b, atol=1e-9)
 
 def test_calc_gf_tau_trivial():
-    from dcore.anacont_spm import _calc_gf_tau
+    from dcore.anacont.spm import _calc_gf_tau
     a = 1.123
     b = 2.234
     beta = 40
@@ -60,7 +60,7 @@ def test_calc_gf_tau_trivial():
     assert np.allclose(gf_tau, gf_tau_expected, atol=1e-10)
 
 def test_calc_gf_tau_nontrivial():
-    from dcore.anacont_spm import _calc_gf_tau, _find_sum_rule_const
+    from dcore.anacont.spm import _calc_gf_tau, _find_sum_rule_const
     beta = 40
     n_matsubara = 1000
     n_matsubara_tail = 100
@@ -87,7 +87,7 @@ def test_calc_gf_tau_nontrivial():
     assert np.allclose(gf_tau, gf_tau_expected, atol=1e-10)
 
 def test_calc_gf_tau_from_gf_matsubara():
-    from dcore.anacont_spm import calc_gf_tau_from_gf_matsubara
+    from dcore.anacont.spm import calc_gf_tau_from_gf_matsubara
     ntau = 11
     n_matsubara = 1000
     n_matsubara_tail = 100
@@ -111,7 +111,7 @@ def test_calc_gf_tau_from_gf_matsubara():
     assert np.allclose(gf_tau, gf_tau_expected, atol=1e-10)
 
 def test_get_kernel_matrix():
-    from dcore.anacont_spm import _get_kernel_matrix
+    from dcore.anacont.spm import _get_kernel_matrix
     beta = 40
     energies = np.linspace(-3, 3, num=5)
     delta_energy = energies[1] - energies[0]
@@ -121,7 +121,7 @@ def test_get_kernel_matrix():
     assert np.allclose(kernel, kernel_expected, atol=1e-10)
 
 def test_get_kernel_matrix_extreme_energies():
-    from dcore.anacont_spm import _get_kernel_matrix
+    from dcore.anacont.spm import _get_kernel_matrix
     beta = 40
     energies = np.linspace(-30, 30, num=5)
     delta_energy = energies[1] - energies[0]
@@ -134,7 +134,7 @@ def test_get_kernel_matrix_extreme_energies():
     assert np.allclose(kernel, kernel_expected, atol=1e-10)
 
 def test_getSVD():
-    from dcore.anacont_spm import _getSVD, _get_kernel_matrix
+    from dcore.anacont.spm import _getSVD, _get_kernel_matrix
     beta = 40
     nsv = 2
     ntau = 3
@@ -150,7 +150,7 @@ def test_getSVD():
     assert Vt.shape == (nsv, nenergies)
 
 def test_get_svd_for_continuation():
-    from dcore.anacont_spm import _get_svd_for_continuation
+    from dcore.anacont.spm import _get_svd_for_continuation
     beta = 40
     nsv = 2
     ntau = 5
@@ -171,7 +171,7 @@ def test_get_svd_for_continuation():
     assert np.allclose(energies_extract, energies_extract_expected, atol=1e-7)
 
 def test_solveProblem():
-    from dcore.anacont_spm import _solveProblem, _find_sum_rule_const, _calc_gf_tau, _get_svd_for_continuation
+    from dcore.anacont.spm import _solveProblem, _find_sum_rule_const, _calc_gf_tau, _get_svd_for_continuation
     beta = 40
     nsv = 100
     n_matsubara = 1000
@@ -202,7 +202,7 @@ def test_solveProblem():
     assert np.allclose(chi2, 3.4763690885418575e-06, atol=1e-5)
 
 def test_integral_kramers_kronig():
-    from dcore.anacont_spm import _integral_kramers_kronig, dos_to_gf_imag
+    from dcore.anacont.spm import _integral_kramers_kronig, dos_to_gf_imag
     from scipy.interpolate import interp1d
 
     energies = np.linspace(-5, 5, num=1000)
@@ -225,7 +225,7 @@ def test_integral_kramers_kronig():
     assert np.allclose(results, expected_results, atol=1e-7)
 
 def test_get_kramers_kronig_realpart():
-    from dcore.anacont_spm import get_kramers_kronig_realpart, dos_to_gf_imag
+    from dcore.anacont.spm import get_kramers_kronig_realpart, dos_to_gf_imag
     energies = np.linspace(-5, 5, num=10)
     dos = _get_dos_semicircular(energies, 4)
     gf_imag = dos_to_gf_imag(dos)
@@ -237,7 +237,7 @@ def test_get_kramers_kronig_realpart():
     assert np.allclose(gf_real_result, gf_real_expected, atol=1e-7)
 
 def test_get_single_continuation():
-    from dcore.anacont_spm import get_single_continuation, calc_gf_tau_from_gf_matsubara, calc_gf_tau_from_gf_matsubara
+    from dcore.anacont.spm import get_single_continuation, calc_gf_tau_from_gf_matsubara, calc_gf_tau_from_gf_matsubara
     beta = 40
     nsv = 24
     emin = -10
@@ -269,7 +269,7 @@ def test_get_single_continuation():
 
 
 def test_get_multiple_continuations():
-    from dcore.anacont_spm import get_multiple_continuations, calc_gf_tau_from_gf_matsubara, calc_gf_tau_from_gf_matsubara
+    from dcore.anacont.spm import get_multiple_continuations, calc_gf_tau_from_gf_matsubara, calc_gf_tau_from_gf_matsubara
     beta = 40
     nsv = 24
     emin = -10