diff --git a/img_formula_70.png b/img_formula_70.png
new file mode 100644
index 0000000..5b8eed4
Binary files /dev/null and b/img_formula_70.png differ
diff --git a/main.py b/main.py
index b9dc3ee..9bdb8e0 100644
--- a/main.py
+++ b/main.py
@@ -17,101 +17,497 @@ poles = vf.poles
 residues = vf.residues
 y = vf.network.y
 
-fig, ax = plt.subplots(2, 2)
-fig.set_size_inches(12, 8)
-vf.plot("mag",0,0,freqs,ax=ax[0][0],parameter="y")
-rms_error11 = vf.get_rms_error(0,0,"y")
-print("rms_error11",rms_error11)
-vf.plot("mag",1,0,freqs,ax=ax[1][0],parameter="y")
-rms_error21 = vf.get_rms_error(1,0,"y")
-print("rms_error21",rms_error21)
-vf.plot("mag",0,1,freqs,ax=ax[0][1],parameter="y")
-rms_error12 = vf.get_rms_error(0,1,"y")
-print("rms_error12",rms_error12)
-vf.plot("mag",1,1,freqs,ax=ax[1][1],parameter="y")
-rms_error22 = vf.get_rms_error(1,1,"y")
-print("rms_error22",rms_error22)
-fig.tight_layout()
-plt.show()
-plt.savefig(f"img.png")
+# fig, ax = plt.subplots(2, 2)
+# fig.set_size_inches(12, 8)
+# vf.plot("mag",0,0,freqs,ax=ax[0][0],parameter="y")
+# rms_error11 = vf.get_rms_error(0,0,"y")
+# print("rms_error11",rms_error11)
+# vf.plot("mag",1,0,freqs,ax=ax[1][0],parameter="y")
+# rms_error21 = vf.get_rms_error(1,0,"y")
+# print("rms_error21",rms_error21)
+# vf.plot("mag",0,1,freqs,ax=ax[0][1],parameter="y")
+# rms_error12 = vf.get_rms_error(0,1,"y")
+# print("rms_error12",rms_error12)
+# vf.plot("mag",1,1,freqs,ax=ax[1][1],parameter="y")
+# rms_error22 = vf.get_rms_error(1,1,"y")
+# print("rms_error22",rms_error22)
+# fig.tight_layout()
+# plt.show()
+# plt.savefig(f"img.png")
 
-diag_values = [y[i][0][0] for i in range(len(y))]
-H = np.diag(diag_values)
-P = 4
-start_poles = generate_starting_poles(P,beta_min=1e8,beta_max=freqs[-1])
-basis = generate_laguerre_basis(start_poles,s).T
-print("start_poles",start_poles)
-print("basis",basis)
+def formula_67(s,y):
+    diag_values = [y[i][0][0] for i in range(len(y))]
+    H = np.diag(diag_values)
+    P = 2
+    start_poles = generate_starting_poles(P,beta_min=1e8,beta_max=freqs[-1])
+    basis = generate_laguerre_basis(start_poles,s).T
+    print("start_poles",start_poles)
+    print("basis",basis)
 
-# first step iteration
-# A*x = b
+    # first step iteration
+    # A*x = b
 
