diff --git a/src/OID_analysis/data_li22b/derived/pod/energy.json b/src/OID_analysis/data_li22b/derived/pod/energy.json
new file mode 100644
index 0000000..83dd1b8
--- /dev/null
+++ b/src/OID_analysis/data_li22b/derived/pod/energy.json
@@ -0,0 +1,27 @@
+{
+  "n_commands": 8,
+  "n_snapshots": 1600,
+  "dof": 240000,
+  "roi": [
+    800,
+    1400,
+    200,
+    400
+  ],
+  "r99": 4,
+  "e10": 0.999856945066893,
+  "e2": 0.8630403049012719,
+  "e1": 0.532196928160917,
+  "top10": [
+    0.532196928160917,
+    0.3308433767403549,
+    0.11744563908241336,
+    0.00979119184601673,
+    0.0034631638505219966,
+    0.002136456704568166,
+    0.0018127711252703537,
+    0.0011799071249095198,
+    0.000795472213865046,
+    0.0001920382180559607
+  ]
+}
\ No newline at end of file
diff --git a/src/OID_analysis/li22b/collect_openloop_db.py b/src/OID_analysis/li22b/collect_openloop_db.py
new file mode 100644
index 0000000..3362f71
--- /dev/null
+++ b/src/OID_analysis/li22b/collect_openloop_db.py
@@ -0,0 +1,119 @@
+# OID_analysis/li22b/collect_openloop_db.py
+"""Generate Li22b open-loop database: 50 steady control commands on fluidic pinball.
+
+Uses new CelerisLab Simulation API, 2000x600 uniform-inlet config (no disturbance cylinder).
+Saves ux/uy fields, 9 sensors (3x3 grid), and forces per command.
+
+Usage:
+  conda run -n pycuda_3_10 python3 src/OID_analysis/li22b/collect_openloop_db.py --device 2
+"""
+import os, sys, json, time, argparse
+import numpy as np
+
+_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
+if _REPO not in sys.path:
+    sys.path.insert(0, _REPO)
+
+from CelerisLab import Simulation
+
+U0 = 0.01; L0 = 20.0; R_CYL = 10.0
+NX, NY = 2000, 600
+CY = (NY - 1) / 2.0
+CFG_PATH = os.path.join(_REPO, "configs", "config_lbm_karman_2000x600.json")
+OUT_BASE = os.path.join(os.path.dirname(__file__), "..", "data_li22b")
+WARMUP_STEPS = 4000
+RECORD_STEPS = 200
+DEVICE = 2
+
+PINBALL_FRONT_X = 1000.0
+PINBALL_REAR_X  = 1026.0
+PINBALL_Y_SPAN  = 15.0
+SENSOR_XS = [PINBALL_REAR_X + dx for dx in [100.0, 130.0, 160.0]]
+SENSOR_YS = [CY + 25.0, CY, CY - 25.0]
+
+
+def latin_hypercube(n, d=3, low=-2.0, high=2.0, seed=42):
+    rng = np.random.RandomState(seed)
+    samples = np.zeros((n, d))
+    for j in range(d):
+        perm = rng.permutation(n)
+        samples[:, j] = low + (high - low) * (perm + rng.rand(n)) / n
+    return samples
+
+
+def run_one(cmd_id, b):
+    out_dir = os.path.join(OUT_BASE, f"{cmd_id:03d}")
+    fields_fp = os.path.join(out_dir, "fields.npz")
+    if os.path.isfile(fields_fp):
+        print(f"  [{cmd_id:03d}] SKIP (exists)")
+        return
+    os.makedirs(out_dir, exist_ok=True)
+
+    sim = Simulation(CFG_PATH, device_id=DEVICE)
+    sim.add_body("circle", center=(PINBALL_FRONT_X, CY, 0.0), radius=R_CYL)
+    sim.add_body("circle", center=(PINBALL_REAR_X,  CY + PINBALL_Y_SPAN, 0.0), radius=R_CYL)
+    sim.add_body("circle", center=(PINBALL_REAR_X,  CY - PINBALL_Y_SPAN, 0.0), radius=R_CYL)
+    for sx in SENSOR_XS:
+        for sy in SENSOR_YS:
+            sim.add_body("sensor", center=(sx, sy, 0.0), radius=5.0)
+    sim.initialize()
+
+    for i in range(3):
+        sim.set_body(i, omega=float(b[i]) * U0 / R_CYL)
+
+    sim.run(WARMUP_STEPS)
+
+    ux_list, uy_list = [], []
+    sens_list, force_list = [], []
+    for step in range(RECORD_STEPS):
+        sim.run(1)
+        macro = sim.get_macroscopic()
+        ux_list.append(macro["ux"].copy())
+        uy_list.append(macro["uy"].copy())
+        ss, fs = [], []
+        for si in range(9):
+            sv = sim.read_sensor(si, normalize=False)
+            ss.extend([float(sv[0]), float(sv[1])])
+        sens_list.append(ss)
+        for fi in range(3):
+            fv = sim.read_force(fi)
+            fs.extend([float(fv[0]), float(fv[1])])
+        force_list.append(fs)
+        if (step + 1) % 100 == 0:
+            print(f"    {step+1}/{RECORD_STEPS}")
+
+    sim.close()
+    np.savez_compressed(fields_fp,
+                        ux=np.array(ux_list, dtype=np.float32),
+                        uy=np.array(uy_list, dtype=np.float32))
+    np.savez_compressed(os.path.join(out_dir, "sensors.npz"),
+                        sensors=np.array(sens_list, dtype=np.float32))
+    np.savez_compressed(os.path.join(out_dir, "forces.npz"),
+                        forces=np.array(force_list, dtype=np.float32))
+    with open(os.path.join(out_dir, "config.json"), "w") as f:
+        json.dump({"b": b.tolist(), "nu": 0.004, "U0": U0,
+                   "Re_D": 50, "Re_code": 100}, f, indent=2)
+    print(f"  [{cmd_id:03d}] b={b[0]:+.3f},{b[1]:+.3f},{b[2]:+.3f}  saved")
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--device", type=int, default=2)
+    ap.add_argument("--seed", type=int, default=42)
+    ap.add_argument("--num", type=int, default=10)
+    ap.add_argument("--start", type=int, default=0)
+    args = ap.parse_args()
+    global DEVICE; DEVICE = args.device
+
+    b_all = latin_hypercube(args.num, seed=args.seed)
+    print(f"{args.num} commands, device={DEVICE}, seed={args.seed}")
+    print(f"Output: {OUT_BASE}")
+    t0 = time.time()
+    for i in range(args.start, args.num):
+        run_one(i, b_all[i])
+        print(f"  [ETA: {(time.time()-t0)*(args.num-i-1)/((i-args.start+1)*60):.0f} min]")
+    print(f"Done in {(time.time()-t0)/60:.1f} min")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/OID_analysis/li22b/collect_reference.py b/src/OID_analysis/li22b/collect_reference.py
new file mode 100644
index 0000000..5418c34
--- /dev/null
+++ b/src/OID_analysis/li22b/collect_reference.py
@@ -0,0 +1,82 @@
+# OID_analysis/li22b/collect_reference.py
+"""Collect reference fields: q_in (empty channel) and q_blk (uncontrolled pinball).
+
+Usage:
+  conda run -n pycuda_3_10 python3 src/OID_analysis/li22b/collect_reference.py --device 2
+"""
+import os, sys, json, time, argparse
+import numpy as np
+
+_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
+if _REPO not in sys.path:
+    sys.path.insert(0, _REPO)
+from CelerisLab import Simulation
+
+U0 = 0.01; R_CYL = 10.0
+NX, NY = 2000, 600
+CY = (NY - 1) / 2.0
+CFG_PATH = os.path.join(_REPO, "configs", "config_lbm_karman_2000x600.json")
+OUT_BASE = os.path.join(os.path.dirname(__file__), "..", "data_li22b")
+STEPS = 200
+
+
+def collect_empty_channel(device):
+    """q_in: no cylinders, just uniform flow."""
+    out_dir = os.path.join(OUT_BASE, "q_in")
+    if os.path.isfile(os.path.join(out_dir, "fields.npz")):
+        print("q_in: SKIP (exists)")
+        return
+    os.makedirs(out_dir, exist_ok=True)
+    sim = Simulation(CFG_PATH, device_id=device)
+    sim.initialize()
+    sim.run(4000)
+    ux_list, uy_list = [], []
+    for _ in range(STEPS):
+        sim.run(1)
+        macro = sim.get_macroscopic()
+        ux_list.append(macro["ux"].copy())
+        uy_list.append(macro["uy"].copy())
+    sim.close()
+    np.savez_compressed(os.path.join(out_dir, "fields.npz"),
+                        ux=np.array(ux_list, dtype=np.float32),
+                        uy=np.array(uy_list, dtype=np.float32))
+    print(f"q_in: saved {STEPS} snaps")
+
+
+def collect_uncontrolled_pinball(device):
+    """q_blk: pinball with b=[0,0,0]."""
+    out_dir = os.path.join(OUT_BASE, "q_blk")
+    if os.path.isfile(os.path.join(out_dir, "fields.npz")):
+        print("q_blk: SKIP (exists)")
+        return
+    os.makedirs(out_dir, exist_ok=True)
+    sim = Simulation(CFG_PATH, device_id=device)
+    sim.add_body("circle", center=(1000.0, CY, 0.0), radius=R_CYL)
+    sim.add_body("circle", center=(1026.0, CY + 15.0, 0.0), radius=R_CYL)
+    sim.add_body("circle", center=(1026.0, CY - 15.0, 0.0), radius=R_CYL)
+    sim.initialize()
+    sim.run(4000)
+    ux_list, uy_list = [], []
+    for _ in range(STEPS):
+        sim.run(1)
+        macro = sim.get_macroscopic()
+        ux_list.append(macro["ux"].copy())
+        uy_list.append(macro["uy"].copy())
+    sim.close()
+    np.savez_compressed(os.path.join(out_dir, "fields.npz"),
+                        ux=np.array(ux_list, dtype=np.float32),
+                        uy=np.array(uy_list, dtype=np.float32))
+    print(f"q_blk: saved {STEPS} snaps")
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--device", type=int, default=2)
+    args = ap.parse_args()
+    collect_empty_channel(args.device)
+    collect_uncontrolled_pinball(args.device)
+    print("Done.")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/OID_analysis/li22b/phase_a2_pod.py b/src/OID_analysis/li22b/phase_a2_pod.py
new file mode 100644
index 0000000..7bb5356
--- /dev/null
+++ b/src/OID_analysis/li22b/phase_a2_pod.py
@@ -0,0 +1,99 @@
+# OID_analysis/li22b/phase_a2_pod.py
+"""Phase A.2: Snapshot POD on Li22b open-loop database.
+
+Replicates Li22b 4.1:
+  - Load all commands, compute ensemble mean over b,t
+  - POD via method-of-snapshots on ROI (x=[800:1400], y=[200:400])
+  - Truncate to 99% energy
+  - Visualize mode shapes + physical interpretation
+
+Usage (after data collected):
+  PYTHONPATH="src:$PYTHONPATH" python3 src/OID_analysis/li22b/phase_a2_pod.py
+"""
+import os, sys, json, glob
+import numpy as np
+
+_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
+if _REPO not in sys.path:
+    sys.path.insert(0, _REPO)
+from OID_analysis.utils.analysis import compute_pod
+
+DATA_BASE = os.path.join(os.path.dirname(__file__), "..", "data_li22b")
+DERIVED = os.path.join(DATA_BASE, "derived", "pod")
+os.makedirs(DERIVED, exist_ok=True)
+
+NX, NY = 2000, 600
+# ROI: pinball wake region
+ROI_X0, ROI_X1 = 800, 1400
+ROI_Y0, ROI_Y1 = 200, 400
+NX_ROI = ROI_X1 - ROI_X0; NY_ROI = ROI_Y1 - ROI_Y0
+
+
+def load_all_fields():
+    """Load all ux, uy snapshots from all commands, plus b vectors."""
+    cmd_dirs = sorted(glob.glob(os.path.join(DATA_BASE, "[0-9][0-9][0-9]")))
+    all_ux, all_uy, all_b = [], [], []
+    for d in cmd_dirs:
+        fp = os.path.join(d, "fields.npz")
+        if not os.path.isfile(fp): continue
+        try:
+            data = np.load(fp)
+        except Exception as e:
+            print(f"  SKIP {os.path.basename(d)}: {e}")
+            continue
+        ux = data["ux"][:, ROI_Y0:ROI_Y1, ROI_X0:ROI_X1]
+        uy = data["uy"][:, ROI_Y0:ROI_Y1, ROI_X0:ROI_X1]
+        all_ux.append(ux); all_uy.append(uy)
+        with open(os.path.join(d, "config.json")) as f:
+            cfg = json.load(f)
+        all_b.append(cfg["b"])
+        print(f"  Loaded {os.path.basename(d)}: {ux.shape[0]} snaps, b={cfg['b']}")
+    return (np.concatenate(all_ux, axis=0), np.concatenate(all_uy, axis=0),
+            np.array(all_b), len(cmd_dirs))
+
+
+def build_snapshot_matrix(ux, uy):
+    N = ux.shape[0]
+    DOF = NX_ROI * NY_ROI * 2
+    Q = np.zeros((N, DOF), dtype=np.float64)
+    for t in range(N):
+        Q[t] = np.concatenate([ux[t].ravel(), uy[t].ravel()])
+    return Q
+
+
+def main():
+    ux_all, uy_all, b_all, n_cmds = load_all_fields()
+    N = ux_all.shape[0]
+    print(f"\nTotal: {n_cmds} commands, {N} snapshots, ROI={NX_ROI}x{NY_ROI}")
+
+    # POD
+    Q = build_snapshot_matrix(ux_all, uy_all)
+    print(f"Snapshot matrix: {Q.shape} ({Q.nbytes/1e9:.2f} GB)")
+    pod = compute_pod(Q, rank=min(100, N-1))
+    S = pod["S"]; energy = pod["energy"]; cum = pod["cum_energy"]
+    print(f"POD: S[0]={S[0]:.2e}, cum10={cum[9]:.4f}, cum2={cum[1]:.4f}")
+
+    # Truncation to 99%
+    r99 = int(np.searchsorted(cum, 0.99)) + 1
+    print(f"99% energy: {r99} modes (Li22b: 78 modes)")
+
+    # Save
+    np.savez_compressed(os.path.join(DERIVED, "pod_modes.npz"),
+                        modes=pod["modes"], mean=pod["mean"])
+    np.savez(os.path.join(DERIVED, "pod_coefs.npy"),
+             coefs=pod["coefs"], S=S, energy=energy, cum_energy=cum)
+    json.dump({"n_commands": n_cmds, "n_snapshots": N,
+               "dof": Q.shape[1], "roi": [ROI_X0,ROI_X1,ROI_Y0,ROI_Y1],
+               "r99": r99, "e10": float(cum[9]), "e2": float(cum[1]),
+               "e1": float(energy[0]), "top10": energy[:10].tolist()},
+              open(os.path.join(DERIVED, "energy.json"), "w"), indent=2)
+
+    # Comparison table
+    print(f"\n=== Energy Comparison ===")
+    print(f"Li22b: e1+2=44.9%, e1-10=78.9%, r99=78")
+    print(f"Ours:  e1+2={cum[1]*100:.1f}%, e1-10={cum[9]*100:.1f}%, r99={r99}")
+    print(f"\nResults saved to {DERIVED}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/OID_analysis/li22b/phase_a3_lse.py b/src/OID_analysis/li22b/phase_a3_lse.py
new file mode 100644
index 0000000..39f30e6
--- /dev/null
+++ b/src/OID_analysis/li22b/phase_a3_lse.py
@@ -0,0 +1,134 @@
+# OID_analysis/li22b/phase_a3_lse.py
+"""Phase A.3: LSE — [sensors, b] → POD coefficients.
+
+Li22b Eq 2.7-2.8: T_ij solves a_i = sum_j T_ij * q_j
+where q = [sensors (18 channels), b (3)] = 21 inputs.
+
+Usage:
+  PYTHONPATH="src:$PYTHONPATH" python3 src/OID_analysis/li22b/phase_a3_lse.py
+"""
+import os, sys, json, glob
+import numpy as np
+
+_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
+if _REPO not in sys.path:
+    sys.path.insert(0, _REPO)
+
+DATA_BASE = os.path.join(os.path.dirname(__file__), "..", "data_li22b")
+DERIVED = os.path.join(DATA_BASE, "derived")
+os.makedirs(os.path.join(DERIVED, "lse"), exist_ok=True)
+
+
+def load_all():
+    cmd_dirs = sorted(glob.glob(os.path.join(DATA_BASE, "[0-9][0-9][0-9]")))
+    all_sens, all_forces, all_b_snap, all_configs = [], [], [], []
+    for d in cmd_dirs:
+        sfp = os.path.join(d, "sensors.npz")
+        ffp = os.path.join(d, "forces.npz")
+        if not os.path.isfile(sfp): continue
+        sens = np.load(sfp)["sensors"]   # (200, 18)
+        force = np.load(ffp)["forces"]   # (200, 6)
+        all_sens.append(sens); all_forces.append(force)
+        with open(os.path.join(d, "config.json")) as f:
+            cfg = json.load(f)
+        all_configs.append(cfg)
+        # Repeat b for each time step
+        b_rep = np.tile(cfg["b"], (sens.shape[0], 1))  # (200, 3)
+        all_b_snap.append(b_rep)
+    S = np.concatenate(all_sens, axis=0)
+    B = np.concatenate(all_b_snap, axis=0)
+    N = S.shape[0]
+    print(f"Loaded {N} snapshots from {len(cmd_dirs)} commands")
+    return S, B, all_configs, len(cmd_dirs)
+
+
+def main():
+    S, B, configs, n_cmds = load_all()
+    N = S.shape[0]
+
+    # Load POD coefs (handles .npy or .npy.npz naming)
+    coef_base = os.path.join(DERIVED, "pod", "pod_coefs.npy")
+    for ext in [".npy", ".npy.npz", ".npz"]:
+        fp = coef_base + ext
+        if os.path.isfile(fp):
+            pod = np.load(fp)
+            break
+    else:
+        raise FileNotFoundError(f"No pod_coefs found")
+    A = pod["coefs"]  # (N, r)
+    r = A.shape[1]
+    print(f"POD coefs: {A.shape}, rank={r}")
+
+    # Build input matrix Q = [S, B] (21-dim)
+    # Standardize sensors per-channel
+    S_std = np.zeros_like(S)
+    S_mean, S_std_scale = [], []
+    for j in range(18):
+        m = np.mean(S[:, j]); s = np.std(S[:, j]) or 1.0
+        S_std[:, j] = (S[:, j] - m) / s
+        S_mean.append(m); S_std_scale.append(s)
+    Q = np.hstack([S_std, B])  # (N, 21)
+
+    # Train/test split by commands (80/20 per Li22b)
+    rng = np.random.RandomState(42)
+    cmd_indices = np.arange(n_cmds)
+    rng.shuffle(cmd_indices)
+    n_train = int(0.8 * n_cmds)
+    # Map back to snapshot indices
+    snaps_per_cmd = N // n_cmds
+    train_mask = np.zeros(N, dtype=bool)
+    for ci in cmd_indices[:n_train]:
+        train_mask[ci*snaps_per_cmd:(ci+1)*snaps_per_cmd] = True
+    test_mask = ~train_mask
+
+    # Solve LSE: T = (Q^T Q)^{-1} Q^T A
+    Q_train = Q[train_mask]; A_train = A[train_mask]
+    QTQ = Q_train.T @ Q_train
+    T = np.linalg.solve(QTQ + 1e-6 * np.eye(21), Q_train.T @ A_train)  # (21, r)
+    print(f"T matrix shape: {T.shape}")
+
+    # Predict test set
+    Q_test = Q[test_mask]; A_test = A[test_mask]
+    A_pred = Q_test @ T
+
+    # Error metrics (Li22b Eq 3.3, 3.5)
+    # Per-command error ε_a(b)
+    eps_a = []
+    for ci in cmd_indices[n_train:]:
+        i0 = ci * snaps_per_cmd; i1 = (ci + 1) * snaps_per_cmd
+        mask = np.zeros(N, dtype=bool); mask[i0:i1] = True
+        mask_test = mask[test_mask]
+        if np.sum(mask_test) == 0: continue
+        a_true = A_test[mask_test]; a_pred = A_pred[mask_test]
+        num = np.mean(np.sum((a_true - a_pred)**2, axis=1))
+        den = np.mean(np.sum(a_true**2, axis=1))
+        eps_a.append(float(np.sqrt(num / den)) if den > 1e-30 else 0.0)
+    E = float(np.sqrt(np.mean(np.array(eps_a)**2)))
+
+    print(f"\n=== LSE Results ===")
+    print(f"Train commands: {n_train}, Test: {n_cmds - n_train}")
+    print(f"Mean ε_a: {np.mean(eps_a):.4f} ± {np.std(eps_a):.4f}")
+    print(f"Overall E: {E:.4f} (Li22b ~0.15-0.25 for periodic)")
+    print(f"Best ε_a: {np.min(eps_a):.4f}, Worst: {np.max(eps_a):.4f}")
+
+    # Per-mode estimation error
+    per_mode_mse = np.mean((A_test - A_pred)**2, axis=0)
+    per_mode_energy = np.mean(A_test**2, axis=0)
+    per_mode_err = np.sqrt(per_mode_mse / (per_mode_energy + 1e-30))
+    print(f"\nTop 10 mode errors: {per_mode_err[:10].round(4).tolist()}")
+
+    # Save
+    out_dir = os.path.join(DERIVED, "lse")
+    np.savez(os.path.join(out_dir, "lse_T.npz"), T=T,
+             S_mean=S_mean, S_std=S_std_scale)
+    json.dump({"n_train": n_train, "n_test": n_cmds - n_train,
+               "E": E, "mean_eps_a": float(np.mean(eps_a)),
+               "std_eps_a": float(np.std(eps_a)),
+               "eps_a_all": eps_a,
+               "per_mode_err": per_mode_err[:20].tolist()},
+              open(os.path.join(out_dir, "lse_results.json"), "w"), indent=2)
+    print(f"\nSaved to {out_dir}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/OID_analysis/li22b/phase_b.py b/src/OID_analysis/li22b/phase_b.py
new file mode 100644
index 0000000..1eacaef
--- /dev/null
+++ b/src/OID_analysis/li22b/phase_b.py
@@ -0,0 +1,139 @@
+# OID_analysis/li22b/phase_b.py
+"""Phase B: OID on Li22b DB + Cross-mapping + Joint-input OID.
+
+B.1: Correction-field OID: Delta_q_ctl(b) = q_ctl(b) - q_blk,
+     Force-OID with Y = force per command.
+B.2: Li22b POD modes (phase A.2) vs OID modes cross-mapping via subspace overlap.
+B.3: Joint-input OID: Y = [forces, sensors] concatenated.
+
+Usage:
+  PYTHONPATH="src:$PYTHONPATH" python3 src/OID_analysis/li22b/phase_b.py
+"""
+import os, sys, json, glob
+import numpy as np
+
+_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
+if _REPO not in sys.path:
+    sys.path.insert(0, _REPO)
+from OID_analysis.utils.analysis import compute_pod, compute_force_oid, standardize, reconstruct_oid_modes
+
+DATA_BASE = os.path.join(os.path.dirname(__file__), "..", "data_li22b")
+DERIVED = os.path.join(DATA_BASE, "derived")
+
+ROI_X0, ROI_X1 = 800, 1400
+ROI_Y0, ROI_Y1 = 200, 400
+NY_ROI = ROI_Y1 - ROI_Y0; NX_ROI = ROI_X1 - ROI_X0
+
+
+def load_qblk():
+    fp = os.path.join(DATA_BASE, "q_blk", "fields.npz")
+    d = np.load(fp)
+    return d["ux"][:, ROI_Y0:ROI_Y1, ROI_X0:ROI_X1], d["uy"][:, ROI_Y0:ROI_Y1, ROI_X0:ROI_X1]
+
+
+def load_cmd(cmd_dir):
+    d = np.load(os.path.join(cmd_dir, "fields.npz"))
+    sens = np.load(os.path.join(cmd_dir, "sensors.npz"))["sensors"]
+    force = np.load(os.path.join(cmd_dir, "forces.npz"))["forces"]
+    with open(os.path.join(cmd_dir, "config.json")) as f:
+        b = json.load(f)["b"]
+    return (d["ux"][:, ROI_Y0:ROI_Y1, ROI_X0:ROI_X1],
+            d["uy"][:, ROI_Y0:ROI_Y1, ROI_X0:ROI_X1], sens, force, b)
+
+
+def snapshot_matrix(ux, uy):
+    N = ux.shape[0]; DOF = NY_ROI * NX_ROI * 2
+    Q = np.zeros((N, DOF), dtype=np.float64)
+    for t in range(N):
+        Q[t] = np.concatenate([ux[t].ravel(), uy[t].ravel()])
+    return Q
+
+
+def main():
+    out_dir = os.path.join(DERIVED, "oid_li22b")
+    os.makedirs(out_dir, exist_ok=True)
+
+    ux_blk, uy_blk = load_qblk()
+    print(f"q_blk: {ux_blk.shape}")
+
+    cmd_dirs = sorted(glob.glob(os.path.join(DATA_BASE, "[0-9][0-9][0-9]")))
+    valid_dirs = [d for d in cmd_dirs if os.path.isfile(os.path.join(d, "forces.npz"))]
+    print(f"Commands: {len(valid_dirs)}")
+
+    # --- B.1: Correction-field OID per command ---
+    all_delta = []; all_force = []; all_sens = []; all_b = []
+    for d in valid_dirs:
+        ux_ctl, uy_ctl, sens, force, b = load_cmd(d)
+        N = min(ux_ctl.shape[0], ux_blk.shape[0])
+        dux = ux_ctl[:N] - ux_blk[:N]; duy = uy_ctl[:N] - uy_blk[:N]
+        all_delta.append((dux, duy))
+        # Aggregate force: mean over time
+        all_force.append(np.mean(force[:N], axis=0))
+        all_sens.append(np.mean(sens[:N], axis=0))
+        all_b.append(b)
+
+    # Concatenate all Delta q_ctl
+    ux_all = np.concatenate([d[0] for d in all_delta], axis=0)
+    uy_all = np.concatenate([d[1] for d in all_delta], axis=0)
+    Q = snapshot_matrix(ux_all, uy_all)
+    print(f"Delta-q_ctl snapshot: {Q.shape}")
+
+    # POD on Delta-q_ctl
+    pod_li22b = compute_pod(Q, rank=min(20, Q.shape[0]-1))
+    print(f"Delta-q POD: cum5={pod_li22b['cum_energy'][4]:.4f}")
+
+    # Force-OID: Y = force per command, expanded to all snapshots
+    snaps_per_cmd = ux_all.shape[0] // len(valid_dirs)
+    Y_force = np.repeat(np.array(all_force), snaps_per_cmd, axis=0)[:ux_all.shape[0]]
+    Yf_std, _, _ = standardize(Y_force.astype(np.float64))
+    A_std, _, _ = standardize(pod_li22b["coefs"].astype(np.float64))
+    oid_force = compute_force_oid(A_std, Yf_std)
+    print(f"Force-OID on Li22b DB: S[0]={oid_force['S'][0]:.4f}")
+
+    # --- B.2: Cross-mapping Li22b POD ↔ OID modes ---
+    # Load Li22b full POD modes (from phase A.2)
+    li22b_pod_fp = os.path.join(DERIVED, "pod", "pod_modes.npz")
+    if os.path.isfile(li22b_pod_fp):
+        li22b_modes = np.load(li22b_pod_fp)["modes"]
+        r_min = min(li22b_modes.shape[1], pod_li22b["modes"].shape[1], 10)
+        # Subspace overlap: O(Phi_Li22b, Phi_OID)
+        overlap_matrix = np.zeros((r_min, r_min))
+        for i in range(r_min):
+            for j in range(r_min):
+                overlap_matrix[i,j] = abs(np.dot(li22b_modes[:,i], pod_li22b["modes"][:,j]))
+        print(f"\nLi22b POD x OID subspace overlap (first {r_min} modes):")
+        print("    " + "  ".join([f"OID{j}" for j in range(r_min)]))
+        for i in range(r_min):
+            row = "  ".join([f"{overlap_matrix[i,j]:.3f}" for j in range(r_min)])
+            print(f"Li{i:>2d} {row}")
+        np.savez(os.path.join(out_dir, "crossmap.npz"), overlap=overlap_matrix)
+        # Top matches
+        for i in range(min(3, r_min)):
+            top_j = int(np.argmax(overlap_matrix[i]))
+            print(f"  Li22b mode {i} best matches OID mode {top_j} (overlap={overlap_matrix[i,top_j]:.3f})")
+
+    # --- B.3: Joint-input OID ---
+    Y_sens = np.repeat(np.array(all_sens), snaps_per_cmd, axis=0)[:ux_all.shape[0]]
+    Ys_std, _, _ = standardize(Y_sens.astype(np.float64))
+    Y_joint = np.hstack([Yf_std, Ys_std])
+    oid_joint = compute_force_oid(A_std, Y_joint)
+    print(f"\nJoint-OID: S[0]={oid_joint['S'][0]:.4f}, cum2={oid_joint['cum_corr'][1]:.4f}")
+
+    # Compare: does joint combine both channels?
+    print("\n=== Comparison ===")
+    print(f"Force-OID S[0]: {oid_force['S'][0]:.4f}")
+    print(f"Joint-OID S[0]: {oid_joint['S'][0]:.4f}")
+
+    # Save
+    np.savez(os.path.join(out_dir, "oid_force.npz"),
+             U=oid_force["U"], S=oid_force["S"], z=oid_force["z"])
+    np.savez(os.path.join(out_dir, "oid_joint.npz"),
+             U=oid_joint["U"], S=oid_joint["S"], z=oid_joint["z"])
+    np.savez(os.path.join(out_dir, "pod_delta_q.npz"),
+             modes=pod_li22b["modes"], mean=pod_li22b["mean"],
+             energy=pod_li22b["energy"])
+    print(f"\nSaved to {out_dir}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/OID_analysis/li22b/phase_c_synthesis.py b/src/OID_analysis/li22b/phase_c_synthesis.py
new file mode 100644
index 0000000..9ec6b91
--- /dev/null
+++ b/src/OID_analysis/li22b/phase_c_synthesis.py
@@ -0,0 +1,139 @@
+# OID_analysis/li22b/phase_c_synthesis.py
+"""Phase C: Three-framework synthesis — SR, Li22b, OID.
+
+C.1: Unified POD (Li22b DB + PPO DB global POD)
+C.2: Cross-framework comparison table + interpretation
+
+Usage:
+  PYTHONPATH="src:$PYTHONPATH" python3 src/OID_analysis/li22b/phase_c_synthesis.py
+"""
+import os, sys, json, glob, numpy as np
+
+_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
+if _REPO not in sys.path:
+    sys.path.insert(0, _REPO)
+
+DATA_BASE = os.path.join(os.path.dirname(__file__), "..", "data_li22b")
+DERIVED = os.path.join(DATA_BASE, "derived")
+OID_DERIVED = os.path.join(os.path.dirname(__file__), "..", "data", "derived")
+
+
+def load_li22b_lse():
+    fp = os.path.join(DERIVED, "lse", "lse_results.json")
+    if not os.path.isfile(fp): return None
+    return json.load(open(fp))
+
+
+def load_oid_results():
+    """Load existing OID results from old PPO data."""
+    fp = os.path.join(OID_DERIVED, "master", "master_table.json")
+    if not os.path.isfile(fp): return None
+    return json.load(open(fp))
+
+
+def load_li22b_oid():
+    fp = os.path.join(DERIVED, "oid_li22b", "crossmap.npz")
+    if not os.path.isfile(fp): return None
+    return np.load(fp)["overlap"]
+
+
+def load_li22b_pod_energy():
+    fp = os.path.join(DERIVED, "pod", "energy.json")
+    if not os.path.isfile(fp): return None
+    return json.load(open(fp))
+
+
+def main():
+    lse = load_li22b_lse()
+    oid_res = load_oid_results()
+    crossmap = load_li22b_oid()
+    pod_en = load_li22b_pod_energy()
+
+    print("=" * 70)
+    print("Three-Framework Synthesis: SR + Li22b + OID")
+    print("=" * 70)
+
+    print("\n--- 1. Li22b LSE ---")
+    if lse:
+        print(f"  Overall error E: {lse['E']:.4f}")
+        print(f"  Mean per-command eps_a: {lse['mean_eps_a']:.4f} +/- {lse['std_eps_a']:.4f}")
+    else:
+        print("  (not yet computed — run phase_a3_lse.py)")
+
+    print("\n--- 2. Li22b POD Energy ---")
+    if pod_en:
+        print(f"  Commands: {pod_en['n_commands']}, Snapshots: {pod_en['n_snapshots']}")
+        print(f"  E2={pod_en['e2']*100:.1f}%, E10={pod_en['e10']*100:.1f}%")
+        print(f"  99% truncation: {pod_en['r99']} modes")
+        print(f"  Li22b reference: E2=44.9%, E10=78.9%, r99=78")
+
+    print("\n--- 3. Li22b OID ---")
+    if crossmap is not None:
+        r = crossmap.shape[0]
+        # Top cross-mappings
+        print(f"  Mode cross-mapping ({r}x{r}):")
+        for i in range(min(5, r)):
+            top_j = int(np.argmax(crossmap[i]))
+            print(f"    Li22b mode {i} → OID mode {top_j} (overlap={crossmap[i,top_j]:.3f})")
+
+    print("\n--- 4. PPO OID Results (from previous work) ---")
+    if oid_res:
+        print(f"  (loaded from {OID_DERIVED}/master/master_table.json)")
+    else:
+        print("  (not found)")
+
+    print("\n--- 5. Synthesis ---")
+    interpretations = []
+
+    # Is Li22b POD energy distribution consistent with PPO POD?
+    if pod_en and pod_en['e10'] > 0.7:
+        interpretations.append(
+            "Li22b POD confirms low-dimensionality of pinball flows (E10~{:.0f}%). "
+            "Supports OID's default r=10 POD truncation.".format(pod_en['e10']*100))
+
+    # Does LSE error explain why SR works?
+    if lse and lse['E'] < 0.5:
+        interpretations.append(
+            "LSE achieves moderate estimation accuracy (E={:.3f}), suggesting "
+            "a linear component EXISTS in [s,b]→field mapping. This explains why "
+            "SR can find clean formulas: SR's obs→act chain operates on the same "
+            "sensor channels, and act is lower-dimensional than full field.".format(lse['E']))
+
+    # Cross-mapping interpretation
+    if crossmap is not None:
+        # Check if any Li22b mode strongly overlaps with OID modes
+        max_overlap = float(np.max(crossmap))
+        interpretations.append(
+            "Max Li22b-OID mode overlap = {:.3f}. ".format(max_overlap) +
+            ("The two POD bases share significant structure — steady-control "
+             "and dynamic-control correction fields lie in overlapping subspaces."
+             if max_overlap > 0.7 else
+             "Steady and dynamic control engage notably different POD structures, "
+             "confirming that control temporal dynamics fundamentally change the "
+             "correction space."))
+
+    # Key difference: Li22b uses steady open-loop, we use dynamic PPO
+    interpretations.append(
+        "CRITICAL DISTINCTION: Li22b estimates FULL fields from [s,b] under STEADY "
+        "controls. Our OID diagnoses CORRECTION structures in Δq_ctl under DYNAMIC "
+        "PPO. The bridge between them quantifies how much of the dynamic control "
+        "maneuver can be captured by the steady-control POD basis.")
+
+    for it in interpretations:
+        print(f"  • {it}")
+
+    # Save synthesis
+    out = {
+        "li22b_lse_E": lse['E'] if lse else None,
+        "li22b_pod_e10": pod_en['e10'] if pod_en else None,
+        "li22b_pod_r99": pod_en['r99'] if pod_en else None,
+        "max_crossmap_overlap": float(np.max(crossmap)) if crossmap is not None else None,
+        "interpretations": interpretations,
+    }
+    with open(os.path.join(DERIVED, "synthesis.json"), "w") as f:
+        json.dump(out, f, indent=2)
+    print(f"\nSaved synthesis to {DERIVED}/synthesis.json")
+
+
+if __name__ == "__main__":
+    main()