第4章: ラマン分光法 (Raman Spectroscopy)

4.3.2 励起波長の選択

4.4 表面増強ラマン分光法（SERS）

4.4.1 SERSの原理

表面増強ラマン分光法（Surface-Enhanced Raman Spectroscopy, SERS）は、金や銀などの金属ナノ構造表面で生じるラマン信号の巨大な増強効果を利用した分析手法です。増強因子は10⁶から10¹⁴に達し、単一分子検出も可能です。

4.4.2 SERS基板の種類

コード例3: SERS増強因子のシミュレーション

import numpy as np
import matplotlib.pyplot as plt

def sers_enhancement_factor(distance, nanoparticle_radius=50, wavelength=532):
    """
    ナノ粒子表面からの距離に対するSERS増強因子を計算

    Parameters:
    -----------
    distance : array
        表面からの距離 (nm)
    nanoparticle_radius : float
        ナノ粒子半径 (nm)
    wavelength : float
        励起波長 (nm)

    Returns:
    --------
    ef : array
        電磁気的増強因子
    """
    # 簡略化モデル: |E/E0|^4 = [(R/(R+d))^3]^4
    # R: 粒子半径, d: 表面からの距離
    R = nanoparticle_radius

    # 電場増強（双極子近似）
    field_enhancement = (R / (R + distance))**3

    # SERS増強因子（|E|^4に比例）
    ef = field_enhancement**4

    return ef

def simulate_sers_spectrum():
    """
    通常ラマンとSERSスペクトルの比較
    """
    wavenumber = np.linspace(400, 2000, 1600)

    # 4-アミノチオフェノール（4-ATP）の代表的ピーク
    peak_positions = [1077, 1142, 1390, 1435, 1575]  # cm^-1
    peak_widths = [15, 12, 10, 12, 15]
    peak_intensities = [0.6, 0.4, 0.3, 0.5, 1.0]

    # 通常ラマンスペクトル
    normal_raman = np.zeros_like(wavenumber)
    for pos, width, intensity in zip(peak_positions, peak_widths, peak_intensities):
        normal_raman += lorentzian_peak(wavenumber, intensity, pos, width)

    # ノイズ追加（通常ラマンは信号が弱いためノイズが目立つ）
    normal_raman += np.random.normal(0, 0.05, len(wavenumber))
    normal_raman = np.maximum(normal_raman, 0)

    # SERSスペクトル（10^6倍増強、ノイズ比が大幅に改善）
    enhancement = 1e6
    sers_spectrum = normal_raman * np.sqrt(enhancement)  # 可視化のためスケール調整
    sers_spectrum += np.random.normal(0, 0.5, len(wavenumber))

    # プロット
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))

    # 通常ラマンスペクトル
    axes[0, 0].plot(wavenumber, normal_raman, 'b-', linewidth=1.5)
    axes[0, 0].set_xlabel('ラマンシフト (cm$^{-1}$)', fontsize=12)
    axes[0, 0].set_ylabel('ラマン強度 (a.u.)', fontsize=12)
    axes[0, 0].set_title('通常ラマンスペクトル（4-ATP）', fontsize=14, fontweight='bold')
    axes[0, 0].grid(alpha=0.3)

    # SERSスペクトル
    axes[0, 1].plot(wavenumber, sers_spectrum, 'r-', linewidth=1.5)
    axes[0, 1].set_xlabel('ラマンシフト (cm$^{-1}$)', fontsize=12)
    axes[0, 1].set_ylabel('SERS強度 (a.u.)', fontsize=12)
    axes[0, 1].set_title('SERSスペクトル（4-ATP on Ag）', fontsize=14, fontweight='bold')
    axes[0, 1].grid(alpha=0.3)

    # 距離依存性
    distances = np.linspace(0.1, 30, 100)  # nm
    radii = [20, 40, 60, 80]  # nm

    for R in radii:
        ef = sers_enhancement_factor(distances, nanoparticle_radius=R)
        axes[1, 0].plot(distances, ef, linewidth=2, label=f'R = {R} nm')

    axes[1, 0].set_xlabel('表面からの距離 (nm)', fontsize=12)
    axes[1, 0].set_ylabel('増強因子 EF', fontsize=12)
    axes[1, 0].set_title('SERS増強因子の距離依存性', fontsize=14, fontweight='bold')
    axes[1, 0].legend()
    axes[1, 0].set_yscale('log')
    axes[1, 0].grid(alpha=0.3)

    # ホットスポットの模式図（2粒子間ギャップ）
    theta = np.linspace(0, 2*np.pi, 100)
    R = 1
    gap = 0.1

    # 左粒子
    x1 = R * np.cos(theta) - R - gap/2
    y1 = R * np.sin(theta)

    # 右粒子
    x2 = R * np.cos(theta) + R + gap/2
    y2 = R * np.sin(theta)

    axes[1, 1].fill(x1, y1, color='gold', alpha=0.8, label='Au nanoparticle')
    axes[1, 1].fill(x2, y2, color='gold', alpha=0.8)
    axes[1, 1].scatter([0], [0], s=200, c='red', marker='*', zorder=5, label='Hot spot')
    axes[1, 1].set_xlim(-3, 3)
    axes[1, 1].set_ylim(-2, 2)
    axes[1, 1].set_aspect('equal')
    axes[1, 1].set_xlabel('x (a.u.)', fontsize=12)
    axes[1, 1].set_ylabel('y (a.u.)', fontsize=12)
    axes[1, 1].set_title('ナノギャップ ホットスポット', fontsize=14, fontweight='bold')
    axes[1, 1].legend(loc='upper right')
    axes[1, 1].grid(alpha=0.3)

    plt.tight_layout()
    plt.show()

