text:
# Meep Materials Library
# https://github.com/NanoComp/meep/blob/master/python/materials.py
# default unit length is 1 um
um_scale = 1.0

# conversion factor for eV to 1/um [=1/hc]
eV_um_scale = um_scale / 1.23984193

#------------------------------------------------------------------
# crystalline silicon (c-Si) from A. Deinega et al., J. Optical Society of America A, Vol. 28, No. 5, pp. 770-77 (2011)
# based on experimental data for intrinsic silicon at T=300K from M.A. Green and M. Keevers, Progress in Photovoltaics, Vol. 3, pp. 189-92 (1995)
# wavelength range: 0.4 - 1.0 um

cSi_range = FreqRange(min = um_scale, max = um_scale / 0.4)

cSi_frq1 = 3.64 / um_scale
cSi_gam1 = 0
cSi_sig1 = 8
cSi_frq2 = 2.76 / um_scale
cSi_gam2 = 2 * 0.063 / um_scale
cSi_sig2 = 2.85
cSi_frq3 = 1.73 / um_scale
cSi_gam3 = 2 * 2.5 / um_scale
cSi_sig3 = -0.107

cSi_susc = [LorentzianSusceptibility(frequency = cSi_frq1, gamma = cSi_gam1, sigma = cSi_sig1),
    LorentzianSusceptibility(frequency = cSi_frq2, gamma = cSi_gam2, sigma = cSi_sig2),
    LorentzianSusceptibility(frequency = cSi_frq3, gamma = cSi_gam3, sigma = cSi_sig3)]

cSi = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = cSi_susc, valid_freq_range = cSi_range))

#------------------------------------------------------------------
# amorphous silicon (a-Si) from Horiba Technical Note 08: Lorentz Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Lorentz_Dispersion_Model.pdf
# wavelength range: 0.21 - 0.83 um

aSi_range = FreqRange(min = um_scale / 0.83, max = um_scale / 0.21)

aSi_frq1 = 1 / (0.315481407124682 * um_scale)
aSi_gam1 = 1 / (0.645751005208333 * um_scale)
aSi_sig1 = 14.571

aSi_susc = [LorentzianSusceptibility(frequency = aSi_frq1, gamma = aSi_gam1, sigma = aSi_sig1)]

aSi = Medium(epsilon_model = DispersionModel(inf = 3.109, susceptibilities = aSi_susc, valid_freq_range = aSi_range))

#------------------------------------------------------------------
# hydrogenated amorphous silicon (a-Si:H) from Horiba Technical Note 08: Lorentz Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Lorentz_Dispersion_Model.pdf
# wavelength range: 0.21 - 0.83 um

aSi_H_range = FreqRange(min = um_scale / 0.83, max = um_scale / 0.21)

aSi_H_frq1 = 1 / (0.334189199460916 * um_scale)
aSi_H_gam1 = 1 / (0.579365387850467 * um_scale)
aSi_H_sig1 = 12.31

aSi_H_susc = [LorentzianSusceptibility(frequency = aSi_H_frq1, gamma = aSi_H_gam1, sigma = aSi_H_sig1)]

aSi_H = Medium(epsilon_model = DispersionModel(inf = 3.22, susceptibilities = aSi_H_susc, valid_freq_range = aSi_H_range))

#------------------------------------------------------------------
# indium tin oxide (ITO) from Horiba Technical Note 08: Lorentz Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Lorentz_Dispersion_Model.pdf
# wavelength range: 0.21 - 0.83 um

ITO_range = FreqRange(min = um_scale / 0.83, max = um_scale / 0.21)

ITO_frq1 = 1 / (0.182329695588235 * um_scale)
ITO_gam1 = 1 / (1.94637665620094 * um_scale)
ITO_sig1 = 2.5

ITO_susc = [LorentzianSusceptibility(frequency = ITO_frq1, gamma = ITO_gam1, sigma = ITO_sig1)]

ITO = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = ITO_susc, valid_freq_range = ITO_range))

#------------------------------------------------------------------
# alumina (Al2O3) from Horiba Technical Note 08: Lorentz Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Lorentz_Dispersion_Model.pdf
# wavelength range: 0.21 - 2.07 um

Al2O3_range = FreqRange(min = um_scale / 2.07, max = um_scale / 0.21)

Al2O3_frq1 = 1 / (0.101476668030774 * um_scale)
Al2O3_gam1 = 0
Al2O3_sig1 = 1.52

Al2O3_susc = [LorentzianSusceptibility(frequency = Al2O3_frq1, gamma = Al2O3_gam1, sigma = Al2O3_sig1)]

Al2O3 = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Al2O3_susc, valid_freq_range = Al2O3_range))

#------------------------------------------------------------------
# aluminum nitride (AlN) from Horiba Technical Note 08: Lorentz Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Lorentz_Dispersion_Model.pdf
# wavelength range: 0.26 - 1.65 um

AlN_range = FreqRange(min = um_scale / 1.65, max = um_scale / 0.26)

AlN_frq1 = 1 / (0.139058089950651 * um_scale)
AlN_gam1 = 0
AlN_sig1 = 3.306

AlN_susc = [LorentzianSusceptibility(frequency = AlN_frq1, gamma = AlN_gam1, sigma = AlN_sig1)]

AlN = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = AlN_susc, valid_freq_range = AlN_range))

#------------------------------------------------------------------
# aluminum arsenide (AlAs) from R.E. Fern and A. Onton, J. Applied Physics, Vol. 42, pp. 3499-500, 1971
# ref: https://refractiveindex.info/?shelf=main&book=AlAs&page=Fern
# wavelength range: 0.56 - 2.2 um

AlAs_range = FreqRange(min = um_scale / 2.2, max = um_scale / 0.56)

AlAs_frq1 = 1 / (0.2822 * um_scale)
AlAs_gam1 = 0
AlAs_sig1 = 6.0840
AlAs_frq2 = 1 / (27.62 * um_scale)
AlAs_gam2 = 0
AlAs_sig2 = 1.900

AlAs_susc = [LorentzianSusceptibility(frequency = AlAs_frq1, gamma = AlAs_gam1, sigma = AlAs_sig1),
    LorentzianSusceptibility(frequency = AlAs_frq2, gamma = AlAs_gam2, sigma = AlAs_sig2)]

AlAs = Medium(epsilon_model = DispersionModel(inf = 2.0792, susceptibilities = AlAs_susc, valid_freq_range = AlAs_range))

#------------------------------------------------------------------
# borosilicate glass (BK7) from SCHOTT Zemax catalog 2017-01-20b
# ref: https://refractiveindex.info/?shelf=glass&book=BK7&page=SCHOTT
# wavelength range: 0.3 - 2.5 um

BK7_range = FreqRange(min = um_scale / 2.5, max = um_scale / 0.3)

BK7_frq1 = 1 / (0.07746417668832478 * um_scale)
BK7_gam1 = 0
BK7_sig1 = 1.03961212
BK7_frq2 = 1 / (0.14148467902921502 * um_scale)
BK7_gam2 = 0
BK7_sig2 = 0.231792344
BK7_frq3 = 1 / (10.176475470417055 * um_scale)
BK7_gam3 = 0
BK7_sig3 = 1.01046945

BK7_susc = [LorentzianSusceptibility(frequency = BK7_frq1, gamma = BK7_gam1, sigma = BK7_sig1),
    LorentzianSusceptibility(frequency = BK7_frq2, gamma = BK7_gam2, sigma = BK7_sig2),
    LorentzianSusceptibility(frequency = BK7_frq3, gamma = BK7_gam3, sigma = BK7_sig3)]

BK7 = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = BK7_susc, valid_freq_range = BK7_range))

#------------------------------------------------------------------
# fused quartz (silica) from I.H. Malitson, J. Optical Society of America, Vol. 55, pp. 1205-9, 1965
# ref: https://refractiveindex.info/?shelf=glass&book=fused_silica&page=Malitson
# wavelength range: 0.21 - 6.7 um

fused_quartz_range = FreqRange(min = um_scale / 6.7, max = um_scale / 0.21)

fused_quartz_frq1 = 1 / (0.0684043 * um_scale)
fused_quartz_gam1 = 0
fused_quartz_sig1 = 0.696166300
fused_quartz_frq2 = 1 / (0.1162414 * um_scale)
fused_quartz_gam2 = 0
fused_quartz_sig2 = 0.407942600
fused_quartz_frq3 = 1 / (9.896161 * um_scale)
fused_quartz_gam3 = 0
fused_quartz_sig3 = 0.897479400

fused_quartz_susc = [LorentzianSusceptibility(frequency = fused_quartz_frq1, gamma = fused_quartz_gam1, sigma = fused_quartz_sig1),
    LorentzianSusceptibility(frequency = fused_quartz_frq2, gamma = fused_quartz_gam2, sigma = fused_quartz_sig2),
    LorentzianSusceptibility(frequency = fused_quartz_frq3, gamma = fused_quartz_gam3, sigma = fused_quartz_sig3)]

fused_quartz = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = fused_quartz_susc, valid_freq_range = fused_quartz_range))

#------------------------------------------------------------------
# gallium arsenide (GaAs) from T. Skauli et al., J. Applied Physics, Vol. 94, pp. 6447-55, 2003
# ref: https://refractiveindex.info/?shelf=main&book=GaAs&page=Skauli
# wavelength range: 0.97 - 17 um

GaAs_range = FreqRange(min = um_scale / 17, max = um_scale / 0.97)

GaAs_frq1 = 1 / (0.4431307 * um_scale)
GaAs_gam1 = 0
GaAs_sig1 = 5.466742
GaAs_frq2 = 1 / (0.8746453 * um_scale)
GaAs_gam2 = 0
GaAs_sig2 = 0.02429960
GaAs_frq3 = 1 / (36.9166 * um_scale)
GaAs_gam3 = 0
GaAs_sig3 = 1.957522

GaAs_susc = [LorentzianSusceptibility(frequency = GaAs_frq1, gamma = GaAs_gam1, sigma = GaAs_sig1),
    LorentzianSusceptibility(frequency = GaAs_frq2, gamma = GaAs_gam2, sigma = GaAs_sig2),
    LorentzianSusceptibility(frequency = GaAs_frq3, gamma = GaAs_gam3, sigma = GaAs_sig3)]

GaAs = Medium(epsilon_model = DispersionModel(inf = 5.372514, susceptibilities = GaAs_susc, valid_freq_range = GaAs_range))

#------------------------------------------------------------------
# silicon nitride (Si3N4) from H. R. Philipp, J. Electrochemical Society 120, 295-300, 1973
# ref: https://refractiveindex.info/?shelf=main&book=Si3N4&page=Philipp
# wavelength range: 0.207 - 1.24 um

Si3N4_VISNIR_range = FreqRange(min = um_scale / 1.24, max = um_scale / 0.207)

Si3N4_VISNIR_frq1 = 1 / (0.13967 * um_scale)
Si3N4_VISNIR_gam1 = 0
Si3N4_VISNIR_sig1 = 2.8939

Si3N4_VISNIR_susc = [LorentzianSusceptibility(frequency = Si3N4_VISNIR_frq1, gamma = Si3N4_VISNIR_gam1, sigma = Si3N4_VISNIR_sig1)]

Si3N4_VISNIR = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Si3N4_VISNIR_susc, valid_freq_range = Si3N4_VISNIR_range))

#------------------------------------------------------------------
# silicon nitride (Si3N4) from K. Luke, et. al., Optics Letters, Vol. 40, pp. 4823-26, 2015
# ref: https://refractiveindex.info/?shelf=main&book=Si3N4&page=Luke
# wavelength range: 0.310 - 5.504 um

Si3N4_NIR_range = FreqRange(min = um_scale / 5.504, max = um_scale / 0.310)

Si3N4_NIR_frq1 = 1 / (0.1353406 * um_scale)
Si3N4_NIR_gam1 = 0
Si3N4_NIR_sig1 = 3.0249
Si3N4_NIR_frq2 = 1 / (1239.842 * um_scale)
Si3N4_NIR_gam2 = 0
Si3N4_NIR_sig2 = 40314

Si3N4_NIR_susc = [LorentzianSusceptibility(frequency = Si3N4_NIR_frq1, gamma = Si3N4_NIR_gam1, sigma = Si3N4_NIR_sig1),
    LorentzianSusceptibility(frequency = Si3N4_NIR_frq2, gamma = Si3N4_NIR_gam2, sigma = Si3N4_NIR_sig2)]

Si3N4_NIR = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Si3N4_NIR_susc, valid_freq_range = Si3N4_NIR_range))

#------------------------------------------------------------------
# elemental metals from A.D. Rakic et al., Applied Optics, Vol. 37, No. 22, pp. 5271-83 (1998)
# wavelength range: 0.2 - 12.4 um

metal_range = FreqRange(min = um_scale / 12.398, max = um_scale / 0.24797)

# silver (Ag)

Ag_plasma_frq = 9.01 * eV_um_scale
Ag_f0 = 0.845
Ag_frq0 = 1e-10
Ag_gam0 = 0.048 * eV_um_scale
Ag_sig0 = Ag_f0 * Ag_plasma_frq^2 / Ag_frq0^2
Ag_f1 = 0.065
Ag_frq1 = 0.816 * eV_um_scale      # 1.519 um
Ag_gam1 = 3.886 * eV_um_scale
Ag_sig1 = Ag_f1 * Ag_plasma_frq^2 / Ag_frq1^2
Ag_f2 = 0.124
Ag_frq2 = 4.481 * eV_um_scale      # 0.273 um
Ag_gam2 = 0.452 * eV_um_scale
Ag_sig2 = Ag_f2 * Ag_plasma_frq^2 / Ag_frq2^2
Ag_f3 = 0.011
Ag_frq3 = 8.185 * eV_um_scale      # 0.152 um
Ag_gam3 = 0.065 * eV_um_scale
Ag_sig3 = Ag_f3 * Ag_plasma_frq^2 / Ag_frq3^2
Ag_f4 = 0.840
Ag_frq4 = 9.083 * eV_um_scale      # 0.137 um
Ag_gam4 = 0.916 * eV_um_scale
Ag_sig4 = Ag_f4 * Ag_plasma_frq^2 / Ag_frq4^2
Ag_f5 = 5.646
Ag_frq5 = 20.29 * eV_um_scale      # 0.061 um
Ag_gam5 = 2.419 * eV_um_scale
Ag_sig5 = Ag_f5 * Ag_plasma_frq^2 / Ag_frq5^2

Ag_susc = [DrudeSusceptibility(frequency = Ag_frq0, gamma = Ag_gam0, sigma = Ag_sig0),
    LorentzianSusceptibility(frequency = Ag_frq1, gamma = Ag_gam1, sigma = Ag_sig1),
    LorentzianSusceptibility(frequency = Ag_frq2, gamma = Ag_gam2, sigma = Ag_sig2),
    LorentzianSusceptibility(frequency = Ag_frq3, gamma = Ag_gam3, sigma = Ag_sig3),
    LorentzianSusceptibility(frequency = Ag_frq4, gamma = Ag_gam4, sigma = Ag_sig4),
    LorentzianSusceptibility(frequency = Ag_frq5, gamma = Ag_gam5, sigma = Ag_sig5)]

Ag = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Ag_susc, valid_freq_range = metal_range))

#------------------------------------------------------------------
# gold (Au)

metal_range = FreqRange(min = um_scale / 6.1992, max = um_scale / 0.24797)

