sample script for fitting arbitrary complex refractive index over bro…

…ad bandwidth to sum of Lorenztians (#2005)

doc/docs/FAQ.md

@@ -306,10 +306,6 @@ In nonlinear media involving e.g. $\chi^{(2)}$, it is possible for an odd-symmet
 Usage: Materials
 ----------------
 
-### How do I import n and k values into Meep?
-
-You can import any arbitrary complex permittivity profile via n and k values into Meep by fitting the wavelength- or frequency-dependent data to a sum of Drude-Lorentz polarizability terms as described in [Materials](Materials.md#material-dispersion). In general, you have to use nonlinear optimization to do the fit (e.g., to minimize the sum-of-squares errors or whatever error criterion you prefer). Enough Lorentzians should form a complete basis, so you should be able to fit any function given enough Lorentzians. A wavelength-dependent, purely-real permittivity (i.e., with no loss) which can be represented using the [Sellmeier equation](https://en.wikipedia.org/wiki/Sellmeier_equation) can be directly [transferred to the Lorentz model using a simple substitution of variables](Materials.md#sellmeier-coefficients). Note that Meep only does [subpixel averaging of the nondispersive part of $\varepsilon$ (and $\mu$)](Subpixel_Smoothing.md#what-about-dispersive-materials).
-
 ### Is there a materials library?
 
 Yes. A [materials library](https://github.com/NanoComp/meep/blob/master/python/materials.py) is available containing [crystalline silicon](https://en.wikipedia.org/wiki/Crystalline_silicon) (c-Si), [amorphous silicon](https://en.wikipedia.org/wiki/Amorphous_silicon) (a-Si) including the hydrogenated form, [silicon dioxide](https://en.wikipedia.org/wiki/Silicon_dioxide) (SiO<sub>2</sub>), [indium tin oxide](https://en.wikipedia.org/wiki/Indium_tin_oxide) (ITO), [alumina](https://en.wikipedia.org/wiki/Aluminium_oxide) (Al<sub>2</sub>O<sub>3</sub>), [gallium arsenide](https://en.wikipedia.org/wiki/Gallium_arsenide) (GaAs), [gallium nitride](https://en.wikipedia.org/wiki/Gallium_nitride) (GaN), [aluminum arsenide](https://en.wikipedia.org/wiki/Aluminium_arsenide) (AlAs), [aluminum nitride](https://en.wikipedia.org/wiki/Aluminium_nitride) (AlN), [borosilicate glass](https://en.wikipedia.org/wiki/Borosilicate_glass) (BK7), [fused quartz](https://en.wikipedia.org/wiki/Fused_quartz), [silicon nitride](https://en.wikipedia.org/wiki/Silicon_nitride) (Si<sub>3</sub>N<sub>4</sub>), [germanium](https://en.wikipedia.org/wiki/Germanium) (Ge), [indium phosphide](https://en.wikipedia.org/wiki/Indium_phosphide) (InP), [lithium niobate](https://en.wikipedia.org/wiki/Lithium_niobate) (LiNbO<sub>3</sub>), as well as 11 elemental metals: [silver](https://en.wikipedia.org/wiki/Silver) (Ag), [gold](https://en.wikipedia.org/wiki/Gold) (Au), [copper](https://en.wikipedia.org/wiki/Copper) (Cu), [aluminum](https://en.wikipedia.org/wiki/Aluminium) (Al), [berylium](https://en.wikipedia.org/wiki/Beryllium) (Be), [chromium](https://en.wikipedia.org/wiki/Chromium) (Cr), [nickel](https://en.wikipedia.org/wiki/Nickel) (Ni), [palladium](https://en.wikipedia.org/wiki/Palladium) (Pd), [platinum](https://en.wikipedia.org/wiki/Platinum) (Pt), [titanium](https://en.wikipedia.org/wiki/Titanium) (Ti), and [tungsten](https://en.wikipedia.org/wiki/Tungsten) (W). Additional information is provided in [Materials](Materials.md#materials-library).

doc/docs/Materials.md

@@ -51,6 +51,14 @@ $$\frac{i \sigma_n(\mathbf{x}) \cdot \omega_n^2 }{\omega (\gamma_n - i\omega)}$$
 
 which is equivalent to the Lorentzian model except that the $\omega_n^2$ term has been omitted from the denominator, and asymptotes to a conductivity $\sigma_n \omega_n^2 / \gamma_n$ as $\omega\to 0$. In this case, $\omega_n^2$ is just a dimensional scale factor and has no interpretation as a resonance frequency.
 
+### How do I import n and k values into Meep?
+
+You can import into Meep any arbitrary complex permittivity profile via $n$ and $k$ values obtained via e.g. [ellipsometry](https://en.wikipedia.org/wiki/Ellipsometry) by fitting the wavelength- or frequency-dependent data to a sum of Lorentzian polarizability terms. In general, you have to use nonlinear optimization to perform the fit (e.g., to minimize the sum-of-squares errors or whatever error criterion you prefer). Enough Lorentzians should form a complete basis, so you should be able to fit any function given enough Lorentzians. Unfortunately, the fitting process requires some trial and error to specify the number of fitting parameters and their initial values. For a demonstration, see this [script](../../python/examples/eps_fit_lorentzian.py).
+
+A wavelength-dependent, purely-real permittivity (i.e., with no loss) which can be represented using the [Sellmeier equation](https://en.wikipedia.org/wiki/Sellmeier_equation) can be directly [transferred to the Lorentz model using a simple substitution of variables](#sellmeier-coefficients).
+
+Note: Meep only does [subpixel averaging of the nondispersive part of $\varepsilon$ (and $\mu$)](Subpixel_Smoothing.md#what-about-dispersive-materials).
+
 ### Sellmeier Coefficients
 
 For a wavelength-dependent, purely-real permittivity (i.e., with no loss) which can be represented via the [Sellmeier equation](https://en.wikipedia.org/wiki/Sellmeier_equation):

python/examples/eps_fit_lorentzian.py

@@ -0,0 +1,106 @@
+# fit complex refractive index profile over broad bandwidth imported from
+# a file to a sum of Lorentzian polarizability terms using gradient-based
+# optimization via NLopt (nlopt.readthedocs.io)
+
+import nlopt
+import numpy as np
+import meep as mp
+
+
+def lorentzfunc(p, x):
+    """Return the complex ε profile given a set of Lorentzian parameters p
+       (σ_0, ω_0, γ_0, σ_1, ω_1, γ_1, ...) for a set of frequencies x.
+    """
+    N = int(len(p)/3)
+    y = np.zeros(len(x))
+    for n in range(N):
+        A_n = p[3*n+0]
+        x_n = p[3*n+1]
+        g_n = p[3*n+2]
+        y = y + A_n/(np.square(x_n) - np.square(x) - 1j*x*g_n)
+    return y
+
+def lorentzerr(p, x, y, grad):
+    """Return the error (or residual or loss funnction) as the L2 norm
+       of the difference of ε(p,x) and y over a set of frequencies x as
+       well as the gradient of this error with respect to each Lorentzian
+       polarizability parameter in p and saving the result in grad.
+    """
+    N = int(len(p)/3)
+    yp = lorentzfunc(p, x)
+    val = np.sum(np.square(abs(y - yp)))
+    for n in range(N):
+        A_n = p[3*n+0]
+        x_n = p[3*n+1]
+        g_n = p[3*n+2]
+        d = 1 / (np.square(x_n) - np.square(x) - 1j*x*g_n)
+        if grad.size > 0:
+            grad[3*n+0] = 2 * np.real(np.dot(np.conj(yp - y), d))
+            grad[3*n+1] = -4 * x_n * A_n * np.real(np.dot(np.conj(yp - y), np.square(d)))
+            grad[3*n+2] = -2 * A_n * np.imag(np.dot(np.conj(yp - y), x * np.square(d)))
+    return val
+
+def lorentzfit(p0, x, y, alg=nlopt.LD_LBFGS, tol=1e-25, maxeval=10000):
+    """Return the optimal Lorentzian polarizability parameters and error
+       which minimize the error in ε(p0,x) relative to y for an initial
+       set of Lorentzian polarizability parameters p0 over a set of
+       frequencies x using the NLopt algorithm alg for a relative
+       tolerance tol and a maximum number of iterations maxeval.
+    """
+    opt = nlopt.opt(alg, len(p0))
+    opt.set_ftol_rel(tol)
+    opt.set_maxeval(maxeval)
+    opt.set_lower_bounds(np.zeros(len(p0)))
+    opt.set_upper_bounds(float('inf')*np.ones(len(p0)))
+    opt.set_min_objective(lambda p, grad: lorentzerr(p, x, y, grad))
+    local_opt = nlopt.opt(nlopt.LD_LBFGS, len(p0))
+    local_opt.set_ftol_rel(1e-10)
+    local_opt.set_xtol_rel(1e-8)
+    opt.set_local_optimizer(local_opt)
+    popt = opt.optimize(p0);
+    minf = opt.last_optimum_value()
+    return popt, minf
+
+
+# format of input data is three columns: wavelength (nm), real(n), complex(n)
+# each column is separated by a comma
+mydata = np.genfromtxt('mymaterial.csv',delimiter=',')
+n = mydata[:,1] + 1j*mydata[:,2]
+eps_inf = 1.2 # should be > 1.0 for stability; choose this value such that np.amin(np.real(eps)) is ~1.0
+eps = np.square(n) - eps_inf
+wl = mydata[:,0]
+wl_min = 399 # minimum wavelength (units of nm)
+wl_max = 701 # maximum wavelength (units of nm)
+start_idx = np.where(wl > wl_min)
+idx_start = start_idx[0][0]
+end_idx = np.where(wl < wl_max)
+idx_end = end_idx[0][-1]+1
+freqs = 1000/wl # units of 1/μm
+freqs_reduced = freqs[idx_start:idx_end]
+wl_reduced = wl[idx_start:idx_end]
+eps_reduced = eps[idx_start:idx_end]
+
+num_repeat = 30 # number of times to repeat local optimization with random initial values
+ps = np.zeros((num_repeat,3))
+mins = np.zeros(num_repeat)
+for m in range(num_repeat):
+    # specify one Lorentzian polarizability term consisting of three parameters (σ_0, ω_0, γ_0)
+    # with random initial values
+    # note: for the case of no absorption, γ (third parameter) is zero
+    p_rand = [10**(np.random.random()), 10**(np.random.random()), 10**(np.random.random())]
+    ps[m,:], mins[m] = lorentzfit(p_rand, freqs_reduced, eps_reduced, nlopt.LD_MMA, 1e-25, 50000)
+    print(f"iteration{m}:, {ps[m,:]}, {mins[m]}")
+
+# find the best performing set of parameters
+idx_min = np.where(np.min(mins) == mins)
+print(f"optimal:, {ps[idx_min]}, {mins[idx_min]}")
+
+# define a Meep `Medium` class object with the fitting parameters
+mymaterial_eps_inf = eps_inf
+mymaterial_freq = ps[idx_min][1]
+mymaterial_gamma = ps[idx_min][2]
+mymaterial_sigma = ps[idx_min][0] / mymaterial_freq**2
+mymaterial = mp.Medium(epsilon=mymaterial_eps_inf,
+                       E_susceptibilities=[mp.LorentzianSusceptibility(frequency=mymaterial_freq,
+                                                                       gamma=mymaterial_gamma,
+                                                                       sigma=mymaterial_sigma)])