# 実行
simulate_sers_spectrum()

4.5 炭素材料のラマン分光分析

4.5.1 炭素材料の特徴的ラマンバンド

炭素材料（グラファイト、グラフェン、カーボンナノチューブ、ダイヤモンドなど）はラマン分光法で特徴的なスペクトルを示し、構造評価に広く利用されています。

4.5.2 D/G比による欠陥評価

コード例4: グラフェンのラマンスペクトル解析とD/G比計算

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.signal import find_peaks

def gaussian_peak(x, amplitude, center, width):
    """ガウス関数によるピーク形状"""
    return amplitude * np.exp(-((x - center)**2) / (2 * width**2))

def multi_peak_model(x, *params):
    """
    複数ピークモデル（Dバンド、Gバンド、D'バンド、2Dバンド）
    params: [A_D, c_D, w_D, A_G, c_G, w_G, A_D', c_D', w_D', A_2D, c_2D, w_2D]
    """
    n_peaks = len(params) // 3
    result = np.zeros_like(x)
    for i in range(n_peaks):
        A, c, w = params[3*i:3*i+3]
        result += gaussian_peak(x, A, c, w)
    return result

def analyze_graphene_raman(wavenumber, intensity, peak_names=['D', 'G', "D'", '2D']):
    """
    グラフェンラマンスペクトルの解析

    Parameters:
    -----------
    wavenumber : array
        ラマンシフト (cm^-1)
    intensity : array
        ラマン強度
    peak_names : list
        フィットするピーク名

    Returns:
    --------
    results : dict
        各ピークのパラメータと評価指標
    """
    # 初期推定値（グラフェンの典型値）
    initial_params = {
        'D': [100, 1350, 30],      # [振幅, 中心, 幅]
        'G': [500, 1580, 15],
        "D'": [50, 1620, 10],
        '2D': [400, 2700, 25]
    }

    # フィッティング用パラメータ
    p0 = []
    for name in peak_names:
        p0.extend(initial_params[name])

    # バウンド設定
    lower_bounds = []
    upper_bounds = []
    for name in peak_names:
        center = initial_params[name][1]
        lower_bounds.extend([0, center - 50, 5])
        upper_bounds.extend([np.max(intensity)*2, center + 50, 100])

    try:
        popt, pcov = curve_fit(multi_peak_model, wavenumber, intensity,
                               p0=p0, bounds=(lower_bounds, upper_bounds), maxfev=10000)
    except RuntimeError:
        print("フィッティングが収束しませんでした。")
        return None

    # 結果の整理
    results = {'peaks': {}}
    for i, name in enumerate(peak_names):
        A, c, w = popt[3*i:3*i+3]
        # ピーク面積（ガウスの場合: A * sqrt(2*pi) * w）
        area = A * np.sqrt(2 * np.pi) * w
        results['peaks'][name] = {
            'amplitude': A,
            'center': c,
            'width_fwhm': 2.355 * w,  # FWHM = 2.355 * sigma
            'area': area
        }

    # D/G比の計算
    if 'D' in results['peaks'] and 'G' in results['peaks']:
        results['D_G_ratio'] = results['peaks']['D']['area'] / results['peaks']['G']['area']
        results['D_G_intensity_ratio'] = results['peaks']['D']['amplitude'] / results['peaks']['G']['amplitude']

    # 2D/G比の計算
    if '2D' in results['peaks'] and 'G' in results['peaks']:
        results['2D_G_ratio'] = results['peaks']['2D']['area'] / results['peaks']['G']['area']
        results['2D_G_intensity_ratio'] = results['peaks']['2D']['amplitude'] / results['peaks']['G']['amplitude']

    # フィッティング曲線
    results['fitted_curve'] = multi_peak_model(wavenumber, *popt)
    results['popt'] = popt

    return results