-A11 = np.real(basis)
-print("A11 shape:",A11.shape)
-A12 = np.real(- H @ basis[:,1:])
-print("A12 shape:",A12.shape)
-# print("A11",A11)
-A21 = np.imag(basis)
-print("A21 shape:",A21.shape)
-A22 = np.imag(- H @ basis[:,1:])
-print("A22 shape:",A22.shape)
-A1 = np.hstack([A11,A12])
-A2 = np.hstack([A21,A22])
-A = np.vstack([A1,A2])
-print ("A shape:",A.shape)
-
-b1 = np.real(H @ basis[:,0])
-b2 = np.imag(H @ basis[:,0])
-b = np.hstack([b1,b2])
-Q, R = np.linalg.qr(A, mode='reduced')
-print("Q :",Q)
-print("R :",R)
-print("b shape:",b.shape)
-# x = np.linalg.solve(R, Q.T @ b)
-x_star, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
-
-# print("x_qr",x_star)
-print("residuals",residuals)
-# print("rank",rank)
-# print("s",s)
-
-x_ne = np.linalg.inv(A.T @ A) @ A.T @ b
-print("x_ne",x_ne)
-
-# sk iteration
-target_vec = np.array([1+0j] + list(x_ne[len(x_ne)//2+1:]))
-D_t = basis @ target_vec
-print("D_t",D_t)
-K = 25
-
-for i in range(K):
-    print(f"Iteration {i+1}/{K}")
-    A11 = np.real(basis / D_t[:, np.newaxis])
-    A12 = np.real(- H @ basis[:,1:] / D_t[:, np.newaxis])
-    A21 = np.imag(basis / D_t[:, np.newaxis])
-    A22 = np.imag(- H @ basis[:,1:] / D_t[:, np.newaxis])
+    A11 = np.real(basis)
+    print("A11 shape:",A11.shape)
+    A12 = np.real(- H @ basis[:,1:])
+    print("A12 shape:",A12.shape)
+    # print("A11",A11)
+    A21 = np.imag(basis)
+    print("A21 shape:",A21.shape)
+    A22 = np.imag(- H @ basis[:,1:])
+    print("A22 shape:",A22.shape)
     A1 = np.hstack([A11,A12])
     A2 = np.hstack([A21,A22])
     A = np.vstack([A1,A2])
+    print ("A shape:",A.shape)
+
     b1 = np.real(H @ basis[:,0])
     b2 = np.imag(H @ basis[:,0])
     b = np.hstack([b1,b2])
+    Q, R = np.linalg.qr(A, mode='reduced')
+    print("Q :",Q)
+    print("R :",R)
+    print("b shape:",b.shape)
+    # x = np.linalg.solve(R, Q.T @ b)
     x_star, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
-    # print("x_lstsq",x_star)
-    # print("residuals",residuals)
+
+    # print("x_qr",x_star)
+    print("residuals",residuals)
     # print("rank",rank)
     # print("s",s)
+
     x_ne = np.linalg.inv(A.T @ A) @ A.T @ b
-    # print("x_ne",x_ne)
+    print("x_ne",x_ne)
+
+    # sk iteration
     target_vec = np.array([1+0j] + list(x_ne[len(x_ne)//2+1:]))
-    D_t_pre = D_t
     D_t = basis @ target_vec
-    print("D_t/D_t_pre",np.linalg.norm(D_t)/np.linalg.norm(D_t_pre))
-    # print("D_t",D_t)
-    n_target_vec = np.array(list(x_ne[:len(x_ne)//2+1]))
-    N_t = basis @ n_target_vec
-    # print("H = N_t / D_t",np.linalg.norm(H - N_t / D_t)/np.linalg.norm(H))
+    print("D_t",D_t)
+    K = 25
+
+    for i in range(K):
+        print(f"Iteration {i+1}/{K}")
+        A11 = np.real(basis / D_t[:, np.newaxis])
+        A12 = np.real(- H @ basis[:,1:] / D_t[:, np.newaxis])
+        A21 = np.imag(basis / D_t[:, np.newaxis])
+        A22 = np.imag(- H @ basis[:,1:] / D_t[:, np.newaxis])
+        A1 = np.hstack([A11,A12])
+        A2 = np.hstack([A21,A22])
+        A = np.vstack([A1,A2])
+        b1 = np.real(H @ basis[:,0])
+        b2 = np.imag(H @ basis[:,0])
+        b = np.hstack([b1,b2])
+        x_star, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
+        # print("x_lstsq",x_star)
+        # print("residuals",residuals)
+        # print("rank",rank)
+        # print("s",s)
+        x_ne = np.linalg.inv(A.T @ A) @ A.T @ b
+        # print("x_ne",x_ne)
+        target_vec = np.array([1+0j] + list(x_ne[len(x_ne)//2+1:]))
+        D_t_pre = D_t
+        D_t = basis @ target_vec
+        print("Dt", D_t)
+        # print("D_t/D_t_pre",D_t/D_t_pre)
+        # print("D_t",D_t)
+        n_target_vec = np.array(list(x_ne[:len(x_ne)//2+1]))
+        N_t = basis @ n_target_vec
+        # print("H = N_t / D_t",np.abs(N_t / D_t))
+
+class formula_70:
+
+    # ---------- helpers ----------
+    def _psi(self,s, z_prev):
+        """Psi[k,p] = 1/(s_k + z_prev[p])  (shape: Kf x P)"""
+        s = np.asarray(s, dtype=np.complex128).reshape(-1)
+        z_prev = np.asarray(z_prev, dtype=np.complex128).reshape(-1)
+        return 1.0 / (s[:, None] + z_prev[None, :])
+
+    def _enforce_lhp_conjugate(self, poles):
+        """Reflect any RHP poles across the imag axis and pair conjugates loosely."""
+        z = np.array(poles, dtype=np.complex128)
+        for i in range(len(z)):
+            if z[i].real > 0:
+                z[i] = -np.conj(z[i])
+        return z
+
+    # ---------- build Ax=b for (70) ----------1
+    def build_vf70_system(self,s, H_list, z_prev, d0=1.0):
+        """
+        s:      (Kf,) complex sample points (e.g. j*2πf)
+        H_list: list of (Kf,) complex responses, one per MIMO element you want to fit
+        z_prev: (P,) complex poles from previous iteration
+        d0:     fixed scalar, usually 1.0
+
+        Returns: A (2*Kf*M, (M+1)*P), b (2*Kf*M,), plus Psi for convenience
+        Unknown vector x = [c (P), r^(1) (P), ..., r^(M) (P)]
+        """
+        H_list = [np.asarray(h, dtype=np.complex128).reshape(-1) for h in H_list]
+        Kf = len(s)
+        P = len(z_prev)
+        M = len(H_list)
+
+        Psi = self._psi(s, z_prev)              # Kf x P
+        zerosP = np.zeros((Kf, P))         # convenience
+
+        rows = []
+        rhs  = []
+
+        for m, H in enumerate(H_list):
+            Hp = H[:, None]                # Kf x 1
+
+            # Left block (shared c):  H * Psi
+            L_re = np.real(Hp * Psi)       # Kf x P
+            L_im = np.imag(Hp * Psi)       # Kf x P
+
+            # Right block for this response: -Psi
+            R_re = -np.real(Psi)           # Kf x P
+            R_im = -np.imag(Psi)           # Kf x P
+
+            # Assemble columns: [c | r^(1) | r^(2) | ... | r^(M)]
+            # Real rows
+            cols_re = [L_re]
+            for j in range(M):
+                cols_re.append(R_re if j == m else zerosP)
+            rows.append(np.hstack(cols_re))
+
+            # Imag rows
+            cols_im = [L_im]
+            for j in range(M):
+                cols_im.append(R_im if j == m else zerosP)
+            rows.append(np.hstack(cols_im))
+
+            # RHS (move d0*H to the right)
+            rhs.append(-np.real(d0 * H))
+            rhs.append(-np.imag(d0 * H))
+
+        A = np.vstack(rows)                # (2*Kf*M) x ((M+1)*P)
+        b = np.concatenate(rhs)            # (2*Kf*M,)
+        return A, b, Psi
+
+    # ---------- one relocation step ----------
+    def vf70_step_once(self,s, H_list, z_prev, d0=1.0, column_scale=True):
+        """
+        Returns: z_new (P,), c (P,), residues (list of M arrays length P)
+        """
+        P = len(z_prev)
+        M = len(H_list)
+
+        A, b, Psi = self.build_vf70_system(s, H_list, z_prev, d0=d0)
+
+        # optional simple column scaling for conditioning
+        if column_scale:
+            col_norm = np.maximum(np.linalg.norm(A, axis=0), 1e-12)
+            x, *_ = np.linalg.lstsq(A / col_norm, b, rcond=None)
+            x = x / col_norm
+        else:
+            x, *_ = np.linalg.lstsq(A, b, rcond=None)
+
+        # unpack unknowns
+        c = x[:P]
+        residues = []
+        off = P
+        for _ in range(M):
+            residues.append(x[off:off+P])
+            off += P
+
+        # pole relocation (zeros of d0 + Psi @ c)
+        Sigma = np.diag(-np.asarray(z_prev, dtype=np.complex128))
+        T = Sigma - (np.ones((P, 1), dtype=np.complex128) @ (c.reshape(1, -1) / d0))
+        z_new = -np.linalg.eigvals(T)
+        z_new = self._enforce_lhp_conjugate(z_new)
+
+        return z_new, c, residues
+
+    # ---------- iterate K times (K == n_iter) ----------
+    def vf70_iterate(self,s, H_list, z0, n_iter=15, d0=1.0, verbose=True):
+        z = np.array(z0, dtype=np.complex128)
+        c = None
+        residues = None
+        for it in range(n_iter):
+            z_next, c, residues = self.vf70_step_once(s, H_list, z, d0=d0)
+            if verbose:
+                rel = np.linalg.norm(z_next) / max(1.0, np.linalg.norm(z))
+                print(f"[VF-70] iter {it+1:02d}/{n_iter:02d}  Δz_rel = {rel}")
+            z = z_next
+        return z, c, residues
+
+    def eval_Hfit_from_70(self,s, z_prev, c, residues_m, d0=1.0):
+        """
+        Evaluate H_fit for one response m using (70):
+        H_fit = (Psi @ residues_m) / (d0 + Psi @ c)
+        """
+        Psi = 1.0 / (s[:, None] + np.asarray(z_prev)[None, :])
+        D_ratio = d0 + Psi @ c
+        N_part  = Psi @ residues_m
+        return N_part / D_ratio
+
+    def rms_abs(self,y):
+        return np.sqrt(np.mean(np.abs(y)**2))
+
+    def rms_error(self,H, Hfit):
+        """
+        Returns (abs_rms, relative_rms). Relative RMS is normalized by RMS of H.
+        """
+        err = Hfit - H
+        abs_rms = self.rms_abs(err)
+        ref = max(self.rms_abs(H), 1e-16)
+        rel_rms = abs_rms / ref
+        return abs_rms, rel_rms
+
+    # --- run K (== n_iter) VF-(70) steps and record history ---
+    def run_vf70_with_history(self,s, H, z0, K=25, d0=1.0, verbose=True):
+        """
+        H: 1D complex array (single response). If you want multi-response,
+        pass a list [H11, H21, ...] and adapt the evaluation below.
+        """
+        # Wrap single response as a list for the step function
+        H_list = [np.asarray(H, dtype=np.complex128)]
+        z = np.array(z0, dtype=np.complex128)
+        hist = {
+            "z": [z.copy()],
+            "c": [],
+            "r": [],
+            "Hfit": [],
+            "rms_abs": [],
+            "rms_rel": [],
+        }
+        for it in range(K):
+            z_new, c, residues = self.vf70_step_once(s, H_list, z, d0=d0, column_scale=True)
+            # Evaluate fitted response after this iteration
+            Hfit = self.eval_Hfit_from_70(s, z, c, residues[0], d0=d0)
+            abs_r, rel_r = self.rms_error(H, Hfit)
+
+            if verbose:
+                drel = np.linalg.norm(z_new - z) / max(1.0, np.linalg.norm(z))
+                print(f"[VF-70] iter {it+1:02d}/{K:02d}  Δz_rel={drel:.3e}  RMS={abs_r:.3e}  RMSrel={rel_r:.3e}")
+
+            # save
+            hist["c"].append(c.copy())
+            hist["r"].append(residues[0].copy())
+            hist["Hfit"].append(Hfit.copy())
+            hist["rms_abs"].append(abs_r)
+            hist["rms_rel"].append(rel_r)
+            z = z_new
+            hist["z"].append(z.copy())
+
+        return hist
+    
+    def evaluate(self, s_eval, z_ref, c, residues, d0=1.0):
+        """
+        Evaluate the (70) model on any frequency grid.
+
+        Parameters
+        ----------
+        s_eval : (K_eval,) complex
+            Points where to evaluate (e.g., 1j*2π*freq_eval).
+        z_ref  : (P,) complex
+            The *reference* poles used in the LS system that produced (c, residues).
+            For the last iteration in run_vf70_with_history, this is hist['z'][-2].
+        c      : (P,) complex
+            Denominator-update coefficients from (70).
+        residues : (P,) complex  OR  list of (P,) for multi-response
+            Residues for the response(s) from (70).
+        d0     : scalar, default 1.0
+
+        Returns
+        -------
+        H_eval : (K_eval,) complex  OR  (K_eval, M) for M responses
+        """
+        s_eval = np.asarray(s_eval, dtype=np.complex128).reshape(-1)
+        z_ref = np.asarray(z_ref, dtype=np.complex128).reshape(-1)
+        Psi_eval = 1.0 / (s_eval[:, None] + z_ref[None, :])  # K_eval x P
+        Den = d0 + Psi_eval @ np.asarray(c, dtype=np.complex128)
+
+        # Single response
+        if not isinstance(residues, (list, tuple)):
+            r = np.asarray(residues, dtype=np.complex128).reshape(-1)
+            Num = Psi_eval @ r
+            return Num / Den
+
+        # Multi-response
+        M = len(residues)
+        H_eval = np.empty((s_eval.size, M), dtype=np.complex128)
+        for m, r_m in enumerate(residues):
+            r_m = np.asarray(r_m, dtype=np.complex128).reshape(-1)
+            H_eval[:, m] = (Psi_eval @ r_m) / Den
+        return H_eval
+    
+def auto_select(H, freq,
+                n_baseline=64,          # log-spaced backbone points
+                peak_prominence=0.05,   # fraction of |H| dB dynamic range for peak detection
+                peak_window=5,          # take ±peak_window samples around each peak
+                topgrad_q=0.98,         # keep top 2% largest slope/phase-change points
+                max_points=25,         # final cap on selected samples (None = no cap)
+                ensure_ends=True):
+    """
+    Select several significant sample points for vector fitting.
+
+    Strategy:
+      1) Always keep endpoints (optional).
+      2) Add a log-spaced baseline over the band.
+      3) Detect resonance peaks in |H| (on a log scale) and keep small windows around them.
+      4) Add points with the largest magnitude slope and phase-change (w.r.t log-f).
+      5) De-duplicate, sort, and optionally thin to 'max_points' with priority
+         to endpoints and detected peaks.
+
+    Parameters
+    ----------
+    H : (N,) complex array
+        Frequency response samples.
+    freq : (N,) float array
+        Frequency axis [Hz], strictly increasing.
+    n_baseline : int
+        Count of log-spaced baseline samples across the band.
+    peak_prominence : float
+        Peak prominence threshold as a fraction of the dynamic range in log|H|.
+        0.05 ≈ keep peaks ≥ 5% of the range.
+    peak_window : int
+        Number of neighbor indices to include on each side of every detected peak.
+    topgrad_q : float in (0,1)
+        Quantile for selecting strong slope/phase points.
+        0.98 ⇒ keep the top 2% largest derivatives.
+    max_points : int or None
+        If not None, cap the total number of selected indices to this value.
+    ensure_ends : bool
+        Always include the first and last samples.
+
+    Returns
+    -------
+    H_sel : (K,) complex array
+    freq_sel : (K,) float array
+    """
+    H = np.asarray(H).reshape(-1)
+    f = np.asarray(freq).reshape(-1)
+    if H.size != f.size:
+        raise ValueError("H and freq must have the same length.")
+    N = f.size
+    if N < 4:
+        return H.copy(), f.copy()
+
+    eps = 1e-16
+    mag = np.abs(H)
+    logmag = np.log10(mag + eps)
+    phase = np.unwrap(np.angle(H))
+
+    # log-frequency axis (scale-invariant derivatives)
+    # keep it linear if any non-positive freq sneaks in
+    if np.all(f > 0):
+        lf = np.log(f)
+    else:
+        lf = f.copy()
+
+    dlf = np.gradient(lf)
+    d_logmag = np.gradient(logmag) / (dlf + 1e-16)
+    d_phase  = np.gradient(phase)  / (dlf + 1e-16)
+
+    idx = set()
+    if ensure_ends:
+        idx.update([0, N-1])
+
+    # 1) log-spaced baseline
+    if n_baseline > 0:
+        # map a log grid to nearest indices
+        grid = np.linspace(lf.min(), lf.max(), n_baseline)
+        base_idx = np.clip(np.searchsorted(lf, grid), 0, N-1)
+        idx.update(np.unique(base_idx).tolist())
+
+    # 2) peaks in |H|
+    try:
+        from scipy.signal import find_peaks
+        dyn = logmag.max() - logmag.min()
+        prom = peak_prominence * (dyn + 1e-12)
+        peaks, _ = find_peaks(logmag, prominence=prom)
+    except Exception:
+        # simple fallback: strict local maxima
+        peaks = np.where((mag[1:-1] > mag[:-2]) & (mag[1:-1] > mag[2:]))[0] + 1
+
+    for p in peaks:
+        lo = max(0, p - peak_window)
+        hi = min(N, p + peak_window + 1)
+        idx.update(range(lo, hi))
+
+    # 3) strongest slope / phase-change points
+    thr_slope = np.quantile(np.abs(d_logmag), topgrad_q)
+    thr_phase = np.quantile(np.abs(d_phase),  topgrad_q)
+    idx.update(np.where(np.abs(d_logmag) >= thr_slope)[0].tolist())
+    idx.update(np.where(np.abs(d_phase)  >= thr_phase)[0].tolist())
+
+    # 4) finalize set
+    sel = np.array(sorted(idx), dtype=int)
+
+    # 5) optional thinning with priority to endpoints and peaks
+    if max_points is not None and sel.size > max_points:
+        priority = np.zeros(sel.size, dtype=int)
+        if ensure_ends:
+            priority[(sel == 0) | (sel == N-1)] = 3
+        if peaks.size:
+            priority[np.isin(sel, peaks)] = np.maximum(priority[np.isin(sel, peaks)], 2)
+
+        keep = []
+        budget = max_points
+        # keep highest-priority first
+        for lev in (3, 2, 1, 0):
+            cand = sel[priority == lev]
+            if cand.size == 0:
+                continue
+            if cand.size <= budget:
+                keep.extend(cand.tolist())
+                budget -= cand.size
+            else:
+                step = max(1, int(np.ceil(cand.size / budget)))
+                keep.extend(cand[::step][:budget].tolist())
+                budget = 0
+            if budget == 0:
+                break
+        sel = np.array(sorted(set(keep)), dtype=int)
+
+    return H[sel], f[sel]
 
 
+if __name__ == "__main__":
+    # formula_67(s,y)
+    H11 = np.array([y[i,0,0] for i in range(len(y))])
+    H11_slice,freqs_slice = auto_select(H11,freqs,max_points=20)
+    s_slice = freqs_slice * 2j * np.pi
+    P_pairs = 2
+    z0 = generate_starting_poles(P_pairs, beta_min=1e8, beta_max=freqs_slice[-1])
+    z0 = np.array(z0, dtype=np.complex128)
+
+    f70 = formula_70()
+    K = 25
+    hist = f70.run_vf70_with_history(s_slice, H11_slice, z0, K=K, d0=1.0, verbose=True)
+
+    it_show = K-1
+    Hfit_final = hist["Hfit"][it_show]
+
+    s_dense = 1j * 2*np.pi * freqs
+    c_final = hist["c"][it_show]
+    r_final = hist["r"][it_show]          # single response
+    z_ref   = hist["z"][it_show]
+    Hfit_dense = f70.evaluate(s_dense, z_ref, c_final, r_final, d0=1.0)
+
+    fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=False)
+
+    # (1) Magnitude plot like your screenshot
+    ax0 = axes[0]
+    ax0.plot(freqs, np.abs(H11), 'o', ms=4, color='red', label='Samples')
+    ax0.plot(freqs, np.abs(Hfit_dense), '-', lw=2, color='k', label='Fit')
+    ax0.plot(freqs_slice, np.abs(H11_slice), 'x', ms=4, color='blue', label='Input Samples')
+    ax0.set_title("Response i=0, j=0")
+    ax0.set_ylabel("Magnitude")
+    ax0.legend(loc="best")
+
+    # (2) RMS error vs iteration
+    ax1 = axes[1]
+    its = np.arange(1, K+1)
+    ax1.plot(its, hist["rms_rel"], '-o', lw=2)
+    ax1.set_xlabel("Iteration")
+    ax1.set_ylabel("RMS error (relative)")
+    ax1.grid(True, alpha=0.3)
+    ax1.set_title(f"RMS(final) = {hist['rms_rel'][-1]:.3e}")
+
+    fig.tight_layout()
+    plt.savefig(f"img_formula_70.png")