半導体信頼性物理

import numpy as np
import matplotlib.pyplot as plt

# プログラム名: reliability_analysis.py
# Program Name: reliability_analysis.py
# 内容: 信頼性関数、故障密度関数、故障率、MTTF、およびいくつかの分布モデルの計算と可視化
# Purpose: Calculation and visualization of reliability function, failure density, failure rate, MTTF, and distribution models

# --- 信頼度関数 / Reliability Function R(t) ---
def reliability_function(N0, C_t, t):
    return (N0 - C_t) / N0

# --- 故障密度関数 / Failure Density Function f(t) ---
def failure_density_function(R_t):
    return -np.gradient(R_t)

# --- 故障率 / Failure Rate λ(t) ---
def failure_rate(R_t, t):
    dR_t = np.gradient(R_t, t)
    return -dR_t / R_t

# --- 寿命 / MTTF (Mean Time To Failure) ---
def mttf(failure_density, time):
    return np.trapz(time * failure_density, time)

# --- 指数分布 / Exponential Distribution ---
def exponential_distribution(lmbda, t):
    R_t = np.exp(-lmbda * t)
    f_t = lmbda * np.exp(-lmbda * t)
    return R_t, f_t

# --- ワイブル分布 / Weibull Distribution ---
def weibull_distribution(m, t, t0=1, gamma=0):
    R_t = np.exp(-((t - gamma) / t0) ** m)
    f_t = m / t0 * ((t - gamma) / t0) ** (m - 1) * np.exp(-((t - gamma) / t0) ** m)
    return R_t, f_t

# --- グラフ描画 / Plotting Functions ---
def plot_graphs():
    t = np.linspace(0.1, 100, 500)
    N0 = 1000
    C_t = 100 * np.sin(t / 20) + 100
    R_t = reliability_function(N0, C_t, t)
    f_t = failure_density_function(R_t)
    λ_t = failure_rate(R_t, t)
    
    plt.figure(figsize=(12, 8))
    plt.subplot(3, 1, 1)
    plt.plot(t, R_t)
    plt.title("Reliability Function R(t)")
    plt.xlabel("Time t")
    plt.ylabel("R(t)")
    plt.grid(True)
    
    plt.subplot(3, 1, 2)
    plt.plot(t, f_t, color='red')
    plt.title("Failure Density Function f(t)")
    plt.xlabel("Time t")
    plt.ylabel("f(t)")
    plt.grid(True)
    
    plt.subplot(3, 1, 3)
    plt.plot(t, λ_t, color='green')
    plt.title("Failure Rate λ(t)")
    plt.xlabel("Time t")
    plt.ylabel("λ(t)")
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()

# 実行
plot_graphs()

# reliability_analysis.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import comb, factorial

# 信頼度関数 R(t)
def reliability_function(N0, C_t, t):
    """
    N0: 使用開始時の良品数
    C_t: 時間 t までに故障した製品の数
    t: 時間
    """
    return (N0 - C_t) / N0

# 故障密度関数 f(t)
def failure_density_function(R_t):
    """
    R_t: 信頼度関数
    """
    return -np.gradient(R_t)

# 故障率 λ(t)
def failure_rate(R_t, t):
    """
    R_t: 信頼度関数
    t: 時間
    """
    dR_t = np.gradient(R_t, t)
    return -dR_t / R_t

# 寿命 MTTF (Mean Time To Failure)
def mttf(failure_density, time):
    """
    failure_density: 故障密度関数
    time: 時間
    """
    return np.trapz(time * failure_density, time)

# --- 幾何分布 Geometric Distribution ---
def geometric_distribution(p, k):
    """
    p: 故障確率
    k: 時刻(離散)
    """
    q = 1 - p
    f_k = p * q**(k - 1)
    R_k = q**k
    mu = 1 / p
    return f_k, R_k, mu

# --- 負の二項分布 Negative Binomial Distribution ---
def negative_binomial_distribution(p, m, k):
    """
    p: 故障確率
    m: 成功（耐性）回数
    k: 試行回数 (k >= m)
    """
    q = 1 - p
    f_k = comb(k - 1, m - 1) * p**m * q**(k - m)
    mu = m * q / p
    return f_k, mu

# --- ポアソン分布 Poisson Distribution ---
def poisson_distribution(lmbda, r):
    """
    lmbda: 故障回数の平均
    r: 故障回数
    """
    P_r = (lmbda**r * np.exp(-lmbda)) / factorial(r)
    return P_r

# --- 二項分布 Binomial Distribution ---
def binomial_distribution(n, p, r):
    """
    n: 試行回数
    p: 故障確率
    r: 故障回数
    """
    return comb(n, r) * p**r * (1 - p)**(n - r)

