/*
  Permanent dipoles and orientational polarization
  Maxima 5.49-compatible batch worksheet

  Run with:
      maxima --very-quiet -b permanent-dipoles-orientational-polarization.mac
*/

kill(all)$
display2d:false$
fpprintprec:12$

print("PERMANENT DIPOLES AND ORIENTATIONAL POLARIZATION")$

/* 1. Langevin function and its limiting forms */
L(x) := coth(x) - 1/x$

L_series : taylor(L(x), x, 0, 9)$
print("Langevin series through x^9:", L_series)$

L_prime : trigsimp(diff(L(x), x))$
print("dL/dx:", L_prime)$
print("lim_(x->0) L(x)/x =", limit(L(x)/x, x, 0))$
print("lim_(x->0) dL/dx =", limit(L_prime, x, 0))$
print("lim_(x->infinity) L(x) =", limit(L(x), x, inf))$
print("lim_(x->infinity) x*(1-L(x)) =",
      limit(x*(1-L(x)), x, inf))$

/* Check the cubic nonlinear-polarization coefficient. */
P_series : N*mu*subst(mu*Eloc/(kB*T), x, L_series)$
P_series : expand(P_series)$
print("Series for P = N mu L(mu Eloc/kBT):", P_series)$
print("Linear coefficient dP/dEloc at Eloc=0:",
      limit(diff(N*mu*L(mu*Eloc/(kB*T)), Eloc), Eloc, 0))$

/* 2. Lorentz local field and Debye--Clausius--Mossotti relation */
Aeff : alpha_ind + mu^2/(3*kB*T)$
local_field_equation :
    eps0*(er-1) = N*Aeff*(1 + (er-1)/3)$
er_solution : solve(local_field_equation, er)$
print("Solution for relative permittivity:", er_solution)$

CM_check : ratsimp((er-1)/(er+2) - N*Aeff/(3*eps0))$
CM_check_after_solution :
    ratsimp(subst(rhs(first(er_solution)), er, CM_check))$
print("Clausius--Mossotti residual after substitution (must be 0):",
      CM_check_after_solution)$

/* 3. Numerical example in SI units */
kb_num   : 1.380649e-23$
eps0_num : 8.8541878128e-12$
debye    : 3.33564e-30$
N_num    : 2.50e25$
mu_num   : 2.00*debye$
T_num    : 300.0$
E_num    : 1.00e6$

x_num : mu_num*E_num/(kb_num*T_num)$
P_exact : N_num*mu_num*L(x_num)$
P_linear : N_num*mu_num^2*E_num/(3*kb_num*T_num)$
chi_or : P_linear/(eps0_num*E_num)$
E_thermal : kb_num*T_num/mu_num$
P_sat : N_num*mu_num$
relative_error : (P_exact-P_linear)/P_linear$
series_relative_error : -x_num^2/15$

print("x = mu E/(kB T) =", float(x_num))$
print("Exact Langevin polarization (C/m^2) =", float(P_exact))$
print("Linear Debye polarization (C/m^2) =", float(P_linear))$
print("Orientational susceptibility =", float(chi_or))$
print("Thermal alignment field kB T/mu (V/m) =", float(E_thermal))$
print("Saturation polarization N mu (C/m^2) =", float(P_sat))$
print("Exact relative correction =", float(relative_error))$
print("Leading series correction -x^2/15 =", float(series_relative_error))$

/* 4. Debye rotational relaxation */
eps_real(y, eps_s, eps_inf) :=
    eps_inf + (eps_s-eps_inf)/(1+y^2)$
eps_loss(y, eps_s, eps_inf) :=
    (eps_s-eps_inf)*y/(1+y^2)$

loss_derivative : factor(diff(eps_loss(y, eps_s, eps_inf), y))$
print("d(epsilon loss)/d(omega tau) =", loss_derivative)$
print("Stationary points of the loss curve:",
      solve((y-1)*(y+1)=0, y))$
print("Loss at omega*tau=1:", eps_loss(1, eps_s, eps_inf))$
print("Low-frequency real permittivity:",
      limit(eps_real(y, eps_s, eps_inf), y, 0))$
print("High-frequency real permittivity:",
      limit(eps_real(y, eps_s, eps_inf), y, inf))$

/* Example relaxation spectrum: eps_s=20, eps_inf=4, tau=1 ps. */
eps_s_num : 20.0$
eps_inf_num : 4.0$
tau_num : 1.0e-12$
omega_peak : 1/tau_num$
f_peak : omega_peak/(2*%pi)$
print("For tau=1 ps, loss-peak angular frequency (rad/s) =",
      float(omega_peak))$
print("For tau=1 ps, loss-peak frequency (Hz) =", float(f_peak))$
print("Peak epsilon double-prime =",
      float(eps_loss(1, eps_s_num, eps_inf_num)))$

quit();