def simulate_graphene_spectrum(defect_level='low'):
    """
    グラフェンラマンスペクトルのシミュレーション

    Parameters:
    -----------
    defect_level : str
        'pristine', 'low', 'medium', 'high'
    """
    wavenumber = np.linspace(1100, 3100, 2000)

    # 欠陥レベルに応じたD/G比
    defect_params = {
        'pristine': {'D': 0, 'G': 500, 'D_prime': 0, '2D': 1200},
        'low': {'D': 50, 'G': 500, 'D_prime': 20, '2D': 1000},
        'medium': {'D': 200, 'G': 500, 'D_prime': 80, '2D': 600},
        'high': {'D': 400, 'G': 500, 'D_prime': 150, '2D': 300}
    }

    params = defect_params[defect_level]

    # スペクトル生成
    D_band = gaussian_peak(wavenumber, params['D'], 1350, 25)
    G_band = gaussian_peak(wavenumber, params['G'], 1580, 12)
    D_prime = gaussian_peak(wavenumber, params['D_prime'], 1620, 10)
    two_D = gaussian_peak(wavenumber, params['2D'], 2700, 25)

    spectrum = D_band + G_band + D_prime + two_D
    noise = np.random.normal(0, 10, len(wavenumber))
    spectrum_noisy = spectrum + noise

    return wavenumber, spectrum_noisy

# 異なる欠陥レベルのグラフェンスペクトルを比較
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
defect_levels = ['pristine', 'low', 'medium', 'high']
titles = ['高品質グラフェン', '低欠陥グラフェン', '中欠陥グラフェン', '高欠陥グラフェン']

for ax, level, title in zip(axes.flatten(), defect_levels, titles):
    wavenumber, spectrum = simulate_graphene_spectrum(level)

    ax.plot(wavenumber, spectrum, 'b-', linewidth=1.5)
    ax.axvline(1350, color='red', linestyle='--', alpha=0.5, label='D band')
    ax.axvline(1580, color='green', linestyle='--', alpha=0.5, label='G band')
    ax.axvline(2700, color='purple', linestyle='--', alpha=0.5, label='2D band')
    ax.set_xlabel('ラマンシフト (cm$^{-1}$)', fontsize=11)
    ax.set_ylabel('ラマン強度 (a.u.)', fontsize=11)
    ax.set_title(title, fontsize=12, fontweight='bold')
    ax.legend(loc='upper right', fontsize=9)
    ax.grid(alpha=0.3)

plt.tight_layout()
plt.show()

# 定量解析の例
print("グラフェンラマンスペクトル解析結果:")
print("-" * 60)
for level, title in zip(defect_levels, titles):
    wavenumber, spectrum = simulate_graphene_spectrum(level)
    results = analyze_graphene_raman(wavenumber, spectrum)
    if results:
        print(f"\n{title}:")
        print(f"  D/G強度比: {results.get('D_G_intensity_ratio', 0):.3f}")
        print(f"  2D/G強度比: {results.get('2D_G_intensity_ratio', 0):.3f}")

4.5.3 カーボンナノチューブのRBMモード

単層カーボンナノチューブ（SWCNT）は、径方向呼吸振動（Radial Breathing Mode, RBM）と呼ばれる特徴的な低波数モードを示します。RBM振動数からナノチューブの直径を決定できます。

コード例5: CNTのRBMモード解析と直径決定

import numpy as np
import matplotlib.pyplot as plt

def rbm_to_diameter(omega_rbm, A=234, B=10):
    """
    RBM振動数からCNT直径を計算

    Parameters:
    -----------
    omega_rbm : float or array
        RBM振動数 (cm^-1)
    A : float
        定数 (cm^-1 nm)、デフォルト 234
    B : float
        定数 (cm^-1)、デフォルト 10

    Returns:
    --------
    diameter : float or array
        CNT直径 (nm)
    """
    return A / (omega_rbm - B)

def diameter_to_rbm(diameter, A=234, B=10):
    """
    CNT直径からRBM振動数を計算
    """
    return A / diameter + B

def chirality_to_diameter(n, m, a_cc=0.142):
    """
    カイラリティ(n, m)からCNT直径を計算

    Parameters:
    -----------
    n, m : int
        カイラル指数
    a_cc : float
        C-C結合長 (nm)、デフォルト 0.142 nm

    Returns:
    --------
    diameter : float
        CNT直径 (nm)
    """
    a = np.sqrt(3) * a_cc  # 格子定数
    diameter = a * np.sqrt(n**2 + n*m + m**2) / np.pi
    return diameter