# --- グラフを描画するためのサンプルデータ ---
def plot_graphs():
    k = np.arange(1, 20)
    r = np.arange(0, 20)

    # 使用開始時の良品数
    N0 = 1000
    
    # 時間 t までに故障した製品の数 (累積故障)
    C_t = 100 * np.sin(k / 100) + 100  # 仮の故障データ
    
    # 信頼度関数 R(t)
    R_t = reliability_function(N0, C_t, k)
    
    # 故障密度関数 f(t)
    f_t = failure_density_function(R_t)
    
    # 故障率 λ(t)
    λ_t = failure_rate(R_t, k)
    
    # 幾何分布
    p_geom = 0.2
    geom_f = [geometric_distribution(p_geom, ki)[0] for ki in k]

    # 負の二項分布
    p_nb = 0.3
    m = 3
    neg_bin_f = [negative_binomial_distribution(p_nb, m, ki)[0] for ki in k]

    # ポアソン分布
    lam = 5
    pois_f = [poisson_distribution(lam, ri) for ri in r]

    # 二項分布
    n = 20
    p_bin = 0.2
    binom_f = [binomial_distribution(n, p_bin, ri) for ri in r]

    # プロット
    fig, axs = plt.subplots(3, 2, figsize=(12, 8))

    # 信頼度関数 R(t)
    axs[0, 0].plot(k, R_t, label="Reliability Function R(t)")
    axs[0, 0].set_title("Reliability Function R(t)")
    axs[0, 0].set_xlabel("Time (t)")
    axs[0, 0].set_ylabel("R(t)")
    axs[0, 0].grid(True)

    # 故障密度関数 f(t)
    axs[0, 1].plot(k, f_t, label="Failure Density Function f(t)", color='red')
    axs[0, 1].set_title("Failure Density Function f(t)")
    axs[0, 1].set_xlabel("Time (t)")
    axs[0, 1].set_ylabel("f(t)")
    axs[0, 1].grid(True)

    # 故障率 λ(t)
    axs[1, 0].plot(k, λ_t, label="Failure Rate λ(t)", color='green')
    axs[1, 0].set_title("Failure Rate λ(t)")
    axs[1, 0].set_xlabel("Time (t)")
    axs[1, 0].set_ylabel("λ(t)")
    axs[1, 0].grid(True)

    # 幾何分布
    axs[1, 1].bar(k, geom_f, color='skyblue')
    axs[1, 1].set_title("Geometric Distribution")
    axs[1, 1].set_xlabel("k")
    axs[1, 1].set_ylabel("f(k)")

    # 負の二項分布
    axs[2, 0].bar(k, neg_bin_f, color='salmon')
    axs[2, 0].set_title("Negative Binomial Distribution")
    axs[2, 0].set_xlabel("k")
    axs[2, 0].set_ylabel("f(k)")

    # ポアソン分布
    axs[2, 1].bar(r, pois_f, color='limegreen')
    axs[2, 1].set_title("Poisson Distribution")
    axs[2, 1].set_xlabel("r")
    axs[2, 1].set_ylabel("P(r)")

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    plot_graphs()

import numpy as np
import matplotlib.pyplot as plt
#transition_delay_test.py
# ゼロタイム不良生成方法 (Zero-Time Defect Generation)
def zero_time_defect(T, Y, alpha, beta):
    """
    T: テストパターンの故障検出率 (Test Pattern Failure Detection Rate)
    Y: 歩留 (Yield)
    alpha: 修整パラメータ
    beta: 系数
    """
    # ゼロタイム不良生成関数
    DL = 1 - Y * (1 - T)
    FDL = alpha * (1 - Y * (1 - T)) + beta
    return DL, FDL

# データのプロット
def plot_defect_generation():
    # パラメータ設定
    T = np.linspace(0, 1, 500)  # テストパターンの故障検出率 (0から1の範囲)
    Y = 0.8  # 歩留 (仮設定)
    alpha = 0.5  # 修整パラメータ (仮設定)
    beta = 0.1  # 系数 (仮設定)
    
    # ゼロタイム不良生成
    DL, FDL = zero_time_defect(T, Y, alpha, beta)
    
    # プロット
    plt.figure(figsize=(8, 6))
    
    # ゼロタイム不良生成
    plt.plot(T, DL, label="Zero-Time Defect Generation (DL)", color='blue')
    plt.plot(T, FDL, label="Failure Detection Rate (FDL)", color='red')
    
    plt.xlabel("Test Pattern Failure Detection Rate (T)")
    plt.ylabel("Defect Generation / Failure Detection Rate")
    plt.title("Zero-Time Defect Generation and Failure Detection Rate")
    plt.legend()
    plt.grid(True)
    plt.show()

# 実行
if __name__ == "__main__":
    plot_defect_generation()

import numpy as np
import matplotlib.pyplot as plt

# バスタブカーブのシミュレーション関数
def bathtub_curve(t):
    """
    バスタブカーブをシミュレートする関数。
    時間tにおける故障率を計算します。
    
    t: 時間（使用時間）
    
    初期故障期、偶発故障期、摩耗故障期に分けて、故障率を計算します。
    """
   # bathtub_curve_simulation.py


    # 初期故障期 (早期故障) の減少
    # 初期故障期は指数関数的に故障率が減少する
    early_failure_rate = np.exp(-0.05 * t)
    
    # 偶発故障期 (安定期) の一定故障率
    # 安定期では故障率は一定
    normal_failure_rate = 0.02
    
    # 摩耗故障期 (後期故障) の増加
    # 使用時間が100を過ぎると故障率が増加する
    wear_out_failure_rate = 0.02 * np.exp(0.01 * (t - 100))
    
    # バスタブカーブの故障率
    # t < 100では初期故障期、100を超えると摩耗故障期に入る
    failure_rate = np.where(t < 100, early_failure_rate, wear_out_failure_rate)
    
    # 偶発故障期は時間が100を超えても一定の故障率を維持
    failure_rate = np.where((t >= 100) & (failure_rate < normal_failure_rate), normal_failure_rate, failure_rate)
    
    return failure_rate

# グラフを描画する関数
def plot_bathtub_curve():
    """
    バスタブカーブの故障率を時間tに対してプロットする関数。
    """
    
    # 時間（使用時間）の範囲を設定
    t = np.linspace(0, 200, 1000)  # 時間 t: 0から200まで1000点
    
    # バスタブカーブの故障率を計算
    failure_rate = bathtub_curve(t)
    
    # バスタブカーブをプロット
    plt.figure(figsize=(10, 6))
    plt.plot(t, failure_rate, label="Failure Rate", color='blue')
    
    # グラフの設定
    plt.title("Bathtub Curve - Reliability over Time")
    plt.xlabel("Time (t) - Usage Time")
    plt.ylabel("Failure Rate")
    plt.grid(True)
    plt.legend()
    
    # グラフを表示
    plt.show()

# 実行部分
if __name__ == "__main__":
    plot_bathtub_curve()

# zero_delay_defect_model.py
import numpy as np
import matplotlib.pyplot as plt

# --- 不良混入率 DLを計算する関数 / Defect Level Calculation ---
def calculate_DL(Y, T):
    """
    Parameters:
    Y: 歩留 (Yield)
    T: テストカバレッジ (Test coverage)
    """
    return 1 - Y * (1 - T)

