qcode_example_v6_enhanced.py

"""
QCoDeS Lab Tools - SQUID CPR Analysis Script - Refactored Version
===================================================================

This script implements a staged fitting approach for SQUID CPR analysis.

The workflow includes:
1.  Defining a modular SQUIDAnalyzer class for core physics and fitting.
2.  Defining helper functions for data loading, preprocessing, visualization, and interpretation.
3.  Executing the analysis workflow in a clear, step-by-step main block:
    a. Loading a specific QCoDeS dataset.
    b. Preprocessing data: unit conversion and centering the magnetic field.
    c. Extracting critical current vs. magnetic field (CPR) data.
    d. Running a staged analysis:
        i.  Pre-fitting with a simple sine model.
        ii. Fitting with a single-transparency model.
        iii.Fitting with a dual-transparency model.
    e. Visualizing the results with comprehensive plots, including a dV/dI heatmap.
    f. Interpreting and printing the final physical parameters.
"""
import os
import json
from pathlib import Path
import qcodes as qc
from qcodes.dataset.data_set import DataSet
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
from scipy.signal import find_peaks, detrend
import holoviews as hv
from holoviews import opts
import panel as pn
from scipy.interpolate import griddata
import scipy.constants as sc
from scipy.fft import fft, fftfreq
from lmfit import minimize, Parameters, report_fit, Model
from lmfit.model import ModelResult
from typing import Dict, Tuple, Optional, Any, List
import itertools
from rich.progress import track, Progress
from joblib import Parallel, delayed
from matplotlib import lines
import time
import psutil # For system resource monitoring
import argparse

# Import qcodes_optimized for enhanced parameter table
# from qcodes_optimized import ParameterTableBuilder

# --- MCMC functionality removed ---
# JAX and BlackJAX setup has been removed to simplify the codebase
# MCMC sampling capabilities have been disabled

def parse_run_id(run_id_str: str) -> list:
    """
    Parse RUN_ID parameter. Supports single ID (e.g., "500") or range (e.g., "500-511").
    Returns a list of integers.
    """
    if '-' in run_id_str:
        start, end = map(int, run_id_str.split('-'))
        return list(range(start, end + 1))
    else:
        return [int(run_id_str)]

def setup_output_directories(base_dir: Path, run_ids: list) -> dict:
    """
    Set up output directories based on RUN_ID configuration.
    Returns a dictionary with output paths for different file types.
    """
    if len(run_ids) == 1:
        # Single RUN_ID: fit_500
        output_base = base_dir / f"fit_{run_ids[0]}"
    else:
        # Multiple RUN_IDs: fit_500_511
        output_base = base_dir / f"fit_{run_ids[0]}_{run_ids[-1]}"

    # Create subdirectories
    subdirs = {
        'plot': output_base / 'plot',
        'csv': output_base / 'csv',
        'json': output_base / 'json',
        'html': output_base / 'html'
    }

    for subdir in subdirs.values():
        subdir.mkdir(parents=True, exist_ok=True)

    return {
        'base': output_base,
        'plot': subdirs['plot'],
        'csv': subdirs['csv'],
        'json': subdirs['json'],
        'html': subdirs['html']
    }

# --- 1. SQUID Analyzer Class ---

