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

# 物理定数
KB = 1.380649e-23
QE = 1.60217663e-19

def calculate_material_params(args):
    T = args.temp
    nc = 2.51e19 * (args.men**1.5) * ((T/300)**1.5)
    nv = 2.51e19 * (args.mhp**1.5) * ((T/300)**1.5)
    
    efn = args.ecn - (KB * T / QE) * math.log(args.ndn / nc)
    efp = args.evp + (KB * T / QE) * math.log(args.nap / nv)
    v_bi = efp - efn

    dn = args.mun * (KB * T / QE)
    dp = args.mup * (KB * T / QE)
    ln, lp = math.sqrt(dn * args.taun), math.sqrt(dp * args.taup)

    term_exp = math.exp(-QE * v_bi / (KB * T))
    js = QE * (dn / (ln)) * args.ndn * term_exp + QE * (dp / (lp)) * args.nap * term_exp
    return js, v_bi

def calculate_diode_current(v_target, temp, js, n_diode, rs, area, v_initial):
    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

    for _ in range(max_iter):
        get_i = lambda vd: js * (math.exp(vd / vt) - 1.0) * area
        im, ip = get_i(v_diode - dv), get_i(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, js * (math.exp(v_diode / vt) - 1.0) * area

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--temp', type=float, default=300.0)
    parser.add_argument('--area', type=float, default=0.01)
    parser.add_argument('--rs', type=float, default=10.0, help='Series Resistance [ohm]')
    parser.add_argument('--ndn', type=float, default=1.0e16)
    parser.add_argument('--ecn', type=float, default=4.05)
    parser.add_argument('--men', type=float, default=0.19)
    parser.add_argument('--mun', type=float, default=1500.0)
    parser.add_argument('--taun', type=float, default=1e-5)
    parser.add_argument('--nap', type=float, default=1.0e16)
    parser.add_argument('--evp', type=float, default=5.17)
    parser.add_argument('--mhp', type=float, default=0.16)
    parser.add_argument('--mup', type=float, default=500.0)
    parser.add_argument('--taup', type=float, default=1e-5)
    parser.add_argument('--v0', type=float, default=-1.0)
    parser.add_argument('--v1', type=float, default=1.0)
    parser.add_argument('--step', type=float, default=0.02)
    parser.add_argument('--noplot', action='store_true', help='Disable plot display')
    args = parser.parse_args()

    js, v_bi = calculate_material_params(args)
    
    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:
        _, i = calculate_diode_current(v_app, args.temp, js, 1.0, args.rs, args.area, v_prev)
        currents.append(i)
        v_prev = _

    abs_currents = np.abs(currents)

    # グラフ描画
    if not args.noplot:
        plt.figure(figsize=(8, 6))
        plt.semilogy(v_apps, abs_currents, 'b-', label=f'Rs = {args.rs} ohm')
        plt.title("PN Junction IV Characteristics")
        plt.xlabel("Applied Voltage (V)")
        plt.ylabel("|Current| (A)")
        plt.grid(True, which="both", ls="-", alpha=0.5)
        plt.legend()
        
        # 0Vの線を表示
        plt.axvline(0, color='black', lw=1)
        
        print(f"Calculated Vbi: {v_bi:.4f} V")
        print(f"Calculated Js : {js:.6e} A/cm2")
        plt.show()

if __name__ == "__main__":
    main()