def simulate_swcnt_raman():
    """
    SWCNTのRBMスペクトルをシミュレート
    """
    wavenumber_rbm = np.linspace(100, 350, 500)
    wavenumber_high = np.linspace(1200, 1700, 500)

    # 混合サンプル中の異なる直径のSWCNT
    # (n, m)カイラリティの例
    chiralities = [(10, 10), (8, 7), (6, 5), (11, 0)]

    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    # RBMスペクトル
    rbm_spectrum = np.zeros_like(wavenumber_rbm)
    cnt_info = []

    for n, m in chiralities:
        d = chirality_to_diameter(n, m)
        omega = diameter_to_rbm(d)
        # 半導体か金属かを判定
        mod_value = (n - m) % 3
        if mod_value == 0:
            cnt_type = '金属'
            color = 'red'
        else:
            cnt_type = '半導体'
            color = 'blue'

        amplitude = np.random.uniform(0.5, 1.0)
        rbm_spectrum += gaussian_peak(wavenumber_rbm, amplitude, omega, 3)
        cnt_info.append({
            'chirality': (n, m),
            'diameter': d,
            'rbm': omega,
            'type': cnt_type
        })

    noise = np.random.normal(0, 0.02, len(wavenumber_rbm))
    rbm_spectrum += noise

    axes[0].plot(wavenumber_rbm, rbm_spectrum, 'b-', linewidth=1.5)
    for info in cnt_info:
        axes[0].axvline(info['rbm'], color='red' if info['type']=='金属' else 'blue',
                       linestyle='--', alpha=0.7)
        axes[0].text(info['rbm'], 1.1, f"({info['chirality'][0]},{info['chirality'][1]})",
                    rotation=90, fontsize=9, ha='center')

    axes[0].set_xlabel('ラマンシフト (cm$^{-1}$)', fontsize=12)
    axes[0].set_ylabel('ラマン強度 (a.u.)', fontsize=12)
    axes[0].set_title('SWCNT RBMスペクトル', fontsize=14, fontweight='bold')
    axes[0].grid(alpha=0.3)
    axes[0].set_ylim(0, 1.5)

    # 高波数領域（G+, G-）
    G_plus = gaussian_peak(wavenumber_high, 1.0, 1590, 10)
    G_minus = gaussian_peak(wavenumber_high, 0.4, 1570, 15)  # 金属CNTでブロード
    D_band = gaussian_peak(wavenumber_high, 0.15, 1350, 25)
    high_spectrum = G_plus + G_minus + D_band + np.random.normal(0, 0.02, len(wavenumber_high))

    axes[1].plot(wavenumber_high, high_spectrum, 'g-', linewidth=1.5)
    axes[1].axvline(1590, color='red', linestyle='--', alpha=0.7, label='G$^+$')
    axes[1].axvline(1570, color='blue', linestyle='--', alpha=0.7, label='G$^-$')
    axes[1].axvline(1350, color='purple', linestyle='--', alpha=0.7, label='D band')
    axes[1].set_xlabel('ラマンシフト (cm$^{-1}$)', fontsize=12)
    axes[1].set_ylabel('ラマン強度 (a.u.)', fontsize=12)
    axes[1].set_title('SWCNT D/Gバンド', fontsize=14, fontweight='bold')
    axes[1].legend()
    axes[1].grid(alpha=0.3)

    # RBM-直径関係
    diameters = np.linspace(0.5, 2.5, 100)
    rbm_frequencies = diameter_to_rbm(diameters)

    axes[2].plot(diameters, rbm_frequencies, 'b-', linewidth=2)
    for info in cnt_info:
        color = 'red' if info['type'] == '金属' else 'blue'
        axes[2].scatter(info['diameter'], info['rbm'], s=100, c=color,
                       edgecolor='black', linewidth=1.5, zorder=5)
        axes[2].annotate(f"({info['chirality'][0]},{info['chirality'][1]})",
                        (info['diameter'], info['rbm']),
                        textcoords="offset points", xytext=(10, 5), fontsize=9)

    axes[2].set_xlabel('CNT直径 (nm)', fontsize=12)
    axes[2].set_ylabel('RBM振動数 (cm$^{-1}$)', fontsize=12)
    axes[2].set_title('RBM振動数と直径の関係', fontsize=14, fontweight='bold')
    axes[2].grid(alpha=0.3)

    plt.tight_layout()
    plt.show()

    # 結果表示
    print("検出されたSWCNT:")
    print("-" * 60)
    print(f"{'カイラリティ':<15} {'直径 (nm)':<12} {'RBM (cm^-1)':<12} {'タイプ':<10}")
    print("-" * 60)
    for info in cnt_info:
        print(f"({info['chirality'][0]:2d},{info['chirality'][1]:2d}){'':<10} "
              f"{info['diameter']:<12.3f} {info['rbm']:<12.1f} {info['type']:<10}")