Au_plasma_frq = 9.03 * eV_um_scale
Au_f0 = 0.760
Au_frq0 = 1e-10
Au_gam0 = 0.053 * eV_um_scale
Au_sig0 = Au_f0 * Au_plasma_frq^2 / Au_frq0^2
Au_f1 = 0.024
Au_frq1 = 0.415 * eV_um_scale      # 2.988 um
Au_gam1 = 0.241 * eV_um_scale
Au_sig1 = Au_f1 * Au_plasma_frq^2 / Au_frq1^2
Au_f2 = 0.010
Au_frq2 = 0.830 * eV_um_scale      # 1.494 um
Au_gam2 = 0.345 * eV_um_scale
Au_sig2 = Au_f2 * Au_plasma_frq^2 / Au_frq2^2
Au_f3 = 0.071
Au_frq3 = 2.969 * eV_um_scale      # 0.418 um
Au_gam3 = 0.870 * eV_um_scale
Au_sig3 = Au_f3 * Au_plasma_frq^2 / Au_frq3^2
Au_f4 = 0.601
Au_frq4 = 4.304 * eV_um_scale      # 0.288 um
Au_gam4 = 2.494 * eV_um_scale
Au_sig4 = Au_f4 * Au_plasma_frq^2 / Au_frq4^2
Au_f5 = 4.384
Au_frq5 = 13.32 * eV_um_scale      # 0.093 um
Au_gam5 = 2.214 * eV_um_scale
Au_sig5 = Au_f5 * Au_plasma_frq^2 / Au_frq5^2

Au_susc = [DrudeSusceptibility(frequency = Au_frq0, gamma = Au_gam0, sigma = Au_sig0),
    LorentzianSusceptibility(frequency = Au_frq1, gamma = Au_gam1, sigma = Au_sig1),
    LorentzianSusceptibility(frequency = Au_frq2, gamma = Au_gam2, sigma = Au_sig2),
    LorentzianSusceptibility(frequency = Au_frq3, gamma = Au_gam3, sigma = Au_sig3),
    LorentzianSusceptibility(frequency = Au_frq4, gamma = Au_gam4, sigma = Au_sig4),
    LorentzianSusceptibility(frequency = Au_frq5, gamma = Au_gam5, sigma = Au_sig5)]

Au = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Au_susc, valid_freq_range = metal_range))

#------------------------------------------------------------------
# copper (Cu)

metal_range = FreqRange(min = um_scale / 12.398, max = um_scale / 0.20664)

Cu_plasma_frq = 10.83 * eV_um_scale
Cu_f0 = 0.575
Cu_frq0 = 1e-10
Cu_gam0 = 0.030 * eV_um_scale
Cu_sig0 = Cu_f0 * Cu_plasma_frq^2 / Cu_frq0^2
Cu_f1 = 0.061
Cu_frq1 = 0.291 * eV_um_scale      # 4.261 um
Cu_gam1 = 0.378 * eV_um_scale
Cu_sig1 = Cu_f1 * Cu_plasma_frq^2 / Cu_frq1^2
Cu_f2 = 0.104
Cu_frq2 = 2.957 * eV_um_scale      # 0.419 um
Cu_gam2 = 1.056 * eV_um_scale
Cu_sig2 = Cu_f2 * Cu_plasma_frq^2 / Cu_frq2^2
Cu_f3 = 0.723
Cu_frq3 = 5.300 * eV_um_scale      # 0.234 um
Cu_gam3 = 3.213 * eV_um_scale
Cu_sig3 = Cu_f3 * Cu_plasma_frq^2 / Cu_frq3^2
Cu_f4 = 0.638
Cu_frq4 = 11.18 * eV_um_scale      # 0.111 um
Cu_gam4 = 4.305 * eV_um_scale
Cu_sig4 = Cu_f4 * Cu_plasma_frq^2 / Cu_frq4^2

Cu_susc = [DrudeSusceptibility(frequency = Cu_frq0, gamma = Cu_gam0, sigma = Cu_sig0),
    LorentzianSusceptibility(frequency = Cu_frq1, gamma = Cu_gam1, sigma = Cu_sig1),
    LorentzianSusceptibility(frequency = Cu_frq2, gamma = Cu_gam2, sigma = Cu_sig2),
    LorentzianSusceptibility(frequency = Cu_frq3, gamma = Cu_gam3, sigma = Cu_sig3),
    LorentzianSusceptibility(frequency = Cu_frq4, gamma = Cu_gam4, sigma = Cu_sig4)]

Cu = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Cu_susc, valid_freq_range = metal_range))

#------------------------------------------------------------------
# aluminum (Al)

Al_plasma_frq = 14.98 * eV_um_scale
Al_f0 = 0.523
Al_frq0 = 1e-10
Al_gam0 = 0.047 * eV_um_scale
Al_sig0 = Al_f0 * Al_plasma_frq^2 / Al_frq0^2
Al_f1 = 0.227
Al_frq1 = 0.162 * eV_um_scale      # 7.654 um
Al_gam1 = 0.333 * eV_um_scale
Al_sig1 = Al_f1 * Al_plasma_frq^2 / Al_frq1^2
Al_f2 = 0.050
Al_frq2 = 1.544 * eV_um_scale      # 0.803 um
Al_gam2 = 0.312 * eV_um_scale
Al_sig2 = Al_f2 * Al_plasma_frq^2 / Al_frq2^2
Al_f3 = 0.166
Al_frq3 = 1.808 * eV_um_scale      # 0.686 um
Al_gam3 = 1.351 * eV_um_scale
Al_sig3 = Al_f3 * Al_plasma_frq^2 / Al_frq3^2
Al_f4 = 0.030
Al_frq4 = 3.473 * eV_um_scale      # 0.357 um
Al_gam4 = 3.382 * eV_um_scale
Al_sig4 = Al_f4 * Al_plasma_frq^2 / Al_frq4^2

Al_susc = [DrudeSusceptibility(frequency = Al_frq0, gamma = Al_gam0, sigma = Al_sig0),
    LorentzianSusceptibility(frequency = Al_frq1, gamma = Al_gam1, sigma = Al_sig1),
    LorentzianSusceptibility(frequency = Al_frq2, gamma = Al_gam2, sigma = Al_sig2),
    LorentzianSusceptibility(frequency = Al_frq3, gamma = Al_gam3, sigma = Al_sig3),
    LorentzianSusceptibility(frequency = Al_frq4, gamma = Al_gam4, sigma = Al_sig4)]

Al = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Al_susc, valid_freq_range = metal_range))

#------------------------------------------------------------------
# beryllium (Be)

Be_plasma_frq = 18.51 * eV_um_scale
Be_f0 = 0.084
Be_frq0 = 1e-10
Be_gam0 = 0.035 * eV_um_scale
Be_sig0 = Be_f0 * Be_plasma_frq^2 / Be_frq0^2
Be_f1 = 0.031
Be_frq1 = 0.100 * eV_um_scale     # 12.398 um
Be_gam1 = 1.664 * eV_um_scale
Be_sig1 = Be_f1 * Be_plasma_frq^2 / Be_frq1^2
Be_f2 = 0.140
Be_frq2 = 1.032 * eV_um_scale      # 1.201 um
Be_gam2 = 3.395 * eV_um_scale
Be_sig2 = Be_f2 * Be_plasma_frq^2 / Be_frq2^2
Be_f3 = 0.530
Be_frq3 = 3.183 * eV_um_scale      # 0.390 um
Be_gam3 = 4.454 * eV_um_scale
Be_sig3 = Be_f3 * Be_plasma_frq^2 / Be_frq3^2
Be_f4 = 0.130
Be_frq4 = 4.604 * eV_um_scale      # 0.269 um
Be_gam4 = 1.802 * eV_um_scale
Be_sig4 = Be_f4 * Be_plasma_frq^2 / Be_frq4^2

Be_susc = [DrudeSusceptibility(frequency = Be_frq0, gamma = Be_gam0, sigma = Be_sig0),
    LorentzianSusceptibility(frequency = Be_frq1, gamma = Be_gam1, sigma = Be_sig1),
    LorentzianSusceptibility(frequency = Be_frq2, gamma = Be_gam2, sigma = Be_sig2),
    LorentzianSusceptibility(frequency = Be_frq3, gamma = Be_gam3, sigma = Be_sig3),
    LorentzianSusceptibility(frequency = Be_frq4, gamma = Be_gam4, sigma = Be_sig4)]

Be = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Be_susc, valid_freq_range = metal_range))

#------------------------------------------------------------------
# chromium (Cr)

Cr_plasma_frq = 10.75 * eV_um_scale
Cr_f0 = 0.168
Cr_frq0 = 1e-10
Cr_gam0 = 0.047 * eV_um_scale
Cr_sig0 = Cr_f0 * Cr_plasma_frq^2 / Cr_frq0^2
Cr_f1 = 0.151
Cr_frq1 = 0.121 * eV_um_scale     # 10.247 um
Cr_gam1 = 3.175 * eV_um_scale
Cr_sig1 = Cr_f1 * Cr_plasma_frq^2 / Cr_frq1^2
Cr_f2 = 0.150
Cr_frq2 = 0.543 * eV_um_scale      # 2.283 um
Cr_gam2 = 1.305 * eV_um_scale
Cr_sig2 = Cr_f2 * Cr_plasma_frq^2 / Cr_frq2^2
Cr_f3 = 1.149
Cr_frq3 = 1.970 * eV_um_scale      # 0.629 um
Cr_gam3 = 2.676 * eV_um_scale
Cr_sig3 = Cr_f3 * Cr_plasma_frq^2 / Cr_frq3^2
Cr_f4 = 0.825
Cr_frq4 = 8.775 * eV_um_scale      # 0.141 um
Cr_gam4 = 1.335 * eV_um_scale
Cr_sig4 = Cr_f4 * Cr_plasma_frq^2 / Cr_frq4^2

Cr_susc = [DrudeSusceptibility(frequency = Cr_frq0, gamma = Cr_gam0, sigma = Cr_sig0),
    LorentzianSusceptibility(frequency = Cr_frq1, gamma = Cr_gam1, sigma = Cr_sig1),
    LorentzianSusceptibility(frequency = Cr_frq2, gamma = Cr_gam2, sigma = Cr_sig2),
    LorentzianSusceptibility(frequency = Cr_frq3, gamma = Cr_gam3, sigma = Cr_sig3),
    LorentzianSusceptibility(frequency = Cr_frq4, gamma = Cr_gam4, sigma = Cr_sig4)]

Cr = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Cr_susc, valid_freq_range = metal_range))

#------------------------------------------------------------------
# nickel (Ni)

Ni_plasma_frq = 15.92 * eV_um_scale
Ni_f0 = 0.096
Ni_frq0 = 1e-10
Ni_gam0 = 0.048 * eV_um_scale
Ni_sig0 = Ni_f0 * Ni_plasma_frq^2 / Ni_frq0^2
Ni_f1 = 0.100
Ni_frq1 = 0.174 * eV_um_scale      # 7.126 um
Ni_gam1 = 4.511 * eV_um_scale
Ni_sig1 = Ni_f1 * Ni_plasma_frq^2 / Ni_frq1^2
Ni_f2 = 0.135
Ni_frq2 = 0.582 * eV_um_scale      # 2.130 um
Ni_gam2 = 1.334 * eV_um_scale
Ni_sig2 = Ni_f2 * Ni_plasma_frq^2 / Ni_frq2^2
Ni_f3 = 0.106
Ni_frq3 = 1.597 * eV_um_scale      # 0.776 um
Ni_gam3 = 2.178 * eV_um_scale
Ni_sig3 = Ni_f3 * Ni_plasma_frq^2 / Ni_frq3^2
Ni_f4 = 0.729
Ni_frq4 = 6.089 * eV_um_scale      # 0.204 um
Ni_gam4 = 6.292 * eV_um_scale
Ni_sig4 = Ni_f4 * Ni_plasma_frq^2 / Ni_frq4^2

Ni_susc = [DrudeSusceptibility(frequency = Ni_frq0, gamma = Ni_gam0, sigma = Ni_sig0),
    LorentzianSusceptibility(frequency = Ni_frq1, gamma = Ni_gam1, sigma = Ni_sig1),
    LorentzianSusceptibility(frequency = Ni_frq2, gamma = Ni_gam2, sigma = Ni_sig2),
    LorentzianSusceptibility(frequency = Ni_frq3, gamma = Ni_gam3, sigma = Ni_sig3),
    LorentzianSusceptibility(frequency = Ni_frq4, gamma = Ni_gam4, sigma = Ni_sig4)]

Ni = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Ni_susc, valid_freq_range = metal_range))

#------------------------------------------------------------------
# palladium (Pd)

Pd_plasma_frq = 9.72 * eV_um_scale
Pd_f0 = 0.330
Pd_frq0 = 1e-10
Pd_gam0 = 0.008 * eV_um_scale
Pd_sig0 = Pd_f0 * Pd_plasma_frq^2 / Pd_frq0^2
Pd_f1 = 0.649
Pd_frq1 = 0.336 * eV_um_scale      # 3.690 um
Pd_gam1 = 2.950 * eV_um_scale
Pd_sig1 = Pd_f1 * Pd_plasma_frq^2 / Pd_frq1^2
Pd_f2 = 0.121
Pd_frq2 = 0.501 * eV_um_scale      # 2.475 um
Pd_gam2 = 0.555 * eV_um_scale
Pd_sig2 = Pd_f2 * Pd_plasma_frq^2 / Pd_frq2^2
Pd_f3 = 0.638
Pd_frq3 = 1.659 * eV_um_scale      # 0.747 um
Pd_gam3 = 4.621 * eV_um_scale
Pd_sig3 = Pd_f3 * Pd_plasma_frq^2 / Pd_frq3^2
Pd_f4 = 0.453
Pd_frq4 = 5.715 * eV_um_scale      # 0.217 um
Pd_gam4 = 3.236 * eV_um_scale
Pd_sig4 = Pd_f4 * Pd_plasma_frq^2 / Pd_frq4^2

Pd_susc = [DrudeSusceptibility(frequency = Pd_frq0, gamma = Pd_gam0, sigma = Pd_sig0),
    LorentzianSusceptibility(frequency = Pd_frq1, gamma = Pd_gam1, sigma = Pd_sig1),
    LorentzianSusceptibility(frequency = Pd_frq2, gamma = Pd_gam2, sigma = Pd_sig2),
    LorentzianSusceptibility(frequency = Pd_frq3, gamma = Pd_gam3, sigma = Pd_sig3),
    LorentzianSusceptibility(frequency = Pd_frq4, gamma = Pd_gam4, sigma = Pd_sig4)]

Pd = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Pd_susc, valid_freq_range = metal_range))

#------------------------------------------------------------------
# platinum (Pt)

Pt_plasma_frq = 9.59 * eV_um_scale
Pt_f0 = 0.333
Pt_frq0 = 1e-10
Pt_gam0 = 0.080 * eV_um_scale
Pt_sig0 = Pt_f0 * Pt_plasma_frq^2 / Pt_frq0^2
Pt_f1 = 0.191
Pt_frq1 = 0.780 * eV_um_scale      # 1.590 um
Pt_gam1 = 0.517 * eV_um_scale
Pt_sig1 = Pt_f1 * Pt_plasma_frq^2 / Pt_frq1^2
Pt_f2 = 0.659
Pt_frq2 = 1.314 * eV_um_scale      # 0.944 um
Pt_gam2 = 1.838 * eV_um_scale
Pt_sig2 = Pt_f2 * Pt_plasma_frq^2 / Pt_frq2^2
Pt_f3 = 0.547
Pt_frq3 = 3.141 * eV_um_scale      # 0.395 um
Pt_gam3 = 3.668 * eV_um_scale
Pt_sig3 = Pt_f3 * Pt_plasma_frq^2 / Pt_frq3^2
Pt_f4 = 3.576
Pt_frq4 = 9.249 * eV_um_scale      # 0.134 um
Pt_gam4 = 8.517 * eV_um_scale
Pt_sig4 = Pt_f4 * Pt_plasma_frq^2 / Pt_frq4^2

