/* Doping, carrier concentration, and mobility. Every printed residual/check is 0. */
kill(all)$
assume(N>0, P>0, ni>0, kT>0, m>0, tau>0, mul>0, mui>0)$

/* Exact compensated-semiconductor roots. */
nroot : (N+sqrt(N^2+4*ni^2))/2$
proot_from_n : ni^2/nroot$
proot : (P+sqrt(P^2+4*ni^2))/2$
nroot_from_p : ni^2/proot$
print("n-type neutrality residual =", radcan(nroot-proot_from_n-N))$
print("n-type mass-action residual =", ratsimp(nroot*proot_from_n-ni^2))$
print("p-type neutrality residual =", radcan(proot-nroot_from_p-P))$
print("p-type mass-action residual =", ratsimp(proot*nroot_from_p-ni^2))$

/* Fermi-level displacement signs. */
Ef_n : kT*log(n/ni)$
Ef_p : -kT*log(p/ni)$
print("n-type Fermi relation residual =", ratsimp(exp(Ef_n/kT)-n/ni))$
print("p-type Fermi relation residual =", ratsimp(exp(-Ef_p/kT)-p/ni))$

/* Intrinsic temperature power when Nc and Nv are proportional to T^(3/2). */
NcT : Cc*T^(3/2)$
NvT : Cv*T^(3/2)$
niT : sqrt(NcT*NvT)*exp(-Eg/(2*k*T))$
niT_expected : sqrt(Cc*Cv)*T^(3/2)*exp(-Eg/(2*k*T))$
assume(Cc>0,Cv>0,T>0,k>0)$
print("intrinsic-temperature factor residual =", radcan(niT-niT_expected))$

/* Relaxation-time mobility and Matthiessen's rule. */
vd : qc*tau*E/m$
print("steady-drift equation residual =", ratsimp(qc*E-m*vd/tau))$
mu_total : mul*mui/(mul+mui)$
print("Matthiessen residual =", ratsimp(1/mu_total-1/mul-1/mui))$

q0 : 1.602176634e-19$
kb_ev : 8.617333262145e-5$
T0 : 300$
kT0 : kb_ev*T0$

/* Solved problem 1. */
ND1 : 8.00e21$
NA1 : 2.00e21$
ni1 : 1.00e16$
N1 : ND1-NA1$
n1 : (N1+sqrt(N1^2+4*ni1^2))/2$
p1 : ni1^2/n1$
shift1 : kT0*log(n1/ni1)$
sigma1 : q0*(n1*0.135+p1*0.0480)$
rho1 : 1/sigma1$
print("solved-1 electron-density check =", round(1e11*n1/1e21)-600000000002)$
print("solved-1 hole-density check =", round(1000*p1/1e10)-1667)$
print("solved-1 Fermi-shift check =", round(1e4*shift1)-3440)$
print("solved-1 conductivity check =", round(10*sigma1)-1298)$
print("solved-1 resistivity check =", round(1e6*rho1)-7706)$

/* Solved problem 2. */
mu2 : 1/(1/0.180+1/0.450)$
sigma2 : q0*3.00e21*mu2$
rho2 : 1/sigma2$
print("solved-2 mobility check =", round(1e4*mu2)-1286)$
print("solved-2 conductivity check =", round(100*sigma2)-6180)$
print("solved-2 resistivity check =", round(1e5*rho2)-1618)$

/* Numerical problem 1: infer compensation from measured carriers. */
nn1 : 500000000000000000000$
pn1 : 200000000000$
Nsum_n1 : 1400000000000000000000$
Nnet_n1 : nn1-pn1$
ND_n1 : (Nsum_n1+Nnet_n1)/2$
NA_n1 : (Nsum_n1-Nnet_n1)/2$
comp_n1 : 100*NA_n1/ND_n1$
print("numerical-1 donor-density check =", round(1000*ND_n1/1e20)-9500)$
print("numerical-1 acceptor-density check =", round(1000*NA_n1/1e20)-4500)$
print("numerical-1 compensation-percent check =", round(100*comp_n1)-4737)$
print("numerical-1 neutrality reconstruction residual =", ratsimp(ND_n1-NA_n1-(nn1-pn1)))$
print("numerical-1 implant-total reconstruction residual =", ratsimp(ND_n1+NA_n1-Nsum_n1))$

/* Numerical problem 2: intrinsic-density temperature scaling. */
kb_exact_n2 : 8617333262145/10^17$
Tlow_n2 : 300$
Thigh_n2 : 400$
Eg_n2 : 28/25$
ni_low_n2 : 10^16$
dos_ratio_n2 : (Thigh_n2/Tlow_n2)^(3/2)$
thermal_exponent_n2 : Eg_n2*(1/Tlow_n2-1/Thigh_n2)/(2*kb_exact_n2)$
ni_ratio_n2 : dos_ratio_n2*exp(thermal_exponent_n2)$
ni_high_n2 : ni_low_n2*ni_ratio_n2$
print("numerical-2 density-ratio check =", round(10*ni_ratio_n2)-3462)$
print("numerical-2 high-temperature density check =", round(1000*ni_high_n2/10^18)-3462)$
print("numerical-2 DOS-scaling residual =", radcan(dos_ratio_n2^2-(Thigh_n2/Tlow_n2)^3))$
print("numerical-2 thermal-factor residual =", ratsimp(ni_ratio_n2/dos_ratio_n2-exp(thermal_exponent_n2)))$

/* Numerical problem 3. */
m0 : 9.1093837015e-31$
mu_n3 : q0*2.00e-13/(0.26*m0)$
print("numerical-3 mobility check =", round(1e4*mu_n3)-1353)$

/* Numerical problem 4. */
mu_n4 : 0.140*(400/300)^(-3/2)$
print("numerical-4 mobility check =", round(1e5*mu_n4)-9093)$

/* Numerical problem 5: extract impurity-limited mobility. */
mu_meas_n5 : 12/100$
mu_lattice_n5 : 20/100$
mu_impurity_n5 : 1/(1/mu_meas_n5-1/mu_lattice_n5)$
print("numerical-5 impurity-mobility check =", ratsimp(mu_impurity_n5-3/10))$
print("numerical-5 Matthiessen reconstruction residual =", ratsimp(1/mu_meas_n5-1/mu_lattice_n5-1/mu_impurity_n5))$

/* Numerical problem 6: lattice/impurity scattering crossover. */
Tref_n6 : 300$
muL_ref_n6 : 2/5$
muI_ref_n6 : 1/20$
x_n6 : (muL_ref_n6/muI_ref_n6)^(1/3)$
Tc_n6 : Tref_n6*x_n6$
muL_n6 : muL_ref_n6*x_n6^(-3/2)$
muI_n6 : muI_ref_n6*x_n6^(3/2)$
mu_total_n6 : 1/(1/muL_n6+1/muI_n6)$
print("numerical-6 crossover-temperature check =", radcan(Tc_n6-600))$
print("numerical-6 equal-limit residual =", radcan(muL_n6-muI_n6))$
print("numerical-6 total-mobility check =", round(10^5*mu_total_n6)-7071)$
print("numerical-6 exact-total-mobility residual =", radcan(mu_total_n6-sqrt(2)/20))$
print("numerical-6 Matthiessen residual =", ratsimp(1/mu_total_n6-1/muL_n6-1/muI_n6))$