# 実行
simulate_swcnt_raman()

4.6 ラマンスペクトルの前処理とピークフィッティング

4.6.1 ベースライン補正

ラマンスペクトルには、蛍光バックグラウンドや装置由来のベースラインが含まれることがあります。定量解析の前にベースライン補正が必要です。

コード例6: ベースライン補正とスペクトル前処理

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
from scipy import sparse
from scipy.sparse.linalg import spsolve

def baseline_als(y, lam=1e5, p=0.01, niter=10):
    """
    非対称最小二乗法によるベースライン推定（ALS法）

    Parameters:
    -----------
    y : array
        入力スペクトル
    lam : float
        平滑化パラメータ（大きいほど滑らか）
    p : float
        非対称パラメータ（0 < p < 1）
    niter : int
        反復回数

    Returns:
    --------
    baseline : array
        推定されたベースライン
    """
    L = len(y)
    D = sparse.diags([1, -2, 1], [0, -1, -2], shape=(L, L-2))
    D = lam * D.dot(D.T)
    w = np.ones(L)
    W = sparse.spdiags(w, 0, L, L)

    for i in range(niter):
        W.setdiag(w)
        Z = W + D
        baseline = spsolve(Z, w * y)
        w = p * (y > baseline) + (1 - p) * (y < baseline)

    return baseline

def preprocess_raman_spectrum(wavenumber, raw_spectrum, smooth_window=11, smooth_order=3):
    """
    ラマンスペクトルの前処理パイプライン

    Parameters:
    -----------
    wavenumber : array
        ラマンシフト
    raw_spectrum : array
        生スペクトル
    smooth_window : int
        Savitzky-Golayフィルタの窓幅
    smooth_order : int
        多項式次数

    Returns:
    --------
    processed : dict
        前処理後のスペクトルと中間結果
    """
    # 1. スムージング
    smoothed = savgol_filter(raw_spectrum, smooth_window, smooth_order)

    # 2. ベースライン推定（ALS法）
    baseline = baseline_als(smoothed, lam=1e6, p=0.001)

    # 3. ベースライン補正
    corrected = smoothed - baseline

    # 4. 正規化（最大値で規格化）
    normalized = corrected / np.max(corrected)

    return {
        'wavenumber': wavenumber,
        'raw': raw_spectrum,
        'smoothed': smoothed,
        'baseline': baseline,
        'corrected': corrected,
        'normalized': normalized
    }

def demonstrate_preprocessing():
    """
    前処理の効果をデモンストレーション
    """
    # 蛍光バックグラウンドを含むシミュレーションデータ
    wavenumber = np.linspace(800, 2000, 1200)

    # 真のラマンピーク（グラファイト様）
    D_band = gaussian_peak(wavenumber, 100, 1350, 30)
    G_band = gaussian_peak(wavenumber, 300, 1580, 15)

    # 蛍光バックグラウンド（広い帯域の曲線）
    fluorescence = 50 * np.exp(-((wavenumber - 1400)**2) / (2 * 300**2))

    # 傾斜ベースライン
    slope_baseline = 0.05 * (wavenumber - 800)

    # 合成スペクトル
    raw = D_band + G_band + fluorescence + slope_baseline

    # ノイズ追加
    noise = np.random.normal(0, 5, len(wavenumber))
    raw_noisy = raw + noise

    # 前処理実行
    result = preprocess_raman_spectrum(wavenumber, raw_noisy)

    # プロット
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))

    # 生データとスムージング
    axes[0, 0].plot(wavenumber, raw_noisy, 'gray', alpha=0.5, linewidth=0.5, label='生データ')
    axes[0, 0].plot(wavenumber, result['smoothed'], 'b-', linewidth=1.5, label='スムージング後')
    axes[0, 0].set_xlabel('ラマンシフト (cm$^{-1}$)', fontsize=12)
    axes[0, 0].set_ylabel('強度 (a.u.)', fontsize=12)
    axes[0, 0].set_title('Step 1: ノイズ除去（Savitzky-Golay）', fontsize=14, fontweight='bold')
    axes[0, 0].legend()
    axes[0, 0].grid(alpha=0.3)

    # ベースライン推定
    axes[0, 1].plot(wavenumber, result['smoothed'], 'b-', linewidth=1.5, label='スムージング後')
    axes[0, 1].plot(wavenumber, result['baseline'], 'r-', linewidth=2, label='推定ベースライン')
    axes[0, 1].set_xlabel('ラマンシフト (cm$^{-1}$)', fontsize=12)
    axes[0, 1].set_ylabel('強度 (a.u.)', fontsize=12)
    axes[0, 1].set_title('Step 2: ベースライン推定（ALS法）', fontsize=14, fontweight='bold')
    axes[0, 1].legend()
    axes[0, 1].grid(alpha=0.3)

    # ベースライン補正後
    axes[1, 0].plot(wavenumber, result['corrected'], 'g-', linewidth=1.5)
    axes[1, 0].axhline(0, color='black', linestyle='--', alpha=0.5)
    axes[1, 0].set_xlabel('ラマンシフト (cm$^{-1}$)', fontsize=12)
    axes[1, 0].set_ylabel('強度 (a.u.)', fontsize=12)
    axes[1, 0].set_title('Step 3: ベースライン補正後', fontsize=14, fontweight='bold')
    axes[1, 0].grid(alpha=0.3)

    # 正規化後
    axes[1, 1].plot(wavenumber, result['normalized'], 'purple', linewidth=1.5)
    axes[1, 1].axvline(1350, color='red', linestyle='--', alpha=0.7, label='D band')
    axes[1, 1].axvline(1580, color='green', linestyle='--', alpha=0.7, label='G band')
    axes[1, 1].set_xlabel('ラマンシフト (cm$^{-1}$)', fontsize=12)
    axes[1, 1].set_ylabel('正規化強度', fontsize=12)
    axes[1, 1].set_title('Step 4: 正規化（最大値 = 1）', fontsize=14, fontweight='bold')
    axes[1, 1].legend()
    axes[1, 1].grid(alpha=0.3)

    plt.tight_layout()
    plt.show()