Pt_susc = [DrudeSusceptibility(frequency = Pt_frq0, gamma = Pt_gam0, sigma = Pt_sig0),
    LorentzianSusceptibility(frequency = Pt_frq1, gamma = Pt_gam1, sigma = Pt_sig1),
    LorentzianSusceptibility(frequency = Pt_frq2, gamma = Pt_gam2, sigma = Pt_sig2),
    LorentzianSusceptibility(frequency = Pt_frq3, gamma = Pt_gam3, sigma = Pt_sig3),
    LorentzianSusceptibility(frequency = Pt_frq4, gamma = Pt_gam4, sigma = Pt_sig4)]

Pt = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Pt_susc, valid_freq_range = metal_range))

#------------------------------------------------------------------
# titanium (Ti)

Ti_plasma_frq = 7.29 * eV_um_scale
Ti_f0 = 0.148
Ti_frq0 = 1e-10
Ti_gam0 = 0.082 * eV_um_scale
Ti_sig0 = Ti_f0 * Ti_plasma_frq^2 / Ti_frq0^2
Ti_f1 = 0.899
Ti_frq1 = 0.777 * eV_um_scale      # 1.596 um
Ti_gam1 = 2.276 * eV_um_scale
Ti_sig1 = Ti_f1 * Ti_plasma_frq^2 / Ti_frq1^2
Ti_f2 = 0.393
Ti_frq2 = 1.545 * eV_um_scale      # 0.802 um
Ti_gam2 = 2.518 * eV_um_scale
Ti_sig2 = Ti_f2 * Ti_plasma_frq^2 / Ti_frq2^2
Ti_f3 = 0.187
Ti_frq3 = 2.509 * eV_um_scale      # 0.494 um
Ti_gam3 = 1.663 * eV_um_scale
Ti_sig3 = Ti_f3 * Ti_plasma_frq^2 / Ti_frq3^2
Ti_f4 = 0.001
Ti_frq4 = 19.43 * eV_um_scale      # 0.064 um
Ti_gam4 = 1.762 * eV_um_scale
Ti_sig4 = Ti_f4 * Ti_plasma_frq^2 / Ti_frq4^2

Ti_susc = [DrudeSusceptibility(frequency = Ti_frq0, gamma = Ti_gam0, sigma = Ti_sig0),
    LorentzianSusceptibility(frequency = Ti_frq1, gamma = Ti_gam1, sigma = Ti_sig1),
    LorentzianSusceptibility(frequency = Ti_frq2, gamma = Ti_gam2, sigma = Ti_sig2),
    LorentzianSusceptibility(frequency = Ti_frq3, gamma = Ti_gam3, sigma = Ti_sig3),
    LorentzianSusceptibility(frequency = Ti_frq4, gamma = Ti_gam4, sigma = Ti_sig4)]

Ti = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Ti_susc, valid_freq_range = metal_range))

#------------------------------------------------------------------
# tungsten (W)

W_plasma_frq = 13.22 * eV_um_scale
W_f0 = 0.206
W_frq0 = 1e-10
W_gam0 = 0.064 * eV_um_scale
W_sig0 = W_f0 * W_plasma_frq^2 / W_frq0^2
W_f1 = 0.054
W_frq1 = 1.004 * eV_um_scale      # 1.235 um
W_gam1 = 0.530 * eV_um_scale
W_sig1 = W_f1 * W_plasma_frq^2 / W_frq1^2
W_f2 = 0.166
W_frq2 = 1.917 * eV_um_scale      # 0.647 um
W_gam2 = 1.281 * eV_um_scale
W_sig2 = W_f2 * W_plasma_frq^2 / W_frq2^2
W_f3 = 0.706
W_frq3 = 3.580 * eV_um_scale      # 0.346 um
W_gam3 = 3.332 * eV_um_scale
W_sig3 = W_f3 * W_plasma_frq^2 / W_frq3^2
W_f4 = 2.590
W_frq4 = 7.498 * eV_um_scale      # 0.165 um
W_gam4 = 5.836 * eV_um_scale
W_sig4 = W_f4 * W_plasma_frq^2 / W_frq4^2

W_susc = [DrudeSusceptibility(frequency = W_frq0, gamma = W_gam0, sigma = W_sig0),
    LorentzianSusceptibility(frequency = W_frq1, gamma = W_gam1, sigma = W_sig1),
    LorentzianSusceptibility(frequency = W_frq2, gamma = W_gam2, sigma = W_sig2),
    LorentzianSusceptibility(frequency = W_frq3, gamma = W_gam3, sigma = W_sig3),
    LorentzianSusceptibility(frequency = W_frq4, gamma = W_gam4, sigma = W_sig4)]

W = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = W_susc, valid_freq_range = metal_range))

#------------------------------------------------------------------
# metals from D. Barchiesi and T. Grosges, J. Nanophotonics, Vol. 8, 08996, 2015
# wavelength range: 0.4 - 0.8 um

metal_visible_range = FreqRange(min = um_scale / 0.8, max = um_scale / 0.4)

# gold (Au)
# fit to P.B. Johnson and R.W. Christy, Physical Review B, Vol. 6, pp. 4370-9, 1972

Au_JC_visible_frq0 = 1 / (0.139779231751333 * um_scale)
Au_JC_visible_gam0 = 1 / (26.1269913352870 * um_scale)
Au_JC_visible_sig0 = 1

Au_JC_visible_frq1 = 1 / (0.404064525036786 * um_scale)
Au_JC_visible_gam1 = 1 / (1.12834046202759 * um_scale)
Au_JC_visible_sig1 = 2.07118534879440

Au_JC_visible_susc = [DrudeSusceptibility(frequency = Au_JC_visible_frq0, gamma = Au_JC_visible_gam0, sigma = Au_JC_visible_sig0),
    LorentzianSusceptibility(frequency = Au_JC_visible_frq1, gamma = Au_JC_visible_gam1, sigma = Au_JC_visible_sig1)]

Au_JC_visible = Medium(epsilon_model = DispersionModel(inf = 6.1599, susceptibilities = Au_JC_visible_susc, valid_freq_range = metal_visible_range))

#------------------------------------------------------------------
# gold (Au)
# fit to E.D. Palik, Handbook of Optical Constants, Academic Press, 1985

Au_visible_frq0 = 1 / (0.0473629248511456 * um_scale)
Au_visible_gam0 = 1 / (0.255476199605166 * um_scale)
Au_visible_sig0 = 1

Au_visible_frq1 = 1 / (0.800619321082804 * um_scale)
Au_visible_gam1 = 1 / (0.381870287531951 * um_scale)
Au_visible_sig1 = -169.060953137985

Au_visible_susc = [DrudeSusceptibility(frequency = Au_visible_frq0, gamma = Au_visible_gam0, sigma = Au_visible_sig0),
    LorentzianSusceptibility(frequency = Au_visible_frq1, gamma = Au_visible_gam1, sigma = Au_visible_sig1)]

Au_visible = Medium(epsilon_model = DispersionModel(inf = 0.6888, susceptibilities = Au_visible_susc, valid_freq_range = metal_visible_range))

#------------------------------------------------------------------
## WARNING: unstable; field divergence may occur

# silver (Au)
# fit to E.D. Palik, Handbook of Optical Constants, Academic Press, 1985

Ag_visible_frq0 = 1 / (0.142050162130618 * um_scale)
Ag_visible_gam0 = 1 / (18.0357292925015 * um_scale)
Ag_visible_sig0 = 1

Ag_visible_frq1 = 1 / (0.115692151792108 * um_scale)
Ag_visible_gam1 = 1 / (0.257794324096575 * um_scale)
Ag_visible_sig1 = 3.74465275944019

Ag_visible_susc = [DrudeSusceptibility(frequency = Ag_visible_frq0, gamma = Ag_visible_gam0, sigma = Ag_visible_sig0),
    LorentzianSusceptibility(frequency = Ag_visible_frq1, gamma = Ag_visible_gam1, sigma = Ag_visible_sig1)]

Ag_visible = Medium(epsilon_model = DispersionModel(inf = 0.0067526, susceptibilities = Ag_visible_susc, valid_freq_range = metal_visible_range))

#------------------------------------------------------------------
## WARNING: unstable; field divergence may occur

# aluminum (Al)
# fit to E.D. Palik, Handbook of Optical Constants, Academic Press, 1985

Al_visible_frq0 = 1 / (0.0625841659042985 * um_scale)
Al_visible_gam0 = 1 / (0.606007002962666 * um_scale)
Al_visible_sig0 = 1

Al_visible_frq1 = 1 / (0.528191199577075 * um_scale)
Al_visible_gam1 = 1 / (0.291862527666814 * um_scale)
Al_visible_sig1 = -44.4456675577921

Al_visible_susc = [DrudeSusceptibility(frequency = Al_visible_frq0, gamma = Al_visible_gam0, sigma = Al_visible_sig0),
    LorentzianSusceptibility(frequency = Al_visible_frq1, gamma = Al_visible_gam1, sigma = Al_visible_sig1)]

Al_visible = Medium(epsilon_model = DispersionModel(inf = 0.13313, susceptibilities = Al_visible_susc, valid_freq_range = metal_visible_range))

#------------------------------------------------------------------
# chroimium (Cr)
# fit to E.D. Palik, Handbook of Optical Constants, Academic Press, 1985

Cr_visible_frq0 = 1 / (0.118410119507342 * um_scale)
Cr_visible_gam0 = 1 / (0.628596264869804 * um_scale)
Cr_visible_sig0 = 1

Cr_visible_frq1 = 1 / (0.565709598452496 * um_scale)
Cr_visible_gam1 = 1 / (0.731117670900812 * um_scale)
Cr_visible_sig1 = 13.2912419951294

Cr_visible_susc = [DrudeSusceptibility(frequency = Cr_visible_frq0, gamma = Cr_visible_gam0, sigma = Cr_visible_sig0),
    LorentzianSusceptibility(frequency = Cr_visible_frq1, gamma = Cr_visible_gam1, sigma = Cr_visible_sig1)]

Cr_visible = Medium(epsilon_model = DispersionModel(inf = 2.7767, susceptibilities = Cr_visible_susc, valid_freq_range = metal_visible_range))

#------------------------------------------------------------------
## WARNING: unstable; field divergence may occur

# titanium (Ti)
# fit to E.D. Palik, Handbook of Optical Constants, Academic Press, 1985

Ti_visible_frq0 = 1 / (0.101331651921602 * um_scale)
Ti_visible_gam0 = 1 / (0.365743382258719 * um_scale)
Ti_visible_sig0 = 1

Ti_visible_frq1 = 1 / (4.56839173979216e-09 * um_scale)
Ti_visible_gam1 = 1 / (5.86441957443603e-10 * um_scale)
Ti_visible_sig1 = 54742662.1963414

Ti_visible_susc = [DrudeSusceptibility(frequency = Ti_visible_frq0, gamma = Ti_visible_gam0, sigma = Ti_visible_sig0),
    LorentzianSusceptibility(frequency = Ti_visible_frq1, gamma = Ti_visible_gam1, sigma = Ti_visible_sig1)]

Ti_visible = Medium(epsilon_model = DispersionModel(inf = -5.4742e7, susceptibilities = Ti_visible_susc, valid_freq_range = metal_visible_range))

#------------------------------------------------------------------
# aluminum (Al) from Horiba Technical Note 09: Drude Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Drude_Dispersion_Model.pdf
# wavelength range: 0.19 - 0.83 um

Al_drude_range = FreqRange(min = um_scale / 0.83, max = um_scale / 0.19)

Al_drude_frq = 1 / (0.0789607648707171 * um_scale)
Al_drude_gam = 1 / (1.78138208333333 * um_scale)
Al_drude_sig = 1

Al_drude_susc = [DrudeSusceptibility(frequency = Al_drude_frq, gamma = Al_drude_gam, sigma = Al_drude_sig)]

Al_drude = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Al_drude_susc, valid_freq_range = Al_drude_range))

#------------------------------------------------------------------
# cobalt (Co) from Horiba Technical Note 09: Drude Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Drude_Dispersion_Model.pdf
# wavelength range: 0.26 - 1.65 um

Co_range = FreqRange(min = um_scale / 1.65, max = um_scale / 0.26)

Co_frq = 1 / (0.0789607648707171 * um_scale)
Co_gam = 1 / (0.213802712536644 * um_scale)
Co_sig = 1

Co_susc = [DrudeSusceptibility(frequency = Co_frq, gamma = Co_gam, sigma = Co_sig)]

Co = Medium(epsilon_model = DispersionModel(inf = 3.694, susceptibilities = Co_susc, valid_freq_range = Co_range))

#------------------------------------------------------------------
## WARNING: unstable; field divergence may occur

# molybdenum (Mo) from Horiba Technical Note 09: Drude Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Drude_Dispersion_Model.pdf
# wavelength range: 0.25 - 0.83 um

Mo_range = FreqRange(min = um_scale / 0.83, max = um_scale / 0.25)

Mo_frq = 1 / (0.0620790071099539 * um_scale)
Mo_gam = 1 / (0.148359690080172 * um_scale)
Mo_sig = 1

Mo_susc = [DrudeSusceptibility(frequency = Mo_frq, gamma = Mo_gam, sigma = Mo_sig)]

Mo = Medium(epsilon_model = DispersionModel(inf = -1.366, susceptibilities = Mo_susc, valid_freq_range = Mo_range))

#------------------------------------------------------------------
# nickel chrome (NiCr) from Horiba Technical Note 09: Drude Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Drude_Dispersion_Model.pdf
# wavelength range: 0.25 - 0.83 um

NiCr_range = FreqRange(min = um_scale / 0.83, max = um_scale / 0.25)

NiCr_frq = 1 / (0.0868845080588648 * um_scale)
NiCr_gam = 1 / (0.308418390547264 * um_scale)
NiCr_sig = 1

NiCr_susc = [DrudeSusceptibility(frequency = NiCr_frq, gamma = NiCr_gam, sigma = NiCr_sig)]

NiCr = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = NiCr_susc, valid_freq_range = NiCr_range))

#------------------------------------------------------------------
# nickel iron (NiFe) from Horiba Technical Note 09: Drude Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Drude_Dispersion_Model.pdf
# wavelength range: 0.25 - 0.83 um

NiFe_range = FreqRange(min = um_scale / 0.83, max = um_scale / 0.25)

NiFe_frq = 1 / (0.0838297450980392 * um_scale)
NiFe_gam = 1 / (0.259381156903766 * um_scale)
NiFe_sig = 1

NiFe_susc = [DrudeSusceptibility(frequency = NiFe_frq, gamma = NiFe_gam, sigma = NiFe_sig)]

NiFe = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = NiFe_susc, valid_freq_range = NiFe_range))

#------------------------------------------------------------------
# titanium (Ti) from Horiba Technical Note 09: Drude Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Drude_Dispersion_Model.pdf
# wavelength range: 0.21 - 1.24 um

Ti_drude_range = FreqRange(min = um_scale / 1.24, max = um_scale / 0.21)

Ti_drude_frq = 1 / (0.113746966055046 * um_scale)
Ti_drude_gam = 1 / (0.490056098814229 * um_scale)
Ti_drude_sig = 1

Ti_drude_susc = [DrudeSusceptibility(frequency = Ti_drude_frq, gamma = Ti_drude_gam, sigma = Ti_drude_sig)]

Ti_drude = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Ti_drude_susc, valid_freq_range = Ti_drude_range))

#------------------------------------------------------------------
# silicon nitride (SiN), non-stoichiometric, from Horiba Technical Note 08: Lorentz Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Lorentz_Dispersion_Model.pdf
# wavelength range: 0.21 - 2.07 um

