import numpy as np
import math
import argparse
import matplotlib.pyplot as plt

# 物理定数
KB = 1.380649e-23
QE = 1.60217663e-19
EPS0 = 8.854187e-12
ME = 9.109383e-31

def calculate_params(args):
    T = args.temp
    # リチャードソン定数の計算
    reff = 1.20173e6 * args.men # A/m2/K2
    am_cm2 = reff * 1e-4        # A/cm2/K2
    
    # 半導体側の統計
    nc = 2.51e19 * (args.men**1.5) * ((T/300)**1.5) # cm-3
    efn = args.ecn - (KB * T / QE) * math.log(args.ndn / nc)
    
    # 障壁と内蔵電位
    phi_b = -(args.ecn - args.efm)
    v_bi = -(efn - args.efm)
    
    # 抵抗と飽和電流
    sigma_n = QE * args.ndn * args.mun
    rn = 1.0 / sigma_n * args.dn * 1e-7 # ohm cm2
    rs = rn
    
    js = am_cm2 * T**2 * math.exp(-QE * phi_b / (KB * T))
    
    return js, v_bi, phi_b, nc, rs

def calculate_schottky_current(model_type, v_target, temp, mus, ncs, nds, epss, js, phi_b, v_bi, n_diode, rs, area, v_initial, men, n_tunnel):
    """
    argsオブジェクトを直接参照せず、必要な値(men, n_tunnel)を引数で受け取るように修正
    """
    if v_target == 0: return 0.0, 0.0
    v_diode = v_initial
    dv, eps, max_iter = 1e-4, 1e-7, 100
    vt = (n_diode * KB * temp) / QE
    
    # 単位換算 (cm -> m)
    mus_m = mus * 1e-4
    ncs_m = ncs * 1e6
    nds_m = nds * 1e6

    # get_current をループの外で定義できるよう vd を引数に取る形式に
    def get_current(vd):
        if v_target >= 0.0 or model_type == 'Simple' or (v_bi - vd) < 0:
            return js * (math.exp(vd / vt) - 1.0) * area
        
        elif model_type == 'Diffusion':
            e_max = math.sqrt(2.0 * QE * nds_m * (v_bi - vd) / epss)
            k_exp = QE * mus_m * ncs_m * e_max * math.exp(-QE * phi_b / (KB * temp))
            return k_exp * (math.exp(vd / vt) - 1.0) * 1e-4 * area
        
        elif model_type == 'Schottky':
            k1 = (QE / epss)**1.5 * math.sqrt(2) / (4.0 * math.pi) * math.sqrt(nds_m) * math.sqrt(v_bi - vd)
            k_exp = math.exp(QE / (KB * temp) * math.sqrt(k1))
            return js * k_exp * (math.exp(vd / vt) - 1.0) * area

        elif model_type == 'Tunneling':
            # 簡易的な指数関数モデル
            k_tunnel = math.exp(vd / (n_tunnel * 0.026)) 
            return js * k_tunnel * (math.exp(vd / vt) - 1.0) * area

        elif model_type == 'TFE':
            hbar = 1.0545718e-34
            m_eff = men * ME
            # E00の計算
            e00_j = (QE * hbar / 2.0) * math.sqrt(nds_m / (epss * m_eff))
            e00 = e00_j / QE # [eV]
            kt_q = (KB * temp) / QE
            e0 = e00 * (1.0 / math.tanh(e00 / kt_q))
            
            if vd >= 0:
                # 順方向：傾きが 1/e0 になる
                return js * (math.exp(vd / e0) - 1.0) * area
            else:
                # 逆方向：電界による障壁の薄層化を考慮（簡易WKB）
                # 電界が強くなる(vdがマイナスに大きくなる)ほど指数関数的に増大させる
                e_max = math.sqrt(2.0 * QE * nds_m * (v_bi - vd) / epss)
                # 補正項：E00が大きくなる（ドーピングが高い）ほど透過率が爆発的に増える
                tunnel_factor = math.exp(abs(vd) * (e00 / kt_q)**2) 
                return js * tunnel_factor * (math.exp(vd / e0) - 1.0) * area        
        return 0

    for _ in range(max_iter):
        im = get_current(v_diode - dv)
        ip = get_current(v_diode + dv)
        
        v_total = 0.5 * ((v_diode - dv + im * rs) + (v_diode + dv + ip * rs))
        if abs(v_total - v_target) < eps:
            return v_diode, 0.5 * (im + ip)
        
        dv_dv = ((v_diode + dv + ip * rs) - (v_diode - dv + im * rs)) / (2.0 * dv)
        v_diode += (v_target - v_total) / dv_dv
        
    return v_diode, get_current(v_diode)

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--model', choices=['Simple', 'Diffusion', 'Schottky', 'Tunneling', 'TFE'], default='Simple')
    parser.add_argument('--temp', type=float, default=300.0)
    parser.add_argument('--area', type=float, default=0.01)
    parser.add_argument('--efm', type=float, default=4.4, help='Metal Fermi level [eV]')
    parser.add_argument('--ecn', type=float, default=4.05, help='SC Electron Affinity [eV]')
    parser.add_argument('--ndn', type=float, default=1.0e16)
    parser.add_argument('--mun', type=float, default=1500.0)
    parser.add_argument('--eps_r', type=float, default=11.9, help='Relative permittivity')
    parser.add_argument('--dn', type=float, default=1000.0, help='SC thickness [nm]')
    parser.add_argument('--men', type=float, default=0.19)
    # トンネルモデル用のパラメータを追加
    parser.add_argument('--n_tunnel', type=float, default=2.0, help='Tunneling ideality factor')
    parser.add_argument('--v0', type=float, default=-2.0)
    parser.add_argument('--v1', type=float, default=1.0)
    parser.add_argument('--step', type=float, default=0.02)
    args = parser.parse_args()

    # パラメータ計算
    js, v_bi, phi_b, ncs, rs_unit = calculate_params(args)
    epss = args.eps_r * EPS0
    rs = rs_unit / args.area 

    v_apps = np.arange(args.v0, args.v1 + args.step, args.step)
    currents = []
    v_prev = args.v0 * 0.2
    
    for v_app in v_apps:
        v_d, i = calculate_schottky_current(
            args.model, v_app, args.temp, args.mun, ncs, args.ndn, epss, 
            js, phi_b, v_bi, 1.0, rs, args.area, v_prev,
            args.men, args.n_tunnel  # 不足していた引数を追加
        )
        currents.append(i)
        v_prev = v_d

    # 結果表示
    print(f"Model      : {args.model}")
    print(f"Barrier phiB: {phi_b:.4f} eV")
    print(f"V_bi       : {v_bi:.4f} V")
    print(f"Js         : {js:.6e} A/cm2")

    # グラフ表示
    plt.figure(figsize=(8, 6))
    plt.semilogy(v_apps, np.abs(currents), label=f'{args.model} model')
    plt.xlabel("Applied Voltage [V]")
    plt.ylabel("|Current| [A]")
    plt.grid(True, which='both', alpha=0.3)
    plt.legend()
    plt.title(f"Schottky Diode IV Characteristics ({args.model})")
    plt.show()

if __name__ == "__main__":
    main()