# 実行
demonstrate_preprocessing()

4.6.2 ピークデコンボリューション

重なり合ったピークを分離するため、複数のピーク関数（ガウス、ローレンツ、Voigt）の和でスペクトルをフィッティングします。

コード例7: D/Gバンドの詳細フィッティング

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit, minimize
from scipy.special import wofz

def voigt_peak(x, amplitude, center, sigma, gamma):
    """
    Voigt関数（ガウスとローレンツの畳み込み）

    Parameters:
    -----------
    x : array
        波数
    amplitude : float
        ピーク振幅
    center : float
        ピーク中心
    sigma : float
        ガウス幅
    gamma : float
        ローレンツ幅

    Returns:
    --------
    y : array
        Voigt関数の値
    """
    z = ((x - center) + 1j*gamma) / (sigma * np.sqrt(2))
    return amplitude * np.real(wofz(z)) / (sigma * np.sqrt(2*np.pi))

def five_peak_model(x, *params):
    """
    5ピークモデル（D1, D2, D3, D4, G）
    炭素材料の詳細解析用

    params: [A1, c1, w1, A2, c2, w2, ...]
    """
    n_peaks = len(params) // 3
    result = np.zeros_like(x)
    for i in range(n_peaks):
        A, c, w = params[3*i:3*i+3]
        result += gaussian_peak(x, A, c, w)
    return result

def fit_carbon_dg_bands(wavenumber, intensity, model='five_peak'):
    """
    炭素材料のD/Gバンド領域のフィッティング

    Parameters:
    -----------
    wavenumber : array
        ラマンシフト (1000-1800 cm^-1 の範囲が望ましい)
    intensity : array
        ラマン強度
    model : str
        'two_peak' (D, G only) または 'five_peak' (D1, D2, D3, D4, G)

    Returns:
    --------
    result : dict
        フィッティング結果
    """
    if model == 'two_peak':
        # 2ピークモデル（D, G）
        peak_names = ['D', 'G']
        p0 = [np.max(intensity)*0.3, 1350, 50, np.max(intensity), 1580, 20]
        bounds = ([0, 1300, 10, 0, 1550, 5],
                  [np.max(intensity)*2, 1400, 100, np.max(intensity)*2, 1620, 50])
    else:
        # 5ピークモデル（D1, D2, D3, D4, G）
        # D1: ~1350 (主欠陥ピーク)
        # D2/D': ~1620 (欠陥ピーク)
        # D3: ~1500 (アモルファス炭素)
        # D4: ~1200 (sp3炭素/多重フォノン)
        # G: ~1580 (sp2 in-plane)
        peak_names = ['D4', 'D1', 'D3', 'G', 'D2']
        I_max = np.max(intensity)
        p0 = [I_max*0.1, 1200, 50,   # D4
              I_max*0.4, 1350, 40,   # D1
              I_max*0.1, 1500, 80,   # D3
              I_max*0.8, 1580, 20,   # G
              I_max*0.2, 1620, 15]   # D2

        bounds = ([0, 1150, 20, 0, 1300, 20, 0, 1450, 30, 0, 1560, 10, 0, 1600, 5],
                  [I_max*2]*15)
        bounds[1][2] = 1250   # D4 center max
        bounds[1][5] = 1400   # D1 center max
        bounds[1][8] = 1560   # D3 center max
        bounds[1][11] = 1600  # G center max
        bounds[1][14] = 1650  # D2 center max

    n_peaks = len(peak_names)

    try:
        popt, pcov = curve_fit(
            lambda x, *p: five_peak_model(x, *p) if len(p) == 15
                          else multi_peak_model(x, *p),
            wavenumber, intensity, p0=p0, bounds=bounds, maxfev=20000
        )
        perr = np.sqrt(np.diag(pcov))
    except RuntimeError:
        print("Warning: フィッティングが収束しませんでした。")
        return None

    # 結果の整理
    result = {'peaks': {}, 'total_fit': None}
    total_area = 0

    for i, name in enumerate(peak_names):
        A, c, w = popt[3*i:3*i+3]
        area = A * np.sqrt(2*np.pi) * w
        total_area += area
        result['peaks'][name] = {
            'amplitude': A,
            'center': c,
            'width': w,
            'fwhm': 2.355 * w,
            'area': area
        }

    # 面積比の計算
    D_total = sum([result['peaks'][k]['area'] for k in result['peaks'] if k.startswith('D')])
    G_area = result['peaks']['G']['area']
    result['D_G_area_ratio'] = D_total / G_area
    result['D1_G_ratio'] = result['peaks'].get('D1', result['peaks'].get('D', {})).get('area', 0) / G_area

    # フィット曲線
    result['total_fit'] = five_peak_model(wavenumber, *popt) if len(popt) == 15 else multi_peak_model(wavenumber, *popt)
    result['popt'] = popt
    result['peak_names'] = peak_names

    return result