SiN_range = FreqRange(min = um_scale / 2.07, max = um_scale / 0.21)

SiN_frq1 = 1 / (0.190891752117013 * um_scale)
SiN_gam1 = 1 / (3.11518072864322 * um_scale)
SiN_sig1 = 1.2650

SiN_susc = [LorentzianSusceptibility(frequency = SiN_frq1, gamma = SiN_gam1, sigma = SiN_sig1)]

SiN = Medium(epsilon_model = DispersionModel(inf = 2.320, susceptibilities = SiN_susc, valid_freq_range = SiN_range))

#------------------------------------------------------------------
# silicon nitride (Si3N4), stoichiometric, from Horiba Technical Note 08: Lorentz Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Lorentz_Dispersion_Model.pdf
# wavelength range: 0.23 - 0.83 um

Si3N4_range = FreqRange(min = um_scale / 0.83, max = um_scale / 0.23)

Si3N4_frq1 = 1 / (0.389153148148148 * um_scale)
Si3N4_gam1 = 1 / (0.693811936205932 * um_scale)
Si3N4_sig1 = 4.377

Si3N4_susc = [LorentzianSusceptibility(frequency = Si3N4_frq1, gamma = Si3N4_gam1, sigma = Si3N4_sig1)]

Si3N4 = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Si3N4_susc, valid_freq_range = Si3N4_range))

#------------------------------------------------------------------
# silicon dioxide (SiO2) from Horiba Technical Note 08: Lorentz Dispersion Model
# ref: http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Lorentz_Dispersion_Model.pdf
# wavelength range: 0.25 - 1.77 um

SiO2_range = FreqRange(min = um_scale / 1.77, max = um_scale / 0.25)

SiO2_frq1 = 1 / (0.103320160833333 * um_scale)
SiO2_gam1 = 1 / (12.3984193000000 * um_scale)
SiO2_sig1 = 1.12

SiO2_susc = [LorentzianSusceptibility(frequency = SiO2_frq1, gamma = SiO2_gam1, sigma = SiO2_sig1)]

SiO2 = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = SiO2_susc, valid_freq_range = SiO2_range))

#------------------------------------------------------------------
# indium phosphide (InP) from Handbook of Optics, 2nd edition, Vol. 2, McGraw-Hill, 1994
# ref: https://refractiveindex.info/?shelf=main&book=InP&page=Pettit
# wavelength range: 0.95 - 10 um

InP_range = FreqRange(min = um_scale / 10, max = um_scale / 0.95)

InP_frq1 = 1 / (0.6263 * um_scale)
InP_gam1 = 0
InP_sig1 = 2.316
InP_frq2 = 1 / (32.935 * um_scale)
InP_gam2 = 0
InP_sig2 = 2.765

InP_susc = [LorentzianSusceptibility(frequency = InP_frq1, gamma = InP_gam1, sigma = InP_sig1),
    LorentzianSusceptibility(frequency = InP_frq2, gamma = InP_gam2, sigma = InP_sig2)]

InP = Medium(epsilon_model = DispersionModel(inf = 7.255, susceptibilities = InP_susc, valid_freq_range = InP_range))

#------------------------------------------------------------------
# germanium (Ge) from N. P. Barnes and M. S. Piltch, J. Optical Society America, Vol. 69, pp. 178-180 (1979)
# ref: https://refractiveindex.info/?shelf=main&book=Ge&page=Icenogle
# wavelength range: 2.5 - 12 um

Ge_range = FreqRange(min = um_scale / 12, max = um_scale / 2.5)

Ge_frq1 = 1 / (0.6641159 * um_scale)
Ge_gam1 = 0
Ge_sig1 = 6.7288
Ge_frq2 = 1 / (62.210127 * um_scale)
Ge_gam2 = 0
Ge_sig2 = 0.21307

Ge_susc = [LorentzianSusceptibility(frequency = Ge_frq1, gamma = Ge_gam1, sigma = Ge_sig1),
    LorentzianSusceptibility(frequency = Ge_frq2, gamma = Ge_gam2, sigma = Ge_sig2)]

Ge = Medium(epsilon_model = DispersionModel(inf = 9.28156, susceptibilities = Ge_susc, valid_freq_range = Ge_range))

#------------------------------------------------------------------
# silicon (Si) from C. D. Salzberg and J. J. Villa, , J. Optical Society America, Vol. 47, pp. 244-246 (1957)
# ref: https://refractiveindex.info/?shelf=main&book=Si&page=Salzberg
# wavelength range: 1.36 - 11 um

Si_range = FreqRange(min = um_scale / 11, max = um_scale / 1.36)

Si_frq1 = 1 / (0.301516485 * um_scale)
Si_gam1 = 0
Si_sig1 = 10.6684293
Si_frq2 = 1 / (1.13475115 * um_scale)
Si_gam2 = 0
Si_sig2 = 0.0030434748
Si_frq3 = 1 / (1104 * um_scale)
Si_gam3 = 0
Si_sig3 = 1.54133408

Si_susc = [LorentzianSusceptibility(frequency = Si_frq1, gamma = Si_gam1, sigma = Si_sig1),
    LorentzianSusceptibility(frequency = Si_frq2, gamma = Si_gam2, sigma = Si_sig2),
    LorentzianSusceptibility(frequency = Si_frq3, gamma = Si_gam3, sigma = Si_sig3)]

Si = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Si_susc, valid_freq_range = Si_range))

#------------------------------------------------------------------
# poly(methyl methacrylate) (PMMA) from N. Sultanova et al., Acta Physica Polonica A, Vol. 116, pp. 585-7 (2009)
# ref: https://refractiveindex.info/?shelf=organic&book=poly%28methyl_methacrylate%29&page=Sultanova
# wavelength range: 0.437 - 1.052 um

PMMA_range = FreqRange(min = um_scale / 1.052, max = um_scale / 0.437)

PMMA_frq1 = 1 / (0.106362587407415 * um_scale)
PMMA_gam1 = 0
PMMA_sig1 = 1.1819

PMMA_susc = [LorentzianSusceptibility(frequency = PMMA_frq1, gamma = PMMA_gam1, sigma = PMMA_sig1)]

PMMA = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = PMMA_susc, valid_freq_range = PMMA_range))

#------------------------------------------------------------------
# polycarbonate (PC) from N. Sultanova et al., Acta Physica Polonica A, Vol. 116, pp. 585-7 (2009)
# ref: https://refractiveindex.info/?shelf=organic&book=polycarbonate&page=Sultanova
# wavelength range: 0.437 - 1.052 um

PC_range = FreqRange(min = um_scale / 1.052, max = um_scale / 0.437)

PC_frq1 = 1 / (0.145958898324152 * um_scale)
PC_gam1 = 0
PC_sig1 = 1.4182

PC_susc = [LorentzianSusceptibility(frequency = PC_frq1, gamma = PC_gam1, sigma = PC_sig1)]

PC = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = PC_susc, valid_freq_range = PC_range))

#------------------------------------------------------------------
# polystyrene (PS) from N. Sultanova et al., Acta Physica Polonica A, Vol. 116, pp. 585-7 (2009)
# ref: https://refractiveindex.info/?shelf=organic&book=polystyren&page=Sultanova
# wavelength range: 0.437 - 1.052 um

PS_range = FreqRange(min = um_scale / 1.052, max = um_scale / 0.437)

PS_frq1 = 1 / (0.142182980697410 * um_scale)
PS_gam1 = 0
PS_sig1 = 1.4435

PS_susc = [LorentzianSusceptibility(frequency = PS_frq1, gamma = PS_gam1, sigma = PS_sig1)]

PS = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = PS_susc, valid_freq_range = PS_range))

#------------------------------------------------------------------
# cellulose (CLS) from N. Sultanova et al., Acta Physica Polonica A, Vol. 116, pp. 585-7 (2009)
# ref: https://refractiveindex.info/?shelf=organic&book=cellulose&page=Sultanova
# wavelength range: 0.437 - 1.052 um

CLS_range = FreqRange(min = um_scale / 1.052, max = um_scale / 0.437)

CLS_frq1 = 1 / (0.105294824184287 * um_scale)
CLS_gam1 = 0
CLS_sig1 = 1.124

CLS_susc = [LorentzianSusceptibility(frequency = CLS_frq1, gamma = CLS_gam1, sigma = CLS_sig1)]

CLS = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = CLS_susc, valid_freq_range = CLS_range))

#------------------------------------------------------------------
# barium borate (BaB2O4), beta phase, from G. Tamosauskas et al., Optical Materials Express, Vol. 8, pp. 1410-18 (2018)
# ref: https://refractiveindex.info/?shelf=main&book=BaB2O4&page=Tamosauskas-o
# ref: https://refractiveindex.info/?shelf=main&book=BaB2O4&page=Tamosauskas-e
# wavelength range: 0.188 - 5.2 um

## NOTE: ordinary (o) axes in X and Y, extraordinary (e) axis in Z

BaB2O4_range = FreqRange(min = um_scale / 5.2, max = um_scale / 0.188)

BaB2O4_frq1 = 1 / (0.06265780079128216 * um_scale)
BaB2O4_gam1 = 0
BaB2O4_sig1 = 0.90291
BaB2O4_frq2 = 1 / (0.13706202975295528 * um_scale)
BaB2O4_gam2 = 0
BaB2O4_sig2 = 0.83155
BaB2O4_frq3 = 1 / (7.746612162745725 * um_scale)
BaB2O4_gam3 = 0
BaB2O4_sig3 = 0.76536

BaB2O4_susc_o = [LorentzianSusceptibility(frequency = BaB2O4_frq1, gamma = BaB2O4_gam1, sigma_diag = BaB2O4_sig1 * [1, 1, 0]),
    LorentzianSusceptibility(frequency = BaB2O4_frq2, gamma = BaB2O4_gam2, sigma_diag = BaB2O4_sig2 * [1, 1, 0]),
    LorentzianSusceptibility(frequency = BaB2O4_frq3, gamma = BaB2O4_gam3, sigma_diag = BaB2O4_sig3 * [1, 1, 0])]

BaB2O4_frq1 = 1 / (0.0845103543951864 * um_scale)
BaB2O4_gam1 = 0
BaB2O4_sig1 = 1.151075
BaB2O4_frq2 = 1 / (0.15029970059850417 * um_scale)
BaB2O4_gam2 = 0
BaB2O4_sig2 = 0.21803
BaB2O4_frq3 = 1 / (16.217274740226856 * um_scale)
BaB2O4_gam3 = 0
BaB2O4_sig3 = 0.656

BaB2O4_susc_e = [
    LorentzianSusceptibility(frequency = BaB2O4_frq1, gamma = BaB2O4_gam1, sigma_diag = BaB2O4_sig1 * [0, 0, 1]),
    LorentzianSusceptibility(frequency = BaB2O4_frq2, gamma = BaB2O4_gam2, sigma_diag = BaB2O4_sig2 * [0, 0, 1]),
    LorentzianSusceptibility(frequency = BaB2O4_frq3, gamma = BaB2O4_gam3, sigma_diag = BaB2O4_sig3 * [0, 0, 1])
]

BaB2O4 = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = [BaB2O4_susc_o; BaB2O4_susc_e], valid_freq_range = BaB2O4_range))

#------------------------------------------------------------------
# lithium niobate (LiNbO3) from D.E. Zelmon et al., J. Optical Society of America B, Vol. 14, pp. 3319-22 (1997)
# ref: https://refractiveindex.info/?shelf=main&book=LiNbO3&page=Zelmon-o
# ref: https://refractiveindex.info/?shelf=main&book=LiNbO3&page=Zelmon-e
# wavelength range: 0.4 - 5.0 um

## NOTE: ordinary (o) axes in X and Y, extraordinary (e) axis in Z

LiNbO3_range = FreqRange(min = um_scale / 5.0, max = um_scale / 0.4)

LiNbO3_frq1 = 1 / (0.13281566172707193 * um_scale)
LiNbO3_gam1 = 0
LiNbO3_sig1 = 2.6734
LiNbO3_frq2 = 1 / (0.24318717071424636 * um_scale)
LiNbO3_gam2 = 0
LiNbO3_sig2 = 1.2290
LiNbO3_frq3 = 1 / (21.78531615561271 * um_scale)
LiNbO3_gam3 = 0
LiNbO3_sig3 = 12.614

LiNbO3_susc_o = [LorentzianSusceptibility(frequency = LiNbO3_frq1, gamma = LiNbO3_gam1, sigma_diag = LiNbO3_sig1 * [1, 1, 0]),
    LorentzianSusceptibility(frequency = LiNbO3_frq2, gamma = LiNbO3_gam2, sigma_diag = LiNbO3_sig2 * [1, 1, 0]),
    LorentzianSusceptibility(frequency = LiNbO3_frq3, gamma = LiNbO3_gam3, sigma_diag = LiNbO3_sig3 * [1, 1, 0])]

LiNbO3_frq1 = 1 / (0.14307340773183533 * um_scale)
LiNbO3_gam1 = 0
LiNbO3_sig1 = 2.9804
LiNbO3_frq2 = 1 / (0.2580697580112788 * um_scale)
LiNbO3_gam2 = 0
LiNbO3_sig2 = 0.5981
LiNbO3_frq3 = 1 / (20.39803912144498 * um_scale)
LiNbO3_gam3 = 0
LiNbO3_sig3 = 8.9543

LiNbO3_susc_e = [LorentzianSusceptibility(frequency = LiNbO3_frq1, gamma = LiNbO3_gam1, sigma_diag = LiNbO3_sig1 * [0, 0, 1]),
    LorentzianSusceptibility(frequency = LiNbO3_frq2, gamma = LiNbO3_gam2, sigma_diag = LiNbO3_sig2 * [0, 0, 1]),
    LorentzianSusceptibility(frequency = LiNbO3_frq3, gamma = LiNbO3_gam3, sigma_diag = LiNbO3_sig3 * [0, 0, 1])]

LiNbO3 = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = [LiNbO3_susc_o; LiNbO3_susc_e], valid_freq_range = LiNbO3_range))

#------------------------------------------------------------------
# calcium tungstate (CaWO4) from W.L. Bond, J. Applied Physics, Vol. 36, pp. 1674-77 (1965)
# ref: https://refractiveindex.info/?shelf=main&book=CaWO4&page=Bond-o
# ref: https://refractiveindex.info/?shelf=main&book=CaWO4&page=Bond-e
# wavelength range: 0.45 - 4.0 um

## NOTE: ordinary (o) axes in X and Y, extraordinary (e) axis in Z

CaWO4_range = FreqRange(min = um_scale / 4.0, max = um_scale / 0.45)

CaWO4_frq1 = 1 / (0.1347 * um_scale)
CaWO4_gam1 = 0
CaWO4_sig1 = 2.5493
CaWO4_frq2 = 1 / (10.815 * um_scale)
CaWO4_gam2 = 0
CaWO4_sig2 = 0.9200

CaWO4_susc_o = [LorentzianSusceptibility(frequency = CaWO4_frq1, gamma = CaWO4_gam1, sigma_diag = CaWO4_sig1 * [1, 1, 0]),
    LorentzianSusceptibility(frequency = CaWO4_frq2, gamma = CaWO4_gam2, sigma_diag = CaWO4_sig2 * [1, 1, 0])]

CaWO4_frq1 = 1 / (0.1379 * um_scale)
CaWO4_gam1 = 0
CaWO4_sig1 = 2.6041
CaWO4_frq2 = 1 / (21.371 * um_scale)
CaWO4_gam2 = 0
CaWO4_sig2 = 4.1237