# --- ゼロタイム不良率 FDLを計算する関数 / Zero-time Defect Rate Calculation ---
def calculate_FDL(Y, T, alpha, beta):
    """
    Parameters:
    Y: 歩留 (Yield)
    T: テストカバレッジ (Test coverage)
    alpha: 製造起因係数
    beta: テスト準備・設計起因係数
    """
    return alpha * (1 - Y * (1 - T)) + beta

# --- グラフ描画 / Visualization ---
def plot_defect_models():
    Y = np.linspace(0.8, 1.0, 100)       # 歩留
    T = 0.95                             # テストカバレッジ（固定）
    alpha = 0.6                          # 係数 α
    beta = 0.002                         # 係数 β

    DL_vals = calculate_DL(Y, T)
    FDL_vals = calculate_FDL(Y, T, alpha, beta)

    plt.figure(figsize=(10, 6))
    plt.plot(Y, DL_vals, label="Defect Level DL", color='blue')
    plt.plot(Y, FDL_vals, label="Zero-time Defect Rate FDL", color='red')
    plt.xlabel("Yield (Y)")
    plt.ylabel("Defect Rate")
    plt.title("DL and FDL vs Yield")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

# --- 実行部 / Main Execution ---
if __name__ == "__main__":
    plot_defect_models()

import numpy as np
import matplotlib.pyplot as plt

# --- 安全な指数関数（オーバーフロー対策） ---
def safe_exp(x):
    return np.exp(np.clip(x, -700, 700))

# --- 各モデルの TTF 関数 ---
def TTF_E_model(A, Ea, k, T, beta):
    return A * safe_exp(Ea / k) * safe_exp(-beta * T)

def TTF_1overE_model(A, Ea, k, T, E):
    return A * safe_exp(Ea / k) * safe_exp(1 / E)

def TTF_VG_model(A, Ea, k, T, gamma, Vg):
    return A * safe_exp(Ea / k) * safe_exp(-gamma * Vg)

def TTF_PowerLaw_model(A, Ea, k, T, eta, Vg):
    return A * safe_exp(Ea / k) * Vg ** (-eta)

# --- Electromigration MTF 関数 ---
def MTF_model(A, J, n, Ea, k, T):
    return A * J ** (-n) * safe_exp(Ea / (k * T))

# --- ストレスマイグレーション速度 R ---
def R_SM_model(C, T0, TY, N, Ea, k, T):
    return C * (T0 - TY) ** N * safe_exp(-Ea / (k * T))

# --- プロット ---
def plot_models():
    A = 1e-3
    Ea = 0.7  # eV
    k = 8.617e-5  # eV/K

    T = np.linspace(250, 400, 100)  # 温度 [K]
    E = np.linspace(1, 5, 100)      # 電界 [V/nm]
    Vg = np.linspace(1, 5, 100)     # ゲート電圧 [V]

    beta = 0.01
    gamma = 0.5
    eta = 3

    # --- 各モデルの計算 ---
    ttf_e = TTF_E_model(A, Ea, k, T, beta)
    ttf_1e = TTF_1overE_model(A, Ea, k, T, E)
    ttf_vg = TTF_VG_model(A, Ea, k, T, gamma, Vg)
    ttf_power = TTF_PowerLaw_model(A, Ea, k, T, eta, Vg)

    # --- Electromigration ---
    J = 1e10  # 電流密度 [A/m^2]
    n = 2
    mtf = MTF_model(A, J, n, Ea, k, T)

    # --- ストレスマイグレーション ---
    C = 1e-3
    T0 = 473  # K
    TY = 398  # K
    N = 2
    r_sm = R_SM_model(C, T0, TY, N, Ea, k, T)

    # --- プロット ---
    plt.figure(figsize=(14, 10))

    plt.subplot(2, 2, 1)
    plt.plot(T, ttf_e, label="E-model")
    plt.plot(T, ttf_1e, label="1/E-model")
    plt.plot(T, ttf_vg, label="VG-model")
    plt.plot(T, ttf_power, label="Power Law-model")
    plt.title("TTF Models vs Temperature")
    plt.xlabel("Temperature (K)")
    plt.ylabel("TTF")
    plt.yscale('log')
    plt.grid(True)
    plt.legend()

    plt.subplot(2, 2, 2)
    plt.plot(T, mtf, color='green', label="Electromigration MTF")
    plt.title("Electromigration MTF vs Temperature")
    plt.xlabel("Temperature (K)")
    plt.ylabel("MTF")
    plt.yscale('log')
    plt.grid(True)
    plt.legend()

    plt.subplot(2, 2, 3)
    plt.plot(T, r_sm, color='red', label="Stress Migration Rate")
    plt.title("Stress Migration Rate vs Temperature")
    plt.xlabel("Temperature (K)")
    plt.ylabel("Rate R")
    plt.yscale('log')
    plt.grid(True)
    plt.legend()

    plt.tight_layout()
    plt.show()

# 実行
plot_models()

import numpy as np
import matplotlib.pyplot as plt

# プログラム名: epoxy_stress_effect.py
# エポキシモールドによるピエゾ抵抗効果と応力変化をモデル化

# --- ピエゾ抵抗係数の定義（例として{100}面） ---
pi_11 = -102.2
pi_12 = 53.4
pi_44 = 138.1

# --- 応力変化率δρを定義（例データ） ---
delta_rho = np.array([
    [4312, 11078, 16076],  # Mount
    [1168, 2455, 3722],    # Mold
    [2280, 3886, 5508]     # Cure 2h
])

labels = ['Mount', 'Mold', 'Cure 2h']