def demonstrate_dg_fitting():
    """
    D/Gバンドフィッティングのデモンストレーション
    """
    # シミュレーションデータ（炭素材料）
    wavenumber = np.linspace(1000, 1800, 800)

    # 5成分スペクトル
    D4 = gaussian_peak(wavenumber, 30, 1210, 50)
    D1 = gaussian_peak(wavenumber, 150, 1350, 45)
    D3 = gaussian_peak(wavenumber, 40, 1510, 70)
    G = gaussian_peak(wavenumber, 300, 1585, 18)
    D2 = gaussian_peak(wavenumber, 80, 1620, 15)

    spectrum = D4 + D1 + D3 + G + D2
    noise = np.random.normal(0, 8, len(wavenumber))
    spectrum_noisy = spectrum + noise
    spectrum_noisy = np.maximum(spectrum_noisy, 0)

    # フィッティング
    result = fit_carbon_dg_bands(wavenumber, spectrum_noisy, model='five_peak')

    if result is None:
        print("フィッティング失敗")
        return

    # プロット
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # フィッティング結果
    axes[0].plot(wavenumber, spectrum_noisy, 'ko', markersize=2, alpha=0.3, label='実測データ')
    axes[0].plot(wavenumber, result['total_fit'], 'r-', linewidth=2, label='フィット（合計）')

    # 各成分を表示
    colors = ['blue', 'green', 'cyan', 'magenta', 'orange']
    for i, name in enumerate(result['peak_names']):
        peak = result['peaks'][name]
        component = gaussian_peak(wavenumber, peak['amplitude'], peak['center'], peak['width'])
        axes[0].fill_between(wavenumber, component, alpha=0.3, color=colors[i], label=f'{name} ({peak["center"]:.0f} cm$^{{-1}}$)')

    axes[0].set_xlabel('ラマンシフト (cm$^{-1}$)', fontsize=12)
    axes[0].set_ylabel('ラマン強度 (a.u.)', fontsize=12)
    axes[0].set_title('D/Gバンド 5ピークフィッティング', fontsize=14, fontweight='bold')
    axes[0].legend(loc='upper left', fontsize=9)
    axes[0].grid(alpha=0.3)

    # 残差
    residual = spectrum_noisy - result['total_fit']
    axes[1].plot(wavenumber, residual, 'b-', linewidth=1)
    axes[1].axhline(0, color='red', linestyle='--')
    axes[1].set_xlabel('ラマンシフト (cm$^{-1}$)', fontsize=12)
    axes[1].set_ylabel('残差', fontsize=12)
    axes[1].set_title('フィッティング残差', fontsize=14, fontweight='bold')
    axes[1].grid(alpha=0.3)

    plt.tight_layout()
    plt.show()

    # 結果表示
    print("\nフィッティング結果:")
    print("-" * 70)
    print(f"{'ピーク':<8} {'中心 (cm^-1)':<15} {'FWHM (cm^-1)':<15} {'面積 (%)':<15}")
    print("-" * 70)

    total_area = sum([p['area'] for p in result['peaks'].values()])
    for name, peak in result['peaks'].items():
        area_pct = peak['area'] / total_area * 100
        print(f"{name:<8} {peak['center']:<15.1f} {peak['fwhm']:<15.1f} {area_pct:<15.1f}")

    print("-" * 70)
    print(f"D/G面積比 (D1/G): {result['D1_G_ratio']:.3f}")
    print(f"D/G面積比 (全D/G): {result['D_G_area_ratio']:.3f}")