CaWO4_susc_e = [LorentzianSusceptibility(frequency = CaWO4_frq1, gamma = CaWO4_gam1, sigma_diag = CaWO4_sig1 * [0, 0, 1]),
    LorentzianSusceptibility(frequency = CaWO4_frq2, gamma = CaWO4_gam2, sigma_diag = CaWO4_sig2 * [0, 0, 1])]

CaWO4 = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = [CaWO4_susc_o; CaWO4_susc_e], valid_freq_range = CaWO4_range))

#------------------------------------------------------------------
# calcium carbonate (CaCO3) from G. Ghosh, Optics Communication, Vol. 163, pp. 95-102 (1999)
# ref: https://refractiveindex.info/?shelf=main&book=CaCO3&page=Ghosh-o
# ref: https://refractiveindex.info/?shelf=main&book=CaCO3&page=Ghosh-e
# wavelength range: 0.204 - 2.172 um

## NOTE: ordinary (o) axes in X and Y, extraordinary (e) axis in Z

CaCO3_range = FreqRange(min = um_scale / 2.172, max = um_scale / 0.204)

CaCO3_frq1 = 1 / (0.13940057496294625 * um_scale)
CaCO3_gam1 = 0
CaCO3_sig1 = 0.96464345
CaCO3_frq2 = 1 / (10.954451150103322 * um_scale)
CaCO3_gam2 = 0
CaCO3_sig2 = 1.82831454

CaCO3_susc_o = [LorentzianSusceptibility(frequency = CaCO3_frq1, gamma = CaCO3_gam1, sigma_diag = CaCO3_sig1 * [1, 1, 0]),
    LorentzianSusceptibility(frequency = CaCO3_frq2, gamma = CaCO3_gam2, sigma_diag = CaCO3_sig2 * [1, 1, 0])]

CaCO3_frq1 = 1 / (0.1032906302623815 * um_scale)
CaCO3_gam1 = 0
CaCO3_sig1 = 0.82427830
CaCO3_frq2 = 1 / (10.954451150103322 * um_scale)
CaCO3_gam2 = 0
CaCO3_sig2 = 0.14429128

CaCO3_susc_e = [LorentzianSusceptibility(frequency = CaCO3_frq1, gamma = CaCO3_gam1, sigma_diag = CaCO3_sig1 * [0, 0, 1]),
    LorentzianSusceptibility(frequency = CaCO3_frq2, gamma = CaCO3_gam2, sigma_diag = CaCO3_sig2 * [0, 0, 1])]

CaCO3 = Medium(epsilon_model = DispersionModel(inf_diag = [1.73358749, 1.73358749, 1.35859695], susceptibilities = [CaCO3_susc_o; CaCO3_susc_e], valid_freq_range = CaCO3_range))

#------------------------------------------------------------------
# silicon dioxide (SiO2) from G. Ghosh, Optics Communication, Vol. 163, pp. 95-102 (1999)
# ref: https://refractiveindex.info/?shelf=main&book=SiO2&page=Ghosh-o
# ref: https://refractiveindex.info/?shelf=main&book=SiO2&page=Ghosh-e
# wavelength range: 0.198 - 2.0531 um

## NOTE: ordinary (o) axes in X and Y, extraordinary (e) axis in Z

SiO2_range = FreqRange(min = um_scale / 2.0531, max = um_scale / 0.198)

SiO2_frq1 = 1 / (0.10029257051247614 * um_scale)
SiO2_gam1 = 0
SiO2_sig1 = 1.07044083
SiO2_frq2 = 1 / (10 * um_scale)
SiO2_gam2 = 0
SiO2_sig2 = 1.10202242

SiO2_susc_o = [LorentzianSusceptibility(frequency = SiO2_frq1, gamma = SiO2_gam1, sigma_diag = SiO2_sig1 * [1, 1, 0]),
    LorentzianSusceptibility(frequency = SiO2_frq2, gamma = SiO2_gam2, sigma_diag = SiO2_sig2 * [1, 1, 0])]

SiO2_frq1 = 1 / (0.10104546699382412 * um_scale)
SiO2_gam1 = 0
SiO2_sig1 = 1.09509924
SiO2_frq2 = 1 / (10 * um_scale)
SiO2_gam2 = 0
SiO2_sig2 = 1.15662475

SiO2_susc_e = [LorentzianSusceptibility(frequency = SiO2_frq1, gamma = SiO2_gam1, sigma_diag = SiO2_sig1 * [0, 0, 1]),
    LorentzianSusceptibility(frequency = SiO2_frq2, gamma = SiO2_gam2, sigma_diag = SiO2_sig2 * [0, 0, 1])]

SiO2_aniso = Medium(epsilon_model = DispersionModel(inf_diag = [1.28604141, 1.28604141, 1.28851804], susceptibilities = [SiO2_susc_o; SiO2_susc_e], valid_freq_range = SiO2_range))

#------------------------------------------------------------------
# gallium nitride (GaN), alpha phase (wurtzite), from A.S. Barker Jr. and M. Ilegems, Physical Review B, Vol. 7, pp. 743-50 (1973)
# ref: https://refractiveindex.info/?shelf=main&book=GaN&page=Barker-o
# ref: https://refractiveindex.info/?shelf=main&book=GaN&page=Barker-e
# wavelength range: 0.35 - 10 um

## NOTE: ordinary (o) axes in X and Y, extraordinary (e) axis in Z

GaN_range = FreqRange(min = um_scale / 10.0, max = um_scale / 0.35)

GaN_frq1 = 1 / (0.256 * um_scale)
GaN_gam1 = 0
GaN_sig1 = 1.75
GaN_frq2 = 1 / (17.86 * um_scale)
GaN_gam2 = 0
GaN_sig2 = 4.1

GaN_susc_o = [LorentzianSusceptibility(frequency = GaN_frq1, gamma = GaN_gam1, sigma_diag = GaN_sig1 * [1, 1, 0]),
    LorentzianSusceptibility(frequency = GaN_frq2, gamma = GaN_gam2, sigma_diag = GaN_sig2 * [1, 1, 0])]

GaN_frq1 = 1 / (18.76 * um_scale)
GaN_gam1 = 0
GaN_sig1 = 5.08

GaN_susc_e = [LorentzianSusceptibility(frequency = GaN_frq1, gamma = GaN_gam1, sigma_diag = GaN_sig1 * [0, 0, 1])]

GaN = Medium(epsilon_model = DispersionModel(inf_diag = [3.6, 3.6, 5.35], susceptibilities = [GaN_susc_o; GaN_susc_e], valid_freq_range = GaN_range))

#------------------------------------------------------------------
# aluminum nitride (AlN) from J. Pastrnak and L. Roskovcova, Physica Status Solidi, Vol. 14, K5-8 (1966)
# ref: https://refractiveindex.info/?shelf=main&book=AlN&page=Pastrnak-o
# ref: https://refractiveindex.info/?shelf=main&book=AlN&page=Pastrnak-e
# wavelength range: 0.22 - 5 um

## NOTE: ordinary (o) axes in X and Y, extraordinary (e) axis in Z

AlN_range = FreqRange(min = um_scale / 5.0, max = um_scale / 0.22)

AlN_frq1 = 1 / (0.1715 * um_scale)
AlN_gam1 = 0
AlN_sig1 = 1.3786
AlN_frq2 = 1 / (15.03 * um_scale)
AlN_gam2 = 0
AlN_sig2 = 3.861

AlN_susc_o = [
    LorentzianSusceptibility(frequency = AlN_frq1, gamma = AlN_gam1, sigma_diag = AlN_sig1 * [1, 1, 0]),
    LorentzianSusceptibility(frequency = AlN_frq2, gamma = AlN_gam2, sigma_diag = AlN_sig2 * [1, 1, 0])]

AlN_frq1 = 1 / (0.1746 * um_scale)
AlN_gam1 = 0
AlN_sig1 = 1.6173
AlN_frq2 = 1 / (15.03 * um_scale)
AlN_gam2 = 0
AlN_sig2 = 4.139

AlN_susc_e = [LorentzianSusceptibility(frequency = AlN_frq1, gamma = AlN_gam1, sigma_diag = AlN_sig1 * [0, 0, 1]),
    LorentzianSusceptibility(frequency = AlN_frq2, gamma = AlN_gam2, sigma_diag = AlN_sig2 * [0, 0, 1])]

AlN_aniso = Medium(epsilon_model = DispersionModel(inf_diag = [3.1399, 3.1399, 3.0729], susceptibilities = [AlN_susc_o; AlN_susc_e], valid_freq_range = AlN_range))

#------------------------------------------------------------------
# alumina/sapphire (Al2O3) from I.H. Malitson and M.J. Dodge, J. Optical Society of America, Vol. 62, pp. 1405 (1972)
# ref: https://refractiveindex.info/?shelf=main&book=Al2O3&page=Malitson-o
# ref: https://refractiveindex.info/?shelf=main&book=Al2O3&page=Malitson-e
# wavelength range: 0.2 - 5 um

## NOTE: ordinary (o) axes in X and Y, extraordinary (e) axis in Z

Al2O3_range = FreqRange(min = um_scale / 5.0, max = um_scale / 0.2)

Al2O3_frq1 = 1 / (0.0726631 * um_scale)
Al2O3_gam1 = 0
Al2O3_sig1 = 1.4313493
Al2O3_frq2 = 1 / (0.1193242 * um_scale)
Al2O3_gam2 = 0
Al2O3_sig2 = 0.65054713
Al2O3_frq3 = 1 / (18.02825 * um_scale)
Al2O3_gam3 = 0
Al2O3_sig3 = 5.3414021

Al2O3_susc_o = [LorentzianSusceptibility(frequency = Al2O3_frq1, gamma = Al2O3_gam1, sigma_diag = Al2O3_sig1 * [1, 1, 0]),
    LorentzianSusceptibility(frequency = Al2O3_frq2, gamma = Al2O3_gam2, sigma_diag = Al2O3_sig2 * [1, 1, 0]),
    LorentzianSusceptibility(frequency = Al2O3_frq3, gamma = Al2O3_gam3, sigma_diag = Al2O3_sig3 * [1, 1, 0])]

Al2O3_frq1 = 1 / (0.0740288 * um_scale)
Al2O3_gam1 = 0
Al2O3_sig1 = 1.5039759
Al2O3_frq2 = 1 / (0.1216529 * um_scale)
Al2O3_gam2 = 0
Al2O3_sig2 = 0.55069141
Al2O3_frq3 = 1 / (20.072248 * um_scale)
Al2O3_gam3 = 0
Al2O3_sig3 = 6.5927379

Al2O3_susc_e = [
    LorentzianSusceptibility(frequency = Al2O3_frq1, gamma = Al2O3_gam1, sigma_diag = Al2O3_sig1 * [0, 0, 1]),
    LorentzianSusceptibility(frequency = Al2O3_frq2, gamma = Al2O3_gam2, sigma_diag = Al2O3_sig2 * [0, 0, 1]),
    LorentzianSusceptibility(frequency = Al2O3_frq3, gamma = Al2O3_gam3, sigma_diag = Al2O3_sig3 * [0, 0, 1])
]

Al2O3_aniso = Medium(epsilon_model = DispersionModel(inf = 1, susceptibilities = [Al2O3_susc_o; Al2O3_susc_e], valid_freq_range = Al2O3_range))

#------------------------------------------------------------------
# yttrium oxide (Y2O3) from Y. Nigara, Japanese J. of Applied Physics, Vol. 7, pp. 404-8 (1968)
# ref: https://refractiveindex.info/?shelf=main&book=Y2O3&page=Nigara
# wavelength range: 0.25 - 9.6 um

Y2O3_range = FreqRange(min = um_scale / 9.6, max = um_scale / 0.25)

Y2O3_frq1 = 1 / (0.1387 * um_scale)
Y2O3_gam1 = 0
Y2O3_sig1 = 2.578
Y2O3_frq2 = 1 / (22.936 * um_scale)
Y2O3_gam2 = 0
Y2O3_sig2 = 3.935

Y2O3_susc = [LorentzianSusceptibility(frequency = Y2O3_frq1, gamma = Y2O3_gam1, sigma = Y2O3_sig1),
    LorentzianSusceptibility(frequency = Y2O3_frq2, gamma = Y2O3_gam2, sigma = Y2O3_sig2)]

Y2O3 = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = Y2O3_susc, valid_freq_range = Y2O3_range))

#------------------------------------------------------------------
# undoped yttrium aluminum garnet (YAG) from D.E. Zelmon et al., Applied Optics, Vol. 37, 4933-5 (1998)
# ref: https://refractiveindex.info/?shelf=main&book=Y3Al5O12&page=Zelmon
# wavelength range: 0.4 - 5.0 um

YAG_range = FreqRange(min = um_scale / 5.0, max = um_scale / 0.4)

YAG_frq1 = 1 / (0.1088577052853862 * um_scale)
YAG_gam1 = 0
YAG_sig1 = 2.28200
YAG_frq2 = 1 / (16.814695953242804 * um_scale)
YAG_gam2 = 0
YAG_sig2 = 3.27644

YAG_susc = [LorentzianSusceptibility(frequency = YAG_frq1, gamma = YAG_gam1, sigma = YAG_sig1),
    LorentzianSusceptibility(frequency = YAG_frq2, gamma = YAG_gam2, sigma = YAG_sig2)]

YAG = Medium(epsilon_model = DispersionModel(inf = 1.0, susceptibilities = YAG_susc, valid_freq_range = YAG_range))

#------------------------------------------------------------------
# cadmium telluride (CdTe) from D.T.F. Marple, J. Applied Physics, Vol. 35, pp. 539-42 (1964)
# ref: https://refractiveindex.info/?shelf=main&book=CdTe&page=Marple
# wavelength range: 0.86 - 2.5 um

CdTe_range = FreqRange(min = um_scale / 2.5, max = um_scale / 0.86)

CdTe_frq1 = 1 / (0.6049793384901669 * um_scale)
CdTe_gam1 = 0
CdTe_sig1 = 1.53

CdTe_susc = [LorentzianSusceptibility(frequency = CdTe_frq1, gamma = CdTe_gam1, sigma = CdTe_sig1)]

CdTe = Medium(epsilon_model = DispersionModel(inf = 5.68, susceptibilities = CdTe_susc, valid_freq_range = CdTe_range))


text:
# ------------------------------------------------------------------------------
# Make SymPy work better
# -------------------------------------------------------------------------------
import SymPy: simplify

# If a number is passed to simplify, then simplify return it. SymPy.simplify now only takes symbolic
# expressions.
# Will be removed when SymPy.jl update.
simplify(val) = val

ext: jl


text:
include("../coin-or-common.jl")

name = "CSDP"
version = v"6.2.0"