# --- 各工程に対する応力成分の計算 ---
def compute_stress_components(delta_rho, pi_11, pi_12, pi_44, crystal='100'):
    """
    応力計算式を使ってσx, σy, τxyを求める。
    """
    if crystal == '100':
        A = pi_11 + pi_12
        B = -pi_44
        C = -2 * (pi_11 - pi_12)
    elif crystal == '111':
        A = pi_11 + pi_12
        B = -pi_11 + pi_12
        C = -2 * (pi_11 - pi_12)
    else:
        raise ValueError("Crystal orientation not supported")

    sigma_x = 0.5 * ((1 / A + 1 / B) * delta_rho[:, 0] + (1 / A + 1 / B) * delta_rho[:, 2])
    sigma_y = 0.5 * ((1 / A + 1 / B) * delta_rho[:, 0] + (1 / A + 1 / B) * delta_rho[:, 2])
    tau_xy = 1 / C * (delta_rho[:, 1] - delta_rho[:, 0])

    return sigma_x, sigma_y, tau_xy

# --- 応力の計算 ---
sigma_x, sigma_y, tau_xy = compute_stress_components(delta_rho, pi_11, pi_12, pi_44, crystal='100')

# --- 結果のプロット ---
x = np.arange(len(labels))

plt.figure(figsize=(10, 6))
plt.plot(x, sigma_x, marker='o', label='Sigma_x')
plt.plot(x, sigma_y, marker='s', label='Sigma_y')
plt.plot(x, tau_xy, marker='^', label='Tau_xy')
plt.xticks(x, labels)
plt.title('Stress Components by Processing Step ({100})')
plt.ylabel('Stress [N/cm²]')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# -*- coding: utf-8 -*-
# Program Name: OC_curve_sampling_plan_simulator.py
# 内容：抜取検査のOC曲線を描画する / Plotting the OC (Operating Characteristic) curve for acceptance sampling

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom

# --- パラメータ設定 / Set parameters ---
n = 20       # サンプルサイズ / Sample size
c = 1        # 合格判定の最大不良数 / Acceptance number
p_def = np.linspace(0, 0.2, 100)  # 不良率の範囲 / Defective rate (from 0% to 20%)

# --- 合格確率を計算 / Calculate acceptance probability ---
P_accept = binom.cdf(c, n, p_def)  # c以下の不良数が出る確率 / Probability of x ≤ c

# --- グラフ描画 / Plot OC Curve ---
plt.figure(figsize=(8, 6))
plt.plot(p_def * 100, P_accept, label='n=20, c=1', color='blue')
plt.title("OC Curve (Operating Characteristic Curve)")
plt.xlabel("Defective Rate [%]")
plt.ylabel("Probability of Acceptance")
plt.grid(True)
plt.axvline(x=1.0, color='gray', linestyle='--', label="AQL (例:1%)")
plt.axvline(x=7.0, color='red', linestyle='--', label="LTPD (例:7%)")
plt.legend()
plt.tight_layout()
plt.show()

# プログラム名: full_reliability_inspection_module.py
# 概要: 抜き取り検査、統計表、信頼性理論、確率モデルの全項目を含むPython統合コード

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

# --- A. 抜き取り検査 / Sampling Inspection ---
def aql_sampling_plan(N, AQL, inspection_level):
    n = min(max(int(np.sqrt(N) * inspection_level), 2), N)
    accept = int(n * AQL / 100)
    print(f"[AQL Sampling] Sample size: {n}, Acceptance number: {accept}")

def ltpd_plan(N, LTPD, beta):
    """簡易LTPDサンプリング計画 / Simplified LTPD plan"""
    n = int(np.log(beta) / np.log(1 - LTPD))
    print(f"[LTPD Sampling] Sample size for LTPD={LTPD*100:.1f}%, β={beta*100:.0f}%: n = {n}")

# --- B. 付表 (統計分布グラフ) / Statistical Tables ---
def plot_normal_distribution():
    x = np.linspace(-4, 4, 500)
    plt.figure()
    plt.plot(x, stats.norm.pdf(x), label='PDF')
    plt.plot(x, stats.norm.cdf(x), label='CDF')
    plt.title("Normal Distribution (PDF & CDF)")
    plt.legend()
    plt.grid(True)
    plt.show()

def plot_poisson(mu=5):
    x = np.arange(0, 20)
    pmf = stats.poisson.pmf(x, mu)
    plt.figure()
    plt.stem(x, pmf, basefmt=" ")
    plt.title(f"Poisson Distribution (mu={mu})")
    plt.grid(True)
    plt.show()

# --- C. 信頼性理論 / Reliability Theory ---
def reliability_exponential(t, lamb):
    return np.exp(-lamb * t)

def plot_reliability_curve():
    t = np.linspace(0, 20, 200)
    plt.figure()
    for lamb in [0.1, 0.2, 0.5]:
        R = reliability_exponential(t, lamb)
        plt.plot(t, R, label=f"λ={lamb}")
    plt.title("Exponential Reliability Function R(t)")
    plt.xlabel("Time t")
    plt.ylabel("Reliability R(t)")
    plt.legend()
    plt.grid(True)
    plt.show()

def series_parallel_system():
    R1 = 0.9
    R2 = 0.95
    R_series = R1 * R2
    R_parallel = 1 - (1 - R1) * (1 - R2)
    print(f"[System Reliability] Series: {R_series:.4f}, Parallel: {R_parallel:.4f}")

# --- D. 確率モデル / Probabilistic Models ---
def plot_discrete_distributions():
    x = np.arange(0, 10)

    # Bernoulli (1 trial)
    p = 0.3
    bern = stats.bernoulli.pmf([0,1], p)
    plt.figure()
    plt.stem([0,1], bern)
    plt.title("Bernoulli Distribution (p=0.3)")
    plt.grid(True)
    plt.show()

    # Binomial
    binom = stats.binom.pmf(x, n=10, p=0.5)
    plt.figure()
    plt.stem(x, binom)
    plt.title("Binomial Distribution (n=10, p=0.5)")
    plt.grid(True)
    plt.show()

    # Geometric
    geom = stats.geom.pmf(x+1, p=0.3)
    plt.figure()
    plt.stem(x+1, geom)
    plt.title("Geometric Distribution (p=0.3)")
    plt.grid(True)
    plt.show()

    # Hypergeometric
    hyper = stats.hypergeom.pmf(x, M=20, n=7, N=12)
    plt.figure()
    plt.stem(x, hyper)
    plt.title("Hypergeometric Distribution (M=20, n=7, N=12)")
    plt.grid(True)
    plt.show()