# 実行
demonstrate_dg_fitting()

4.7 演習問題

lambda_excitation = 532  # nm
lambda_scattered = 545   # nm

nu_excitation = 1e7 / lambda_excitation  # cm^-1
nu_scattered = 1e7 / lambda_scattered    # cm^-1

raman_shift = nu_excitation - nu_scattered
print(f"ラマンシフト: {raman_shift:.0f} cm^-1")

import numpy as np

h = 6.626e-34  # J*s
c = 3e10       # cm/s
k_B = 1.381e-23  # J/K
T = 300        # K
nu = 1000      # cm^-1

ratio = np.exp(-h * c * nu / (k_B * T))
print(f"I_AS / I_S = {ratio:.4f} ({ratio*100:.2f}%)")

A = 234  # cm^-1 nm
B = 10   # cm^-1
omega_rbm = 180  # cm^-1

diameter = A / (omega_rbm - B)
print(f"CNT直径: {diameter:.2f} nm")

lambda_L = 514  # nm
ID_IG = 0.5

La = 2.4e-10 * (lambda_L**4) * (1/ID_IG)
print(f"結晶子サイズ La: {La:.1f} nm")

C_normal = 1e-3  # M (1 mM)
C_sers = 1e-9    # M (1 nM)

EF = C_normal / C_sers
print(f"SERS増強因子: {EF:.0e}")

import numpy as np

h = 6.626e-34
c = 3e10
k_B = 1.381e-23
nu = 520  # cm^-1
IS_IAS = 8.5

IAS_IS = 1 / IS_IAS
T = -h * c * nu / (k_B * np.log(IAS_IS))
print(f"推定温度: {T:.0f} K ({T-273:.0f} C)")

import numpy as np

# RBMから直径
A, B = 234, 10
rbm_peaks = [150, 180, 220, 260]
for rbm in rbm_peaks:
    d = A / (rbm - B)
    print(f"RBM {rbm} cm^-1 -> d = {d:.2f} nm")

# カイラリティから直径と金属/半導体判定
a_cc = 0.142  # nm
a = np.sqrt(3) * a_cc

chiralities = [(10, 10), (8, 7), (9, 4), (6, 5)]
for n, m in chiralities:
    d = a * np.sqrt(n**2 + n*m + m**2) / np.pi
    mod = (n - m) % 3
    cnt_type = "金属" if mod == 0 else "半導体"
    print(f"({n},{m}): d = {d:.2f} nm, {cnt_type}")

import numpy as np
import matplotlib.pyplot as plt

# シミュレーションデータ（10x10マッピング）
np.random.seed(42)
dg_ratio_map = np.random.exponential(0.3, (10, 10))

# ホットスポット（局所的に高D/G領域を作成）
dg_ratio_map[2:4, 7:9] = np.random.uniform(0.6, 1.0, (2, 2))
dg_ratio_map[6:8, 3:5] = np.random.uniform(0.8, 1.2, (2, 2))

# ホットスポット判定（D/G > 0.5）
threshold = 0.5
hotspot_mask = dg_ratio_map > threshold
hotspot_ratio = np.sum(hotspot_mask) / hotspot_mask.size * 100

# プロット
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# D/Gマップ
im1 = axes[0].imshow(dg_ratio_map, cmap='hot', origin='lower')
axes[0].set_title('D/G比マッピング', fontsize=14, fontweight='bold')
axes[0].set_xlabel('X (点)', fontsize=12)
axes[0].set_ylabel('Y (点)', fontsize=12)
plt.colorbar(im1, ax=axes[0], label='D/G ratio')

# ホットスポットマップ
im2 = axes[1].imshow(hotspot_mask, cmap='Reds', origin='lower')
axes[1].set_title(f'ホットスポット（D/G > {threshold}）', fontsize=14, fontweight='bold')
axes[1].set_xlabel('X (点)', fontsize=12)
axes[1].set_ylabel('Y (点)', fontsize=12)
plt.colorbar(im2, ax=axes[1], label='Hotspot')

plt.tight_layout()
plt.show()

print(f"ホットスポット面積割合: {hotspot_ratio:.1f}%")
print(f"ホットスポット数: {np.sum(hotspot_mask)} / {hotspot_mask.size} ピクセル")
print(f"D/G比: 最小 {dg_ratio_map.min():.2f}, 最大 {dg_ratio_map.max():.2f}, 平均 {dg_ratio_map.mean():.2f}")