# Collection of sources required to build Clp
sources = [
    GitSource("https://github.com/coin-or/Csdp.git",
              "0dcf187a159c365b6d4e4e0ed5849f7b706da408"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/Csdp*

atomic_patch -p1 ${WORKSPACE}/srcdir/patches/blegat.patch

CFLAGS="-O2 -fPIC -fopenmp -ansi -Wall -DUSEOPENMP -DSETNUMTHREADS -DUSEGETTIME -I../include"
LIBS="-L../lib -lsdp -lopenblas -lm"

if [[ "${nbits}" == 64 ]] && [[ "${target}" != *aarch64* ]]; then
    CFLAGS="$CFLAGS -m64 -DBIT64"
fi

if [[ "${target}" == *-freebsd* ]] || [[ "${target}" == *-apple-* ]]; then
    CC=gcc
fi

make CFLAGS="$CFLAGS" LIBS="$LIBS" CC="$CC"
make install

cd lib

all_load="--whole-archive"
noall_load="--no-whole-archive"
if [[ "${target}" == *-apple-* ]]; then
  all_load="-all_load"
  noall_load="-noall_load"
fi

mkdir -p ${libdir}
${CC} -fopenmp -fPIC -shared -Wl,${all_load} libsdp.a -Wl,${noall_load} -o ${libdir}/libcsdp.${dlext} -lgomp -lopenblas -lm
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libcsdp", :libcsdp),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("OpenBLAS32_jll"),
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=gcc_version)

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

function ex_inv(x)

# ----------------------------------------------------------------------------------------------------------------------
# Custom inv function
# ----------------------------------------------------------------------------------------------------------------------

    xInv = x\eye(size(x)[1], size(x)[2]);

    return xInv;
end


text:
#!/usr/bin/env julia

using SEM
using LinearAlgebra,Plots, UnPack
#----------------------------------#
kx=1.0
ky=1.0
ux=1.00
uy=0.0

function utrue(x,y,t)
    xx = @. x - ux*t
    yy = @. y - uy*t
    ut = @. sin(kx*pi*xx)*sin.(ky*pi*yy)
    return ut
end

function set0!(u,x,y,t)
    u .= utrue(x,y,t)
    return
end

function set∂!(ub,x,y,t)
    ub .= 0.0
    return
end

function setF!(f,x,y,t)
    f .= 0.0
    return
end

function setν!(ν,x,y,t)
    ν .= 0.0
    return
end

function callback!(cdn::ConvectionDiffusion)
    @unpack fld,mshV = cdn
    @unpack time, istep = cdn.tstep

    ut = utrue(mshV.x,mshV.y,time[1])
    u  = fld.u
    er = norm(ut-u,Inf)
    println("$(cdn.name), Step $(istep[1]), Time=$(time[1]), er=$er")
    plt = meshplt(u,mshV)
    plt = plot!(zlims=(-1,1))
    display(plt)
    #frame(anim)
    return
end

#----------------------------------#
Ex = 5; nr1 = 8; nrd = 8;
Ey = 5; ns1 = 8; nsd = 8;

ifperiodic=[true,false]
m1 = Mesh(nr1,ns1,Ex,Ey,ifperiodic)
md = Mesh(nrd,nsd,Ex,Ey,ifperiodic)
bc = ['N','N','D','D']

fld = Field(bc,m1)
vx = @. 0*m1.x + ux
vy = @. 0*m1.x + uy
tstep = TimeStepper(0.,1.,5e-3)
ps = ConvectionDiffusion("ps",fld,vx,vy,tstep,md,set0!,set∂!,setF!,setν!)

#anim = Animation()
simulate!(ps,callback!)

#gif(anim,"ConvectionDiffusion.gif",fps=30)
#----------------------------------#
nothing


text:
## CONFUSION MATRIX OBJECT

"""
    ConfusionMatrixObject{C}

Confusion matrix with `C ≥ 2` classes. Rows correspond to predicted values
and columns to the ground truth.
"""
struct ConfusionMatrixObject{C}
    mat::Matrix
    labels::Vector{String}
end

"""
    ConfusionMatrixObject(m, labels)

Instantiates a confusion matrix out of a square integer matrix `m`.
Rows are the predicted class, columns the ground truth. See also the
[wikipedia article](https://en.wikipedia.org/wiki/Confusion_matrix).

"""
function ConfusionMatrixObject(m::Matrix{Int}, labels::Vector{String})
    s = size(m)
    s[1] == s[2] || throw(ArgumentError("Expected a square matrix."))
    s[1] > 1 || throw(ArgumentError("Expected a matrix of size ≥ 2x2."))
    length(labels) == s[1] ||
        throw(ArgumentError("As many labels as classes must be provided."))
    ConfusionMatrixObject{s[1]}(m, labels)
end

# allow to access cm[i,j] but not set (it's immutable)
Base.getindex(cm::ConfusionMatrixObject, inds...) = getindex(cm.mat, inds...)

"""
    _confmat(ŷ, y; rev=false)

A private method. General users should use `confmat` or other instance
of the measure type [`ConfusionMatrix`](@ref).

Computes the confusion matrix given a predicted `ŷ` with categorical elements
and the actual `y`. Rows are the predicted class, columns the ground truth.
The ordering follows that of `levels(y)`.

## Keywords

* `rev=false`: in the binary case, this keyword allows to swap the ordering of
               classes.
* `perm=[]`:   in the general case, this keyword allows to specify a permutation
               re-ordering the classes.
* `warn=true`: whether to show a warning in case `y` does not have scientific
               type `OrderedFactor{2}` (see note below).

## Note

To decrease the risk of unexpected errors, if `y` does not have
scientific type `OrderedFactor{2}` (and so does not have a "natural
ordering" negative-positive), a warning is shown indicating the
current order unless the user explicitly specifies either `rev` or
`perm` in which case it's assumed the user is aware of the class
ordering.

The `confusion_matrix` is a measure (although neither a score nor a
loss) and so may be specified as such in calls to `evaluate`,
`evaluate!`, although not in `TunedModel`s.  In this case, however,
there no way to specify an ordering different from `levels(y)`, where
`y` is the target.

"""
function _confmat(ŷ::Vec{<:CategoricalValue}, y::Vec{<:CategoricalValue};
                          rev::Union{Nothing,Bool}=nothing,
                          perm::Union{Nothing,Vector{<:Integer}}=nothing,
                          warn::Bool=true)

    check_dimensions(ŷ, y)
    levels_ = levels(y)
    nc = length(levels_)
    if rev !== nothing && rev && nc > 2
        throw(ArgumentError("Keyword `rev` can only be used in binary case."))
    end
    if perm !== nothing && !isempty(perm)
        length(perm) == nc ||
            throw(ArgumentError("`perm` must be of length matching the "*
                                "number of classes."))
        Set(perm) == Set(collect(1:nc)) ||
            throw(ArgumentError("`perm` must specify a valid permutation of "*
                                "`[1, 2, ..., c]`, where `c` is "*
                                "number of classes."))
    end

    # warning
    if rev === nothing && perm === nothing
        if warn &&
            if nc==2 && !(scitype_union(y) >: OrderedFactor{2})
                @warn "The classes are un-ordered,\n" *
                      "using: negative='$(levels_[1])' and positive='$(levels_[2])'.\n" *
                      "To suppress this warning, consider coercing to OrderedFactor."
            elseif !(scitype_union(y) >: OrderedFactor{nc})
                @warn "The classes are un-ordered,\n" *
                      "using order: $([l for l in levels_]).\n" *
                      "To suppress this warning, consider coercing to OrderedFactor."
            end
        end
        rev  = false
        perm = Int[]
    elseif rev !== nothing && nc == 2
        # rev takes precedence in binary case
        if rev
            perm = [2, 1]
        else
            perm = Int[]
        end
    end

    # No permutation
    if isempty(perm)
        cmat = zeros(Int, nc, nc)
        @inbounds for i in eachindex(y)
            cmat[int(ŷ[i]), int(y[i])] += 1
        end
        return ConfusionMatrixObject(cmat, string.(levels_))
    end

    # With permutation
    cmat = zeros(Int, nc, nc)
    iperm = invperm(perm)
    @inbounds for i in eachindex(y)
        cmat[iperm[int(ŷ[i])], iperm[int(y[i])]] += 1
    end
    return ConfusionMatrixObject(cmat, string.(levels_[perm]))
end

# Machinery to display the confusion matrix in a non-confusing way
# (provided the REPL is wide enough)

splitw(w::Int) = (sp1 = div(w, 2); sp2 = w - sp1; (sp1, sp2))

function Base.show(stream::IO, m::MIME"text/plain", cm::ConfusionMatrixObject{C}
                   ) where C
    width    = displaysize(stream)[2]
    cw       = 13
    textlim  = 9
    totalwidth = cw * (C+1) + C + 2
    width < totalwidth && (show(stream, m, cm.mat); return)

    iob     = IOBuffer()
    wline   = s -> write(iob, s * "\n")
    splitcw = s -> (w = cw - length(s); splitw(w))
    cropw   = s -> length(s) > textlim ? s[1:prevind(s, textlim)] * "…" : s

    # 1.a top box
    " "^(cw+1) * "┌" * "─"^((cw + 1) * C - 1) * "┐" |> wline
    gt = "Ground Truth"
    w  = (cw + 1) * C - 1 - length(gt)
    sp1, sp2 = splitw(w)
    " "^(cw+1) * "│" * " "^sp1 * gt * " "^sp2 * "│" |> wline
    # 1.b separator
    "┌" * "─"^cw * "┼" * ("─"^cw * "┬")^(C-1) * "─"^cw * "┤" |> wline
    # 2.a description line
    pr = "Predicted"
    sp1, sp2 = splitcw(pr)
    partial = "│" * " "^sp1 * pr * " "^sp2 * "│"
    for c in 1:C
        # max = 10
        s = cm.labels[c] |> cropw
        sp1, sp2 = splitcw(s)
        partial *= " "^sp1 * s * " "^sp2 * "│"
    end
    partial |> wline
    # 2.b separating line
    "├" * "─"^cw * "┼" * ("─"^cw * "┼")^(C-1) * ("─"^cw * "┤") |> wline
    # 2.c line by line
    for c in 1:C
        # line
        s  = cm.labels[c] |> cropw
        sp1, sp2 = splitcw(s)
        partial = "│" * " "^sp1 * s * " "^sp2 * "│"
        for r in 1:C
            e = string(cm[c, r])
            sp1, sp2 = splitcw(e)
            partial *= " "^sp1 * e * " "^sp2 * "│"
        end
        partial |> wline
        # separator
        if c < C
            "├" * "─"^cw * "┼" * ("─"^cw * "┼")^(C-1) * ("─"^cw * "┤") |> wline
        end
    end
    # 2.d final line
    "└" * "─"^cw * "┴" * ("─"^cw * "┴")^(C-1) * ("─"^cw * "┘") |> wline
    write(stream, take!(iob))
end

## CONFUSION MATRIX AS MEASURE

struct ConfusionMatrix <: Measure
    perm::Union{Nothing,Vector{<:Integer}}
end

ConfusionMatrix(; perm=nothing) = ConfusionMatrix(perm)

is_measure(::ConfusionMatrix) = true
is_measure_type(::Type{ConfusionMatrix}) = true
human_name(::Type{<:ConfusionMatrix}) = "confusion matrix"
target_scitype(::Type{ConfusionMatrix}) =
    AbstractVector{<:Finite}
supports_weights(::Type{ConfusionMatrix}) = false
prediction_type(::Type{ConfusionMatrix}) = :deterministic
instances(::Type{<:ConfusionMatrix}) = ["confusion_matrix", "confmat"]
orientation(::Type{ConfusionMatrix}) = :other
reports_each_observation(::Type{ConfusionMatrix}) = false
is_feature_dependent(::Type{ConfusionMatrix}) = false
aggregation(::Type{ConfusionMatrix}) = Sum()

@create_aliases ConfusionMatrix

@create_docs(ConfusionMatrix,
body=
"""
If `r` is the return value, then the raw confusion matrix is `r.mat`,
whose rows correspond to predictions, and columns to ground truth.
The ordering follows that of `levels(y)`.

Use `ConfusionMatrix(perm=[2, 1])` to reverse the class order for binary
data. For more than two classes, specify an appropriate permutation, as in
`ConfusionMatrix(perm=[2, 3, 1])`.
""",
scitype=DOC_ORDERED_FACTOR_BINARY)

# calling behaviour:
(m::ConfusionMatrix)(ŷ::Vec{<:CategoricalValue}, y::Vec{<:CategoricalValue}) =
    _confmat(ŷ, y, perm=m.perm)

# overloading addition to make aggregation work:
Base.round(m::MLJBase.ConfusionMatrixObject; kws...) = m
function Base.:+(m1::ConfusionMatrixObject, m2::ConfusionMatrixObject)
    if m1.labels != m2.labels
        throw(ArgumentError("Confusion matrix labels must agree"))
    end
    ConfusionMatrixObject(m1.mat + m2.mat, m1.labels)
end


text:
"""
BagChain(layers...)
BagChain multiple layers / functions together, so that they are called in sequence
on a given input supported by bags.

"""
struct BagChain{T<:Tuple}
    layers::T
    BagChain(xs...) = new{typeof(xs)}(xs)
end

Flux.@functor BagChain

Flux.@forward BagChain.layers Base.getindex, Base.first, Base.last, Base.lastindex
Flux.@forward BagChain.layers Base.iterate

applychain(::Tuple{}, x, bags) = x
applychain(fs::Tuple, x, bags) = applychain(Base.tail(fs), first(fs)(x, bags), bags)

(c::BagChain)(x, bags) = applychain(c.layers, x, bags)
(c::BagChain)(x::ArrayNode, bags) = ArrayNode(applychain(c.layers, x.data, bags))
(c::BagChain)(x::BagNode) = applychain(c.layers, x.data, x.bags)

Base.getindex(c::BagChain, i::AbstractArray) = BagChain(c.layers[i]...)

activations(c::BagChain, x, bags) = accumulate((x, m) -> m(x, bags), c.layers, init = x)


text:
using StanModels, CSV

df = CSV.read(joinpath(@__DIR__, "..", "..", "data", "UCBadmit.csv"), delim=';');

# Preprocess data

df[!, :male] = [(df[!, :gender][i] == "male" ? 1 : 0) for i in 1:size(df, 1)];
df[!, :dept_id] = [Int(df[!, :dept][i][1])-64 for i in 1:size(df, 1)];

m13_2s = "
  data{
      int N;
      int N_depts;
      int applications[N];
      int admit[N];
      real male[N];
      int dept_id[N];
  }
  parameters{
      vector[N_depts] a_dept;
      real a;
      real bm;
      real<lower=0> sigma_dept;
  }
  model{
      vector[N] p;
      sigma_dept ~ cauchy( 0 , 2 );
      bm ~ normal( 0 , 1 );
      a ~ normal( 0 , 10 );
      a_dept ~ normal( a , sigma_dept );
      for ( i in 1:N ) {
          p[i] = a_dept[dept_id[i]] + bm * male[i];
          p[i] = inv_logit(p[i]);
      }
      admit ~ binomial( applications , p );
  }
";

# Define the Stanmodel and set the output format to :mcmcchains.

sm = SampleModel("m13_2s", m13_2s);

# Input data for cmdstan

m13_2_data = Dict(
  "N" => size(df, 1),
  "N_depts" => maximum(df[!, :dept_id]),
  "applications" => df[!, :applications],
  "admit" => df[!, :admit],
  "male" => df[!, :male],
  "dept_id"=> df[!, :dept_id]
);

# Sample using cmdstan

(sampleFile, log_file) = stan_sample(sm, data=m13_2_data);

# Describe the draws

if !(sample_file == nothing)
  chn = read_samples(sm)
  describe(chn)
end


text:
using VOTables
using Test

@testset "VOTables.jl" begin
    # Write your tests here.
end


text:
# # Classification of penguin species
#
# ## Packages

using AlgebraOfGraphics
using CairoMakie
using CalibrationErrors
using DataFrames
using Distributions
using MLJ
using MLJNaiveBayesInterface
using PalmerPenguins

using Random

## Plotting settings
set_aog_theme!()
CairoMakie.activate!(; type="svg")

# ## Data
#
# In this example we study the calibration of different models that classify three penguin
# species based on measurements of their bill and flipper lengths.
#
# We use the [Palmer penguins dataset](https://allisonhorst.github.io/palmerpenguins/) to
# to train and validate the models.

penguins = dropmissing(DataFrame(PalmerPenguins.load()))

penguins_mapping =
    data(penguins) * mapping(
        :bill_length_mm => "bill length (mm)", :flipper_length_mm => "flipper length (mm)"
    )
draw(penguins_mapping * mapping(; color=:species) * visual(; alpha=0.7))

# We split the data randomly into a training and validation dataset. The training dataset
# contains around 60% of the samples.

Random.seed!(1234)
n = nrow(penguins)
k = floor(Int, 0.7 * n)
Random.seed!(100)
penguins.train = shuffle!(vcat(trues(k), falses(n - k)))

## Plot the training and validation data
dataset = :train => renamer(true => "training", false => "validation") => "Dataset"
plt = penguins_mapping * mapping(; color=:species, col=dataset) * visual(; alpha=0.7)
draw(plt; axis=(height=300,))

# ## Fitting normal distributions
#
# For each species, we fit independent normal distributions to the observations of the bill
# and flipper length in the training data, using maximum likelihood estimation.

y, X = unpack(
    penguins,
    ==(:species),
    x -> x === :bill_length_mm || x === :flipper_length_mm;
    :species => Multiclass,
    :bill_length_mm => MLJ.Continuous,
    :flipper_length_mm => MLJ.Continuous,
)
model = fit!(machine(GaussianNBClassifier(), X, y); rows=penguins.train);

# We plot the estimated normal distributions.

## plot datasets
fg = draw(plt; axis=(height=300,))

## plot Gaussian distributions
xgrid = range(extrema(penguins.bill_length_mm)...; length=100)
ygrid = range(extrema(penguins.flipper_length_mm)...; length=100)
let f = (x, y, dist) -> pdf(dist, [x, y])
    for (class, color) in zip(classes(y), Makie.wong_colors())
        pdfs = f.(xgrid, ygrid', Ref(model.fitresult.gaussians[class]))
        contour!(fg.figure[1, 1], xgrid, ygrid, pdfs; color=color)
        contour!(fg.figure[1, 2], xgrid, ygrid, pdfs; color=color)
    end
end

fg

# ## Naive Bayes classifier
#
# Let us assume that the bill and flipper length are conditionally independent given the
# penguin species. Then Bayes' theorem implies that
# ```math
# \begin{aligned}
# \mathbb{P}(\mathrm{species} \,|\, \mathrm{bill}, \mathrm{flipper})
# &= \frac{\mathbb{P}(\mathrm{species}) \mathbb{P}(\mathrm{bill}, \mathrm{flipper} \,|\, \mathrm{species})}{\mathbb{P}(\mathrm{bill}, \mathrm{flipper})} \\
# &= \frac{\mathbb{P}(\mathrm{species}) \mathbb{P}(\mathrm{bill} \,|\, \mathrm{species}) \mathbb{P}(\mathrm{flipper} \,|\, \mathrm{species})}{\mathbb{P}(\mathrm{bill}, \mathrm{flipper})}.
# \end{aligned}
# ```
# This predictive model is known as
# [naive Bayes classifier](https://en.wikipedia.org/wiki/Naive_Bayes_classifier).
#
# In the section above, we estimated $\mathbb{P}(\mathrm{species})$,
# $\mathbb{P}(\mathrm{bill} \,|\, \mathrm{species})$, and
# $\mathbb{P}(\mathrm{flipper} \,|\, \mathrm{species})$ for each penguin species from
# the training data. For the conditional distributions we used a Gaussian approximation.

predictions = MLJ.predict(model)
train_predict = predictions[penguins.train]
val_predict = predictions[.!penguins.train]

## Plot datasets
fg = draw(plt; axis=(height=300,))

## Plot predictions
predictions_grid = reshape(
    MLJ.predict(model, reduce(hcat, vcat.(xgrid, ygrid'))'), length(xgrid), length(ygrid)
)
for (class, color) in zip(classes(y), Makie.wong_colors())
    p = pdf.(predictions_grid, class)
    contour!(fg.figure[1, 1], xgrid, ygrid, p; color=color)
    contour!(fg.figure[1, 2], xgrid, ygrid, p; color=color)
end

fg

# ## Evaluation
#
# We evaluate the probabilistic predictions of the naive Bayes classifier that we just
# trained.
#
# ### Log-likelihood
#
# We compute the average log-likelihood of the validation data. It is equivalent to the
# negative cross-entropy.

val_y = y[.!penguins.train]
-mean(cross_entropy(val_predict, val_y))

# ### Brier score
#
# The average log-likelihood is also equivalent to the
# [logarithmic score](https://sites.stat.washington.edu/raftery/Research/PDF/Gneiting2007jasa.pdf).
# The Brier score is another strictly proper scoring rule that can be used for evaluating
# probabilistic predictions.

mean(brier_score(val_predict, val_y))

# ### Expected calibration error
#
# As all proper scoring rules, the logarithmic and the Brier score can be [decomposed in
# three terms that quantify the sharpness and calibration of the predictive model and the
# irreducible uncertainty of the targets that is inherent to the prediction
# problem](https://doi.org/10.1002/qj.456). The calibration term in this decomposition is
# the expected calibration error (ECE)
# ```math
# \mathbb{E} d\big(P_X, \mathrm{law}(Y \,|\, P_X)\big)
# ```
# with respect to the score divergence $d$.
#
# Scoring rules, however, include also the sharpness and the uncertainty term. Thus models
# can trade off calibration for sharpness and therefore scoring rules are not suitable for
# specifically evaluating calibration of predictive models.
#
# The score divergence to the logarithmic and the Brier score are the Kullback-Leibler (KL)
# divergence
# ```math
# d\big(P_X, \mathrm{law}(Y \,|\, P_X)\big) = \sum_{y} \mathbb{P}(Y = y \,|\, P_X)
# \log\big(\mathbb{P}(Y = y \,|\, P_X) / P_X(\{y\})\big)
# ```
# and the squared Euclidean distance
# ```math
# d\big(P_X, \mathrm{law}(Y \,|\, P_X)\big) = \sum_{y} \big(P_X - \mathrm{law}(Y \,|\, P_X)\big)^2(\{y\}),
# ```
# respectively. The KL divergence is defined only if $\mathrm{law}(Y \,|\, P_X)$ is
# absolutely continuous with respect to $P_X$, i.e., if $P_X(\{y\}) = 0$ implies
# $\mathbb{P}(Y = y \,|\, P_X) = 0$.

# We estimate the ECE by binning the probability simplex of predictions $P_X$ and computing
# the weighted average of the distances between the mean prediction and the distribution of
# targets in each bin.
#
# One approach is to use bins of uniform size.

ece = ECE(UniformBinning(10), (μ, y) -> kl_divergence(y, μ));

# We have to work with a numerical encoding of the true penguin species and a
# corresponding vector of predictions. We use [`RowVecs`](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/api/#KernelFunctions.RowVecs)
# to indicate that the rows in the matrix of probabilities returned by `pdf`
# are the predictions. If we would provide predictions as columns of a matrix, we would have
# to use [`ColVecs`](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/api/#KernelFunctions.ColVecs).

val_yint = map(MLJ.levelcode, val_y)
val_probs = RowVecs(pdf(val_predict, MLJ.classes(y)));

# We compute the estimate on the validation data:

ece(val_probs, val_yint)

# For the squared Euclidean distance we obtain:

ece = ECE(UniformBinning(10), SqEuclidean())
ece(val_probs, val_yint)

# Alternatively, one can use a data-dependent binning scheme that tries to split the
# predictions in a way that minimizes the variance in each bin.
#
# With the KL divergence we get:

ece = ECE(MedianVarianceBinning(5), (μ, y) -> kl_divergence(y, μ))
ece(val_probs, val_yint)

# For the squared Euclidean distance we obtain:

ece = ECE(MedianVarianceBinning(5), SqEuclidean())
ece(val_probs, val_yint)

# We see that the estimates (of the same theoretical quantity!) are highly dependent on the
# chosen binning scheme.

# ### Kernel calibration error
#
# As an alternative to the ECE, we estimate the kernel calibration error (KCE). We keep it
# simple here, and use the tensor product kernel
# ```math
# k\big((\mu, y), (\mu', y')\big) = \delta_{y,y'} \exp{\bigg(-\frac{{\|\mu - \mu'\|}_2^2}{2\nu^2} \bigg)}
# ```
# with length scale $\nu > 0$ for predictions $\mu,\mu'$ and corresponding targets $y, y'$.
# For simplicity, we estimate length scale $\nu$ with the median heuristic.

distances = pairwise(SqEuclidean(), RowVecs(pdf(train_predict, MLJ.classes(y))))
ν = sqrt(median(distances[i] for i in CartesianIndices(distances) if i[1] < i[2]))
kernel = with_lengthscale(GaussianKernel(), ν) ⊗ WhiteKernel();

# We obtain the following biased estimate of the squared KCE (SKCE):

skce = BiasedSKCE(kernel)
skce(val_probs, val_yint)

# Similar to the biased estimates of the ECE, the biased estimates of the SKCE are always
# non-negative. The unbiased estimates can be negative as well, in particular if the model
# is (close to being) calibrated:

skce = UnbiasedSKCE(kernel)
skce(val_probs, val_yint)

# When the datasets are large, the quadratic sample complexity of the standard biased and
# unbiased estimators of the SKCE can become prohibitive. In these cases, one can resort to
# an estimator that averages estimates of non-overlapping blocks of samples. This estimator
# allows to trade off computational cost for increased variance.
#
# Here we consider the extreme case of blocks with two samples, which yields an estimator
# with linear sample complexity:

skce = BlockUnbiasedSKCE(kernel, 2)
skce(val_probs, val_yint)


text:
isfile("result.md") && Base.rm("result.md")
# isfile("tune.json") && Base.rm("tune.json")
isfile("judgement.md") && Base.rm("judgement.md")

for f in readdir(@__DIR__)
    startswith(f, "result") || continue
    endswith(f, "json") || continue
    Base.rm(f)
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
function ChainRulesCore.rrule(::typeof(projection_on_set), d::DefaultDistance, v::T, s::MOI.EqualTo) where {T}
    vproj = projection_on_set(d, v, s)
    function pullback(Δvproj)
        return (ChainRulesCore.NO_FIELDS, ChainRulesCore.DoesNotExist(), ChainRulesCore.Zero(), ChainRulesCore.Composite{typeof(s)}(value=Δvproj))
    end
    return (vproj, pullback)
end

function ChainRulesCore.rrule(::typeof(projection_on_set), d::DefaultDistance, v::T, s::MOI.LessThan) where {T}
    vproj = projection_on_set(d, v, s)
    function pullback(Δvproj)
        if vproj == v # if exactly equal, then constraint inactive
            return (ChainRulesCore.NO_FIELDS, ChainRulesCore.DoesNotExist(), Δvproj, ChainRulesCore.Composite{typeof(s)}(upper=zero(Δvproj)))
        end
        return (ChainRulesCore.NO_FIELDS, ChainRulesCore.DoesNotExist(), zero(Δvproj), ChainRulesCore.Composite{typeof(s)}(upper=Δvproj))
    end
    return (vproj, pullback)
end

function ChainRulesCore.rrule(::typeof(projection_on_set), d::DefaultDistance, v::T, s::MOI.GreaterThan) where {T}
    vproj = projection_on_set(d, v, s)
    function pullback(Δvproj)
        if vproj == v # if exactly equal, then constraint inactive
            return (ChainRulesCore.NO_FIELDS, ChainRulesCore.DoesNotExist(), Δvproj, ChainRulesCore.Composite{typeof(s)}(lower=zero(Δvproj)))
        end
        return (ChainRulesCore.NO_FIELDS, ChainRulesCore.DoesNotExist(), zero(Δvproj), ChainRulesCore.Composite{typeof(s)}(lower=Δvproj))
    end
    return (vproj, pullback)
end

function ChainRulesCore.rrule(::typeof(projection_on_set), d::DefaultDistance, v::AbstractVector{T}, s::MOI.Reals) where {T}
    vproj = projection_on_set(d, v, s)
    function pullback(Δvproj)
        return (ChainRulesCore.NO_FIELDS, ChainRulesCore.DoesNotExist(), Δvproj, ChainRulesCore.DoesNotExist())
    end
    return (vproj, pullback)
end

function ChainRulesCore.rrule(::typeof(projection_on_set), d::DefaultDistance, v::AbstractVector{T}, s::MOI.Zeros) where {T}
    vproj = projection_on_set(d, v, s)
    function pullback(Δvproj)
        return (ChainRulesCore.NO_FIELDS, ChainRulesCore.DoesNotExist(), FillArrays.Zeros(length(v)), ChainRulesCore.DoesNotExist())
    end
    return (vproj, pullback)
end

function ChainRulesCore.frule((_, _, Δv, _), ::typeof(projection_on_set), d::DefaultDistance, v::AbstractVector{T}, s::MOI.Nonnegatives) where {T}
    vproj = projection_on_set(d, v, s)
    ∂vproj = Δv .* (v .>= 0)
    return vproj, ∂vproj
end

function ChainRulesCore.frule((_, _, Δv, _), ::typeof(projection_on_set), d::DefaultDistance, v::AbstractVector{T}, s::MOI.Nonpositives) where {T}
    vproj = projection_on_set(d, v, s)
    ∂vproj = Δv .* (v .<= 0)
    return vproj, ∂vproj
end

function ChainRulesCore.rrule(::typeof(projection_on_set), d::DefaultDistance, v::AbstractVector{T}, s::S) where {T,S <: Union{MOI.Nonnegatives,MOI.Nonpositives}}
    vproj = projection_on_set(d, v, s)
    function pullback(Δvproj)
        v̄ = zeros(eltype(Δvproj), length(Δvproj))
        for i in eachindex(Δvproj)
            if vproj[i] == v[i]
                v̄[i] = Δvproj[i]
            end
        end
        return (ChainRulesCore.NO_FIELDS, ChainRulesCore.DoesNotExist(), v̄, ChainRulesCore.DoesNotExist())
    end
    return (vproj, pullback)
end

function ChainRulesCore.rrule(::typeof(projection_on_set), d::Union{DefaultDistance,NormedEpigraphDistance}, v::AbstractVector{T}, s::MOI.SecondOrderCone) where {T}
    vproj = projection_on_set(d, v, s)
    t = v[1]
    x = v[2:end]
    norm_x = LinearAlgebra.norm2(x)
    function projection_on_set_pullback(Δv)
        Δt = Δv[1]
        Δx = Δv[2:end]
        v̄ = zeros(eltype(Δv), length(Δv))
        result_tuple = (ChainRulesCore.NO_FIELDS, ChainRulesCore.DoesNotExist(), v̄, ChainRulesCore.DoesNotExist())
        if norm_x ≤ t
            v̄ .= Δv
            return result_tuple
        elseif norm_x ≤ -t
            return result_tuple
        end
        inv_norm = inv(2norm_x)
        v̄[1] = inv_norm * sum(Δx[i] * x[i] for i in eachindex(x)) + Δt / 2
        dot_prod = LinearAlgebra.dot(x, Δx)
        inv_sq = inv(norm_x^2)
        cons_mul_last = inv_norm * inv_sq * t * x ⋅ Δx
        for i in eachindex(x)
            v̄[i + 1] = inv_norm * Δt * x[i]
            v̄[i + 1] += inv_norm * (t + norm_x) * Δx[i]
            v̄[i + 1] -= cons_mul_last * x[i]
        end
        return result_tuple
    end
    return (vproj, projection_on_set_pullback)
end

function ChainRulesCore.frule((_, _, Δv, _), ::typeof(projection_on_set), d::DefaultDistance, v::AbstractVector{T}, ::MOI.PositiveSemidefiniteConeTriangle) where {T}
    dim = isqrt(2 * length(v))
    X = unvec_symm(v, dim)
    (λ, U) = LinearAlgebra.eigen(X)
    λmin, λmax = extrema(λ)
    if λmin >= 0
        return (v, Δv)
    end
    if λmax <= 0
        # v zero vector
        return (0 * v, 0 * Δv)
    end
    λp = max.(0, λ)
    vproj = vec_symm(U * Diagonal(λp) * U')
    k = count(λi < 1e-4 for λi in λ)
    # TODO avoid full matrix materialize
    B = zeros(eltype(λ), dim, dim)
    for i in 1:dim, j in 1:dim
        if i > k && j > k
            B[i,j] = 1
        elseif i > k
            B[i,j] = λp[i] / (λp[i] - min(λ[j], 0))
        elseif j > k
            B[i,j] = λp[j] / (λp[j] - min(λ[i], 0))
        end
    end
    M = U * (B .* (U' * unvec_symm(Δv, dim) * U)) * U'
    Δvproj = vec_symm(M)
    return (vproj, Δvproj)
end

function ChainRulesCore.frule((_, _, Δv, _), ::typeof(projection_on_set), d::DefaultDistance, v::AbstractVector{T}, s::MOI.ExponentialCone) where {T}
    if _in_exp_cone(v; dual=false)
        return (v, Δv)
    end
    if _in_exp_cone(-v; dual=true)
        # if in polar cone Ko = -K*
        return (0v, 0Δv)
    end
    if v[1] <= 0 && v[2] <= 0
        vproj = [v[1], 0, max(v[3], 0)]
        Δvproj = [Δv[1], 0, (v[3] >= 0) * Δv[3]]
        return vproj, Δvproj
    end
    vproj = _exp_cone_proj_case_4(v)
    nu = vproj[3] - v[3]
    rs = vproj[1] / vproj[2]
    exp_rs = exp(rs)
    (z1, z2, z3) = vproj
    mat = [
        1 + nu * exp_rs / z2     -nu * exp_rs * rs / z2       0     exp_rs;
        -nu * exp_rs * rs / z2   1 + nu * exp_rs * rs^2 / z2    0     (1 - rs) * exp_rs;
        0                  0                      1     -1
        exp_rs             (1 - rs) * exp_rs          -1    0
    ]
    lin_sol = mat \ [Δv; 0]
    return (vproj, lin_sol[1:3])
end


text:
using Plots

function f(z, C)
    return z * z + C
end

function convergence_judgment(C, inital_z)
    z = inital_z
    i::UInt8 = 0
    while i < 100
        z = f(z, C)
        if abs2(z) > 4
            return i
        end
        i += 1
    end
    return i
end

function generate_grid(N, inital_z)
    grid = Array{UInt8}(undef, N, N)
    points = range(-2, 2, length = N)
    n_pts = length(points)
    Threads.@threads for i in 1:n_pts
        for j in 1:n_pts
            ReC = points[i]
            ImC = points[j]
            C = ComplexF64(ReC, ImC)
            grid[j,i] = convergence_judgment(C, inital_z)
        end
    end
    return grid
end

function main(N)
    inital_z = 0.0
    @time grid = generate_grid(N, inital_z)
    heatmap(1:N, 1:N, grid, aspect_ratio = 1)
end

main(5000)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module StatsModelsTools

using StatsBase
using StatsModels

include("rms_splines.jl")
include("polynomials.jl")

using .RMSSplines
using .Polynomials

export rcs, rms_generate_knots, rms_quantiles_for_knots
export poly

end # module


text:
#!/usr/bin/julia
# http://julia.readthedocs.io/en/latest/manual/calling-c-and-fortran-code/

println(now())

# --[ ccall ]--
# you need to compile your library with -shared and -fPIC options.

# check libkernel.so first
if stat("./libkernel.so").inode == 0 then
   println(" => ./libkernel.so not found")
else
   println(" => ./libkernel.so found")
end

# handle the LD_LIBRARY_PATH environment
function getenv(name::AbstractString) # getenv(3)
   val = ccall( (:getenv, "libc.so.6"), Cstring, (Cstring,), name )
   if val == C_NULL then
      error("getenv: undefined variable: ", name)
   end
   bytestring(val) # use unsafe_string for julia 5.0
end
function setenv(name::AbstractString, value::AbstractString, overwrite::Int64) # setenv(3)
   val = ccall( (:setenv, "libc.so.6"), Int32, (Cstring, Cstring, Cint),
                name, value, overwrite)
   if val == -1 then
      error("setenv: failed")
   end
   val
end

println(" => SHELL ", getenv("SHELL"))
setenv("LD_LIBRARY_PATH", "./", 1)
println(" => LD_LIBRARY_PATH ", getenv("LD_LIBRARY_PATH"))

a = ccall( (:hello, "libkernel.so"), Int32, () )
println(a)

b = [ 1.0, 2.0, 3.0, 4.0, 5.0, -5.0 ]
println(sizeof(b), " ", length(b))
ccall( (:dvdump, "libkernel.so"), Void, (Ptr{Cdouble}, Csize_t), b, length(b))
a = ccall( (:dasum, "libkernel.so"), Float64, (Ptr{Cdouble}, Csize_t), b, length(b))
println(a)
a = ccall( (:pdasum, "libkernel.so"), Float64, (Ptr{Cdouble}, Csize_t), b, length(b))

b = zeros(1, 5)
ccall( (:boxmuller, "libboxmuller.so"), Ptr{Float64}, (Ptr{Cdouble}, Csize_t), b, length(b))
println(b)
println(length(b))
#for i = 1:length(b)
#   println(b[i])
#end

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
if !@isdefined(BodePlots)
include("BodePlots.jl")
end

module PhaseLock

import Graphics.BoundingBox

using InspectDR
import InspectDR: Plot2D
import Printf: @sprintf
using NumericIO
using Colors

#Colour constants
#-------------------------------------------------------------------------------
const black = RGB24(0, 0, 0); const white = RGB24(1, 1, 1)
const red = RGB24(1, 0, 0); const green = RGB24(0, 1, 0); const blue = RGB24(0, 0, 1)

#Convenience aliases
#-------------------------------------------------------------------------------
mag(x) = abs(x); phase(x) = rad2deg(angle(x))
SI(x) = formatted(x, :SI, ndigits=3)
const j = im

#==Main types
===============================================================================#

#1st order transfer function:
mutable struct Xf1
	K::Float64
	ωz::Float64
	ωp::Float64
end
Xf1(;K=1, ωz=1, ωp=1) = Xf1(K, ωz, ωp)

#2nd order system (standard form):
mutable struct Sys2
	α::Float64
	ζ::Float64
	ωn::Float64
end
Sys2(;α=1, ζ=1, ωn=1) = Sys2(α, ζ, ωn)

function Sys2(G::Xf1)
	ωn = sqrt(G.K*G.ωp)
	ζ = 0.5*(G.ωp/ωn+ωn/G.ωz)
	α = ωn/(2*ζ*G.K)
	return Sys2(α=α, ζ=ζ, ωn=ωn)
end

#==Math: 2nd order transfer functions
===============================================================================#

#Open-loop transfer function, G(s) (including freq->phase integrator):
xf_ol(ω::Vector, G::Xf1) = (s=j*ω;
	return G.K .* (1 .+ s ./ G.ωz)./(1 .+ s ./ G.ωp) ./ s)

#Compute unity gain frequency of G(s) (including freq->phase integrator):
function ω_0(G::Xf1)
	invsq(x) = (result = 1/x; return result*result)
	K² = G.K*G.K
	a = invsq(G.ωp)
	b = (1-K²*invsq(G.ωz))
	c = -K²
	r = Main.BodePlots.roots(a,b,c) #Roots of ω²

	thresh = deg2rad(1) #WANTCONST: 1 degree
	rvalid = Float64[Inf,Inf]
	for i in 1:2
		val = r[i]
		if abs(angle(val)) < thresh #Positive real root
			rvalid[i] = sqrt(abs(val))
		end
	end

	x1 = minimum(rvalid) #Select first crossing of |G|=1
	#Don't return Inf - Likely Inf because roots are complex:
	return isfinite(x1) ? x1 : NaN
end

function step_err(s::Sys2, t::Vector)
	α = s.α #WANTCONST
	ωn = s.ωn #WANTCONST
	ζ = s.ζ #WANTCONST
	csqrt(x) = sqrt(complex(x))
	ζroot = csqrt(ζ^2-1)
	herr(f) = (ζ*(2α-1)/ζroot+f*1)*exp.(-ωn*t*(ζ+f*ζroot))

	if ζroot == 0 # ζ = 1
		return (1+(1-2α)*(ωn*t)).*exp.(-ωn*t)
	else
		return (herr(1)-herr(-1)) / 2
	end
end

step_resp(s::Sys2, t::Vector) = 1 .- step_err(s, t)

#==Main functions
===============================================================================#
function newplot()
	w = 500; h = w/1.6 #WANTCONST Select plot width/height
	yext_step = InspectDR.PExtents1D(0, 1.6) #WANTCONST

	mplot = InspectDR.Multiplot(title="PLL Characteristics")
	Δh = mplot.layout[:valloc_title]
	mplot.layout[:ncolumns] = 2
	#Bode plot looks better with wider aspect ratio:
	mplot.layout[:halloc_plot] = w; mplot.layout[:valloc_plot] = h/2

	plot_step = add(mplot, InspectDR.Plot2D)
	plot_bode = add(mplot, Main.BodePlots.new(InspectDR.Plot))
	plot_stepnorm = add(mplot, InspectDR.Plot2D)

	#Control where plots render:
	plot_step.plotbb = BoundingBox(0, w, Δh, Δh+h/2)
	plot_stepnorm.plotbb = BoundingBox(0, w, Δh+h/2, Δh+h)
	plot_bode.plotbb = BoundingBox(w, 2w, Δh, Δh+h)

	#Add title to bode plot:
	a = plot_bode.annotation
		a.title = "Open-Loop Characteristics"
	plot_bode.strips[1].yext_full = InspectDR.PExtents1D(-40, 100)

	#Step-response plot:
	a = plot_step.annotation
		a.title = "Step Response"
		a.xlabel = "Time (s)"
		a.ylabels = ["Amplitude"]
	plot_step.strips[1].yext_full = yext_step
#		plot_step.layout.grid = grid(vmajor=true, vminor=true, hmajor=false)

	a = plot_stepnorm.annotation; lyt = plot_stepnorm.layout
		a.title = ""; lyt[:valloc_top] = 8 #subtitle/space
		a.xlabel = "ωn*t (rad)"
		a.ylabels = ["Amplitude"]
	plot_stepnorm.strips[1].yext_full = yext_step

	return mplot
end

#Modify previously-generated Bode plot (base algorithm):
#-------------------------------------------------------------------------------
function update_step(plot::InspectDR.Plot2D, t::Vector, y::Vector)
	ldata = line(color=blue, width=3, style=:solid) #WANTCONST

	plot.data = [] #Clear old data

	#Reset extents:
	plot.xext = InspectDR.PExtents1D()
	plot.strips[1].yext = InspectDR.PExtents1D()

	wfrm = add(plot, t, mag.(y))
		wfrm.line = ldata

	return plot
end

#Update annotation on Bode plot (enabled=false to clear annotation):
#-------------------------------------------------------------------------------
function update_bodeannot(plot_bode::Plot2D, G::Xf1, ωn)
	lmarker = line(style=:dash, width=2.5) #WANTCONST
	lmarker_light = line(style=:solid, width=2.5, color=RGB24(.4,.4,.4)) #WANTCONST
	fp = G.ωp/(2pi) #WANTCONST
	fz = G.ωz/(2pi) #WANTCONST
	fn = ωn/(2pi) #WANTCONST
	ω0 = ω_0(G) #WANTCONST
	f0 = ω0/(2pi) #WANTCONST
	phase0 = phase.(xf_ol([ω0], G)[1]) #WANTCONST: Unity gain phase

	#Add specialized annotation
	afont = plot_bode.layout[:font_annotation]

	#Add annotation to Magnitude plot:
#	add(plot_bode, atext("0dB", y=0, xoffset_rel=0.5, yoffset=3, font=afont, align=:bc, strip=1))
	add(plot_bode, hmarker(0, lmarker_light, strip=1))

	polelist = [("fp", fp), ("fz", fz), ("fn", fn), ("f0", f0)]
	#Hack if f0>0: TODO: Fix bug where annotating text @ 0 on log plot causes glitch.
	for (id, f) in polelist
		if isfinite(f) && f > 0
			fstr = id * "=" * SI(f) * "Hz"
			add(plot_bode, atext(fstr, x=f, xoffset=-3, yoffset_rel=0.5, font=afont, angle=-90, align=:bc, strip=1))

			#Add vertical markers to both plots:
			add(plot_bode, vmarker(f, lmarker, strip=0))
		end
	end

	#Add annotation to Phase plot:
	if isfinite(phase0)
		pmargin = 180+phase0
		fstr = "PM=" * @sprintf("%.1f°", pmargin)
		add(plot_bode, atext(fstr, y=(phase0-180)/2, xoffset_rel=0.5, font=afont, align=:cc, strip=2))
		add(plot_bode, hmarker(phase0, lmarker, strip=2))
	end

	return
end

#Modify previously-generated plot:
#-------------------------------------------------------------------------------
#(Frequencies: in Hz; G: Open-loop characteristics)
#ωnperiods: in # of ωn periods
function update(mplot::InspectDR.Multiplot, fmin, fmax, tmax, ωnperiods, G::Xf1, annot=true)
	npts = 500 #WANTCONST: for plots
	plot_step = mplot.subplots[1] #WANTCONST
	plot_bode = mplot.subplots[2] #WANTCONST
	plot_stepnorm = mplot.subplots[3] #WANTCONST
	sys = Sys2(G) #WANTCONST: Closed-loop response

	f = 10 .^ range(log10(fmin), stop=log10(fmax), length=npts)
	y = xf_ol(2pi*f, G)
	Main.BodePlots.update(plot_bode, f, y)

	#Reset extents/annotation:
	plot_bode.xext = InspectDR.PExtents1D()
	plot_bode.strips[1].yext = InspectDR.PExtents1D()
	plot_bode.userannot = [] #Clear old markers/text annotation
	if annot; update_bodeannot(plot_bode, G, sys.ωn); end

	#Absolute step:
	tstep = tmax/npts
	t = collect(0:tstep:tmax)
	y = abs.(step_resp(sys, t))
	update_step(plot_step, t, y)

	#Relative step:
	xstep = ωnperiods/npts
	x = collect(0:xstep:ωnperiods)
	t = x./(sys.ωn/2pi)
	y = abs.(step_resp(sys, t))
	update_step(plot_stepnorm, x, y)

	return mplot
end

#==Show functions
===============================================================================#

function Base.show(io::IO, s::Sys2)
	print(io, "Sys2(α=", SI(s.α), ", ζ=", SI(s.ζ), ", ωn=", SI(s.ωn), ")")
end

end #module


text:
"""
$(SIGNATURES)

Return whether `point` is in `hull`.
"""
function Base.in(point::PointLike, hull::ConvexHull)
    op = orientation_comparator(hull)
    vertices = hull.vertices
    n = length(vertices)
    @inbounds begin
        if n === 0
            return false
        elseif n === 1
            return point == hull.vertices[1]
        elseif n === 2
            p′ = point - vertices[1]
            δ = vertices[2] - vertices[1]
            cross2(p′, δ) == 0 && 0 <= p′ ⋅ δ <= δ ⋅ δ
        else
            δ = vertices[1] - vertices[n]
            for i in Base.OneTo(n - 1)
                op(cross2(point - vertices[i], δ), 0) && return false
                δ = vertices[i + 1] - vertices[i]
            end
            return !op(cross2(point - vertices[n], δ), 0)
        end
    end
end