#########################################################
# file: magnetic_electric_lorentz.lsf
#
# Description:
# 	This script will visualize the values of eps
#          and mu that will be created when a 
#          "Magnetic Electric Lorentz" material
#          is used. These properties cannot be
#          seen in the Materials Explorer. The user
#          simply needs to set the name of the materials
#          and choose a wavelength range.
#	
# Copyright 2013, Lumerical Solutions, Inc.
#########################################################

#########################################################
# choose the material, base material (if there is one),
# and wavelength range for the results
clear;
matname = "negative index";
basematerial = "";  # set to "" if no base material is used
lambda = linspace(500e-9,1500e-9,1000);
f = c/lambda;
axis = 1; # relevant for anisotropic media only, can be 1,2,3
#########################################################

#########################################################
# the remainer of the script should not need to be modified
#########################################################

# show the material properties
?getmaterial(matname);

# get electric properties
chi0_e = getmaterial(matname,"χ0 electric");
deps = getmaterial(matname,"Δε electric");
delta_e = getmaterial(matname,"δ electric");
w_e = getmaterial(matname,"ω0 electric");
exclude_electric = getmaterial(matname,"exclude electric");

# get getmagnetic properties
chi0_m = getmaterial(matname,"χ0 magnetic");
dmu = getmaterial(matname,"Δμ magnetic");
delta_m = getmaterial(matname,"δ magnetic");
w_m = getmaterial(matname,"ω0 magnetic");
exclude_magnetic = getmaterial(matname,"exclude magnetic");

# base material properties
mu =  matrix(length(f)) + 1;
eps = matrix(length(f)) + 1;
if(basematerial != "") {
  eps = getindex(basematerial,f);
}

# create full index
w = 2*pi*f;
eps = eps + chi0_e(axis) + deps(axis)*w_e(axis)^2/(w_e(axis)^2-2i*delta_e(axis)*w-w^2);
mu = mu + chi0_m(axis) + dmu(axis)*w_m(axis)^2/(w_m(axis)^2-2i*delta_m(axis)*w-w^2);

# construct a dataset

epsilon = matrixdataset("epsilon");
epsilon.addparameter("lambda",lambda,"w",w);
epsilon.addattribute("epsilon",eps);
#epsilon.addparameter("f",f);
mum = matrixdataset("mu");
mum.addparameter("lambda",lambda,"w",w);
mum.addattribute("mu",mu);
#roperties.addattribute("mu",mu);

#
properties = matrixdataset("properties");
properties.addparameter("lambda",lambda,"f",f,"w",w);
properties.addattribute("eps",eps);
properties.addattribute("mu",mu);

# visualize results
visualize(epsilon);
visualize(mum);