class SQUIDAnalyzer:
    """
    Encapsulates the physical models and fitting procedures for SQUID CPR analysis.
    """
    def __init__(self):
        self.PHI0 = sc.h / (2 * sc.e)
        self.scaling = {
            'current_factor': 1.0, 'field_factor': 1.0, 'area_factor': 1e12,
            'current_unit': 'μA', 'field_unit': 'μT', 'area_unit': 'μm²',
            'original_current_unit': 'μA', 'original_field_unit': 'μT', 'original_area_unit': 'm²'
        }
        self.results = {}  # Store fitting results
        self.detrend_info = {}  # Store detrend information

    # --- Physical Models ---

    @staticmethod
    def _cpr_single_junction(phi: np.ndarray, Ic: float, T: float) -> np.ndarray:
        """Current-phase relation for a single Josephson junction."""
        denominator = np.sqrt(1 - T * np.sin(phi / 2)**2)
        return Ic * np.sin(phi) / np.where(denominator < 1e-10, 1e-10, denominator)

    def _flux_phase(self, B_centered_uT: np.ndarray, B_0_uT: float, A_eff_um2: float) -> np.ndarray:
        """Calculates the magnetic flux phase."""
        B_T = (B_centered_uT - B_0_uT) * 1e-6
        A_eff_m2 = A_eff_um2 / self.scaling['area_factor']
        return 2 * np.pi * A_eff_m2 * B_T / self.PHI0

    def pre_fit_model(self, B_centered_uT, Ic0, r, A_eff_um2, B_0_uT, Ic1, **kwargs):
        """Simple sine model: Ic = Ic0 + r * φ + Ic1 * sin(φ) where φ = 2π * A_eff * (B - B_0) / PHI0."""
        flux_phase = self._flux_phase(B_centered_uT, B_0_uT, A_eff_um2)
        return Ic0 + r * flux_phase + Ic1 * np.sin(flux_phase)

    def single_junction_model(self, B_centered_uT, Ic0, r, A_eff_um2, B_0_uT, Ic1, **kwargs):
        """Simple sine model: Ic = Ic0 + r * φ + Ic1 * sin(φ) (same as pre_fit_model)."""
        flux_phase = self._flux_phase(B_centered_uT, B_0_uT, A_eff_um2)
        return Ic0 + r * flux_phase + Ic1 * np.sin(flux_phase)

    def single_T_model(self, B_centered_uT, Ic0, r, A_eff_um2, B_0_uT, Ic1, **kwargs):
        """Simple sine model: Ic = Ic0 + r * φ + Ic1 * sin(φ) (same as pre_fit_model)."""
        flux_phase = self._flux_phase(B_centered_uT, B_0_uT, A_eff_um2)
        return Ic0 + r * flux_phase + Ic1 * np.sin(flux_phase)

    def dual_T_model(self, B_centered_uT, Ic0, r, A_eff_um2, B_0_uT, Ic1, **kwargs):
        """Simple sine model: Ic = Ic0 + r * φ + Ic1 * sin(φ) (same as pre_fit_model)."""
        flux_phase = self._flux_phase(B_centered_uT, B_0_uT, A_eff_um2)
        return Ic0 + r * flux_phase + Ic1 * np.sin(flux_phase)

    # --- JAX methods removed ---
    # JAX-compatible methods have been removed since MCMC functionality was disabled

    @staticmethod
    def _set_param_bounds(params, name, center_val, perc_bounds):
        """Helper to set parameter bounds based on a center value and percentage."""
        delta = abs(center_val * perc_bounds)
        params[name].set(min=center_val - delta, max=center_val + delta)

    @staticmethod
    def _set_ic_param_bounds(params, name, center_val, current_type, branch):
        """
        Set bounds for Ic parameters based on current_type to enforce physical constraints.
        
        For Ic+ current_type: min=0 (positive currents)
        For Ic- current_type: max=0 (negative currents)  
        For Ic+ w/ Ic- current_type: branch-specific bounds (pos uses Ic+ rules, neg uses Ic- rules)
        For mixed types: use percentage-based bounds
        """
        if current_type == "Ic+":
            # Positive currents only
            params[name].set(min=0 + 1e-9, max=center_val * 2.0)
        elif current_type == "Ic-":
            # Negative currents only
            params[name].set(min=center_val * 2.0, max=0 - 1e-9)
        elif current_type == "Ic+ w/ Ic-":
            # Branch-specific bounds for mixed current type
            if branch == 'pos':
                # Positive branch uses Ic+ rules
                params[name].set(min=0 + 1e-9, max=center_val * 2.0)
            elif branch == 'neg':
                # Negative branch uses Ic- rules
                params[name].set(min=center_val * 2.0, max=0 - 1e-9)
            else:
                # Fallback to percentage bounds
                delta = abs(center_val * 0.5)
                params[name].set(min=center_val - delta + 1e-9, max=center_val + delta - 1e-9)
        else:
            # Mixed or unknown type - use percentage bounds
            delta = abs(center_val * 0.5)
            params[name].set(min=center_val - delta + 1e-9, max=center_val + delta - 1e-9)

    def reconstruct_original_data_and_fit(self, x_data: np.ndarray, y_data_detrended: np.ndarray, 
                                         x_fine: np.ndarray, branch: str, fit_stage: str):
        """
        Reconstruct original data and fit curves when detrend was applied.
        
        Parameters:
        - x_data: B field coordinates (μT)
        - y_data_detrended: Current data after detrend was applied  
        - x_fine: Fine B field array for smooth curve plotting
        - branch: Branch name ('pos' or 'neg')
        - fit_stage: Fitting stage ('pre_fit', 'single_T', 'dual_T')
        
        Returns:
        - x_plot: X coordinates for plotting
        - y_plot: Y data with trend added back for plotting
        - fit_curve: Fit curve with trend added back
        """
        if not hasattr(self, 'detrend_info') or not self.detrend_info.get('detrend_applied', False):
            # No detrend was applied, return data as is
            # Need to get fit curve from results
            from typing import Dict, Any
            if hasattr(self, 'results'):
                fit_result = self.results.get(f'{fit_stage}_{branch}', None)
                if fit_result is not None:
                    if fit_stage == 'pre_fit':
                        fit_curve = self.pre_fit_model(x_fine, **fit_result.params.valuesdict())
                    elif fit_stage == 'single_T':
                        fit_curve = self.single_T_model(x_fine, **fit_result.params.valuesdict())
                    elif fit_stage == 'dual_T':
                        fit_curve = self.dual_T_model(x_fine, **fit_result.params.valuesdict())
                    else:
                        fit_curve = np.zeros_like(x_fine)
                else:
                    fit_curve = np.zeros_like(x_fine)
            else:
                fit_curve = np.zeros_like(x_fine)
            return x_data, y_data_detrended, fit_curve
            
        # Detrend was applied, reconstruct original data
        if branch not in self.detrend_info:
            # No detrend info for this branch
            return x_data, y_data_detrended, np.zeros_like(x_fine)
            
        trend_slope = self.detrend_info[branch]['slope']
        trend_intercept = self.detrend_info[branch]['intercept']
        
        # Reconstruct original data by adding trend back
        y_original = y_data_detrended + trend_slope * x_data + trend_intercept
        
        # Get fit curve on detrended data
        if hasattr(self, 'results'):
            fit_result = self.results.get(f'{fit_stage}_{branch}', None)
            if fit_result is not None:
                if fit_stage == 'pre_fit':
                    fit_curve_detrended = self.pre_fit_model(x_fine, **fit_result.params.valuesdict())
                elif fit_stage == 'single_T':
                    fit_curve_detrended = self.single_T_model(x_fine, **fit_result.params.valuesdict())
                elif fit_stage == 'dual_T':
                    fit_curve_detrended = self.dual_T_model(x_fine, **fit_result.params.valuesdict())
                else:
                    fit_curve_detrended = np.zeros_like(x_fine)
                    
                # Add trend back to fit curve for visualization
                fit_curve_original = fit_curve_detrended + trend_slope * x_fine + trend_intercept
            else:
                fit_curve_original = np.zeros_like(x_fine)
        else:
            fit_curve_original = np.zeros_like(x_fine)
            
        return x_data, y_original, fit_curve_original
        
    # --- Fitting Execution ---

    def _run_fit(self, model_func, params, x_data, y_data, method='least_squares', fit_kws=None):
        """
        Generic fitting function using lmfit.Model.
        Uses 'soft_l1' loss for robustness against outliers by default when using least_squares.
        """
        model = Model(model_func, independent_vars=['B_centered_uT'])
        
        if fit_kws is None:
            fit_kws = {}
            
        # Set default loss function for least-squares methods
        if 'loss' not in fit_kws and method == 'least_squares':
            fit_kws['loss'] = 'soft_l1'
            
        result = model.fit(y_data, params, B_centered_uT=x_data, method=method, fit_kws=fit_kws)
        return result

    def _run_robust_fit_stage2(self, pre_fit_vals: Parameters, x_clean: np.ndarray, y_clean: np.ndarray, branch: str, current_type: str) -> ModelResult:
        """
        Performs a robust fit for Stage 2 by trying multiple initial values and fitting methods.
        Now fits for the ratio of Ic2/Ic1.
        Uses joblib for parallel processing.
        """
        print("\n[Stage 2] Starting robust fit for Single-Transparency Model...")
        print("   Model: Ic = Ic0 + r * φ + Ic1 * sin(φ)")
        print("   where φ = 2π * A_eff * (B - B_0) / PHI0")
        
        # Check if Ic1 is effectively zero (no modulation detected)
        if abs(pre_fit_vals['Ic1'].value) < 1e-10:
            print(f"   Warning: No modulation detected in {branch} branch (Ic1 ≈ 0). Skipping Stage 2 fit.")
            return None

        methods_to_try = ['least_squares', 'nelder', 'lbfgsb']
        t1_guesses = [0.5]
        phi_21_guesses = np.linspace(-np.pi, np.pi, 24, endpoint=False).tolist()

        best_result = None
        min_chi_sq = np.inf

        # --- Base parameters ---
        base_params = Parameters()
        # Freeze parameters from Stage 1
        base_params.add('Ic0', value=pre_fit_vals['Ic0'].value, vary=False)
        base_params.add('r', value=pre_fit_vals['r'].value, vary=False)
        base_params.add('B_0_uT', value=pre_fit_vals['B_0_uT'].value, vary=False)
        base_params.add('A_eff_um2', value=pre_fit_vals['A_eff_um2'].value, vary=False)
        
        # Sine parameter from Stage 1 (will be varied in the final refinement step)
        base_params.add('Ic1', value=pre_fit_vals['Ic1'].value, vary=False)
        # ---

        search_space = list(itertools.product(methods_to_try, t1_guesses, phi_21_guesses))
        
        print(f"Starting parallel grid search over {len(search_space)} combinations for {branch} branch...")

        def fit_worker(method, t1_init, phi_21_init):
            current_params = base_params.copy()
            current_params.add('T1', value=t1_init, min=0.25, max=0.75)
            current_params.add('phi_21', value=phi_21_init, min=-np.pi, max=np.pi)
            # Skip additional parameters since we now use simplified CPR model
            # All stages use the same basic model: Ic = Ic0 + r * CPR
            try:
                # First pass with frozen global parameters
                result = self._run_fit(self.single_T_model, current_params, x_clean, y_clean, method=method)
                if result.success:
                    return result
            except Exception:
                return None
            return None

        results = Parallel(n_jobs=-1)(delayed(fit_worker)(m, t1, p21) for m, t1, p21 in track(search_space, description="Fitting..."))

        # Filter out None results and find the best one
        successful_results = [res for res in results if res is not None]
        
        if not successful_results:
            print(f"   Warning: Stage 2 robust fitting failed to find any valid result for {branch} branch.")
            return None

        best_result = min(successful_results, key=lambda r: r.chisqr)
        min_chi_sq = best_result.chisqr

        print("\nRobust fit for Stage 2 complete.")
        if best_result is None:
            print(f"   Warning: Stage 2 robust fitting failed to find any valid result for {branch} branch.")
            return None
            
        # --- Second pass: Refine the best fit by unfreezing global parameters ---
        print("Refining best Stage 2 fit with global parameters (A_eff, B_0, Ic0) varying...")
        final_params = best_result.params.copy()
        
        # Unfreeze and set bounds for global parameters
        final_params['Ic0'].set(vary=True)
        self._set_ic_param_bounds(final_params, 'Ic0', pre_fit_vals['Ic0'].value, current_type, branch)
        
        # Unfreeze r parameter
        final_params['r'].set(vary=True, min=-10000, max=10000)
        
        final_params['B_0_uT'].set(vary=True)
        self._set_param_bounds(final_params, 'B_0_uT', pre_fit_vals['B_0_uT'].value, 0.05)

        final_params['A_eff_um2'].set(vary=True)
        self._set_param_bounds(final_params, 'A_eff_um2', pre_fit_vals['A_eff_um2'].value, 0.05)

        # Also unfreeze Ic1 for final refinement with current_type-aware bounds
        final_params['Ic1'].set(vary=True)
        self._set_ic_param_bounds(final_params, 'Ic1', pre_fit_vals['Ic1'].value, current_type, branch)

        final_result = self._run_fit(self.single_T_model, final_params, x_clean, y_clean, method='least_squares')
            
        return final_result

    # --- MCMC functionality removed ---
    # _run_blackjax_and_bands method has been removed

    # --- MCMC functionality removed ---
    # _run_emcee_and_bands method has been removed

    # --- MCMC functionality removed ---
    # _run_emcee_and_bands method has been removed

    def get_junction_components(self, fit_params: Dict[str, float], x_fine: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
        """
        Extract individual junction currents from dual-T model fit parameters.
        Returns (I1, I2) where I1 and I2 are the currents through each junction.
        """
        # Calculate flux phase
        flux_phase = self._flux_phase(x_fine, fit_params['B_0_uT'], fit_params['A_eff_um2'])
        
        # Calculate sine component (simple model)
        I_sine = fit_params['Ic1'] * np.sin(flux_phase)
        
        # Return same component twice for compatibility (all stages use same model now)
        return I_sine, I_sine

    def detect_linear_trend(self, branch_data, field_data):
        """檢測數據中的線性趨勢強度"""
        if len(branch_data) < 10 or len(field_data) < 10:
            return 0, 0, False
            
        # 計算線性擬合
        valid_mask = ~np.isnan(branch_data) & ~np.isnan(field_data)
        if np.sum(valid_mask) < 10:
            return 0, 0, False
            
        slope, intercept = np.polyfit(field_data[valid_mask], branch_data[valid_mask], 1)
        
        # 計算線性趨勢的強度 vs SQUID 振蕩振幅
        data_range = np.ptp(branch_data[valid_mask])
        oscillation_amplitude = np.std(branch_data[valid_mask])
        linear_trend_strength = abs(slope) * np.ptp(field_data[valid_mask])
        
        # 判斷是否需要 detrend（線性趨勢 > 振蕩振幅的 30%）
        needs_detrend = linear_trend_strength > 0.3 * oscillation_amplitude
        
        return slope, linear_trend_strength, needs_detrend
        
    def analyze_data_quality(self, ds_cpr):
        """分析數據質量並提供建議"""
        recommendations = []
        
        # 檢查 positive branch
        if 'crit_current_pos' in ds_cpr:
            y_pos = ds_cpr['crit_current_pos'].values
            x_data = ds_cpr.coords['B_y_centered_μT'].values
            
            valid_mask = np.isfinite(y_pos) & np.isfinite(x_data)
            if np.sum(valid_mask) > 10:
                slope_pos, trend_strength_pos, needs_detrend_pos = self.detect_linear_trend(
                    y_pos[valid_mask], x_data[valid_mask]
                )
                
                if needs_detrend_pos:
                    recommendations.append({
                        'branch': 'positive',
                        'slope': slope_pos,
                        'trend_strength': trend_strength_pos,
                        'message': f'檢測到強線性趨勢 (斜率: {slope_pos:.3f} μA/μT)，建議使用 --use-detrend 選項獲得更好的擬合結果'
                    })
                    
        # 檢查 negative branch  
        if 'crit_current_neg' in ds_cpr:
            y_neg = ds_cpr['crit_current_neg'].values
            x_data = ds_cpr.coords['B_y_centered_μT'].values
            
            valid_mask = np.isfinite(y_neg) & np.isfinite(x_data)
            if np.sum(valid_mask) > 10:
                slope_neg, trend_strength_neg, needs_detrend_neg = self.detect_linear_trend(
                    y_neg[valid_mask], x_data[valid_mask]
                )
                
                if needs_detrend_neg:
                    recommendations.append({
                        'branch': 'negative',
                        'slope': slope_neg, 
                        'trend_strength': trend_strength_neg,
                        'message': f'檢測到強線性趨勢 (斜率: {slope_neg:.3f} μA/μT)，建議使用 --use-detrend 選項獲得更好的擬合結果'
                    })
                    
        return recommendations
        
    def print_data_quality_analysis(self, recommendations):
        """輸出數據質量分析結果"""
        if recommendations:
            print("\n" + "="*60)
            print("📊 數據質量分析 & 建議")
            print("="*60)
            
            for rec in recommendations:
                print(f"⚠️  {rec['branch'].upper()} BRANCH: {rec['message']}")
            
            print("\n💡 使用 detrend 的優勢:")
            print("   • 移除測量系統的線性漂移")
            print("   • 突出 SQUID 干涉振蕩信號")
            print("   • 提高擬合穩定性和物理合理性")
            print("   • 減少參數間的相互干擾")
            print("="*60 + "\n")
        else:
            print("✅ 數據質量良好，無明顯線性趨勢\n")
    
    def estimate_oscillation_parameters(self, branch_data, field_data):
        """智能估計振蕩參數，用於改善非 detrend 模式的初始化"""
        if len(branch_data) < 20 or len(field_data) < 20:
            return None
            
        valid_mask = ~np.isnan(branch_data) & ~np.isnan(field_data)
        if np.sum(valid_mask) < 20:
            return None
            
        y_clean = branch_data[valid_mask]
        x_clean = field_data[valid_mask]
        
        # 臨時線性去趨勢用於振蕩參數估計
        from scipy.signal import detrend
        y_detrended = detrend(y_clean)
        
        # 估計振蕩振幅 (使用標準差的 2.5 倍作為峰峰值估計)
        oscillation_amplitude = np.std(y_detrended) * 2.5
        
        # 使用 FFT 估計主頻率和相位
        try:
            fft_vals = np.fft.fft(y_detrended)
            freqs = np.fft.fftfreq(len(y_detrended), np.mean(np.diff(x_clean)))
            
            # 找到主頻率 (忽略 DC 成分)
            valid_freq_mask = freqs > 0
            if np.any(valid_freq_mask):
                main_freq_idx = np.argmax(np.abs(fft_vals[valid_freq_mask])) 
                main_freq = freqs[valid_freq_mask][main_freq_idx]
                
                # 從主頻率估計 A_eff (基於 φ = 2π * A_eff * (B - B_0) / PHI0)
                PHI0 = sc.h / (2 * sc.e) * 1e6  # Flux quantum in μT⋅μm²
                estimated_A_eff = abs(main_freq) * PHI0 / (2 * np.pi) if main_freq != 0 else 350
                estimated_A_eff = np.clip(estimated_A_eff, 250, 500)  # 物理合理範圍
            else:
                estimated_A_eff = 350  # 默認值
                
        except:
            estimated_A_eff = 350  # 發生錯誤時的默認值
        
        return {
            'ic1_estimate': oscillation_amplitude,
            'a_eff_estimate': estimated_A_eff,
            'data_std': np.std(y_detrended)
        }

    def analyze(self, ds_cpr: xr.Dataset, output_dir: Path, current_type: str = "Ic+ w/ Ic-") -> Dict[str, Any]:
        """
        Performs the full staged fitting procedure for both positive and negative branches.
        current_type: "Ic+", "Ic+ w/ Ic-", or "Ic- w/ Ir+"
        """
        results = {}
        
        # 1. 分析數據質量
        quality_recommendations = self.analyze_data_quality(ds_cpr)
        self.print_data_quality_analysis(quality_recommendations)
        
        # Determine which branches to process based on current_type
        if current_type == "Ic+":
            branches_to_process = ['pos']
        elif current_type == "Ic+ w/ Ic-":
            branches_to_process = ['pos', 'neg']
        elif current_type == "Ic- w/ Ir+":
            branches_to_process = ['neg']  # Assuming neg corresponds to Ic-
        else:
            branches_to_process = ['pos', 'neg']  # Default fallback
        
        for branch in branches_to_process:
            print(f"\n--- Analyzing Branch: Ic{'+' if branch == 'pos' else '-'} ---")
            
            y_data = ds_cpr[f'crit_current_{branch}'].values
            x_data = ds_cpr.coords['B_y_centered_μT'].values
            mask = np.isfinite(y_data) & np.isfinite(x_data)
            x_clean, y_clean = x_data[mask], y_data[mask]
            
            # Check if we have any valid data points
            if len(x_clean) == 0:
                print(f"   Warning: No valid data points found for {branch} branch. Skipping analysis.")
                results[f'pre_fit_{branch}'] = None
                results[f'single_junction_fit_{branch}'] = None
                results[f'single_T_fit_{branch}'] = None
                results[f'dual_T_fit_{branch}'] = None
                continue

            # Stage 1: Pre-fitting with simple sine model
            print("\n[Stage 1] Pre-fitting with simple sine model...")
            print("   Model: Ic = Ic0 + r * φ + Ic1 * sin(φ)")
            print("   where φ = 2π * A_eff * (B - B_0) / PHI0")
            
            # 智能參數估計 (特別針對非 detrend 模式的改善)
            data_range = float(np.max(y_clean) - np.min(y_clean))
            is_detrended = data_range < 20.0
            has_strong_trend = len(quality_recommendations) > 0
            
            smart_params = None
            if not is_detrended and has_strong_trend:
                # 對原始數據使用智能參數估計
                smart_params = self.estimate_oscillation_parameters(y_clean, x_clean)
                if smart_params:
                    print(f"🧠 智能參數估計: Ic1≈{smart_params['ic1_estimate']:.2f}μA, A_eff≈{smart_params['a_eff_estimate']:.0f}μm²")
            
            pre_fit_params = Parameters()
            mean_ic, std_ic = np.mean(y_clean), np.std(y_clean)
            
            # Set Ic0 with current_type-aware bounds
            pre_fit_params.add('Ic0', value=mean_ic)
            self._set_ic_param_bounds(pre_fit_params, 'Ic0', mean_ic, current_type, branch)
            
            # Add r parameter (CPR scaling factor) - improved initialization
            # Estimate linear trend from data
            if len(x_clean) > 2:
                linear_slope, _ = np.polyfit(x_clean, y_clean, 1)
                r_initial = linear_slope * 100  # Scale factor for μT to μA conversion
                r_initial = np.clip(r_initial, -1.0, 1.0)  # Keep within reasonable bounds
            else:
                r_initial = 0.1  # Conservative default
            
            pre_fit_params.add('r', value=r_initial, min=-2.0, max=2.0)  # 允許更大的線性響應範圍
            
            # Set Ic1 with adaptive estimation based on data characteristics
            std_ic = float(np.std(y_clean))  # 數據標準差（總是正值）
            peak_to_peak = float(np.max(y_clean) - np.min(y_clean))  # 峰谷差值
            
            # 根据分支类型和数据特征设定 Ic1 参数
            if (current_type == 'Ic+') or (current_type == 'Ic+ w/ Ic-' and branch == 'pos'):
                if is_detrended:
                    # 去趨勢數據：使用更大的比例
                    ic1_value = 0.6 * peak_to_peak
                    ic1_min = 0.2 * peak_to_peak
                    ic1_max = 1.2 * peak_to_peak
                else:
                    # 原始數據：使用智能估計或保守設定
                    if smart_params and smart_params['ic1_estimate'] > 0:
                        # 使用智能估計，設定合理範圍防止過度壓縮
                        ic1_value = smart_params['ic1_estimate']
                        ic1_min = smart_params['ic1_estimate'] * 0.7  # 防止被壓縮到過小
                        ic1_max = smart_params['ic1_estimate'] * 1.5
                        print(f"   📌 使用智能 Ic1 估計: {ic1_value:.2f}μA (範圍: {ic1_min:.2f}-{ic1_max:.2f})")
                    else:
                        # 備用設定：更積極的初始值
                        ic1_value = max(0.15 * peak_to_peak, 2.5 * std_ic)  # 更大的初始估計
                        ic1_min = max(0.05 * peak_to_peak, 0.8 * std_ic)   # 設定合理下界
                        ic1_max = 0.4 * peak_to_peak
                        print(f"   📌 使用備用 Ic1 設定: {ic1_value:.2f}μA (範圍: {ic1_min:.2f}-{ic1_max:.2f})")
            elif (current_type == 'Ic-') or (current_type == 'Ic+ w/ Ic-' and branch == 'neg'):
                # 負分支：使用負值  
                ic1_value = -2.0 * std_ic
                ic1_min = -5.0 * std_ic
                ic1_max = -0.5 * std_ic
            else:
                # 備用邏輯
                ic1_value = 2.0 * std_ic
                ic1_min = 0.5 * std_ic
                ic1_max = 5.0 * std_ic
            
            pre_fit_params.add('Ic1', value=ic1_value, min=ic1_min, max=ic1_max)
            
            # Physical parameters for magnetic field to phase conversion
            # 根據數據特徵和智能估計調整 A_eff 參數
            if is_detrended:
                # 去趨勢數據：使用優化過的設置
                a_eff_value, a_eff_min, a_eff_max = 330, 300, 400
            else:
                # 原始數據：使用智能估計或合理設定
                if smart_params and 250 <= smart_params['a_eff_estimate'] <= 500:
                    # 使用智能估計
                    a_eff_value = smart_params['a_eff_estimate']
                    a_eff_min = max(280, a_eff_value - 50)
                    a_eff_max = min(450, a_eff_value + 50)
                    print(f"   📌 使用智能 A_eff 估計: {a_eff_value:.0f}μm² (範圍: {a_eff_min:.0f}-{a_eff_max:.0f})")
                else:
                    # 備用設定：基於成功 detrend 結果的經驗值
                    a_eff_value, a_eff_min, a_eff_max = 340, 300, 420
                    print(f"   📌 使用經驗 A_eff 設定: {a_eff_value:.0f}μm² (範圍: {a_eff_min:.0f}-{a_eff_max:.0f})")
            
            pre_fit_params.add('A_eff_um2', value=a_eff_value, min=a_eff_min, max=a_eff_max)
            
            # B_0 參數設置，根據數據特徵和物理合理性調整
            b_range = float(np.max(x_clean) - np.min(x_clean))
            if is_detrended:
                # 去趨勢數據：使用用戶判斷的 -1μT 附近
                b0_initial = -1.0
                b0_range = min(b_range * 0.15, 8.0)
                b0_min, b0_max = b0_initial - b0_range, b0_initial + b0_range
            else:
                # 原始數據：基於成功 detrend 結果，使用更聚焦的範圍
                if has_strong_trend:
                    # 有強線性趨勢時，使用接近 detrend 成功結果的設定
                    b0_initial = -1.2  # 稍微偏移以適應原始數據
                    b0_min, b0_max = -3.5, 0.5  # 更聚焦的範圍，基於 detrend 結果 -0.88μT
                    print(f"   📌 針對強線性趨勢優化 B_0: {b0_initial:.1f}μT (範圍: {b0_min:.1f}-{b0_max:.1f})")
                else:
                    # 無強線性趨勢時使用標準設定
                    b0_initial = -1.0
                    b0_range = min(b_range * 0.2, 12.0)
                    b0_min, b0_max = -b0_range, b0_range
            
            pre_fit_params.add('B_0_uT', value=b0_initial, min=b0_min, max=b0_max)
            
            pre_fit_result = self._run_fit(self.pre_fit_model, pre_fit_params, x_clean, y_clean)
            pre_fit_vals = pre_fit_result.params
            print("Pre-fit results:")
            report_fit(pre_fit_result)
            results[f'pre_fit_{branch}'] = pre_fit_result

            # Stage 1.2: Single Harmonic Transparency Model Fit
            # Check if Ic1 is effectively zero (no modulation detected)
            if abs(pre_fit_vals['Ic1'].value) < 1e-10:
                print(f"   Warning: No modulation detected in {branch} branch (Ic1 ≈ 0). Skipping single junction fit.")
                results[f'single_junction_fit_{branch}'] = None
                results[f'single_T_fit_{branch}'] = None
                results[f'dual_T_fit_{branch}'] = None
                # Skip to next stage or return early
                continue
            
            single_junction_params = Parameters()
            
            # Use pre-fit results as initial values with current_type-aware bounds
            single_junction_params.add('Ic0', value=pre_fit_vals['Ic0'].value)
            self._set_ic_param_bounds(single_junction_params, 'Ic0', pre_fit_vals['Ic0'].value, current_type, branch)
            
            # Add r parameter from pre-fit results (CPR scaling factor)
            single_junction_params.add('r', value=pre_fit_vals['r'].value, min=-2.0, max=2.0)
            
            # Physical parameters for magnetic field to phase conversion
            # 根據 pre-fit 結果和數據特徵設置範圍
            a_eff_value = pre_fit_vals['A_eff_um2'].value
            if a_eff_value < 350:  # 去趨勢數據的情況
                single_junction_params.add('A_eff_um2', value=a_eff_value, min=300, max=450)
            else:  # 原始數據的情況
                single_junction_params.add('A_eff_um2', value=a_eff_value, min=300, max=500)
            # 動態設置 B_0 範圍基於數據特徵
            b_range = float(np.max(x_clean) - np.min(x_clean))
            b0_range = min(b_range * 0.2, 10.0)
            single_junction_params.add('B_0_uT', value=pre_fit_vals['B_0_uT'].value, min=-b0_range, max=b0_range)
            
            # Sine modulation parameter from pre-fit results
            single_junction_params.add('Ic1', value=pre_fit_vals['Ic1'].value)
            self._set_ic_param_bounds(single_junction_params, 'Ic1', pre_fit_vals['Ic1'].value, current_type, branch)
            
            single_junction_result = self._run_fit(self.single_junction_model, single_junction_params, x_clean, y_clean)
            print("Single junction fit results:")
            report_fit(single_junction_result)
            results[f'single_junction_fit_{branch}'] = single_junction_result

            # after pre_fit_result
            x_fine = np.linspace(x_clean.min(), x_clean.max(), 600)
            # MCMC functionality removed - skipping MCMC sampling

            # Stage 2: Single-T Model Robust Fit
            single_T_result = self._run_robust_fit_stage2(pre_fit_vals, x_clean, y_clean, branch, current_type)
            if single_T_result is None:
                print(f"   Warning: Single-T fit failed for {branch} branch. Skipping dual-T fit.")
                results[f'single_T_fit_{branch}'] = None
                results[f'dual_T_fit_{branch}'] = None
                continue
            print("\nBest Single-T fit results from robust search:")
            report_fit(single_T_result)
            results[f'single_T_fit_{branch}'] = single_T_result

            # MCMC functionality removed - skipping MCMC sampling

            # Stage 3: Dual-T Model Fit
            print("\n[Stage 3] Fitting with Dual-Transparency Model...")
            print("   Model: Ic = Ic0 + r * φ + Ic1 * sin(φ)")
            
            if single_T_result is None:
                print(f"   Warning: Single-T fit failed for {branch} branch. Skipping dual-T fit.")
                results[f'dual_T_fit_{branch}'] = None
                continue
            
            # Temporarily skip Stage 3 due to parameter compatibility issues
            print("   Warning: Stage 3 temporarily disabled for simplified CPR model.")
            results[f'dual_T_fit_{branch}'] = single_T_result  # Use Stage 2 result as fallback
            continue
            
            # --- First pass: Freeze global parameters to stabilize fit ---
            print("   [Stage 3a] Freezing global parameters to fit T1, T2, Ic ratios...")
            dual_T_params_pass1 = Parameters()
            p2_vals = single_T_result.params

            # Freeze global parameters from Stage 2 results
            dual_T_params_pass1.add('Ic0', value=p2_vals['Ic0'].value, vary=False)
            dual_T_params_pass1.add('r', value=p2_vals['r'].value, vary=False)
            dual_T_params_pass1.add('B_0_uT', value=p2_vals['B_0_uT'].value, vary=False)
            dual_T_params_pass1.add('A_eff_um2', value=p2_vals['A_eff_um2'].value, vary=False)
            
            # Freeze CPR parameters from Stage 2 results
            dual_T_params_pass1.add('Ic_cpr', value=p2_vals['Ic_cpr'].value, vary=False)
            dual_T_params_pass1.add('T_cpr', value=p2_vals['T_cpr'].value, vary=False)

            result_pass1 = self._run_fit(self.dual_T_model, dual_T_params_pass1, x_clean, y_clean)
            print("   Fit results after Pass 1 (frozen globals):")
            report_fit(result_pass1)

            # --- Second pass: Unfreeze global parameters for final refinement ---
            print("\n   [Stage 3b] Unfreezing global parameters for final refinement...")
            dual_T_params_pass2 = result_pass1.params.copy()

            # Unfreeze and set bounds for global parameters
            dual_T_params_pass2['Ic0'].set(vary=True)
            self._set_ic_param_bounds(dual_T_params_pass2, 'Ic0', p2_vals['Ic0'].value, current_type, branch)
            
            # Unfreeze r parameter (CPR scaling factor)
            dual_T_params_pass2['r'].set(vary=True, min=-10.0, max=10.0)
            
            dual_T_params_pass2['B_0_uT'].set(vary=True)
            self._set_param_bounds(dual_T_params_pass2, 'B_0_uT', p2_vals['B_0_uT'].value, 0.1)

            dual_T_params_pass2['A_eff_um2'].set(vary=True)
            self._set_param_bounds(dual_T_params_pass2, 'A_eff_um2', p2_vals['A_eff_um2'].value, 0.1)
            
            # Unfreeze CPR parameters 
            dual_T_params_pass2['Ic_cpr'].set(vary=True)
            self._set_ic_param_bounds(dual_T_params_pass2, 'Ic_cpr', p2_vals['Ic_cpr'].value, current_type, branch)
            dual_T_params_pass2['T_cpr'].set(vary=True, min=0.1, max=0.9)

            dual_T_result = self._run_fit(self.dual_T_model, dual_T_params_pass2, x_clean, y_clean)
            print("\nFinal Dual-T fit results after refinement:")
            report_fit(dual_T_result)
            results[f'dual_T_fit_{branch}'] = dual_T_result

            # MCMC functionality removed - skipping MCMC sampling and benchmarks
            
        return results

# --- 2. Helper Functions ---

def load_squid_analyzer(script_path: str) -> type:
    """Loads the SQUIDAnalyzer class. In this standalone script, it just returns the class."""
    return SQUIDAnalyzer

def df_print_sci(df: pd.DataFrame, digits: int = 3):
    """Prints a DataFrame with scientific notation."""
    print(df.to_string(
        formatters={
            col: (lambda x: f"{x:.{digits}e}")
            for col, dtype in df.dtypes.items()
            if pd.api.types.is_float_dtype(dtype)
        }
    ))

def get_branches_for_current_type(current_type: str) -> List[str]:
    """根據 current_type 確定需要處理的分支"""
    branch_map = {
        "Ic+": ['pos'],
        "Ic+ w/ Ic-": ['pos', 'neg'],
        "Ic- w/ Ir+": ['neg'],
        "Ic-": ['neg']
    }
    return branch_map.get(current_type, ['pos', 'neg'])

def should_extract_branch(branch: str, current_type: str) -> bool:
    """判斷是否應該提取特定分支的資料"""
    required_branches = get_branches_for_current_type(current_type)
    return branch in required_branches

def prepare_cpr_data_from_qcodes(ds: DataSet, df_raw: pd.DataFrame, current_type: str = "Ic+ w/ Ic-", use_detrend: bool = False) -> Tuple[xr.Dataset, pd.DataFrame]:
    """Extracts CPR data (Ic vs. B) from raw IV curve data with current_type awareness."""
    print(f"\nExtracting CPR data for {current_type}...")
    
    # 根據 current_type 決定要提取的分支
    branches_to_extract = get_branches_for_current_type(current_type)
    
    # 初始化對應的分支列表
    crit_currents = {}
    for branch in ['pos', 'neg']:
        if should_extract_branch(branch, current_type):
            crit_currents[f'crit_current_{branch}'] = []
        else:
            # 對於不需要的分支，建立空列表以維持資料結構一致性
            crit_currents[f'crit_current_{branch}'] = []
    
    unique_B = df_raw['y_field_centered_µT'].unique()
    
    for B in unique_B:
        mask = (df_raw['y_field_centered_µT'] == B)
        i_vals = df_raw.loc[mask, 'appl_current_µA'].values
        v_vals = df_raw.loc[mask, 'meas_voltage_µV'].values
        
        sort_idx = np.argsort(i_vals)
        i_sorted, v_sorted = i_vals[sort_idx], v_vals[sort_idx]
        
        if len(i_sorted) > 1:
            di = np.diff(i_sorted)
            dvdi = np.abs(np.diff(v_sorted) / np.where(di == 0, 1e-9, di))
            i_mid = (i_sorted[:-1] + i_sorted[1:]) / 2

            # 處理所有分支，但只對需要提取的分支進行實際計算
            for branch in ['pos', 'neg']:
                if should_extract_branch(branch, current_type):
                    if branch == 'pos':
                        pos_mask = i_mid > 0
                        if np.any(pos_mask):
                            pos_peaks, _ = find_peaks(dvdi[pos_mask], height=np.max(dvdi[pos_mask]) * 0.5)
                            if len(pos_peaks) > 0:
                                crit_currents['crit_current_pos'].append(i_mid[pos_mask][pos_peaks[0]])
                            else:
                                # 如果沒有找到峰值，嘗試降低閾值
                                pos_peaks_relaxed, _ = find_peaks(dvdi[pos_mask], height=np.max(dvdi[pos_mask]) * 0.3)
                                if len(pos_peaks_relaxed) > 0:
                                    crit_currents['crit_current_pos'].append(i_mid[pos_mask][pos_peaks_relaxed[0]])
                                    print(f"   Debug: Found peak with relaxed threshold at B={B:.1f}µT, Ic={i_mid[pos_mask][pos_peaks_relaxed[0]]:.3f}µA")
                                else:
                                    crit_currents['crit_current_pos'].append(np.nan)
                                    print(f"   Debug: No peaks found in positive branch at B={B:.1f}µT")
                        else:
                            crit_currents['crit_current_pos'].append(np.nan)
                            print(f"   Debug: No positive current data at B={B:.1f}µT")
                            
                    elif branch == 'neg':
                        neg_mask = i_mid < 0
                        if np.any(neg_mask):
                            neg_peaks, _ = find_peaks(dvdi[neg_mask], height=np.max(dvdi[neg_mask]) * 0.5)
                            crit_currents['crit_current_neg'].append(
                                i_mid[neg_mask][neg_peaks[-1]] if len(neg_peaks) > 0 else np.nan
                            )
                        else:
                            crit_currents['crit_current_neg'].append(np.nan)
                else:
                    # 對於不需要提取的分支，加入 NaN 以維持陣列長度一致性
                    crit_currents[f'crit_current_{branch}'].append(np.nan)
        else:
            # 對於所有分支，加入 NaN（因為沒有足夠的資料點）
            for branch in ['pos', 'neg']:
                crit_currents[f'crit_current_{branch}'].append(np.nan)
    
    # 建構 DataFrame，根據 current_type 決定如何處理 NaN 值
    df_data = {'B_y_centered_μT': unique_B}
    
    # 確保所有分支陣列長度一致，如果長度不一致則用 NaN 填充
    max_length = len(unique_B)
    for branch_key, values in crit_currents.items():
        if len(values) < max_length:
            # 用 NaN 填充到相同長度
            values.extend([np.nan] * (max_length - len(values)))
        elif len(values) > max_length:
            # 如果長度超過，截斷到正確長度
            values = values[:max_length]
        df_data[branch_key] = values
    
    # 建立初始 DataFrame
    cpr_df_full = pd.DataFrame(df_data)
    
    # 根據 current_type 決定如何過濾資料
    if current_type == "Ic+":
        # 只保留正分支有有效值的行
        cpr_df = cpr_df_full.dropna(subset=['crit_current_pos']).sort_values('B_y_centered_μT')
    elif current_type == "Ic-":
        # 只保留負分支有有效值的行
        cpr_df = cpr_df_full.dropna(subset=['crit_current_neg']).sort_values('B_y_centered_μT')
    else:
        # 雙分支型別：保留任一分支有有效值的行
        cpr_df = cpr_df_full.dropna(how='all', subset=['crit_current_pos', 'crit_current_neg']).sort_values('B_y_centered_μT')
    
    # Apply detrend if requested
    detrend_info = {}
    if use_detrend and len(cpr_df) > 2:  # Need at least 3 points for detrend
        print(f"Applying linear detrend to CPR data...")
        
        # Store original data for comparison
        original_pos = cpr_df['crit_current_pos'].copy() if 'crit_current_pos' in cpr_df.columns else None
        original_neg = cpr_df['crit_current_neg'].copy() if 'crit_current_neg' in cpr_df.columns else None
        
        # Apply detrend to each branch and save trend information
        for branch in ['pos', 'neg']:
            col_name = f'crit_current_{branch}'
            if col_name in cpr_df.columns:
                # Get valid (non-NaN) data points
                valid_mask = cpr_df[col_name].notna()
                if valid_mask.sum() > 2:  # Need at least 3 valid points
                    valid_data = cpr_df.loc[valid_mask, col_name].values
                    valid_B = cpr_df.loc[valid_mask, 'B_y_centered_μT'].values
                    
                    # Apply detrend and calculate the removed trend
                    detrended_data = detrend(valid_data, type='linear')
                    removed_trend = valid_data - detrended_data
                    
                    # Fit a line to the removed trend to get slope and intercept
                    trend_slope, trend_intercept = np.polyfit(valid_B, removed_trend, 1)
                    
                    # Store trend information
                    detrend_info[branch] = {
                        'slope': trend_slope,
                        'intercept': trend_intercept,
                        'original_data': original_pos if branch == 'pos' else original_neg,
                        'valid_mask': valid_mask.copy(),
                        'B_field': valid_B.copy()
                    }
                    
                    # Update the DataFrame with detrended data
                    cpr_df.loc[valid_mask, col_name] = detrended_data
                    
                    # Calculate and report detrend statistics
                    original_std = np.std(valid_data)
                    detrended_std = np.std(detrended_data)
                    trend_removed = original_std - detrended_std
                    
                    print(f"  {branch} branch: Original std = {original_std:.3f} μA, "
                          f"Detrended std = {detrended_std:.3f} μA, "
                          f"Trend removed = {trend_removed:.3f} μA")
                    print(f"    Removed trend: slope = {trend_slope:.6f} μA/μT, intercept = {trend_intercept:.3f} μA")
                else:
                    print(f"  Warning: Not enough valid data points in {branch} branch for detrend")
    
    # 建立 xarray Dataset
    ds_cpr = xr.Dataset.from_dataframe(cpr_df.set_index('B_y_centered_μT'))
    ds_cpr.attrs['run_id'] = ds.run_id
    ds_cpr.attrs['current_type'] = current_type
    ds_cpr.attrs['use_detrend'] = use_detrend
    
    # Store detrend information if detrend was applied
    if use_detrend and detrend_info:
        ds_cpr.attrs['detrend_info'] = detrend_info
        print(f"Detrend information stored for branches: {list(detrend_info.keys())}")
    
    print(f"Extracted {len(cpr_df)} valid CPR data points for {current_type}")
    print(f"Branches processed: {branches_to_extract}")
    print(f"Data shape - B_field: {len(unique_B)}, Pos branch: {len(crit_currents['crit_current_pos'])}, Neg branch: {len(crit_currents['crit_current_neg'])}")
    print(f"Valid data ratios - Pos: {np.sum(np.isfinite(crit_currents['crit_current_pos']))/len(crit_currents['crit_current_pos']):.1%}, Neg: {np.sum(np.isfinite(crit_currents['crit_current_neg']))/len(crit_currents['crit_current_neg']):.1%}")
    
    # 顯示前幾個有效資料點作為範例
    if len(cpr_df) > 0:
        print("Sample extracted data:")
        sample_data = cpr_df.head(3)
        for idx, row in sample_data.iterrows():
            b_field = row['B_y_centered_μT']
            pos_ic = row.get('crit_current_pos', np.nan)
            neg_ic = row.get('crit_current_neg', np.nan)
            print(f"  B={b_field:.1f}µT: Ic+={pos_ic:.3f}µA, Ic-={neg_ic:.3f}µA")
    
    return ds_cpr, cpr_df

def validate_cpr_data_against_current_type(ds_cpr: xr.Dataset, current_type: str) -> bool:
    """驗證提取的 CPR 資料是否符合 current_type 要求"""
    expected_branches = get_branches_for_current_type(current_type)
    
    validation_passed = True
    for branch in expected_branches:
        if f'crit_current_{branch}' not in ds_cpr.data_vars:
            print(f"⚠️  Missing expected branch '{branch}' for {current_type}")
            validation_passed = False
        else:
            # 檢查資料品質
            data = ds_cpr[f'crit_current_{branch}'].values
            valid_ratio = np.sum(np.isfinite(data)) / len(data)
            if valid_ratio < 0.1:  # 少於 10% 的有效資料
                print(f"⚠️  Low data quality for {branch} branch ({valid_ratio:.1%} valid)")
    
    if validation_passed:
        print(f"✅ CPR data validation passed for {current_type}")
    
    return validation_passed

def get_current_type_from_catalog(run_id: int, catalog_path: str) -> str:
    """Get current_type for a given run_id from Catalog.csv"""
    try:
        df = pd.read_csv(catalog_path)
        row = df[df['run_id'] == run_id]
        if not row.empty:
            return row['current_type'].iloc[0]
        else:
            print(f"Warning: run_id {run_id} not found in catalog, using default 'Ic+ w/ Ic-'")
            return "Ic+ w/ Ic-"
    except Exception as e:
        print(f"Warning: Could not read catalog ({e}), using default 'Ic+ w/ Ic-'")
        return "Ic+ w/ Ic-"

def get_measurement_info_from_catalog(run_id: int, catalog_path: str) -> Tuple[str, str]:
    """Get measurement_type and current_type for a given run_id from Catalog.csv"""
    try:
        df = pd.read_csv(catalog_path)
        row = df[df['run_id'] == run_id]
        if not row.empty:
            measurement_type = row['measurement_type'].iloc[0]
            current_type = row['current_type'].iloc[0]
            return measurement_type, current_type
        else:
            print(f"Warning: run_id {run_id} not found in catalog")
            return None, None
    except Exception as e:
        print(f"Warning: Could not read catalog ({e})")
        return None, None

def get_all_cpr_run_ids(catalog_path: str) -> list:
    """Get all run_ids with measurement_type='CPR' from Catalog.csv"""
    try:
        df = pd.read_csv(catalog_path)
        cpr_df = df[df['measurement_type'] == 'CPR']
        run_ids = sorted(cpr_df['run_id'].tolist())
        print(f"Found {len(run_ids)} CPR measurements in catalog: {run_ids[:10]}{'...' if len(run_ids) > 10 else ''}")
        return run_ids
    except Exception as e:
        print(f"Warning: Could not read catalog ({e})")
        return []

def analyze_cpr_with_squid_analyzer(analyzer: SQUIDAnalyzer, ds_cpr: xr.Dataset, run_id: int, output_dir: Path, current_type: str = "Ic+ w/ Ic-") -> Dict[str, Any]:
    """Wraps the main analysis call and saves results."""
    # Store detrend information in analyzer if available
    if 'detrend_info' in ds_cpr.attrs:
        analyzer.detrend_info = ds_cpr.attrs['detrend_info']
        analyzer.detrend_info['detrend_applied'] = ds_cpr.attrs.get('use_detrend', False)
    
    analysis_results = analyzer.analyze(ds_cpr, output_dir, current_type)
    
    # Store results in analyzer for reconstruction
    analyzer.results = analysis_results
    
    for name, result in analysis_results.items():
        if isinstance(result, ModelResult):
            filename = output_dir / f"fit_results_{name}_run_{run_id}.json"
            result_dict = {
                'params': {k: {'value': v.value, 'stderr': v.stderr} for k, v in result.params.items()},
                'rsquared': result.rsquared
            }
            with open(filename, 'w') as f:
                json.dump(result_dict, f, indent=4, default=lambda x: None) # Handle non-serializable stderr
            print(f"Saved {name} results to {filename}")
        
    return analysis_results

def visualize_stage1_results(analyzer: SQUIDAnalyzer, ds_cpr: xr.Dataset, results: Dict[str, Any], run_id: int, output_dir: Path):
    """
    Creates a plot of the Stage 1 pre-fit results, including Delta Ic,
    and displays the model and fit parameters below the plot.
    """
    print("\nVisualizing Stage 1 pre-fit results...")
    
    # Check if we have any data points
    if len(ds_cpr.coords['B_y_centered_μT']) == 0:
        print("   Warning: No data points available for visualization. Skipping Stage 1 plot.")
        return
    
    # Increase figure height to accommodate the text box
    fig, ax = plt.subplots(figsize=(12, 12)) 

    x_data = ds_cpr.coords['B_y_centered_μT'].values
    x_fine = np.linspace(x_data.min(), x_data.max(), 500)
    
    all_plots = []

    for branch, color in [('pos', 'blue'), ('neg', 'red')]:
        # 檢查分支資料是否存在
        branch_data_var = f'crit_current_{branch}'
        if branch_data_var not in ds_cpr.data_vars:
            print(f"   Skipping {branch} branch visualization (data not available)")
            continue
            
        # Get data for plotting - use reconstruction if detrend was applied
        if hasattr(analyzer, 'detrend_info') and analyzer.detrend_info.get('detrend_applied', False):
            # Use reconstructed original data and fit for visualization
            x_plot, y_plot, fit_curve = analyzer.reconstruct_original_data_and_fit(
                x_data, ds_cpr[branch_data_var].values, x_fine, branch, 'pre_fit'
            )
        else:
            # Use data as is
            x_plot, y_plot = x_data, ds_cpr[branch_data_var].values
            fit_result = results.get(f'pre_fit_{branch}', None)
            if fit_result is None:
                print(f"   Skipping fit visualization for {branch} branch (no valid fit)")
                continue
            fit_curve = analyzer.pre_fit_model(x_fine, **fit_result.params.valuesdict())
            
        # Plot raw data
        p_data, = ax.plot(x_plot, y_plot, marker='o', linestyle='--', label=f'Ic{"+" if branch == "pos" else "-"} data', color=color, alpha=0.7, markersize=4, linewidth=1)
        all_plots.append(p_data)
        
        # Plot fit and error band - 使用 .get() 安全訪問
        fit_result = results.get(f'pre_fit_{branch}', None)
        if fit_result is None:
            print(f"   Skipping fit visualization for {branch} branch (no valid fit)")
            continue
        
        # Plot the fit line
        p_fit, = ax.plot(x_fine, fit_curve, color=color, linestyle='--', lw=2, label=f'Ic{"+" if branch == "pos" else "-"} pre-fit')
        all_plots.append(p_fit)

        # Calculate and plot 1-sigma error band
        try:
            # The 'x' keyword argument must match the name of the independent
            # variable in the objective function `fcn2min`.
            dely = fit_result.eval_uncertainty(sigma=1, B_centered_uT=x_fine)
            ax.fill_between(x_fine, fit_curve - dely, fit_curve + dely, color=color, alpha=0.2)
        except Exception as e:
            print(f"Warning: Could not calculate or plot uncertainty for {branch} branch: {e}")

        # Add B0 and Ic0 annotations
        p = fit_result.params.valuesdict()
        p_b0 = ax.axvline(x=p['B_0_uT'], color=color, linestyle=':', lw=1.5, alpha=0.9, label=f'B₀_{branch} = {p["B_0_uT"]:.2f} µT')
        p_ic0 = ax.axhline(y=p['Ic0'], color=color, linestyle=':', lw=1.5, alpha=0.9, label=f'Ic₀_{branch} = {p["Ic0"]:.2f} μA')
        all_plots.extend([p_b0, p_ic0])

        # --- emcee credible & posterior-predictive bands (Stage 1) ---
        emcee_key = f'pre_fit_{branch}_emcee'
        if emcee_key in results:
            emc = results[emcee_key]
            xx = emc['x_eval']
            mb = emc['post_mean_bands']
            pb = emc['post_pred_bands']

            # 95% credible band (mean model)
            ax.fill_between(xx, mb['p2p5'], mb['p97p5'], alpha=0.12, label=f'95% cred (mean, {branch})', color=color)
            # 68% credible band
            ax.fill_between(xx, mb['p16'], mb['p84'], alpha=0.20, label=f'68% cred (mean, {branch})', color=color)
            # 95% posterior predictive band
            ax.fill_between(xx, pb['p2p5'], pb['p97p5'], alpha=0.10, label=f'95% post. pred. ({branch})', color='gray')

    # --- Delta Ic on secondary axis ---
    ax2 = ax.twinx()
    
    # 檢查是否有足夠的分支資料來計算 Delta Ic
    has_pos_data = 'crit_current_pos' in ds_cpr.data_vars
    has_neg_data = 'crit_current_neg' in ds_cpr.data_vars
    has_both_fits = (results.get('pre_fit_pos') is not None and 
                     results.get('pre_fit_neg') is not None)
    
    if has_pos_data and has_neg_data and has_both_fits:
        # 雙分支情況：正常計算 Delta Ic
        # Get reconstructed data if detrend was applied
        if hasattr(analyzer, 'detrend_info') and analyzer.detrend_info.get('detrend_applied', False):
            # Use reconstructed data for Delta Ic calculation
            _, ic_pos_plot, _ = analyzer.reconstruct_original_data_and_fit(
                x_data, ds_cpr['crit_current_pos'].values, x_fine, 'pos', 'pre_fit'
            )
            _, ic_neg_plot, _ = analyzer.reconstruct_original_data_and_fit(
                x_data, ds_cpr['crit_current_neg'].values, x_fine, 'neg', 'pre_fit'
            )
            ic_pos_data, ic_neg_data = ic_pos_plot, ic_neg_plot
        else:
            ic_pos_data = ds_cpr['crit_current_pos'].values
            ic_neg_data = ds_cpr['crit_current_neg'].values
            
        valid_mask = np.isfinite(ic_pos_data) & np.isfinite(ic_neg_data)
        delta_ic_data = ic_pos_data[valid_mask] - np.abs(ic_neg_data[valid_mask])
        
        fit_pos_result = results.get('pre_fit_pos', None)
        fit_neg_result = results.get('pre_fit_neg', None)
        if fit_pos_result is not None and fit_neg_result is not None:
            # Get reconstructed fit curves if detrend was applied
            if hasattr(analyzer, 'detrend_info') and analyzer.detrend_info.get('detrend_applied', False):
                _, _, fit_pos_curve = analyzer.reconstruct_original_data_and_fit(
                    x_data, ds_cpr['crit_current_pos'].values, x_fine, 'pos', 'pre_fit'
                )
                _, _, fit_neg_curve = analyzer.reconstruct_original_data_and_fit(
                    x_data, ds_cpr['crit_current_neg'].values, x_fine, 'neg', 'pre_fit'
                )
            else:
                fit_pos_curve = analyzer.pre_fit_model(x_fine, **fit_pos_result.values)
                fit_neg_curve = analyzer.pre_fit_model(x_fine, **fit_neg_result.values)
                
            delta_ic_fit = fit_pos_curve - np.abs(fit_neg_curve)
            
            p_delta_data, = ax2.plot(x_data[valid_mask], delta_ic_data, color='green', linestyle='--', marker='.', markersize=5, label='ΔIc (data)')
            p_delta_fit, = ax2.plot(x_fine, delta_ic_fit, 'g-', lw=2, label='ΔIc (pre-fit)')
            all_plots.extend([p_delta_data, p_delta_fit])
        else:
            print("   Warning: Cannot calculate ΔIc pre-fit - missing fit results for one or both branches")
    
    
    ax2.set_ylabel('ΔIc = Ic+ - |Ic-| (μA)', color='g', fontsize=12)
    ax2.tick_params(axis='y', labelcolor='g')
    ax2.grid(True, linestyle=':', color='g', alpha=0.5)

    ax.set_title(f'Stage 1: Pre-fitting with Simple Sine Model (Run {run_id})', fontsize=16)
    ax.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax.set_ylabel('Critical Current (µA)', fontsize=12)
    
    all_labels = [p.get_label() for p in all_plots]
    fig.legend(all_plots, all_labels, loc='upper center', bbox_to_anchor=(0.5, 0.26), ncol=5, fontsize=10)
    ax.grid(True, alpha=0.4)

    # --- Add model and parameter text below plot ---
    model_text = r"Model: $I_c = I_{c0} + r \cdot \phi + I_{c1} \cdot \sin(\phi)$, $\phi = \frac{2\pi A_{eff}(B-B_0)}{\Phi_0}$"
    
    param_details = {
        'Ic0': {'unit': 'μA', 'name': r'$I_{c0}$'},
        'r': {'unit': '', 'name': r'$r$'},
        'Ic1': {'unit': 'μA', 'name': r'$I_{c1}$'},
        'A_eff_um2': {'unit': 'μm²', 'name': r'$A_{eff}$'},
        'B_0_uT': {'unit': 'μT', 'name': r'$B_0$'}
    }

    pos_params = []
    neg_params = []

    # Process Positive Branch
    fit_result_pos = results.get('pre_fit_pos')
    if fit_result_pos is not None:
        # Define parameter order for better readability - unified ordering
        priority_params = ['Ic0', 'Ic1', 'B_0_uT', 'A_eff_um2']
        
        # First, process priority parameters in specified order
        for param_name in priority_params:
            if param_name in fit_result_pos.params:
                param = fit_result_pos.params[param_name]
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')
                if param.stderr is not None:
                    percent = abs(param.stderr / param.value) * 100
                    val_str = f"{param.value:.5g} ± {param.stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{param.value:.5g}"
                init_val_str = f"(init = {param.init_value:.5g})" if hasattr(param, 'init_value') and param.init_value is not None else ""
                if unit:
                    pos_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    pos_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        
        # Then, process remaining parameters in their original order
        processed_params = set(priority_params)
        for name, param in fit_result_pos.params.items():
            if name not in processed_params:
                detail = param_details.get(name, {})
                p_name = detail.get('name', name)
                unit = detail.get('unit', '')
                if param.stderr is not None:
                    percent = abs(param.stderr / param.value) * 100
                    val_str = f"{param.value:.5g} ± {param.stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{param.value:.5g}"
                init_val_str = f"(init = {param.init_value:.5g})" if hasattr(param, 'init_value') and param.init_value is not None else ""
                if unit:
                    pos_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    pos_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        pos_params.append((r'$\chi^2$', f"{fit_result_pos.chisqr:.4g}"))
        pos_params.append((r'$R^2$', f"{fit_result_pos.rsquared:.4f}"))
    else:
        pos_params.append(("No fit", "No valid data for positive branch"))

    # Process Negative Branch
    fit_result_neg = results.get('pre_fit_neg')
    if fit_result_neg is not None:
        # Define parameter order for better readability - unified ordering
        priority_params = ['Ic0', 'Ic1', 'B_0_uT', 'A_eff_um2']
        
        # First, process priority parameters in specified order
        for param_name in priority_params:
            if param_name in fit_result_neg.params:
                param = fit_result_neg.params[param_name]
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')
                if param.stderr is not None:
                    percent = abs(param.stderr / param.value) * 100
                    val_str = f"{param.value:.5g} ± {param.stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{param.value:.5g}"
                init_val_str = f"(init = {param.init_value:.5g})" if hasattr(param, 'init_value') and param.init_value is not None else ""
                if unit:
                    neg_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    neg_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        
        # Then, process remaining parameters in their original order
        processed_params = set(priority_params)
        for name, param in fit_result_neg.params.items():
            if name not in processed_params:
                detail = param_details.get(name, {})
                p_name = detail.get('name', name)
                unit = detail.get('unit', '')
                if param.stderr is not None:
                    percent = abs(param.stderr / param.value) * 100
                    val_str = f"{param.value:.5g} ± {param.stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{param.value:.5g}"
                init_val_str = f"(init = {param.init_value:.5g})" if hasattr(param, 'init_value') and param.init_value is not None else ""
                if unit:
                    neg_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    neg_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        neg_params.append((r'$\chi^2$', f"{fit_result_neg.chisqr:.4g}"))
        neg_params.append((r'$R^2$', f"{fit_result_neg.rsquared:.4f}"))
    else:
        neg_params.append(("No fit", "No valid data for negative branch"))

    # --- New line-by-line layout logic for perfect alignment ---

    # Add text to figure
    y_pos = 0.01
    fig.text(0.5, y_pos + 0.15, f"{model_text}", ha='center', va='bottom', fontsize=25, family='monospace')

    text_y_start = y_pos + 0.13
    line_height = 0.018  # Adjusted line height in figure coordinates

    common_props = {'va': 'top', 'fontsize': 9, 'family': 'monospace'}

    # Column Headers - 動態根據可用分支調整
    header_y = text_y_start
    
    # 檢查哪些分支有有效的參數
    has_pos_params = len(pos_params) > 0 and pos_params[0][0] != "No fit"
    has_neg_params = len(neg_params) > 0 and neg_params[0][0] != "No fit"
    
    if has_pos_params and has_neg_params:
        # 雙欄佈局
        fig.text(0.02, header_y, f'$I_c^{{+}}$ Fit', ha='left', weight='bold', **common_props)
        fig.text(0.52, header_y, f'$I_c^{{-}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.02, 0.08, 0.10
        right_name_x, right_colon_x, right_value_x = 0.52, 0.58, 0.60
        layout_mode = 'dual'
    elif has_pos_params:
        # 只有正分支：單欄置中佈局
        fig.text(0.25, header_y, f'$I_c^{{+}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.25, 0.31, 0.33
        right_name_x, right_colon_x, right_value_x = None, None, None
        layout_mode = 'single_pos'
    elif has_neg_params:
        # 只有負分支：單欄置中佈局
        fig.text(0.25, header_y, f'$I_c^{{-}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.25, 0.31, 0.33
        right_name_x, right_colon_x, right_value_x = None, None, None
        layout_mode = 'single_neg'
    else:
        # 沒有有效參數
        fig.text(0.25, header_y, 'No Valid Fits Available', ha='left', weight='bold', **common_props)
        layout_mode = 'none'
    
    # Loop and plot each parameter line by line
    if layout_mode == 'dual':
        max_len = max(len(pos_params), len(neg_params))
        for i in range(max_len):
            current_y = text_y_start - (i + 1.5) * line_height # +1.5 to leave space after header

            # Left column (Positive)
            if i < len(pos_params):
                name, value = pos_params[i]
                fig.text(left_name_x, current_y, name, ha='left', **common_props)
                fig.text(left_colon_x, current_y, ':', ha='left', **common_props)
                fig.text(left_value_x, current_y, value, ha='left', **common_props)

            # Right column (Negative)
            if i < len(neg_params):
                name, value = neg_params[i]
                fig.text(right_name_x, current_y, name, ha='left', **common_props)
                fig.text(right_colon_x, current_y, ':', ha='left', **common_props)
                fig.text(right_value_x, current_y, value, ha='left', **common_props)
    elif layout_mode in ['single_pos', 'single_neg']:
        # 單欄佈局
        params_to_show = pos_params if layout_mode == 'single_pos' else neg_params
        for i, (name, value) in enumerate(params_to_show):
            current_y = text_y_start - (i + 1.5) * line_height
            fig.text(left_name_x, current_y, name, ha='left', **common_props)
            fig.text(left_colon_x, current_y, ':', ha='left', **common_props)
            fig.text(left_value_x, current_y, value, ha='left', **common_props)

    # Determine max_len if not already set (for single column layouts)
    if 'max_len' not in locals():
        if layout_mode == 'single_pos':
            max_len = len(pos_params)
        elif layout_mode == 'single_neg':
            max_len = len(neg_params)
        else:  # layout_mode == 'none'
            max_len = 0

    # Add a continuous vertical line between the columns (only for dual layout)
    if layout_mode == 'dual':
        num_param_lines = max_len
        line_y_end = text_y_start - line_height # Start below header
        line_y_start = text_y_start - (num_param_lines + 2) * line_height
        line = lines.Line2D([0.5, 0.5], [line_y_start, line_y_end],
                            transform=fig.transFigure, color='grey', linestyle='-', linewidth=1)
        fig.add_artist(line)

    # Adjust rect to make space for legend and text
    fig.tight_layout(rect=(0, 0.25, 1, 0.99))
    
    output_path = output_dir / f"stage1_prefit_run_{run_id}.png"
    plt.savefig(output_path, dpi=300)
    print(f"Stage 1 plot saved to {output_path}")
    plt.close(fig)

def visualize_stage1_2_results(analyzer: SQUIDAnalyzer, ds_cpr: xr.Dataset, results: Dict[str, Any], run_id: int, output_dir: Path):
    """
    Creates a plot of the Stage 1.2 single harmonic transparency model results,
    and displays the model and fit parameters below the plot.
    """
    print("\nVisualizing Stage 1.2 single harmonic transparency model results...")
    
    # Check if we have any data points
    if len(ds_cpr.coords['B_y_centered_μT']) == 0:
        print("   Warning: No data points available for visualization. Skipping Stage 1.2 plot.")
        return
    
    # Increase figure height to accommodate the text box
    fig, ax = plt.subplots(figsize=(12, 14))

    x_data = ds_cpr.coords['B_y_centered_μT'].values
    x_fine = np.linspace(x_data.min(), x_data.max(), 500)

    all_plots = []

    # 檢查哪些分支有數據
    has_pos_data = 'crit_current_pos' in ds_cpr.data_vars and len(ds_cpr['crit_current_pos'].values) > 0
    has_neg_data = 'crit_current_neg' in ds_cpr.data_vars and len(ds_cpr['crit_current_neg'].values) > 0
    
    for branch, color in [('pos', 'blue'), ('neg', 'red')]:
        if branch == 'pos' and not has_pos_data:
            continue
        if branch == 'neg' and not has_neg_data:
            continue
            
        # Plot raw data
        p_data, = ax.plot(x_data, ds_cpr[f'crit_current_{branch}'].values, marker='o', linestyle='--', label=f'Ic{"+" if branch == "pos" else "-"} data', color=color, alpha=0.7, markersize=4, linewidth=1)
        all_plots.append(p_data)

        # Plot fit and error band
        fit_result = results.get(f'single_junction_fit_{branch}', None)
        if fit_result is None:
            print(f"   Skipping fit visualization for {branch} branch (no valid data)")
            continue
            
        fit_curve = analyzer.single_junction_model(x_fine, **fit_result.params.valuesdict())

        # Plot the fit line
        p_fit, = ax.plot(x_fine, fit_curve, color=color, linestyle='-', lw=2, label=f'Ic{"+" if branch == "pos" else "-"} fit')
        all_plots.append(p_fit)

        # Calculate and plot 1-sigma error band
        try:
            dely = fit_result.eval_uncertainty(sigma=1, B_centered_uT=x_fine)
            ax.fill_between(x_fine, fit_curve - dely, fit_curve + dely, color=color, alpha=0.2)
        except Exception as e:
            print(f"Warning: Could not calculate or plot uncertainty for {branch} branch: {e}")

        # Add B0 and Ic0 annotations
        p = fit_result.params.valuesdict()
        p_b0 = ax.axvline(x=p['B_0_uT'], color=color, linestyle=':', lw=1.5, alpha=0.9, label=f'B₀_{branch} = {p["B_0_uT"]:.2f} µT')
        p_ic0 = ax.axhline(y=p['Ic0'], color=color, linestyle=':', lw=1.5, alpha=0.9, label=f'Ic₀_{branch} = {p["Ic0"]:.2f} μA')
        all_plots.extend([p_b0, p_ic0])

    # --- Delta Ic on secondary axis ---
    ax2 = ax.twinx()

    # 檢查是否有足夠的數據來計算 Delta Ic
    if has_pos_data and has_neg_data:
        ic_pos_data = ds_cpr['crit_current_pos'].values
        ic_neg_data = ds_cpr['crit_current_neg'].values
        valid_mask = np.isfinite(ic_pos_data) & np.isfinite(ic_neg_data)
        delta_ic_data = ic_pos_data[valid_mask] - np.abs(ic_neg_data[valid_mask])

        # Check if both fits are available for delta Ic calculation
        fit_pos_result = results.get('single_junction_fit_pos', None)
        fit_neg_result = results.get('single_junction_fit_neg', None)
        if fit_pos_result is None or fit_neg_result is None:
            print("   Warning: Cannot calculate ΔIc fit - missing fit results for one or both branches")
            delta_ic_fit = np.zeros_like(x_fine)  # Use zeros as placeholder
        else:
            fit_pos_curve = analyzer.single_junction_model(x_fine, **fit_pos_result.params.valuesdict())
            fit_neg_curve = analyzer.single_junction_model(x_fine, **fit_neg_result.params.valuesdict())
            delta_ic_fit = fit_pos_curve - np.abs(fit_neg_curve)

        p_delta_data, = ax2.plot(x_data[valid_mask], delta_ic_data, color='green', linestyle='--', marker='.', markersize=5, label='ΔIc (data)')
        p_delta_fit, = ax2.plot(x_fine, delta_ic_fit, 'g-', lw=2, label='ΔIc ( fit)')
        all_plots.extend([p_delta_data, p_delta_fit])

        ax2.set_ylabel('ΔIc = Ic+ - |Ic-| (μA)', color='g', fontsize=12)
        ax2.tick_params(axis='y', labelcolor='g')
        ax2.grid(True, linestyle=':', color='g', alpha=0.5)
    else:
        print("   Warning: Cannot calculate ΔIc - missing data for one or both branches")
        # 隱藏第二個y軸
        ax2.set_visible(False)

    ax.set_title(f'Stage 1.2: Single Harmonic Transparency Model Fit (Run {run_id})', fontsize=16)
    ax.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax.set_ylabel('Critical Current (µA)', fontsize=12)

    all_labels = [p.get_label() for p in all_plots]
    fig.legend(all_plots, all_labels, loc='upper center', bbox_to_anchor=(0.5, 0.30), ncol=5, fontsize=10)
    ax.grid(True, alpha=0.4)

    # --- Add model and parameter text below plot ---
    model_text = r"Model: $I_c = I_{c0} + r \cdot \phi + I_{c1} \cdot \sin(\phi)$, $\phi = \frac{2\pi A_{eff}(B-B_0)}{\Phi_0}$"

    param_details = {
        'Ic0': {'unit': 'μA', 'name': r'$I_{c0}$'},
        'r': {'unit': 'μA/T', 'name': r'$r$'},
        'Ic1': {'unit': 'μA', 'name': r'$I_{c1}$'},
        'A_eff_um2': {'unit': 'μm²', 'name': r'$A_{eff}$'},
        'B_0_uT': {'unit': 'μT', 'name': r'$B_0$'}
    }

    pos_params = []
    neg_params = []

    # Process Positive Branch
    fit_result_pos = results.get('single_junction_fit_pos', None)
    if fit_result_pos is None:
        pos_params.append(("No fit", "No valid data for positive branch"))
    else:
        priority_params = ['Ic0', 'Ic1', 'B_0_uT', 'A_eff_um2']
        
        # First, process priority parameters in specified order
        for param_name in priority_params:
            if param_name in fit_result_pos.params:
                param = fit_result_pos.params[param_name]
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')
                if param.stderr is not None:
                    percent = abs(param.stderr / param.value) * 100
                    val_str = f"{param.value:.5g} ± {param.stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{param.value:.5g}"
                init_val_str = f"(init = {param.init_value:.5g})" if hasattr(param, 'init_value') and param.init_value is not None else ""
                if unit:
                    pos_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    pos_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        
        # Then, process remaining parameters in their original order
        processed_params = set(priority_params)
        for name, param in fit_result_pos.params.items():
            if name not in processed_params:
                detail = param_details.get(name, {})
                p_name = detail.get('name', name)
                unit = detail.get('unit', '')
                if param.stderr is not None:
                    percent = abs(param.stderr / param.value) * 100
                    val_str = f"{param.value:.5g} ± {param.stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{param.value:.5g}"
                init_val_str = f"(init = {param.init_value:.5g})" if hasattr(param, 'init_value') and param.init_value is not None else ""
                if unit:
                    pos_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    pos_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        pos_params.append((r'$\chi^2$', f"{fit_result_pos.chisqr:.4g}"))
        pos_params.append((r'$R^2$', f"{fit_result_pos.rsquared:.4f}"))

    # Process Negative Branch
    fit_result_neg = results.get('single_junction_fit_neg', None)
    if fit_result_neg is None:
        neg_params.append(("No fit", "No valid data for negative branch"))
    else:
        # Define parameter order for better readability - unified ordering
        priority_params = ['Ic0', 'Ic1', 'B_0_uT', 'A_eff_um2']
        
        # First, process priority parameters in specified order
        for param_name in priority_params:
            if param_name in fit_result_neg.params:
                param = fit_result_neg.params[param_name]
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')
                if param.stderr is not None:
                    percent = abs(param.stderr / param.value) * 100
                    val_str = f"{param.value:.5g} ± {param.stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{param.value:.5g}"
                init_val_str = f"(init = {param.init_value:.5g})" if hasattr(param, 'init_value') and param.init_value is not None else ""
                if unit:
                    neg_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    neg_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        
        # Then, process remaining parameters in their original order
        processed_params = set(priority_params)
        for name, param in fit_result_neg.params.items():
            if name not in processed_params:
                detail = param_details.get(name, {})
                p_name = detail.get('name', name)
                unit = detail.get('unit', '')
                if param.stderr is not None:
                    percent = abs(param.stderr / param.value) * 100
                    val_str = f"{param.value:.5g} ± {param.stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{param.value:.5g}"
                init_val_str = f"(init = {param.init_value:.5g})" if hasattr(param, 'init_value') and param.init_value is not None else ""
                if unit:
                    neg_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    neg_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        neg_params.append((r'$\chi^2$', f"{fit_result_neg.chisqr:.4g}"))
        neg_params.append((r'$R^2$', f"{fit_result_neg.rsquared:.4f}"))

    # --- New line-by-line layout logic for perfect alignment ---

    # Add text to figure
    y_pos = 0.01
    fig.text(0.5, y_pos + 0.18, f"{model_text}", ha='center', va='bottom', fontsize=25, family='monospace')
    text_y_start = y_pos + 0.15
    line_height = 0.018  # Adjusted line height in figure coordinates

    common_props = {'va': 'top', 'fontsize': 9, 'family': 'monospace'}

    # Column Headers - 動態根據可用分支調整
    header_y = text_y_start
    
    # 檢查哪些分支有有效的參數
    has_pos_params = len(pos_params) > 0 and pos_params[0][0] != "No fit"
    has_neg_params = len(neg_params) > 0 and neg_params[0][0] != "No fit"
    
    if has_pos_params and has_neg_params:
        # 雙欄佈局
        fig.text(0.02, header_y, f'$I_c^{{+}}$ Fit', ha='left', weight='bold', **common_props)
        fig.text(0.52, header_y, f'$I_c^{{-}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.02, 0.08, 0.10
        right_name_x, right_colon_x, right_value_x = 0.52, 0.58, 0.60
        layout_mode = 'dual'
    elif has_pos_params:
        # 只有正分支：單欄置中佈局
        fig.text(0.25, header_y, f'$I_c^{{+}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.25, 0.31, 0.33
        right_name_x, right_colon_x, right_value_x = None, None, None
        layout_mode = 'single_pos'
    elif has_neg_params:
        # 只有負分支：單欄置中佈局
        fig.text(0.25, header_y, f'$I_c^{{-}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.25, 0.31, 0.33
        right_name_x, right_colon_x, right_value_x = None, None, None
        layout_mode = 'single_neg'
    else:
        # 沒有有效參數
        fig.text(0.25, header_y, 'No Valid Fits Available', ha='left', weight='bold', **common_props)
        layout_mode = 'none'
    
    # Loop and plot each parameter line by line
    if layout_mode == 'dual':
        max_len = max(len(pos_params), len(neg_params))
        for i in range(max_len):
            current_y = text_y_start - (i + 1.5) * line_height # +1.5 to leave space after header

            # Left column (Positive)
            if i < len(pos_params):
                name, value = pos_params[i]
                fig.text(left_name_x, current_y, name, ha='left', **common_props)
                fig.text(left_colon_x, current_y, ':', ha='left', **common_props)
                fig.text(left_value_x, current_y, value, ha='left', **common_props)

            # Right column (Negative)
            if i < len(neg_params):
                name, value = neg_params[i]
                fig.text(right_name_x, current_y, name, ha='left', **common_props)
                fig.text(right_colon_x, current_y, ':', ha='left', **common_props)
                fig.text(right_value_x, current_y, value, ha='left', **common_props)
    elif layout_mode in ['single_pos', 'single_neg']:
        # 單欄佈局
        params_to_show = pos_params if layout_mode == 'single_pos' else neg_params
        for i, (name, value) in enumerate(params_to_show):
            current_y = text_y_start - (i + 1.5) * line_height
            fig.text(left_name_x, current_y, name, ha='left', **common_props)
            fig.text(left_colon_x, current_y, ':', ha='left', **common_props)
            fig.text(left_value_x, current_y, value, ha='left', **common_props)

    # Add a continuous vertical line between the columns (only for dual layout)
    if layout_mode == 'dual':
        num_param_lines = max(len(pos_params), len(neg_params))
        line_y_end = text_y_start - line_height # Start below header
        line_y_start = text_y_start - (num_param_lines + 2) * line_height
        line = lines.Line2D([0.5, 0.5], [line_y_start, line_y_end],
                            transform=fig.transFigure, color='grey', linestyle='-', linewidth=1)
        fig.add_artist(line)

    # Adjust rect to make space for legend and text
    fig.tight_layout(rect=(0, 0.30, 1, 0.99))

    output_path = output_dir / f"stage1_2_single_harmonic_transparency_run_{run_id}.png"
    plt.savefig(output_path, dpi=300)
    print(f"Stage 1.2 plot saved to {output_path}")
    plt.close(fig)

def visualize_stage2_simple_results(analyzer: SQUIDAnalyzer, ds_cpr: xr.Dataset, results: Dict[str, Any], run_id: int, output_dir: Path):
    """
    Creates a simple plot of the Stage 2 single-transparency model results,
    similar to Stage 1 and Stage 1.2 design.
    """
    print("\nVisualizing Stage 2 single-T model results (simple version)...")
    
    # Check if we have any data points
    if len(ds_cpr.coords['B_y_centered_μT']) == 0:
        print("   Warning: No data points available for visualization. Skipping Stage 2 simple plot.")
        return
    
    # Increase figure height to accommodate the text box
    fig, ax = plt.subplots(figsize=(12, 14))

    x_data = ds_cpr.coords['B_y_centered_μT'].values
    x_fine = np.linspace(x_data.min(), x_data.max(), 500)

    all_plots = []

    # 檢查哪些分支有數據
    has_pos_data = 'crit_current_pos' in ds_cpr.data_vars and len(ds_cpr['crit_current_pos'].values) > 0
    has_neg_data = 'crit_current_neg' in ds_cpr.data_vars and len(ds_cpr['crit_current_neg'].values) > 0
    
    for branch, color in [('pos', 'blue'), ('neg', 'red')]:
        if branch == 'pos' and not has_pos_data:
            continue
        if branch == 'neg' and not has_neg_data:
            continue
            
        # Plot raw data
        p_data, = ax.plot(x_data, ds_cpr[f'crit_current_{branch}'].values, marker='o', linestyle='--', label=f'Ic{"+" if branch == "pos" else "-"} data', color=color, alpha=0.7, markersize=4, linewidth=1)
        all_plots.append(p_data)

        # Plot fit and error band
        fit_result = results.get(f'single_T_fit_{branch}', None)
        if fit_result is None:
            print(f"   Skipping fit visualization for {branch} branch (no valid data)")
            continue
        fit_curve = analyzer.single_T_model(x_fine, **fit_result.params.valuesdict())

        # Plot the fit line
        p_fit, = ax.plot(x_fine, fit_curve, color=color, linestyle='-', lw=2, label=f'Ic{"+" if branch == "pos" else "-"} fit')
        all_plots.append(p_fit)

        # Calculate and plot 1-sigma error band
        try:
            dely = fit_result.eval_uncertainty(sigma=1, B_centered_uT=x_fine)
            ax.fill_between(x_fine, fit_curve - dely, fit_curve + dely, color=color, alpha=0.2)
        except Exception as e:
            print(f"Warning: Could not calculate or plot uncertainty for {branch} branch: {e}")

        # Plot junction components (simple sine model)
        p = fit_result.params.valuesdict()
        flux_phase = analyzer._flux_phase(x_fine, p['B_0_uT'], p['A_eff_um2'])
        I_sine = p['Ic1'] * np.sin(flux_phase)
        # For simple model, we have linear and sine components
        
        comp_color = 'darkcyan' if branch == 'pos' else 'darkred'
        p_comp, = ax.plot(x_fine, I_sine + p['Ic0'], color=comp_color, linestyle=':', lw=2, label=f'Ic0+Sine ({branch})')
        all_plots.append(p_comp)

        # Add B0 and Ic0 annotations
        p_b0 = ax.axvline(x=p['B_0_uT'], color=color, linestyle=':', lw=1.5, alpha=0.9, label=f'B₀_{branch} = {p["B_0_uT"]:.2f} µT')
        p_ic0 = ax.axhline(y=p['Ic0'], color=color, linestyle=':', lw=1.5, alpha=0.9, label=f'Ic₀_{branch} = {p["Ic0"]:.2f} μA')
        all_plots.extend([p_b0, p_ic0])

    # --- Delta Ic on secondary axis ---
    ax2 = ax.twinx()

    # 檢查是否有足夠的數據來計算 Delta Ic
    if has_pos_data and has_neg_data:
        ic_pos_data = ds_cpr['crit_current_pos'].values
        ic_neg_data = ds_cpr['crit_current_neg'].values
        valid_mask = np.isfinite(ic_pos_data) & np.isfinite(ic_neg_data)
        delta_ic_data = ic_pos_data[valid_mask] - np.abs(ic_neg_data[valid_mask])

        # Check if both fits are available for delta Ic calculation
        fit_pos_result = results.get('single_T_fit_pos', None)
        fit_neg_result = results.get('single_T_fit_neg', None)
        if fit_pos_result is None or fit_neg_result is None:
            print("   Warning: Cannot calculate ΔIc fit - missing fit results for one or both branches")
            delta_ic_fit = np.zeros_like(x_fine)  # Use zeros as placeholder
        else:
            fit_pos_curve = analyzer.single_T_model(x_fine, **fit_pos_result.params.valuesdict())
            fit_neg_curve = analyzer.single_T_model(x_fine, **fit_neg_result.params.valuesdict())
            delta_ic_fit = fit_pos_curve - np.abs(fit_neg_curve)

        p_delta_data, = ax2.plot(x_data[valid_mask], delta_ic_data, color='green', linestyle='--', marker='.', markersize=5, label='ΔIc (data)')
        p_delta_fit, = ax2.plot(x_fine, delta_ic_fit, 'g-', lw=2, label='ΔIc (fit)')
        all_plots.extend([p_delta_data, p_delta_fit])

        ax2.set_ylabel('ΔIc = Ic+ - |Ic-| (μA)', color='g', fontsize=12)
        ax2.tick_params(axis='y', labelcolor='g')
        ax2.grid(True, linestyle=':', color='g', alpha=0.5)
    else:
        print("   Warning: Cannot calculate ΔIc - missing data for one or both branches")
        # 隱藏第二個y軸
        ax2.set_visible(False)

    ax.set_title(f'Stage 2: Single-Transparency Model Fit (Run {run_id})', fontsize=16)
    ax.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax.set_ylabel('Critical Current (µA)', fontsize=12)

    all_labels = [p.get_label() for p in all_plots]
    fig.legend(all_plots, all_labels, loc='upper center', bbox_to_anchor=(0.5, 0.305), ncol=7, fontsize=10)
    ax.grid(True, alpha=0.4)

    # --- Add model and parameter text below plot ---
    model_text = r"Model: $I_c = I_{c0} + r \cdot CPR(\phi, I_{c\_cpr}, T_{cpr})$, $\phi = \frac{2\pi A_{eff}(B-B_0)}{\Phi_0}$"

    param_details = {
        'Ic0': {'unit': 'μA', 'name': r'$I_{c0}$'},
        'r': {'unit': 'μA/T', 'name': r'$r$'},
        'Ic1': {'unit': 'μA', 'name': r'$I_{c1}$'},
        'Ic2': {'unit': 'μA', 'name': r'$I_{c2}$'},
        'T1': {'unit': '', 'name': r'$T_1$'},
        'phi_21': {'unit': '°', 'name': r'$\phi_{21}$'},
        'ratio_ic2_ic1': {'unit': '', 'name': r'$I_{c2}/I_{c1}$'},
        'A_eff_um2': {'unit': 'μm²', 'name': r'$A_{eff}$'},
        'B_0_uT': {'unit': 'μT', 'name': r'$B_0$'}
    }

    pos_params = []
    neg_params = []

    # Process Positive Branch
    fit_result_pos = results.get('single_T_fit_pos', None)
    if fit_result_pos is None:
        pos_params.append(("No fit", "No valid data for positive branch"))
    else:
        # Define parameter order for better readability - unified ordering
        priority_params = ['Ic0', 'Ic1', 'Ic2', 'ratio_ic2_ic1', 'phi_21', 'T1', 'B_0_uT', 'A_eff_um2']
        
        # First, process priority parameters in specified order
        for param_name in priority_params:
            if param_name in fit_result_pos.params:
                param = fit_result_pos.params[param_name]
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')
                
                # Convert phi_21 from radians to degrees for display
                if param_name == 'phi_21':
                    display_value = param.value * 180 / np.pi
                    display_stderr = param.stderr * 180 / np.pi if param.stderr is not None else None
                    display_init = param.init_value * 180 / np.pi if hasattr(param, 'init_value') and param.init_value is not None else None
                else:
                    display_value = param.value
                    display_stderr = param.stderr
                    display_init = param.init_value if hasattr(param, 'init_value') else None
                    
                if display_stderr is not None:
                    percent = abs(display_stderr / display_value) * 100
                    val_str = f"{display_value:.5g} ± {display_stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{display_value:.5g}"
                init_val_str = f"(init = {display_init:.5g})" if display_init is not None else ""
                if unit:
                    pos_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    pos_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        
        # Then, process remaining parameters in their original order
        processed_params = set(priority_params)
        for param_name, param in fit_result_pos.params.items():
            if param_name not in processed_params:
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')
                
                # Convert phi_21 from radians to degrees for display
                if param_name == 'phi_21':
                    display_value = param.value * 180 / np.pi
                    display_stderr = param.stderr * 180 / np.pi if param.stderr is not None else None
                    display_init = param.init_value * 180 / np.pi if hasattr(param, 'init_value') and param.init_value is not None else None
                else:
                    display_value = param.value
                    display_stderr = param.stderr
                    display_init = param.init_value if hasattr(param, 'init_value') else None
                    
                if display_stderr is not None:
                    percent = abs(display_stderr / display_value) * 100
                    val_str = f"{display_value:.5g} ± {display_stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{display_value:.5g}"
                init_val_str = f"(init = {display_init:.5g})" if display_init is not None else ""
                if unit:
                    pos_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    pos_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        pos_params.append((r'$\chi^2$', f"{fit_result_pos.chisqr:.4g}"))
        pos_params.append((r'$R^2$', f"{fit_result_pos.rsquared:.4f}"))

    # Process Negative Branch
    fit_result_neg = results.get('single_T_fit_neg', None)
    if fit_result_neg is None:
        neg_params.append(("No fit", "No valid data for negative branch"))
    else:
        # Define parameter order for better readability - unified ordering
        priority_params = ['Ic0', 'Ic1', 'Ic2', 'ratio_ic2_ic1', 'phi_21', 'T1', 'B_0_uT', 'A_eff_um2']
        
        # First, process priority parameters in specified order
        for param_name in priority_params:
            if param_name in fit_result_neg.params:
                param = fit_result_neg.params[param_name]
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')
                
                # Convert phi_21 from radians to degrees for display
                if param_name == 'phi_21':
                    display_value = param.value * 180 / np.pi
                    display_stderr = param.stderr * 180 / np.pi if param.stderr is not None else None
                    display_init = param.init_value * 180 / np.pi if hasattr(param, 'init_value') and param.init_value is not None else None
                else:
                    display_value = param.value
                    display_stderr = param.stderr
                    display_init = param.init_value if hasattr(param, 'init_value') else None
                    
                if display_stderr is not None:
                    percent = abs(display_stderr / display_value) * 100
                    val_str = f"{display_value:.5g} ± {display_stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{display_value:.5g}"
                init_val_str = f"(init = {display_init:.5g})" if display_init is not None else ""
                if unit:
                    neg_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    neg_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        
        # Then, process remaining parameters in their original order
        processed_params = set(priority_params)
        for param_name, param in fit_result_neg.params.items():
            if param_name not in processed_params:
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')
                
                # Convert phi_21 from radians to degrees for display
                if param_name == 'phi_21':
                    display_value = param.value * 180 / np.pi
                    display_stderr = param.stderr * 180 / np.pi if param.stderr is not None else None
                    display_init = param.init_value * 180 / np.pi if hasattr(param, 'init_value') and param.init_value is not None else None
                else:
                    display_value = param.value
                    display_stderr = param.stderr
                    display_init = param.init_value if hasattr(param, 'init_value') else None
                    
                if display_stderr is not None:
                    percent = abs(display_stderr / display_value) * 100
                    val_str = f"{display_value:.5g} ± {display_stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{display_value:.5g}"
                init_val_str = f"(init = {display_init:.5g})" if display_init is not None else ""
                if unit:
                    neg_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    neg_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        neg_params.append((r'$\chi^2$', f"{fit_result_neg.chisqr:.4g}"))
        neg_params.append((r'$R^2$', f"{fit_result_neg.rsquared:.4f}"))

    # --- New line-by-line layout logic for perfect alignment ---

    # Add text to figure
    y_pos = 0.01
    fig.text(0.5, y_pos + 0.215, f"{model_text}", ha='center', va='bottom', fontsize=25, family='monospace')
    text_y_start = y_pos + 0.195
    line_height = 0.018  # Adjusted line height in figure coordinates

    common_props = {'va': 'top', 'fontsize': 9, 'family': 'monospace'}

    # Column Headers - 動態根據可用分支調整
    header_y = text_y_start
    
    # 檢查哪些分支有有效的參數
    has_pos_params = len(pos_params) > 0 and pos_params[0][0] != "No fit"
    has_neg_params = len(neg_params) > 0 and neg_params[0][0] != "No fit"
    
    if has_pos_params and has_neg_params:
        # 雙欄佈局
        fig.text(0.02, header_y, f'$I_c^{{+}}$ Fit', ha='left', weight='bold', **common_props)
        fig.text(0.52, header_y, f'$I_c^{{-}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.02, 0.08, 0.10
        right_name_x, right_colon_x, right_value_x = 0.52, 0.58, 0.60
        layout_mode = 'dual'
    elif has_pos_params:
        # 只有正分支：單欄置中佈局
        fig.text(0.25, header_y, f'$I_c^{{+}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.25, 0.31, 0.33
        right_name_x, right_colon_x, right_value_x = None, None, None
        layout_mode = 'single_pos'
    elif has_neg_params:
        # 只有負分支：單欄置中佈局
        fig.text(0.25, header_y, f'$I_c^{{-}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.25, 0.31, 0.33
        right_name_x, right_colon_x, right_value_x = None, None, None
        layout_mode = 'single_neg'
    else:
        # 沒有有效參數
        fig.text(0.25, header_y, 'No Valid Fits Available', ha='left', weight='bold', **common_props)
        layout_mode = 'none'
    
    # Loop and plot each parameter line by line
    if layout_mode == 'dual':
        max_len = max(len(pos_params), len(neg_params))
        for i in range(max_len):
            current_y = text_y_start - (i + 1.5) * line_height # +1.5 to leave space after header

            # Left column (Positive)
            if i < len(pos_params):
                name, value = pos_params[i]
                fig.text(left_name_x, current_y, name, ha='left', **common_props)
                fig.text(left_colon_x, current_y, ':', ha='left', **common_props)
                fig.text(left_value_x, current_y, value, ha='left', **common_props)

            # Right column (Negative)
            if i < len(neg_params):
                name, value = neg_params[i]
                fig.text(right_name_x, current_y, name, ha='left', **common_props)
                fig.text(right_colon_x, current_y, ':', ha='left', **common_props)
                fig.text(right_value_x, current_y, value, ha='left', **common_props)
    elif layout_mode in ['single_pos', 'single_neg']:
        # 單欄佈局
        params_to_show = pos_params if layout_mode == 'single_pos' else neg_params
        for i, (name, value) in enumerate(params_to_show):
            current_y = text_y_start - (i + 1.5) * line_height
            fig.text(left_name_x, current_y, name, ha='left', **common_props)
            fig.text(left_colon_x, current_y, ':', ha='left', **common_props)
            fig.text(left_value_x, current_y, value, ha='left', **common_props)

    # Add a continuous vertical line between the columns (only for dual layout)
    if layout_mode == 'dual':
        num_param_lines = max(len(pos_params), len(neg_params))
        line_y_end = text_y_start - line_height # Start below header
        line_y_start = text_y_start - (num_param_lines + 2) * line_height
        line = lines.Line2D([0.5, 0.5], [line_y_start, line_y_end],
                            transform=fig.transFigure, color='grey', linestyle='-', linewidth=1)
        fig.add_artist(line)

    # Adjust rect to make space for legend and text
    fig.tight_layout(rect=(0, 0.30, 1, 0.99))

    output_path = output_dir / f"stage2_single_T_simple_run_{run_id}.png"
    plt.savefig(output_path, dpi=300)
    print(f"Stage 2 simple plot saved to {output_path}")
    plt.close(fig)

def visualize_stage2_results(analyzer: SQUIDAnalyzer, ds_cpr: xr.Dataset, results: Dict[str, Any], run_id: int, output_dir: Path):
    """Visualizes the results of the single-transparency model fitting (Stage 2)."""
    print("\nVisualizing Stage 2 single-T model results...")
    
    # Check if we have any data points
    if len(ds_cpr.coords['B_y_centered_μT']) == 0:
        print("   Warning: No data points available for visualization. Skipping Stage 2 plot.")
        return
    
    fig, ((ax, ax2), (ax3, ax4), (ax5, ax6), (ax_text, ax_empty)) = plt.subplots(4, 2, figsize=(18, 16), constrained_layout=True)
    ax_empty.axis('off')  # Hide the empty subplot
    ax6.axis('off')  # Hide the unused subplot
    
    x_data = ds_cpr.coords['B_y_centered_μT'].values
    x_fine = np.linspace(x_data.min(), x_data.max(), 500)
    all_plots = []

    for branch, color in [('pos', 'blue'), ('neg', 'red')]:
        y_data = ds_cpr[f'crit_current_{branch}'].values
        fit_result = results.get(f'single_T_fit_{branch}', None)
        if fit_result is None:
            print(f"   Skipping fit visualization for {branch} branch (no valid data)")
            continue
            
        p = fit_result.values

        # Plot data and fit
        p_data, = ax.plot(x_data, y_data, marker='o', linestyle='--', label=f'Ic{branch} data', color=color, alpha=0.6, markersize=4, linewidth=1)
        fit_curve = analyzer.single_T_model(x_fine, **p)
        p_fit, = ax.plot(x_fine, fit_curve, color=color, lw=2, label=f'Ic{branch} fit')
        all_plots.extend([p_data, p_fit])

        # Plot junction components (simple sine model)
        flux_phase = analyzer._flux_phase(x_fine, p['B_0_uT'], p['A_eff_um2'])
        I_sine = p['Ic1'] * np.sin(flux_phase)
        
        comp_color = 'darkcyan' if branch == 'pos' else 'darkred'
        p_comp, = ax.plot(x_fine, I_sine + p['Ic0'], color=comp_color, linestyle=':', lw=2, label=f'Ic0+Sine ({branch})')
        all_plots.append(p_comp)

        # Plot B0 and Ic0 annotations
        p_b0 = ax.axvline(x=p['B_0_uT'], color=color, linestyle='--', lw=1, alpha=0.8, label=f'B₀_{branch} = {p["B_0_uT"]:.2f} µT')
        p_ic0 = ax.axhline(y=p['Ic0'], color=color, linestyle='--', lw=1, alpha=0.8, label=f'Ic₀_{branch} = {p["Ic0"]:.2f} μA')
        all_plots.extend([p_b0, p_ic0])

        # --- emcee credible & posterior-predictive bands (Stage 2) ---
        emcee_key = f'single_T_fit_{branch}_emcee'
        if emcee_key in results:
            emc = results[emcee_key]
            xx = emc['x_eval']
            mb = emc['post_mean_bands']
            pb = emc['post_pred_bands']

            ax.fill_between(xx, mb['p2p5'], mb['p97p5'], alpha=0.12, label=f'95% cred (mean, {branch})', color=color)
            ax.fill_between(xx, mb['p16'], mb['p84'],   alpha=0.20, label=f'68% cred (mean, {branch})', color=color)
            ax.fill_between(xx, pb['p2p5'], pb['p97p5'], alpha=0.10, label=f'95% post. pred. ({branch})', color='gray')

    # --- Delta Ic on secondary axis ---
    ax2 = ax.twinx()
    
    # Calculate Delta Ic from data
    # Ensure that there are no NaNs which could cause issues with subtraction
    ic_pos_data = ds_cpr['crit_current_pos'].values
    ic_neg_data = ds_cpr['crit_current_neg'].values
    valid_mask = np.isfinite(ic_pos_data) & np.isfinite(ic_neg_data)
    delta_ic_data = ic_pos_data[valid_mask] - np.abs(ic_neg_data[valid_mask])
    
    # Calculate Delta Ic from fit - only if both fits are available
    fit_pos_result = results.get('single_T_fit_pos', None)
    fit_neg_result = results.get('single_T_fit_neg', None)
    if fit_pos_result is not None and fit_neg_result is not None:
        fit_pos_curve = analyzer.single_T_model(x_fine, **fit_pos_result.values)
        fit_neg_curve = analyzer.single_T_model(x_fine, **fit_neg_result.values)
        delta_ic_fit = fit_pos_curve - np.abs(fit_neg_curve)
        
        p_delta_data, = ax2.plot(x_data[valid_mask], delta_ic_data, color='green', linestyle='--', marker='.', markersize=5, label='ΔIc (data)')
        p_delta_fit, = ax2.plot(x_fine, delta_ic_fit, 'g-', lw=2, label='ΔIc (fit)')
        all_plots.extend([p_delta_data, p_delta_fit])
    else:
        print("   Warning: Cannot calculate ΔIc fit - missing fit results for one or both branches")
        # Plot only data
        p_delta_data, = ax2.plot(x_data[valid_mask], delta_ic_data, color='green', linestyle='--', marker='.', markersize=5, label='ΔIc (data)')
        all_plots.append(p_delta_data)
    
    ax2.set_ylabel('ΔIc = Ic+ - |Ic-| (μA)', color='g', fontsize=12)
    ax2.tick_params(axis='y', labelcolor='g')
    ax2.grid(True, linestyle=':', color='g', alpha=0.5)

    ax.set_title(f'Stage 2: Single-Transparency Model Fit (Run {run_id})', fontsize=14, pad=20)
    ax.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax.set_ylabel('Critical Current (µA)', fontsize=12)
    
    all_labels = [p.get_label() for p in all_plots]
    ax.legend(all_plots, all_labels, loc='upper center', bbox_to_anchor=(0.5, 1.15), ncol=7, fontsize=8)


    # --- 2. Parameter Text (ax_text) ---
    ax_text.axis('off')
    param_text = "Single-T Fit Parameters:\n\n"
    for branch in ['pos', 'neg']:
        fit_result = results.get(f'single_T_fit_{branch}', None)
        param_text += f'--- Ic{"+" if branch=="pos" else "-"} ---\n'
        if fit_result is None:
            param_text += "No valid fit results\n\n"
            continue
        for name, param in fit_result.params.items():
            param_text += f"{name}: {param.value:.3g} ± {param.stderr:.2g}\n" if param.stderr else f"{name}: {param.value:.3g}\n"
        param_text += f"R²: {fit_result.rsquared:.4f}\n\n"
    ax_text.text(0.05, 0.95, param_text, va='top', family='monospace', fontsize=9)

    # --- 3. Residuals Plot (ax3) ---
    for branch, color in [('pos', 'blue'), ('neg', 'red')]:
        fit_result = results.get(f'single_T_fit_{branch}', None)
        if fit_result is None:
            continue
        y_data = ds_cpr[f'crit_current_{branch}'].values
        y_model_at_data = analyzer.single_T_model(x_data, **fit_result.values)
        residuals = y_data - y_model_at_data
        ax3.plot(x_data, residuals, 'o', color=color, alpha=0.7, label=f'Residuals Ic{"+" if branch=="pos" else "-"}')
    ax3.axhline(0, color='k', linestyle='--')
    ax3.set_title('Fit Residuals', fontsize=14)
    ax3.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax3.set_ylabel('Data - Fit (µA)', fontsize=12)
    ax3.legend()
    ax3.grid(True, alpha=0.4)

    # --- 4. dI/dB Gradient Plot (ax4) ---
    for branch, color in [('pos', 'blue'), ('neg', 'red')]:
        fit_result = results.get(f'single_T_fit_{branch}', None)
        y_data = ds_cpr[f'crit_current_{branch}'].values
        
        # Calculate numerical gradient of data
        grad_data = np.gradient(y_data, x_data)
        ax4.plot(x_data, grad_data, 'o', color=color, alpha=0.5, label=f'dI/dB (data, {branch})')
        
        # Calculate analytical gradient of fit
        if fit_result is not None:
            y_fit_fine = analyzer.single_T_model(x_fine, **fit_result.values)
            grad_fit = np.gradient(y_fit_fine, x_fine)
            ax4.plot(x_fine, grad_fit, '-', color=color, lw=2, label=f'dI/dB (fit, {branch})')

    ax4.set_title(f'CPR Gradient (dI/dB) vs. Magnetic Field', fontsize=14)
    ax4.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax4.set_ylabel('dI/dB (μA/μT)', fontsize=12)
    ax4.legend()
    ax4.grid(True, alpha=0.4)

    # --- 5. Reconstructed Components Plot (ax5) ---
    for branch, color in [('pos', 'blue'), ('neg', 'red')]:
        fit_result = results.get(f'single_T_fit_{branch}', None)
        if fit_result is None:
            continue
        fit_params = fit_result.values
        # For simplified CPR model, we only have one CPR component
        flux_phase = analyzer._flux_phase(x_fine, fit_params['B_0_uT'], fit_params['A_eff_um2'])
        I_sine = fit_params['Ic1'] * np.sin(flux_phase)
        
        linestyle = '-' if branch == 'pos' else '--'
        ax5.plot(x_fine, I_sine, color='purple', linestyle=linestyle, lw=2, label=f'I_Sine ({branch})')
        # For simple model, total current includes linear term and sine term
        total_model = fit_params['Ic0'] + fit_params['r'] * flux_phase + I_sine
        ax5.plot(x_fine, total_model, color=color, linestyle=':', lw=2, alpha=0.7, label=f'Total Fit ({branch})')

    ax5.set_title('Reconstructed Junction Currents from Fit', fontsize=14)
    ax5.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax5.set_ylabel('Junction Current (µA)', fontsize=12)
    ax5.legend(ncol=2)
    ax5.grid(True, alpha=0.4)

    # Layout automatically managed by constrained_layout
    fig.suptitle(f'Stage 2: Single-Transparency Model Analysis - Run {run_id}', fontsize=16)
    output_path = output_dir / f"stage2_single_T_complete_run_{run_id}.png"
    plt.savefig(output_path, dpi=300)
    print(f"Stage 2 plot saved to {output_path}")
    plt.close(fig)

def visualize_stage3_simple_results(analyzer: SQUIDAnalyzer, ds_cpr: xr.Dataset, results: Dict[str, Any], run_id: int, output_dir: Path):
    """
    Creates a simple plot of the Stage 3 dual-transparency model results,
    similar to Stage 1, Stage 1.2, and Stage 2 simple design.
    """
    print("\nVisualizing Stage 3 dual-T model results (simple version)...")
    
    # Check if we have any data points
    if len(ds_cpr.coords['B_y_centered_μT']) == 0:
        print("   Warning: No data points available for visualization. Skipping Stage 3 simple plot.")
        return
    
    # Increase figure height to accommodate the text box
    fig, ax = plt.subplots(figsize=(12, 14))

    x_data = ds_cpr.coords['B_y_centered_μT'].values
    x_fine = np.linspace(x_data.min(), x_data.max(), 500)

    all_plots = []

    # 檢查哪些分支有數據
    has_pos_data = 'crit_current_pos' in ds_cpr.data_vars and len(ds_cpr['crit_current_pos'].values) > 0
    has_neg_data = 'crit_current_neg' in ds_cpr.data_vars and len(ds_cpr['crit_current_neg'].values) > 0
    
    for branch, color in [('pos', 'blue'), ('neg', 'red')]:
        if branch == 'pos' and not has_pos_data:
            continue
        if branch == 'neg' and not has_neg_data:
            continue
            
        # Plot raw data
        p_data, = ax.plot(x_data, ds_cpr[f'crit_current_{branch}'].values, marker='o', linestyle='--', label=f'Ic{"+" if branch == "pos" else "-"} data', color=color, alpha=0.7, markersize=4, linewidth=1)
        all_plots.append(p_data)

        # Plot fit and error band
        fit_result = results.get(f'dual_T_fit_{branch}', None)
        if fit_result is None:
            print(f"   Skipping fit visualization for {branch} branch (no valid data)")
            continue
        fit_curve = analyzer.dual_T_model(x_fine, **fit_result.params.valuesdict())

        # Plot the fit line
        p_fit, = ax.plot(x_fine, fit_curve, color=color, linestyle='-', lw=2, label=f'Ic{"+" if branch == "pos" else "-"} fit')
        all_plots.append(p_fit)

        # Calculate and plot 1-sigma error band
        try:
            dely = fit_result.eval_uncertainty(sigma=1, B_centered_uT=x_fine)
            ax.fill_between(x_fine, fit_curve - dely, fit_curve + dely, color=color, alpha=0.2)
        except Exception as e:
            print(f"Warning: Could not calculate or plot uncertainty for {branch} branch: {e}")

        # Plot junction components (simple sine model)
        p = fit_result.params.valuesdict()
        flux_phase = analyzer._flux_phase(x_fine, p['B_0_uT'], p['A_eff_um2'])
        I_sine = p['Ic1'] * np.sin(flux_phase)
        
        comp_color = 'darkcyan' if branch == 'pos' else 'darkred'
        p_comp, = ax.plot(x_fine, I_sine + p['Ic0'], color=comp_color, linestyle=':', lw=2, label=f'Ic0+Sine ({branch})')
        all_plots.append(p_comp)

        # Add B0 and Ic0 annotations
        p_b0 = ax.axvline(x=p['B_0_uT'], color=color, linestyle=':', lw=1.5, alpha=0.9, label=f'B₀_{branch} = {p["B_0_uT"]:.2f} µT')
        p_ic0 = ax.axhline(y=p['Ic0'], color=color, linestyle=':', lw=1.5, alpha=0.9, label=f'Ic₀_{branch} = {p["Ic0"]:.2f} μA')
        all_plots.extend([p_b0, p_ic0])

    # --- Delta Ic on secondary axis ---
    ax2 = ax.twinx()

    # 檢查是否有足夠的數據來計算 Delta Ic
    if has_pos_data and has_neg_data:
        ic_pos_data = ds_cpr['crit_current_pos'].values
        ic_neg_data = ds_cpr['crit_current_neg'].values
        valid_mask = np.isfinite(ic_pos_data) & np.isfinite(ic_neg_data)
        delta_ic_data = ic_pos_data[valid_mask] - np.abs(ic_neg_data[valid_mask])

        # Calculate Delta Ic from fit - only if both fits are available
        fit_pos_result = results.get('dual_T_fit_pos', None)
        fit_neg_result = results.get('dual_T_fit_neg', None)
        if fit_pos_result is not None and fit_neg_result is not None:
            fit_pos_curve = analyzer.dual_T_model(x_fine, **fit_pos_result.params.valuesdict())
            fit_neg_curve = analyzer.dual_T_model(x_fine, **fit_neg_result.params.valuesdict())
            delta_ic_fit = fit_pos_curve - np.abs(fit_neg_curve)
            
            p_delta_data, = ax2.plot(x_data[valid_mask], delta_ic_data, color='green', linestyle='--', marker='.', markersize=5, label='ΔIc (data)')
            p_delta_fit, = ax2.plot(x_fine, delta_ic_fit, 'g-', lw=2, label='ΔIc (fit)')
            all_plots.extend([p_delta_data, p_delta_fit])
        else:
            print("   Warning: Cannot calculate ΔIc fit - missing fit results for one or both branches")
            # Plot only data
            p_delta_data, = ax2.plot(x_data[valid_mask], delta_ic_data, color='green', linestyle='--', marker='.', markersize=5, label='ΔIc (data)')
            all_plots.append(p_delta_data)

        ax2.set_ylabel('ΔIc = Ic+ - |Ic-| (μA)', color='g', fontsize=12)
        ax2.tick_params(axis='y', labelcolor='g')
        ax2.grid(True, linestyle=':', color='g', alpha=0.5)
    else:
        print("   Warning: Cannot calculate ΔIc - missing data for one or both branches")
        # 隱藏第二個y軸
        ax2.set_visible(False)

    ax.set_title(f'Stage 3: Dual-Transparency Model Fit (Run {run_id})', fontsize=16)
    ax.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax.set_ylabel('Critical Current (µA)', fontsize=12)

    all_labels = [p.get_label() for p in all_plots]
    fig.legend(all_plots, all_labels, loc='upper center', bbox_to_anchor=(0.5, 0.355), ncol=7, fontsize=10)
    ax.grid(True, alpha=0.4)

    # --- Add model and parameter text below plot ---
    model_text = r"Model: $I_c = I_{c0} + r \cdot CPR(\phi, I_{c\_cpr}, T_{cpr})$, $\phi = \frac{2\pi A_{eff}(B-B_0)}{\Phi_0}$"

    param_details = {
        'Ic0': {'unit': 'μA', 'name': r'$I_{c0}$'},
        'r': {'unit': 'μA/T', 'name': r'$r$'},
        'Ic1': {'unit': 'μA', 'name': r'$I_{c1}$'},
        'Ic2': {'unit': 'μA', 'name': r'$I_{c2}$'},
        'T1': {'unit': '', 'name': r'$T_1$'},
        'T2': {'unit': '', 'name': r'$T_2$'},
        'phi_21': {'unit': '°', 'name': r'$\phi_{21}$'},
        'ratio_ic2_ic1': {'unit': '', 'name': r'$I_{c2}/I_{c1}$'},
        'ratio_t2_t1': {'unit': '', 'name': r'$T_2/T_1$'},
        'A_eff_um2': {'unit': 'μm²', 'name': r'$A_{eff}$'},
        'B_0_uT': {'unit': 'μT', 'name': r'$B_0$'}
    }

    pos_params = []
    neg_params = []

    # Process Positive Branch
    fit_result_pos = results.get('dual_T_fit_pos', None)
    if fit_result_pos is None:
        pos_params.append(("No fit", "No valid data for positive branch"))
    else:
        # Define parameter order for better readability - unified ordering
        priority_params = ['Ic0', 'Ic1', 'Ic2', 'ratio_ic2_ic1', 'phi_21', 'T1', 'T2', 'ratio_t2_t1', 'B_0_uT', 'A_eff_um2']

        # First, process priority parameters in specified order
        for param_name in priority_params:
            if param_name in fit_result_pos.params:
                param = fit_result_pos.params[param_name]
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')

                # Convert phi_21 from radians to degrees for display
                if param_name == 'phi_21':
                    display_value = param.value * 180 / np.pi
                    display_stderr = param.stderr * 180 / np.pi if param.stderr is not None else None
                    display_init = param.init_value * 180 / np.pi if hasattr(param, 'init_value') and param.init_value is not None else None
                else:
                    display_value = param.value
                    display_stderr = param.stderr
                    display_init = param.init_value if hasattr(param, 'init_value') else None

                if display_stderr is not None:
                    percent = abs(display_stderr / display_value) * 100
                    val_str = f"{display_value:.5g} ± {display_stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{display_value:.5g}"
                init_val_str = f"(init = {display_init:.5g})" if display_init is not None else ""
                if unit:
                    pos_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    pos_params.append((f"{p_name}", f"{val_str} {init_val_str}"))

        # Then, process remaining parameters in their original order
        processed_params = set(priority_params)
        for param_name, param in fit_result_pos.params.items():
            if param_name not in processed_params:
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')

                # Convert phi_21 from radians to degrees for display
                if param_name == 'phi_21':
                    display_value = param.value * 180 / np.pi
                    display_stderr = param.stderr * 180 / np.pi if param.stderr is not None else None
                    display_init = param.init_value * 180 / np.pi if hasattr(param, 'init_value') and param.init_value is not None else None
                else:
                    display_value = param.value
                    display_stderr = param.stderr
                    display_init = param.init_value if hasattr(param, 'init_value') else None

                if display_stderr is not None:
                    percent = abs(display_stderr / display_value) * 100
                    val_str = f"{display_value:.5g} ± {display_stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{display_value:.5g}"
                init_val_str = f"(init = {display_init:.5g})" if display_init is not None else ""
                if unit:
                    pos_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    pos_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        pos_params.append((r'$\chi^2$', f"{fit_result_pos.chisqr:.4g}"))
        pos_params.append((r'$R^2$', f"{fit_result_pos.rsquared:.4f}"))

    # Process Negative Branch
    fit_result_neg = results.get('dual_T_fit_neg', None)
    if fit_result_neg is None:
        neg_params.append(("No fit", "No valid data for negative branch"))
    else:
        # Define parameter order for better readability - unified ordering
        priority_params = ['Ic0', 'Ic1', 'Ic2', 'ratio_ic2_ic1', 'phi_21', 'T1', 'T2', 'ratio_t2_t1', 'B_0_uT', 'A_eff_um2']

        # First, process priority parameters in specified order
        for param_name in priority_params:
            if param_name in fit_result_neg.params:
                param = fit_result_neg.params[param_name]
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')

                # Convert phi_21 from radians to degrees for display
                if param_name == 'phi_21':
                    display_value = param.value * 180 / np.pi
                    display_stderr = param.stderr * 180 / np.pi if param.stderr is not None else None
                    display_init = param.init_value * 180 / np.pi if hasattr(param, 'init_value') and param.init_value is not None else None
                else:
                    display_value = param.value
                    display_stderr = param.stderr
                    display_init = param.init_value if hasattr(param, 'init_value') else None

                if display_stderr is not None:
                    percent = abs(display_stderr / display_value) * 100
                    val_str = f"{display_value:.5g} ± {display_stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{display_value:.5g}"
                init_val_str = f"(init = {display_init:.5g})" if display_init is not None else ""
                if unit:
                    neg_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    neg_params.append((f"{p_name}", f"{val_str} {init_val_str}"))

        # Then, process remaining parameters in their original order
        processed_params = set(priority_params)
        for param_name, param in fit_result_neg.params.items():
            if param_name not in processed_params:
                detail = param_details.get(param_name, {})
                p_name = detail.get('name', param_name)
                unit = detail.get('unit', '')

                # Convert phi_21 from radians to degrees for display
                if param_name == 'phi_21':
                    display_value = param.value * 180 / np.pi
                    display_stderr = param.stderr * 180 / np.pi if param.stderr is not None else None
                    display_init = param.init_value * 180 / np.pi if hasattr(param, 'init_value') and param.init_value is not None else None
                else:
                    display_value = param.value
                    display_stderr = param.stderr
                    display_init = param.init_value if hasattr(param, 'init_value') else None

                if display_stderr is not None:
                    percent = abs(display_stderr / display_value) * 100
                    val_str = f"{display_value:.5g} ± {display_stderr:.3g} ({percent:.2f}%)"
                else:
                    val_str = f"{display_value:.5g}"
                init_val_str = f"(init = {display_init:.5g})" if display_init is not None else ""
                if unit:
                    neg_params.append((f"{p_name} ({unit})", f"{val_str} {init_val_str}"))
                else:
                    neg_params.append((f"{p_name}", f"{val_str} {init_val_str}"))
        neg_params.append((r'$\chi^2$', f"{fit_result_neg.chisqr:.4g}"))
        neg_params.append((r'$R^2$', f"{fit_result_neg.rsquared:.4f}"))

    # --- New line-by-line layout logic for perfect alignment ---

    # Add text to figure
    y_pos = 0.01
    fig.text(0.5, y_pos + 0.265, f"{model_text}", ha='center', va='bottom', fontsize=25, family='monospace')
    text_y_start = y_pos + 0.245
    line_height = 0.018  # Adjusted line height in figure coordinates

    common_props = {'va': 'top', 'fontsize': 9, 'family': 'monospace'}

    # Column Headers - 動態根據可用分支調整
    header_y = text_y_start
    
    # 檢查哪些分支有有效的參數
    has_pos_params = len(pos_params) > 0 and pos_params[0][0] != "No fit"
    has_neg_params = len(neg_params) > 0 and neg_params[0][0] != "No fit"
    
    if has_pos_params and has_neg_params:
        # 雙欄佈局
        fig.text(0.02, header_y, f'$I_c^{{+}}$ Fit', ha='left', weight='bold', **common_props)
        fig.text(0.52, header_y, f'$I_c^{{-}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.02, 0.08, 0.10
        right_name_x, right_colon_x, right_value_x = 0.52, 0.58, 0.60
        layout_mode = 'dual'
    elif has_pos_params:
        # 只有正分支：單欄置中佈局
        fig.text(0.25, header_y, f'$I_c^{{+}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.25, 0.31, 0.33
        right_name_x, right_colon_x, right_value_x = None, None, None
        layout_mode = 'single_pos'
    elif has_neg_params:
        # 只有負分支：單欄置中佈局
        fig.text(0.25, header_y, f'$I_c^{{-}}$ Fit', ha='left', weight='bold', **common_props)
        left_name_x, left_colon_x, left_value_x = 0.25, 0.31, 0.33
        right_name_x, right_colon_x, right_value_x = None, None, None
        layout_mode = 'single_neg'
    else:
        # 沒有有效參數
        fig.text(0.25, header_y, 'No Valid Fits Available', ha='left', weight='bold', **common_props)
        layout_mode = 'none'
    
    # Loop and plot each parameter line by line
    if layout_mode == 'dual':
        max_len = max(len(pos_params), len(neg_params))
        for i in range(max_len):
            current_y = text_y_start - (i + 1.5) * line_height # +1.5 to leave space after header

            # Left column (Positive)
            if i < len(pos_params):
                name, value = pos_params[i]
                fig.text(left_name_x, current_y, name, ha='left', **common_props)
                fig.text(left_colon_x, current_y, ':', ha='left', **common_props)
                fig.text(left_value_x, current_y, value, ha='left', **common_props)

            # Right column (Negative)
            if i < len(neg_params):
                name, value = neg_params[i]
                fig.text(right_name_x, current_y, name, ha='left', **common_props)
                fig.text(right_colon_x, current_y, ':', ha='left', **common_props)
                fig.text(right_value_x, current_y, value, ha='left', **common_props)
    elif layout_mode in ['single_pos', 'single_neg']:
        # 單欄佈局
        params_to_show = pos_params if layout_mode == 'single_pos' else neg_params
        for i, (name, value) in enumerate(params_to_show):
            current_y = text_y_start - (i + 1.5) * line_height
            fig.text(left_name_x, current_y, name, ha='left', **common_props)
            fig.text(left_colon_x, current_y, ':', ha='left', **common_props)
            fig.text(left_value_x, current_y, value, ha='left', **common_props)

    # Add a continuous vertical line between the columns (only for dual layout)
    if layout_mode == 'dual':
        num_param_lines = max(len(pos_params), len(neg_params))
        line_y_end = text_y_start - line_height # Start below header
        line_y_start = text_y_start - (num_param_lines + 2) * line_height
        line = lines.Line2D([0.5, 0.5], [line_y_start, line_y_end],
                            transform=fig.transFigure, color='grey', linestyle='-', linewidth=1)
        fig.add_artist(line)

    # Adjust rect to make space for legend and text
    fig.tight_layout(rect=(0, 0.35, 1, 0.99))

    output_path = output_dir / f"stage3_dual_T_simple_run_{run_id}.png"
    plt.savefig(output_path, dpi=300)
    print(f"Stage 3 simple plot saved to {output_path}")
    plt.close(fig)

def visualize_stage3_results(analyzer: SQUIDAnalyzer, ds_cpr: xr.Dataset, results: Dict[str, Any], run_id: int, output_dir: Path):
    """Visualizes the results of the dual-transparency model fitting (Stage 3)."""
    print("\nVisualizing Stage 3 dual-T model results...")
    
    # Check if we have any data points
    if len(ds_cpr.coords['B_y_centered_μT']) == 0:
        print("   Warning: No data points available for visualization. Skipping Stage 3 plot.")
        return
    
    fig, ((ax, ax2), (ax3, ax4), (ax5, ax6), (ax_text, ax_empty)) = plt.subplots(4, 2, figsize=(18, 16), constrained_layout=True)
    ax_empty.axis('off')  # Hide the empty subplot
    ax6.axis('off')  # Hide the unused subplot

    x_data = ds_cpr.coords['B_y_centered_μT'].values
    x_fine = np.linspace(x_data.min(), x_data.max(), 500)
    all_plots = []

    for branch, color in [('pos', 'blue'), ('neg', 'red')]:
        y_data = ds_cpr[f'crit_current_{branch}'].values
        fit_result = results.get(f'dual_T_fit_{branch}', None)
        if fit_result is None:
            print(f"   Skipping fit visualization for {branch} branch (no valid data)")
            continue
            
        p = fit_result.values

        # Main plot: Data and fit
        ax.plot(x_data, y_data, marker='o', linestyle='--', label=f'Ic{"+" if branch == "pos" else "-"} data', color=color, alpha=0.7, markersize=4, linewidth=1)
        ax.plot(x_fine, analyzer.dual_T_model(x_fine, **p), color=color, linestyle='-', lw=2, label=f'Ic{"+" if branch == "pos" else "-"} fit')

        # Residuals
        fit_curve = analyzer.dual_T_model(x_data, **p)
        residuals = y_data - fit_curve
        ax2.plot(x_data, residuals, marker='.', linestyle='-', label=f'Residuals {"+" if branch == "pos" else "-"}', color=color, alpha=0.7)

        # Individual junction contributions (simple sine model)
        flux_phase = analyzer._flux_phase(x_fine, p['B_0_uT'], p['A_eff_um2'])
        I_sine = p['Ic1'] * np.sin(flux_phase)

        linestyle = '-' if branch == 'pos' else '--'
        ax3.plot(x_fine, I_sine, color='purple', linestyle=linestyle, lw=2, label=f'I_Sine ({branch})')
        ax3.plot(x_fine, p['r'] * flux_phase, color=color, linestyle=':', lw=2, alpha=0.7, label=f'r * φ ({branch})')

        # Parameter evolution plots (simplified)
        ax4.plot(x_fine, I_sine, color='purple', linestyle=linestyle, lw=2, label=f'I_Sine ({branch})')

        # Linear coefficient comparison (simplified)
        ax5.plot(x_fine, [p['r']] * len(x_fine), color='green', linestyle=linestyle, lw=2, label=f'r coeff ({branch})')

    # Set titles and labels
    ax.set_title('Data and Dual-T Fit', fontsize=14)
    ax.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax.set_ylabel('Critical Current (µA)', fontsize=12)

    ax2.set_title('Fit Residuals', fontsize=14)
    ax2.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax2.set_ylabel('Residuals (µA)', fontsize=12)

    ax3.set_title('Individual Junction Contributions', fontsize=14)
    ax3.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax3.set_ylabel('Junction Current (µA)', fontsize=12)

    ax4.set_title('Junction Current Comparison', fontsize=14)
    ax4.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax4.set_ylabel('Junction Current (µA)', fontsize=12)

    ax5.set_title('Transparency Parameters', fontsize=14)
    ax5.set_xlabel('Centered B-field (µT)', fontsize=12)
    ax5.set_ylabel('Transparency', fontsize=12)

    # Add legends
    for subplot in [ax, ax2, ax3, ax4, ax5]:
        subplot.legend(loc='best', fontsize=10)
        subplot.grid(True, alpha=0.4)

    # Parameter summary text
    param_text = []
    for branch in ['pos', 'neg']:
        fit_result = results.get(f'dual_T_fit_{branch}', None)
        if fit_result is None:
            param_text.append(f"{branch.upper()} Branch:")
            param_text.append("  No valid fit results")
            param_text.append("")
            continue
        p = fit_result.values
        param_text.append(f"{branch.upper()} Branch:")
        param_text.append(f"  Ic0: {p['Ic0']:.3f} μA")
        param_text.append(f"  Ic1: {p['Ic1']:.3f} μA")
        param_text.append(f"  r: {p['r']:.3f}")
        param_text.append(f"  r: {p['r']:.3f}")
        param_text.append(f"  A_eff: {p['A_eff_um2']:.1f} μm²")
        param_text.append("")

    ax_text.text(0.05, 0.95, '\n'.join(param_text), transform=ax_text.transAxes,
                 fontsize=10, verticalalignment='top', fontfamily='monospace')
    ax_text.set_title('Fit Parameters Summary', fontsize=14)
    ax_text.axis('off')

    fig.suptitle(f'Stage 3: Dual-Transparency Model Analysis - Run {run_id}', fontsize=16)
    output_path = output_dir / f"stage3_dual_T_complete_run_{run_id}.png"
    plt.savefig(output_path, dpi=300)
    print(f"Stage 3 complete plot saved to {output_path}")
    plt.close(fig)

def visualize_cpr_heatmap_with_fits(analyzer: SQUIDAnalyzer, df_raw: pd.DataFrame, cpr_df: pd.DataFrame, results: Dict[str, Any], run_id: int, output_dir: Path):
    """Generates an interactive dV/dI heatmap with fit overlays."""
    print("\nGenerating interactive heatmap...")
    
    # Check if we have any raw data
    if df_raw.empty:
        print("   Warning: No raw data available for heatmap. Skipping heatmap generation.")
        return
    
    # Check if we have valid fits - use .get() to handle missing keys
    dual_T_fit_pos = results.get('dual_T_fit_pos')
    dual_T_fit_neg = results.get('dual_T_fit_neg')
    
    if dual_T_fit_pos is None or dual_T_fit_neg is None:
        print("   Warning: No valid dual-T fits available for heatmap overlay. Skipping heatmap generation.")
        return
    
    dvdi_data = []
    for _, group in df_raw.groupby('y_field_centered_µT'):
        group = group.sort_values('appl_current_µA')
        i_vals, v_vals = group['appl_current_µA'].values, group['meas_voltage_µV'].values
        if len(i_vals) > 1:
            di = np.diff(i_vals)
            dvdi = np.divide(np.diff(v_vals), np.where(di == 0, 1e-9, di))
            i_mid = (i_vals[:-1] + i_vals[1:]) / 2
            for i, val in zip(i_mid, dvdi):
                dvdi_data.append({'B': group['y_field_centered_µT'].iloc[0], 'I': i, 'dVdI': val})
    df_dvdi = pd.DataFrame(dvdi_data)

    heatmap = hv.HeatMap(df_dvdi, kdims=['B', 'I'], vdims=['dVdI']).opts(
        title=f"dV/dI Heatmap (Run {run_id})", cmap='hot', colorbar=True,
        tools=['hover', 'pan', 'wheel_zoom'], width=800, height=500
    )

    x_fine = np.linspace(df_raw['y_field_centered_µT'].min(), df_raw['y_field_centered_µT'].max(), 500)
    
    # Only add fit curves if both fits are available
    overlays = [heatmap]
    
    if dual_T_fit_pos is not None:
        fit_pos = analyzer.dual_T_model(x_fine, **dual_T_fit_pos.values)
        curve_pos = hv.Curve((x_fine, fit_pos), 'B', 'I', label='Ic+ fit').opts(color='cyan')
        overlays.append(curve_pos)
    
    if dual_T_fit_neg is not None:
        fit_neg = analyzer.dual_T_model(x_fine, **dual_T_fit_neg.values)
        curve_neg = hv.Curve((x_fine, fit_neg), 'B', 'I', label='Ic- fit').opts(color='lime')
        overlays.append(curve_neg)

    # Create scatter for data points
    scatter_pos = hv.Scatter(cpr_df, 'B_y_centered_μT', 'crit_current_pos', label='Ic+ data').opts(color='blue', size=5)
    scatter_neg = hv.Scatter(cpr_df, 'B_y_centered_μT', 'crit_current_neg', label='Ic- data').opts(color='red', size=5)

    # Create curve for dashed lines connecting the data points
    line_pos = hv.Curve(cpr_df, 'B_y_centered_μT', 'crit_current_pos').opts(color='blue', line_dash='dashed', width=1)
    line_neg = hv.Curve(cpr_df, 'B_y_centered_μT', 'crit_current_neg').opts(color='red', line_dash='dashed', width=1)

    # Combine all overlays
    final_plot = heatmap * scatter_pos * line_pos * scatter_neg * line_neg
    for overlay in overlays[1:]:  # Skip the heatmap itself
        final_plot = final_plot * overlay
        
    final_plot = final_plot.opts(legend_position='top_right')
    
    output_path = output_dir / f"interactive_heatmap_run_{run_id}.html"
    hv.save(final_plot, output_path, backend='bokeh')
    print(f"Interactive heatmap saved to {output_path}")

def interpret_cpr_results(results: Dict[str, Any]):
    """
    Prints a physical interpretation of the CPR fit results, inspired by qcode_example_v2.py.
    Calculates total current, asymmetry, and flux period.
    """
    print("\n" + "="*60)
    print("Physical Interpretation of SQUID Properties (from Dual-T Model)")
    print("="*60)
    
    # Check if we have valid dual-T fits - use .get() to handle missing keys
    dual_T_fit_pos = results.get('dual_T_fit_pos')
    dual_T_fit_neg = results.get('dual_T_fit_neg')
    
    if dual_T_fit_pos is None or dual_T_fit_neg is None:
        print("   Warning: No valid dual-T fits available for interpretation.")
        print("   This may indicate insufficient or invalid CPR data for this run.")
        return

    def print_interpretation(branch_name: str, params: Parameters):
        p = params.valuesdict()
        Ic1, Ic2 = p['Ic1'], p['Ic2']
        A_eff_um2 = p['A_eff_um2']
        
        # Convert area to m^2 for flux calculation
        A_eff_m2 = A_eff_um2 / 1e12
        
        Ic_total = abs(Ic1) + abs(Ic2)
        
        PHI0 = 2.067833848e-15  # Wb (CODATA 2018)
        B_period_uT = (PHI0 / A_eff_m2) * 1e6 if A_eff_m2 > 0 else float('inf')
        
        print(f"\n--- {branch_name} ---")
        print(f"  Total Critical Current (Ic1+Ic2): {Ic_total:.3f} μA")
        print(f"  Effective Area (A_eff): {A_eff_um2:.2f} µm²")
        print(f"  Flux Period (from A_eff): {B_period_uT:.3f} μT/Φ₀")
        
    print_interpretation("Positive Branch (Ic+)", dual_T_fit_pos.params)
    print_interpretation("Negative Branch (Ic-)", dual_T_fit_neg.params)
    
    # Averaged properties
    params_pos = dual_T_fit_pos.params
    params_neg = dual_T_fit_neg.params
    A_eff_avg = (params_pos['A_eff_um2'].value + params_neg['A_eff_um2'].value) / 2
    T1_avg = (params_pos['T1'].value + params_neg['T1'].value) / 2
    T2_avg = (params_pos['T2'].value + params_neg['T2'].value) / 2
    
    print("\n--- Averaged Properties ---")
    print(f"  Average Effective Area (A_eff): {A_eff_avg:.2f} µm²")
    print(f"  Average Transparency 1 (T1): {T1_avg:.3f}")
    print(f"  Average Transparency 2 (T2): {T2_avg:.3f}")
    print("="*60)

def extract_fit_std(analysis_results: Dict[str, Any]) -> Dict[str, float]:
    """
    Extract standard deviations (stderr) from fitting results for key parameters.

    Returns a dictionary mapping parameter names to their standard deviations.
    Uses dual_T_fit results if available, otherwise falls back to single_T_fit.
    """
    std_dict = {}

    # Define the branches and parameter names to extract
    branches = ['pos', 'neg']
    param_names = ['Ic0', 'Ic1', 'A_eff_um2', 'T1', 'T2', 'B_0_uT', 'phi_21', 'ratio_ic2_ic1']

    for branch in branches:
        # Try dual_T_fit first (most complete), then single_T_fit
        fit_result = analysis_results.get(f'dual_T_fit_{branch}')
        if fit_result is None:
            fit_result = analysis_results.get(f'single_T_fit_{branch}')
        if fit_result is None:
            fit_result = analysis_results.get(f'single_junction_fit_{branch}')
        if fit_result is None:
            continue

        # Extract stderr for each parameter
        for param_name in param_names:
            if param_name in fit_result.params:
                param = fit_result.params[param_name]
                if param.stderr is not None:
                    # Use parameter name directly (without branch suffix) for table merging
                    std_dict[param_name] = param.stderr

    return std_dict


def build_cpr_overview_table(cpr_df: pd.DataFrame) -> pd.DataFrame:
    """
    Build CPR data overview table in Parameter Overview format.
    
    Args:
        cpr_df: CPR dataframe with columns like 'B_y_centered_μT', 'crit_current_pos', 'crit_current_neg'
        
    Returns:
        DataFrame with same format as Parameter Overview
    """
    rows = []
    
    # Process each column in CPR dataframe
    for col in cpr_df.columns:
        data = cpr_df[col].dropna()  # Remove NaN values
        
        if len(data) == 0:
            continue
            
        # Determine parameter info based on column name
        if 'B_y' in col or 'field' in col.lower():
            param_name = "B_field_centered"
            unit = "μT" if 'μT' in col or 'uT' in col else "T"
            param_type = "independent"
        elif 'pos' in col and ('current' in col or 'Ic' in col):
            param_name = "Ic_pos"  
            unit = "μA"
            param_type = "dependent"
        elif 'neg' in col and ('current' in col or 'Ic' in col):
            param_name = "Ic_neg"
            unit = "μA" 
            param_type = "dependent"
        else:
            # Skip unrecognized columns
            continue
            
        # Calculate statistics
        data_array = np.array(data)
        
        # Calculate step for independent parameters (like magnetic field)
        if param_type == "independent" and len(data) > 1:
            # Check if data is sorted, if not sort it for step calculation
            sorted_data = np.sort(data_array)
            unique_vals = np.unique(sorted_data)
            if len(unique_vals) > 1:
                diffs = np.diff(unique_vals)
                # Check if steps are approximately equal (within 1% tolerance)
                if np.allclose(diffs, diffs[0], rtol=1e-2):
                    step = float(diffs[0])
                else:
                    step = "varied"
            else:
                step = 0.0
        else:
            step = np.nan
            
        # Build row dict matching Parameter Overview format
        rows.append({
            "parameter": param_name,
            "type": param_type,
            "unit": unit,
            "points": len(data),
            "step": step,
            "start": float(data.iloc[0]),
            "end": float(data.iloc[-1]), 
            "min": float(data.min()),
            "max": float(data.max()),
            "range": float(data.max() - data.min()),
            "median": float(data.median()),
            "std": float(data.std())
        })
    
    df = pd.DataFrame(rows)
    if len(df) > 0:
        # Sort with independent parameters first, then dependent parameters
        type_order = {"independent": 0, "dependent": 1}
        df["type_order"] = df["type"].map(type_order)
        df = df.sort_values(by=["type_order", "parameter"]).drop("type_order", axis=1).reset_index(drop=True)
    
    return df

# --- 3. Main Execution Block ---

if __name__ == "__main__":
    # --- Command Line Argument Parsing ---
    parser = argparse.ArgumentParser(description='QCoDeS SQUID CPR Analysis Script')
    parser.add_argument('--RUN_ID', type=str,
                       help='RUN_ID to process. Can be single ID (e.g., 500), range (e.g., 500-511), or "ALL" for all CPR measurements')
    parser.add_argument('--drop', type=int, default=0,
                       help='Drop the first N data points from CPR data (default: 0)')
    parser.add_argument('--use-detrend', action='store_true',
                       help='Apply linear detrend to CPR data before analysis')
    args = parser.parse_args()

    # --- Configuration ---
    DB_PATH = "/Users/albert-mac/Code/GitHub/Qcodes/analysis/LT_2025-05-18_01.db"
    CATALOG_PATH = "/Users/albert-mac/Code/GitHub/Qcodes/analysis/data_processing/Catalog.csv"
    ANALYZER_SCRIPT_PATH = __file__
    BASE_DIR = Path('.')

    # Handle RUN_ID parsing
    if args.RUN_ID is None:
        print("Error: --RUN_ID is required. Use --RUN_ID ALL to process all CPR measurements, or specify a single ID/range.")
        exit(1)

    if args.RUN_ID.upper() == 'ALL':
        # Get all CPR run_ids from catalog
        run_ids = get_all_cpr_run_ids(CATALOG_PATH)
        if not run_ids:
            print("Error: No CPR measurements found in catalog.")
            exit(1)
    else:
        # Parse RUN_ID as single ID or range
        run_ids = parse_run_id(args.RUN_ID)

    # Setup output directories
    output_dirs = setup_output_directories(BASE_DIR, run_ids)

    print(f"📊 Processing {len(run_ids)} RUN_ID(s): {run_ids}")
    print(f"📁 Output directory: {output_dirs['base'].name}")

    # Process each RUN_ID
    for RUN_ID in run_ids:
        # Get measurement info
        measurement_type, current_type = get_measurement_info_from_catalog(RUN_ID, CATALOG_PATH)

        if measurement_type is None:
            print(f"⚠️  Skipping RUN_ID {RUN_ID}: not found in catalog")
            continue

        if measurement_type != 'CPR':
            print(f"⏭️  Skipping RUN_ID {RUN_ID}: {measurement_type} (only processing CPR)")
            continue

        print(f"\n▶ Processing RUN_ID {RUN_ID} (CPR, {current_type})")

        # --- Setup ---
        plt.style.use('tableau-colorblind10')
        hv.extension('bokeh')

        # Load the custom analyzer
        SQUIDAnalyzerClass = load_squid_analyzer(ANALYZER_SCRIPT_PATH)
        analyzer = SQUIDAnalyzerClass()

        # Load dataset
        db_name = Path(DB_PATH).name
        print(f"🔄 Loading RUN_ID {RUN_ID} from {db_name}")
        if not Path(DB_PATH).exists():
            print(f"❌ Database not found: {db_name}")
            continue
        qc.initialise_or_create_database_at(DB_PATH)
        ds = qc.load_by_id(run_id=RUN_ID)

        # Perform initial analysis and preprocessing
        # param_df = ParameterTableBuilder.build_basic(ds)
        print("\nParameter Overview:")
        # df_print_sci(param_df, digits=3)
        print("Dataset loaded successfully")

        df_dict = ds.to_pandas_dataframe_dict()
        df = df_dict['meas_voltage_K1'].reset_index()
        df.columns = ['y_field', 'appl_current', 'meas_voltage_K1']

        df['y_field_µT'] = df['y_field'] * 1e6
        df['appl_current_µA'] = df['appl_current'] * 1e6
        df['meas_voltage_µV'] = df['meas_voltage_K1'] * 1e6
        df['y_field_centered_µT'] = df['y_field_µT'] - df['y_field_µT'].median()

        # Extract CPR data
        ds_cpr, cpr_df = prepare_cpr_data_from_qcodes(ds, df, current_type, use_detrend=args.use_detrend)
        
        # Apply drop filter if specified
        if args.drop > 0:
            original_points = len(cpr_df)
            if args.drop >= original_points:
                print(f"⚠️  Warning: --drop {args.drop} exceeds total data points ({original_points}). No filtering applied.")
            else:
                print(f"🔄 Dropping first {args.drop} data points from CPR data...")
                print(f"   Original data points: {original_points}")
                
                # Filter pandas DataFrame
                cpr_df = cpr_df.iloc[args.drop:].reset_index(drop=True)
                
                # Filter xarray Dataset
                # Find the main dimension (usually related to B-field or data points)
                main_dims = []
                for dim in ds_cpr.dims:
                    if ds_cpr.sizes[dim] == original_points:
                        main_dims.append(dim)
                
                if main_dims:
                    # Apply drop to all dimensions that match the data point count
                    ds_cpr = ds_cpr.isel({dim: slice(args.drop, None) for dim in main_dims})
                    print(f"   Filtered data points: {len(cpr_df)}")
                else:
                    print(f"   Warning: Could not identify main dimension in ds_cpr for filtering")
        
        output_csv = output_dirs['csv'] / f"CPR_id{ds.run_id}.csv"
        cpr_df.to_csv(output_csv, sep='\t', index=False)
        print(f"📄 CPR data saved: {output_csv.name}")

        # CPR Data Overview
        cpr_overview_df = build_cpr_overview_table(cpr_df)
        if len(cpr_overview_df) > 0:
            print("\nCPR Data Overview:")
            df_print_sci(cpr_overview_df, digits=3)
        else:
            print("\nCPR Data Overview:")
            print("  No valid CPR data found for overview table.")

        # 驗證提取的資料是否符合 current_type 要求
        if not validate_cpr_data_against_current_type(ds_cpr, current_type):
            print(f"⚠️  CPR data validation failed for RUN_ID {RUN_ID}, continuing with caution...")

        # Use the current_type from catalog
        # (already shown in header)

        # Run the main analysis
        analysis_results = analyze_cpr_with_squid_analyzer(analyzer, ds_cpr, ds.run_id, output_dirs['base'], current_type)

        # Visualize and interpret the final results
        visualize_stage1_results(analyzer, ds_cpr, analysis_results, ds.run_id, output_dirs['plot'])
        visualize_stage1_2_results(analyzer, ds_cpr, analysis_results, ds.run_id, output_dirs['plot'])
        visualize_stage2_simple_results(analyzer, ds_cpr, analysis_results, ds.run_id, output_dirs['plot'])
        visualize_stage2_results(analyzer, ds_cpr, analysis_results, ds.run_id, output_dirs['plot'])
        visualize_stage3_simple_results(analyzer, ds_cpr, analysis_results, ds.run_id, output_dirs['plot'])
        visualize_stage3_results(analyzer, ds_cpr, analysis_results, ds.run_id, output_dirs['plot'])
        visualize_cpr_heatmap_with_fits(analyzer, df, cpr_df, analysis_results, ds.run_id, output_dirs['html'])
        interpret_cpr_results(analysis_results)

        # Update parameter table with fitting uncertainties
        fit_std_dict = extract_fit_std(analysis_results)
        if fit_std_dict:
            # Update std column for parameters that have fitting uncertainties
            for param_name, std_val in fit_std_dict.items():
                mask = param_df['parameter'] == param_name
                if mask.any():
                    param_df.loc[mask, 'std'] = std_val
            
            print("\nUpdated Parameter Overview (with fitting uncertainties):")
            df_print_sci(param_df, digits=3)
        else:
            print("\nNo fitting uncertainties available to update Parameter Overview.")

        # Save analysis results to JSON
        results_json = output_dirs['json'] / f"analysis_results_{ds.run_id}.json"
        with open(results_json, 'w') as f:
            json.dump(analysis_results, f, indent=2, default=str)
        print(f"💾 Analysis results saved: {results_json.name}")

        print(f"✅ Analysis complete for RUN_ID {RUN_ID}")

    print(f"\n🎉 Batch processing complete for {len(run_ids)} RUN_ID(s)")
    print(f"📁 All results saved in: {output_dirs['base'].name}")

def visualize_mcmc_corner_plots(results: Dict[str, Any], run_id: int, output_dir: Path):
    """
    MCMC corner plots disabled - MCMC functionality has been removed.
    """
    print("\n--- MCMC Corner Plots Disabled ---")
    print("MCMC functionality has been removed from this version.")
    return

def visualize_dvdi_heatmap(df_raw: pd.DataFrame, run_id: int, output_dir: Path):
    """
    Creates and saves a dV/dI heatmap visualization.
    """
    print("\nGenerating dV/dI heatmap...")
    
    dvdi_data = []
    for _, group in df_raw.groupby('y_field_centered_µT'):
        group = group.sort_values('appl_current_µA')
        i_vals, v_vals = group['appl_current_µA'].values, group['meas_voltage_µV'].values
        if len(i_vals) > 1:
            di = np.diff(i_vals)
            dvdi = np.divide(np.diff(v_vals), np.where(di == 0, 1e-9, di))
            i_mid = (i_vals[:-1] + i_vals[1:]) / 2
            for i, val in zip(i_mid, dvdi):
                dvdi_data.append({'B': group['y_field_centered_µT'].iloc[0], 'I': i, 'dVdI': val})
    df_dvdi = pd.DataFrame(dvdi_data)

    heatmap = hv.HeatMap(df_dvdi, kdims=['B', 'I'], vdims=['dVdI']).opts(
        title=f"dV/dI Heatmap (Run {run_id})", cmap='hot', colorbar=True,
        tools=['hover', 'pan', 'wheel_zoom'], width=800, height=500
    )

    output_path = output_dir / f"dvdi_heatmap_run_{run_id}.html"
    hv.save(heatmap, output_path, backend='bokeh')
    print(f"dV/dI heatmap saved to {output_path}")