def plot_continuous_distributions():
    t = np.linspace(0.01, 10, 500)

    # Exponential
    plt.figure()
    plt.plot(t, stats.expon.pdf(t, scale=1), label="Exponential")
    plt.plot(t, stats.gamma.pdf(t, a=2), label="Gamma a=2")
    plt.plot(t, stats.norm.pdf(t, loc=5, scale=1), label="Normal μ=5")
    plt.plot(t, stats.weibull_min.pdf(t, 1.5), label="Weibull c=1.5")
    plt.title("Continuous Distributions")
    plt.legend()
    plt.grid(True)
    plt.show()

# --- Main 実行部 / Main Routine ---
def main():
    print("\n--- A. Sampling Inspection ---")
    aql_sampling_plan(N=1000, AQL=1.5, inspection_level=1.0)
    ltpd_plan(N=1000, LTPD=0.05, beta=0.1)

    print("\n--- B. Statistical Tables ---")
    plot_normal_distribution()
    plot_poisson()

    print("\n--- C. Reliability Theory ---")
    plot_reliability_curve()
    series_parallel_system()

    print("\n--- D. Probabilistic Models ---")
    plot_discrete_distributions()
    plot_continuous_distributions()

main()

# プログラム名: mosfet_degradation_models.py
# 概要: NBTI, TDDB, ホットキャリアの各種劣化モデルを温度・電界依存で可視化する

import numpy as np
import matplotlib.pyplot as plt

# --- 定数の定義 / Constants ---
k = 8.617e-5  # ボルツマン定数 [eV/K]
Ea = 0.8      # 活性化エネルギー [eV]（例値）
beta = 1.0    # 電界係数 [cm/MV]（仮定値）
A = 1e-5      # 定数（全モデル共通で仮定）

# --- 範囲の設定 / Parameter ranges ---
T = np.linspace(273, 423, 200)      # 温度 [K]
E = np.linspace(1, 10, 200)         # 電界 [MV/cm]
Vds = np.linspace(0.5, 2.0, 200)    # ドレイン電圧 [V]
Isub = np.logspace(-10, -5, 200)    # 基板電流 [A]
m = 3                               # ホットキャリアモデルの指数
B = 150                             # ホットキャリア指数項定数

# --- NBTIおよびTDDB モデル（形式同一） ---
def calc_mttf_exp_field(T, E):
    return A * np.exp(-beta * E) * np.exp(Ea / (k * T))

# --- ホットキャリアモデル式（基板電流・ドレイン電圧） ---
def calc_hot_carrier_current(Isub):
    C = 1e-6
    return C * Isub ** (-m)

def calc_hot_carrier_vds(Vds):
    A_v = 1e-6
    return A_v * np.exp(-B / Vds)

# --- プロット 1: NBTI/TDDB vs 電界強度 E ---
def plot_mttf_vs_field():
    T_fixed = 373  # 100℃ に固定
    mttf = calc_mttf_exp_field(T_fixed, E)
    plt.figure()
    plt.semilogy(E, mttf)
    plt.title("MTTF vs Electric Field (NBTI / TDDB)")
    plt.xlabel("Electric Field [MV/cm]")
    plt.ylabel("MTTF [h]")
    plt.grid(True)
    plt.show()

# --- プロット 2: MTTF vs 温度（TDDB）---
def plot_mttf_vs_temperature():
    E_fixed = 6  # 電界強度を固定
    mttf = calc_mttf_exp_field(T, E_fixed)
    plt.figure()
    plt.semilogy(1000 / T, mttf)
    plt.title("MTTF vs Temperature (TDDB Model)")
    plt.xlabel("1000 / T [1/K]")
    plt.ylabel("MTTF [h]")
    plt.grid(True)
    plt.show()

# --- プロット 3: ホットキャリア vs 基板電流 ---
def plot_hot_carrier_current():
    t = calc_hot_carrier_current(Isub)
    plt.figure()
    plt.loglog(Isub, t)
    plt.title("Hot Carrier Degradation vs Substrate Current")
    plt.xlabel("Isub [A]")
    plt.ylabel("Lifetime t [arb. unit]")
    plt.grid(True)
    plt.show()

# --- プロット 4: ホットキャリア vs ドレイン電圧 ---
def plot_hot_carrier_vds():
    t = calc_hot_carrier_vds(Vds)
    plt.figure()
    plt.plot(Vds, t)
    plt.title("Hot Carrier Degradation vs Vds")
    plt.xlabel("Vds [V]")
    plt.ylabel("Lifetime t [arb. unit]")
    plt.grid(True)
    plt.show()

# --- 全実行 / Run All ---
def main():
    plot_mttf_vs_field()
    plot_mttf_vs_temperature()
    plot_hot_carrier_current()
    plot_hot_carrier_vds()

main()

# プログラム名: reliability_metrics_distribution_models.py
# 概要: 信頼度関数、不信頼度関数、故障率関数、寿命分布（指数分布、正規分布、対数正規分布、ワイブル分布）の可視化

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, lognorm, weibull_min

# --- 時間軸の定義 / Define time axis ---
t = np.linspace(0.01, 10, 500)

# --- λ: 故障率（定数）/ Failure rate (constant) ---
λ = 0.5

# --- 指数分布 / Exponential distribution ---
R_exp = np.exp(-λ * t)                 # 信頼度関数 / Reliability
F_exp = 1 - R_exp                      # 不信頼度関数 / Unreliability
f_exp = λ * np.exp(-λ * t)            # 確率密度関数 / PDF
hazard_exp = f_exp / R_exp            # 故障率関数 / Hazard

# --- 正規分布 / Normal distribution ---
μ, σ = 5, 1
f_norm = norm.pdf(t, loc=μ, scale=σ)
R_norm = 1 - norm.cdf(t, loc=μ, scale=σ)

# --- 対数正規分布 / Log-normal distribution ---
f_lognorm = lognorm.pdf(t, s=0.5, scale=np.exp(μ))
R_lognorm = 1 - lognorm.cdf(t, s=0.5, scale=np.exp(μ))

# --- ワイブル分布 / Weibull distribution ---
m, η = 2.0, 5.0
f_weibull = weibull_min.pdf(t, c=m, scale=η)
R_weibull = weibull_min.sf(t, c=m, scale=η)
hazard_weibull = f_weibull / R_weibull

# --- プロット 1: 指数分布の基本関数 ---
plt.figure()
plt.plot(t, R_exp, label="R(t) Reliability")
plt.plot(t, F_exp, label="F(t) Unreliability")
plt.plot(t, f_exp, label="f(t) PDF")
plt.plot(t, hazard_exp, label="λ(t) Hazard")
plt.title("Exponential Distribution Functions")
plt.xlabel("Time t")
plt.ylabel("Value")
plt.legend()
plt.grid(True)
plt.show()

# --- プロット 2: 正規分布と対数正規分布 ---
plt.figure()
plt.plot(t, f_norm, label="Normal PDF")
plt.plot(t, R_norm, label="Normal R(t)")
plt.plot(t, f_lognorm, label="Log-normal PDF")
plt.plot(t, R_lognorm, label="Log-normal R(t)")
plt.title("Normal and Log-normal Distributions")
plt.xlabel("Time t")
plt.ylabel("Value")
plt.legend()
plt.grid(True)
plt.show()

# --- プロット 3: ワイブル分布 ---
plt.figure()
plt.plot(t, f_weibull, label="Weibull PDF")
plt.plot(t, R_weibull, label="Weibull R(t)")
plt.plot(t, hazard_weibull, label="Weibull λ(t)")
plt.title("Weibull Distribution Functions")
plt.xlabel("Time t")
plt.ylabel("Value")
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import expon, weibull_min, lognorm, norm, gamma, geom, nbinom, binom, poisson
from scipy.integrate import quad

# --- Time axis / 時間軸の定義 ---
t = np.linspace(0.01, 10, 500)

# --- Exponential Distribution / 指数分布 ---
lambda_exp = 1.0
R_exp = np.exp(-lambda_exp * t)
f_exp = lambda_exp * np.exp(-lambda_exp * t)

# --- Weibull Distribution / ワイブル分布 ---
m_weibull = 2.0
t0_weibull = 2.0
R_weibull = np.exp(-(t / t0_weibull) ** m_weibull)
f_weibull = (m_weibull / t0_weibull) * (t / t0_weibull) ** (m_weibull - 1) * np.exp(-(t / t0_weibull) ** m_weibull)

# --- Lognormal Distribution / 対数正規分布 ---
mu_ln, sigma_ln = 1.0, 0.4
f_lognorm = lognorm.pdf(t, s=sigma_ln, scale=np.exp(mu_ln))

# --- Normal Distribution / 正規分布 ---
mu_norm, sigma_norm = 5.0, 1.0
f_norm = norm.pdf(t, loc=mu_norm, scale=sigma_norm)

# --- Gamma Distribution / ガンマ分布 ---
k_gamma, theta_gamma = 2.0, 2.0
f_gamma = gamma.pdf(t, a=k_gamma, scale=theta_gamma)

# --- Plotting continuous distributions ---
plt.figure(figsize=(12, 10))

plt.subplot(3, 2, 1)
plt.plot(t, R_exp, label='R(t) - Exponential')
plt.plot(t, f_exp, label='f(t) - Exponential')
plt.title('Exponential Distribution')
plt.xlabel('Time t')
plt.ylabel('Probability')
plt.legend()

plt.subplot(3, 2, 2)
plt.plot(t, R_weibull, label='R(t) - Weibull')
plt.plot(t, f_weibull, label='f(t) - Weibull')
plt.title('Weibull Distribution')
plt.xlabel('Time t')
plt.ylabel('Probability')
plt.legend()

plt.subplot(3, 2, 3)
plt.plot(t, f_lognorm, label='f(t) - Lognormal')
plt.title('Lognormal Distribution')
plt.xlabel('Time t')
plt.ylabel('Probability')
plt.legend()

plt.subplot(3, 2, 4)
plt.plot(t, f_norm, label='f(t) - Normal')
plt.title('Normal Distribution')
plt.xlabel('Time t')
plt.ylabel('Probability')
plt.legend()

plt.subplot(3, 2, 5)
plt.plot(t, f_gamma, label='f(t) - Gamma')
plt.title('Gamma Distribution')
plt.xlabel('Time t')
plt.ylabel('Probability')
plt.legend()

plt.tight_layout()
plt.show()

# --- Discrete distributions ---
k = np.arange(1, 21)

p_geom = 0.3
pmf_geom = geom.pmf(k, p_geom)

r_nbinom = 5
p_nbinom = 0.4
pmf_nbinom = nbinom.pmf(k - 1, r_nbinom, p_nbinom)

n_binom = 20
p_binom = 0.3
pmf_binom = binom.pmf(k, n_binom, p_binom)

lambda_poisson = 5
pmf_poisson = poisson.pmf(k, mu=lambda_poisson)

plt.figure(figsize=(12, 8))

plt.subplot(2, 2, 1)
plt.stem(k, pmf_geom, basefmt=" ")
plt.title('Geometric Distribution')
plt.xlabel('k')
plt.ylabel('PMF')

plt.subplot(2, 2, 2)
plt.stem(k, pmf_nbinom, basefmt=" ")
plt.title('Negative Binomial Distribution')
plt.xlabel('k')
plt.ylabel('PMF')

plt.subplot(2, 2, 3)
plt.stem(k, pmf_binom, basefmt=" ")
plt.title('Binomial Distribution')
plt.xlabel('k')
plt.ylabel('PMF')

plt.subplot(2, 2, 4)
plt.stem(k, pmf_poisson, basefmt=" ")
plt.title('Poisson Distribution')
plt.xlabel('k')
plt.ylabel('PMF')

plt.tight_layout()
plt.show()

# --- Hazard functions and MTTF ---
mttf_exp, _ = quad(lambda t: t * lambda_exp * np.exp(-lambda_exp * t), 0, np.inf)
haz_exp = lambda_exp * np.ones_like(t)

mttf_weibull, _ = quad(lambda t: t * (m_weibull / t0_weibull) * (t / t0_weibull) ** (m_weibull - 1) * np.exp(-(t / t0_weibull) ** m_weibull), 0, np.inf)
haz_weibull = f_weibull / R_weibull

f_ln_full = lognorm.pdf(t, s=sigma_ln, scale=np.exp(mu_ln))
R_ln_full = 1 - lognorm.cdf(t, s=sigma_ln, scale=np.exp(mu_ln))
mttf_lognorm, _ = quad(lambda t: t * lognorm.pdf(t, s=sigma_ln, scale=np.exp(mu_ln)), 0, np.inf)
haz_lognorm = f_ln_full / R_ln_full

plt.figure(figsize=(10, 6))
plt.plot(t, haz_exp, label=f'Hazard - Exponential (MTTF={mttf_exp:.2f})')
plt.plot(t, haz_weibull, label=f'Hazard - Weibull (MTTF={mttf_weibull:.2f})')
plt.plot(t, haz_lognorm, label=f'Hazard - Lognormal (MTTF={mttf_lognorm:.2f})')
plt.title('Hazard Functions and MTTF')
plt.xlabel('Time t')
plt.ylabel('Hazard Rate λ(t)')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

f_norm_full = norm.pdf(t, loc=mu_norm, scale=sigma_norm)
R_norm_full = 1 - norm.cdf(t, loc=mu_norm, scale=sigma_norm)
mttf_normal, _ = quad(lambda t: t * norm.pdf(t, loc=mu_norm, scale=sigma_norm), 0, np.inf)
haz_normal = f_norm_full / R_norm_full

f_gamma_full = gamma.pdf(t, a=k_gamma, scale=theta_gamma)
R_gamma_full = 1 - gamma.cdf(t, a=k_gamma, scale=theta_gamma)
mttf_gamma, _ = quad(lambda t: t * gamma.pdf(t, a=k_gamma, scale=theta_gamma), 0, np.inf)
haz_gamma = f_gamma_full / R_gamma_full

plt.figure(figsize=(10, 6))
plt.plot(t, haz_normal, label=f'Hazard - Normal (MTTF={mttf_normal:.2f})')
plt.plot(t, haz_gamma, label=f'Hazard - Gamma (MTTF={mttf_gamma:.2f})')
plt.title('Hazard Functions and MTTF (Normal & Gamma)')
plt.xlabel('Time t')
plt.ylabel('Hazard Rate λ(t)')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# プログラム名: accelerated_life_failure_rate.py
# 内容: 加速寿命試験データから故障率を予測（FIT単位）

import numpy as np

# --- 入力値 / Inputs ---
n_devices = 100                 # デバイス数
t_test_hours = 2000            # 試験時間（時間単位）
T_use = 50 + 273.15            # 実使用温度 [K]
T_test = 125 + 273.15          # 試験温度 [K]
Ea = 0.8                       # 活性化エネルギー [eV]
k = 8.617e-5                   # ボルツマン定数 [eV/K]
confidence = 0.60              # 信頼度水準（0.60 または 0.90）
r = 0                          # 故障数

# --- 信頼性水準に対応する仮想故障数 a
alpha_table = {
    60: [0.917, 2.02, 3.00, 3.84, 4.67, 5.48, 6.26, 7.03, 7.79, 8.55],
    90: [2.30, 3.89, 5.32, 6.64, 7.91, 9.15, 10.37, 11.57, 12.76, 13.95]
}
a_lookup = alpha_table[int(confidence * 100)][r] if r < 10 else None

# --- 加速係数（アレニウス式） ---
AT = np.exp((Ea / k) * (1 / T_use - 1 / T_test))  # 温度加速係数のみ（電圧1.0）

# --- 総試験時間（Component Hours） ---
T_total = n_devices * t_test_hours * AT

# --- 故障率λ（FIT = λ × 10^9） ---
lambda_fit = a_lookup / T_total
FIT = lambda_fit * 1e9

# --- 結果出力 ---
print(f"加速係数 AT = {AT:.1f}")
print(f"総試験時間 T = {T_total:.2e} hours")
print(f"信頼性水準{int(confidence*100)}%での仮想故障数 a = {a_lookup}")
print(f"予測故障率 λ = {lambda_fit:.2e} (1/hour)")
print(f"⇒ 予測故障率 = {FIT:.0f} FIT (Failure In Time)")

# プログラム名: lifetime_acceleration_models.py
# Program Name: lifetime_acceleration_models.py
# 内容: 各種加速モデル（アレニウス、Peck、HCI、電界加速）の寿命予測を統合し、可視化
# Description: Unified simulation of lifetime prediction using Arrhenius, Peck, HCI, and Field-Accelerated models

import numpy as np
import matplotlib.pyplot as plt

# --- 定数定義 / Physical constant ---
k = 8.617e-5  # ボルツマン定数 [eV/K] / Boltzmann constant

# --- アレニウスモデル / Arrhenius Model ---
def arrhenius_lifetime(A, Ea, T):
    """アレニウスモデルによる寿命 / Lifetime by Arrhenius model"""
    return A * np.exp(-Ea / (k * T))

# --- Peckモデル（温湿度劣化） / Peck Model (Humidity-Temperature Acceleration) ---
def peck_lifetime(A, RH, n, Ea, T):
    """Peckモデルによる寿命 / Lifetime by Peck model"""
    return A * RH**n * np.exp(Ea / (k * T))

# --- HCIモデル（ホットキャリア） / Hot Carrier Injection Model ---
def hci_lifetime(Isub, W, n, Ea, T):
    """ホットキャリア劣化モデル / Lifetime by Hot Carrier Injection"""
    return (Isub / W) ** (-n) * np.exp(Ea / (k * T))

# --- 電界加速モデル / Field-Accelerated Lifetime Model ---
def field_accelerated_lifetime(C, a, E):
    """電界加速劣化モデル / Lifetime by electric field acceleration"""
    return C * np.exp(-a * E)

# --- パラメータ設定 / Parameter settings ---
T = np.linspace(300, 400, 100)  # 温度範囲 [K] / Temperature range
RH = 0.6  # 湿度 / Relative Humidity
Isub = 1e-7  # 基板電流 [A] / Substrate current
W = 1e-5     # トランジスタ幅 [cm] / Device width
E = np.linspace(0.5, 5, 100)  # 電界 [MV/cm] / Electric field

# --- 各モデルの寿命計算 / Lifetime calculations ---
arr_life = arrhenius_lifetime(A=1e9, Ea=0.8, T=T)
peck_life = peck_lifetime(A=1e6, RH=RH, n=2.0, Ea=0.7, T=T)
hci_life = hci_lifetime(Isub=Isub, W=W, n=1.5, Ea=0.9, T=T)
field_life = field_accelerated_lifetime(C=1e6, a=2.5, E=E)

# --- プロット / Plot results ---
plt.figure(figsize=(12, 6))

# 温度依存モデルのプロット / Temperature-dependent models
plt.subplot(1, 2, 1)
plt.plot(T, arr_life, label='Arrhenius')
plt.plot(T, peck_life, label='Peck')
plt.plot(T, hci_life, label='HCI')
plt.xlabel('Temperature [K]')
plt.ylabel('Lifetime [hours]')
plt.title('Temperature-Dependent Lifetime Models')
plt.grid(True)
plt.legend()

# 電界依存モデルのプロット / Field-dependent model
plt.subplot(1, 2, 2)
plt.plot(E, field_life, label='Field Acceleration', color='purple')
plt.xlabel('Electric Field [MV/cm]')
plt.ylabel('Lifetime [hours]')
plt.title('Field-Accelerated Lifetime Model')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

# プログラム名: reliability_theory_models_summary.py
# Program Name: reliability_theory_models_summary.py
# 内容: 信頼性理論に基づく主要分布（正規、対数正規、ワイブル）と加速モデル（アレニウス、Peck、Miner）の統合可視化
# Description: Integrated visualization of reliability distributions (Normal, Log-normal, Weibull) and acceleration models (Arrhenius, Peck, Miner)

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, lognorm, weibull_min

# --- 時間軸定義 / Time axis ---
t = np.linspace(0.1, 10, 500)

# --- 分布モデル / Reliability distributions ---
μ, σ = 5, 1
shape_ln = 0.5
scale_ln = np.exp(μ)

# Normal distribution
f_norm = norm.pdf(t, loc=μ, scale=σ)
R_norm = 1 - norm.cdf(t, loc=μ, scale=σ)

# Log-normal distribution
f_lognorm = lognorm.pdf(t, s=shape_ln, scale=scale_ln)
R_lognorm = 1 - lognorm.cdf(t, s=shape_ln, scale=scale_ln)

# Weibull distribution
m, η = 2.0, 5.0
f_weibull = weibull_min.pdf(t, c=m, scale=η)
R_weibull = weibull_min.sf(t, c=m, scale=η)

# --- アレニウスモデル / Arrhenius model ---
k = 8.617e-5
T = np.linspace(300, 400, 100)
Ea = 0.8
A = 1e9
arrhenius_life = A * np.exp(-Ea / (k * T))

# --- Peckモデル（湿度加速） / Peck model ---
RH = 0.6
n_peck = 2.0
Ea_peck = 0.7
A_peck = 1e6
peck_life = A_peck * RH**n_peck * np.exp(Ea_peck / (k * T))

# --- Miner則（定数応力蓄積モデル） / Miner rule (cumulative damage) ---
stress_cycles = np.array([0.2, 0.3, 0.5])
life_ratios = 1 / stress_cycles
damage = np.cumsum(stress_cycles / life_ratios)

# --- 描画 / Plotting ---
plt.figure(figsize=(14, 8))

# Reliability distributions
plt.subplot(2, 2, 1)
plt.plot(t, R_norm, label='Normal R(t)')
plt.plot(t, R_lognorm, label='Log-normal R(t)')
plt.plot(t, R_weibull, label='Weibull R(t)')
plt.title('Reliability Functions R(t)')
plt.xlabel('Time t')
plt.ylabel('Reliability R(t)')
plt.grid()
plt.legend()

# PDF comparison
plt.subplot(2, 2, 2)
plt.plot(t, f_norm, label='Normal PDF')
plt.plot(t, f_lognorm, label='Log-normal PDF')
plt.plot(t, f_weibull, label='Weibull PDF')
plt.title('Probability Density Functions f(t)')
plt.xlabel('Time t')
plt.ylabel('Density f(t)')
plt.grid()
plt.legend()

# Acceleration model: Arrhenius and Peck
plt.subplot(2, 2, 3)
plt.plot(T, arrhenius_life, label='Arrhenius')
plt.plot(T, peck_life, label='Peck (RH=60%)')
plt.title('Acceleration Models: Temperature Dependency')
plt.xlabel('Temperature [K]')
plt.ylabel('Lifetime [hours]')
plt.grid()
plt.legend()

# Miner Rule
plt.subplot(2, 2, 4)
plt.plot(np.arange(1, len(damage) + 1), damage, marker='o')
plt.title('Miner’s Rule: Cumulative Damage')
plt.xlabel('Cycle')
plt.ylabel('Accumulated Damage')
plt.grid()

plt.tight_layout()
plt.show()