11
Plasmonic materials

Figures in this chapter that were generated by Mathematica routines contained in the Code section are labelled "Mathematica simulation" in the figure caption.  The reader can regenerate all of these graphs automatically by selecting "Evaluate Notebook" under the "Evaluation" menu bar in Mathematica.  If Mathematica asks "Do you want to automatically evaluate all the initialization cells in the notebook?" it is sufficient to answer "No."  If the reader is interested in generating similar graphs for different parameter ranges or different materials, it is generally straightforward to do this as well.  The Code section should first be expanded by double clicking on the down arrow of the rightmost bracket.  This will expand the notebook to display "Setup" and "Code for generating figures" subsubsections.  If the notebook has not yet been evaluated, the "Setup" subsubsection should be executed first by selecting (clicking once on) its rightmost bracket and entering "Shift-Return" on a standard keyboard or just "Enter" on an extended keyboard.  Then the "Code for generating figures" subsection should be expanded by double clicking on the downward arrow on the rightmost bracket of that subsubsection.  This subsubsection can be further expanded, if necessary, by again double clicking on the right bracket with the downward pointing arrow.  After expanding the subsubsections, eventually a list of figure captions is obtained.  By double-clicking on the up arrow of the figure caption, a text section and a code section for that figure are displayed.  The parameters for generating the graph are described in the text section.  The reader simply needs to modify these parameters in the code section, which is right below the text section, and then replot the figure by selecting (clicking once on) the rightmost bracket that encloses the code section for the figure and hitting "Shift-Return" on the keyboard or just "Enter" on an extended keyboard.  The figure will then be regraphed using the new parameters.  The third line in the Setup subsection of the Code section is normally "option := optC" for generating color figures, but the line can be changed to "option := optB" for generating black and white figures instead if desired.
This notebook also contains a section with all of the figures in the text, a section with the exercises and code for generating relevant figures, and a section with all the chapter references.  Each section can, as usual, be opened by double clicking on the down arrow of the right bracket.

Code

Setup

ClearAll["Global`*"];
Needs["PlotLegends`"];
option:= optC;

Plotting Examples

Graphics Styles

fs=12;(*FontSize for LabelStyle*)
is={350,200};(*ImageSize for one figure*)
gr={400,150};(*ImageSize for GraphicsRow*)
gr2={480,180};(*ImageSize for GraphicsRow*)
gg={400,300};(*ImageSize for GraphicsGrid*)
gs=-30; (*GraphicsSpacing for two figures*)
ggs={-20,20};(*GraphicsSpacing for four figures*)
thinthick=Thickness[0.005];
thick=Thickness[0.007]; (*Thickness of lines*)

Graphics color option: optC

optC={
LabelStyle->{FontSize -> fs},
GridLines->Automatic,
Frame->True,
ImageSize->is,
PlotRange->All,
PlotStyle->{
{Red,Thick},
{Green,Dashing[{0.06,0.02}],Thick},{Blue,Dashing[{0.02,0.01}],Thick},
{Magenta,DotDashed,Thick},{Purple,Thick,Dashing[{0.04,0.01}]},{Orange,Thick,Dashing[{0.01,0.01}]}}};

Graphics black option: optB

optB={
GridLines -> Automatic,
Frame -> True,
ImageSize->is,
LabelStyle -> {FontSize -> fs},
PlotRange->{All,All},
PlotStyle->{
{Black,thick},
{Black,thick,Dashing[{0.06,0.02}]},
{Black,thick,Dashing[{0.02,0.01}]},
{Black,thick,DotDashed},
{Black, thick, Dashing[{0.04, 0.01}]},
{Black,thick,Dashing[{0.04, 0.03}]}}};

General two plots using optX

xmin=-π;xmax=π;
ymin=-1;ymax=1;
Fig01=Plot[
{Sin[x],Sin[2x],Sin[3x],Sin[4x]},
{x,xmin,xmax},
PlotRange->{ymin,ymax},
FrameLabel -> {{"x",""},{"sinx", ""}},
Evaluate[option]]

xmin=-π;xmax=π;
ymin=-1;ymax=1;
Fig02=Plot[
{Sin[x],Sin[2x]},
{x,xmin,xmax},Filling -> Axis,
FillingStyle -> Opacity[0.2],
PlotRange->{ymin,ymax},
FrameLabel -> {{"x",""},{"sinx", ""}},
Evaluate[option]]

GraphicsRow[{Fig01,Fig02},
ImageSize->gr,
Frame->False,
Spacings->gs]

General four plots using optX

GraphicsGrid[{
  {Fig01,Fig02},
  {Fig01,Fig02}},
 ImageSize -> gg,
 Frame -> False,
 Spacings -> gs]

Manipulate plots

amin=1;
amax=4;
Manipulate[
Plot[
 Sin[a θ],
{θ,0,6},
Frame->True,
FrameLabel->{"θ (deg)","sin(a θ)","Function",StringJoin["a = ",ToString[.01Ceiling[100.a]],""]},ImageSize->.75is,
PlotRange->{-1,1},
GridLines->Automatic,
PlotStyle->{
{Red,thinthick},
{Green,thinthick,Dashing[{0.06,0.02}]},{Blue,thinthick,Dashing[{0.02,0.01}]},
{Magenta,thinthick,DotDashed},
{Purple,thinthick,Dashing[{0.04,0.01}]}}],
{{a,amin,"Phase"},amin,amax}]

Tabulated metallic refractive indices

A note on refractive indices:  There are several sources of tabulated refractive indices for silver, gold and other metals in the literature and on the web.  The Handbook of Optical Constants of Solids edited by Palik [Academic, San Diego, 1998] is probably the most frequently cited source in recent years and is one of the more comprehensive sources.  For the noble metals the data of Johnson and Christy [Phys. Rev. B, 6 (1972) 4370-4379] are also frequently cited, especially in literature prior to the publication of the Handbook.  Indices from both of these sources have been included for use in simulations in this chapter.  However, the experimental precision in the Johnson and Christy data is somewhat limited and the data in the Handbook are frequently stitched together from several sources over an extended wavelength range.  The result of this is that calculations of extinction coefficients and peak field intensities often exhibit anomalous peaks at certain wavelengths corresponding to small changes in the refractive index from the limited experimental precision or a data stitch.  In order to avoid these unphysical effects in the modeling, it is generally much better to use an analytical model of the refractive indices that is smooth and continuous and has no precision problems.  Therefore, we have also included the analytical models developed by Rakic et al. [Appl. Opt. 37 (1998) 5271-5283] and generally recommend their use for all modeling.

Refractive indices for gold are taken from D. W. Lynch and W. R. Hunter, "Gold" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 293-295.  Refractive indices for silver are taken from D. W. Lynch and W. R. Hunter, "Silver" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 355-357.  Refractive indices for copper are taken from D. W. Lynch and W. R. Hunter, "Copper" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 284-285.  Refractive indices for aluminum are taken from D. Y. Smith, E. Shiles and M. Inokuti, "The optical properties of metallic aluminum" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 397-400.  Refractive indices for iridium are taken from D. W. Lynch and W. R. Hunter, "Iridium" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 302-303.  Refractive indices for rhodium are taken from D. W. Lynch and W. R. Hunter, "Rhodium" in Handbook of Optical Constants of Solids, ed. E. D. Palik (Academic, San Diego, 1998) 348-349.

GoldIndices = {{0.2883,1.742,1.9},{0.2952,1.776,1.918},{0.3024,1.812,1.92},{0.310, 1.83, 1.916}, {0.3179, 1.84, 1.904}, {0.3263, 1.824, 1.878}, {0.3351, 1.798, 1.86}, {0.3444, 1.766, 1.846}, {0.3542, 1.74, 1.848}, {0.3647, 1.716, 1.862}, {0.3757, 1.696, 1.906}, {0.3875, 1.674, 1.936}, {0.400, 1.658, 1.956}, {0.4133, 1.636, 1.958}, {0.4275, 1.616, 1.94}, {0.4428, 1.562, 1.904}, {0.4592, 1.426, 1.846}, {0.4769, 1.242, 1.796}, {0.4959, 0.916, 1.84}, {0.5166, 0.608, 2.12}, {0.5391, 0.402, 2.54}, {0.5636, 0.306, 2.88}, {0.5904, 0.236, 2.96}, {0.6199, 0.194, 3.05}, {0.6526, 0.166, 3.15}, {0.6888, 0.16, 3.8}, {0.7293, 0.164, 4.35},{0.7749, 0.174, 4.86}, {0.8266, 0.188, 5.39}, {0.8856, 0.21, 5.88}, {0.9537, 0.236, 6.47}, {1.033, 0.272, 7.07}, {1.127, 0.312, 7.93}, {1.240, 0.372, 8.77}, {1.265, 0.389, 8.09}, {1.291, 0.403, 8.25}, {1.319, 0.419, 8.42}, {1.378, 0.454, 8.77}, {1.409, 0.473, 8.96}, {1.442, 0.493, 9.15}, {1.476, 0.515, 9.36}, {1.512, 0.537, 9.58}, {1.550, 0.559, 9.81}, {1.590, 0.583, 10.1}, {1.631, 0.609, 10.3}, {1.675, 0.636, 10.6}, {1.722, 0.665, 10.9}, {1.771, 0.696, 11.2}, {1.823, 0.73, 11.5}, {1.879, 0.767, 11.9}, {1.937, 0.807, 12.2}, {2.000, 0.85, 12.6}, {2.066, 0.896, 13}, {2.138, 0.947, 13.4}, {2.214, 1.002, 13.9}, {2.296, 1.063, 14.4}, {2.384, 1.13, 14.9}, {2.480, 1.205, 15.5}, {2.583, 1.287, 16.1}, {2.695, 1.379, 16.8}, {2.818, 1.482, 17.5}, {2.952, 1.598, 18.3}};

SilverIndices = {{0.200,1.072,1.24},{0.20327869,1.098,1.26},{0.20666667,1.125,1.27},{0.2137931,1.173,1.29},{0.22142857,1.208,1.3},{0.22962963,1.238,1.31},{0.23846154,1.265,1.33},{0.248,1.298,1.35},{0.25306122,1.32,1.35},{0.25833333,1.343,1.35},{0.26382979,1.372,1.35},{0.26956522,1.404,1.33},{0.27555556,1.441,1.31},{0.28181818,1.476,1.26},{0.28837209,1.502,1.19},{0.2952381,1.519,1.08},{0.29879518,1.522,0.992},{0.30243902,1.496,0.882},{0.30617284,1.432,0.766},{0.310,1.323,0.647},{0.31155779,1.246,0.586},{0.31392405,1.149,0.54},{0.31552163,1.044,0.514},{0.31794872,0.932,0.504},{0.31958763,0.815,0.526},{0.32207792,0.708,0.565},{0.32375979,0.616,0.609},{0.32631579,0.526,0.663},{0.33066667,0.371,0.813},{0.33243968,0.321,0.902},{0.33513514,0.294,0.986},{0.33972603,0.259,1.12},{0.34444444,0.238,1.24},{0.35428571,0.209,1.44},{0.36470588,0.186,1.61},{0.37575758,0.2,1.67},{0.3875,0.192,1.81},{0.400,0.173,1.95},{0.41333333,0.173,2.11},{0.42758621,0.16,2.26},{0.44285714,0.157,2.4},{0.45925926,0.144,2.56},{0.47692308,0.132,2.72},{0.496,0.13,2.88},{0.51666667,0.13,3.07},{0.53913043,0.129,3.25},{0.56363636,0.12,3.45},{0.59047619,0.121,3.66},{0.620,0.131,3.88},{0.65263158,0.14,4.15},{0.68888889,0.14,4.44},{0.72941176,0.148,4.74},{0.775,0.143,5.09},{0.82666667,0.145,5.5},{0.88571429,0.163,5.95},{0.95384615,0.198,6.43},{1.03333333,0.226,6.99},{1.12727273,0.251,7.67},{1.240,0.329,8.49},{1.26530612,0.375,7.78},{1.29166667,0.383,7.92},{1.31914894,0.392,8.06},{1.34782609,0.401,8.21},{1.37777778,0.411,8.37},{1.40909091,0.421,8.37},{1.44186047,0.431,8.7},{1.47619048,0.442,8.88},{1.51219512,0.455,9.08},{1.550,0.469,9.32},{1.58974359,0.485,9.57},{1.63157895,0.501,9.84},{1.67567568,0.519,10.1},{1.72222222,0.537,10.4},{1.77142857,0.557,10.7},{1.87878788,0.6,11.4},{1.9375,0.624,11.8},{2.000,0.65,12.2}};

CopperIndices={{0.2952,1.42,1.64},{0.3099,1.34,1.72},{0.3263,1.34,1.81},{0.3444,1.31,1.87},{0.3646,1.27,1.95},{0.3874,1.18,2.21},{0.4133,1.18,2.21},{0.4428,1.17,2.36},{0.4768,1.15,2.5},{0.5166,1.12,2.6},{0.539,1.04,2.59},{0.5635,0.826,2.6},{0.5904,0.468,2.81},{0.6199,0.272,3.24},{0.6525,0.214,3.67},{0.6702,0.215,3.86},{0.688,0.213,4.05},{0.7084,0.214,4.24},{0.7293,0.223,4.43},{0.8265,0.26,5.26},{1.240,0.433,8.46},{1.265,0.496,6.78},{1.291,0.505,6.92},{1.319,0.515,7.06},{1.348,0.525,7.21},{1.378,0.536,7.36},{1.409,0.547,7.53},{1.442,0.559,7.7},{1.476,0.572,7.88},{1.512,0.586,8.06},{1.550,0.606,8.26},{1.589,0.627,8.47},{1.631,0.649,8.69},{1.675,0.672,8.92},{1.722,0.697,9.16},{1.771,0.723,9.41},{1.823,0.752,9.68},{1.878,0.782,9.97},{1.937,0.815,10.3},{2.000,0.85,10.6},{2.066,0.89,11},{2.138,0.933,11.3},{2.214,0.98,11.8},{2.296,1.03,12.2},{2.384,1.09,12.7},{2.480,1.15,13.2},{2.583,1.22,13.7}};

AluminumIndices={{0.300,0.276,3.61},{0.310,0.294,3.74},{0.3179,0.31,3.84},{0.3263,0.326,3.95},{0.3351,0.344,4.06},{0.3444,0.364,4.17},{0.350,0.375,4.24},{0.3542,0.385,4.3},{0.3647,0.407,4.43},{0.3757,0.432,4.56},{0.3875,0.46,4.71},{0.400,0.49,4.86},{0.4133,0.523,5.02},{0.4275,0.558,5.2},{0.4428,0.598,5.38},{0.450,0.618,5.47},{0.4592,0.644,5.58},{0.4769,0.695,5.8},{0.4959,0.755,6.03},{0.500,0.769,6.08},{0.5061,0.789,6.15},{0.5166,0.826,6.28},{0.5276,0.867,6.42},{0.5391,0.912,6.55},{0.550,0.958,6.69},{0.551,0.963,6.7},{0.5636,1.02,6.85},{0.5767,1.08,7},{0.5904,1.15,7.15},{0.600,1.2,7.26},{0.6048,1.22,7.31},{0.6199,1.3,7.48},{0.6358,1.39,7.65},{0.650,1.47,7.79},{0.6526,1.49,7.82},{0.6702,1.6,8.01},{0.6888,1.74,8.21},{0.700,1.83,8.31},{0.7085,1.91,8.39},{0.7293,2.14,8.57},{0.750,2.4,8.62},{0.7514,2.41,8.62},{0.7749,2.63,8.6},{0.7999,2.8,8.45},{0.825,2.75,8.31},{0.8266,2.74,8.31},{0.850,2.61,8.22},{0.8551,2.58,8.21},{0.875,2.38,8.18},{0.8856,2.24,8.21},{0.900,2.06,8.3},{0.9184,1.86,8.44},{0.925,1.77,8.49},{0.950,1.49,8.88},{0.9537,1.47,8.95},{0.9919,1.37,9.49},{1.000,1.35,9.58},{1.033,1.26,10},{1.078,1.21,10.6},{1.127,1.2,11.2},{1.181,1.21,11.8},{1.200,1.21,12},{1.240,1.21,12.5},{1.305,1.23,13.2},{1.378,1.26,14},{1.459,1.33,14.9},{1.500,1.38,15.4},{1.550,1.44,16},{1.653,1.59,17.1},{1.771,1.77,18.3},{1.907,1.99,19.8},{2.000,2.15,20.7}};

IridiumIndices={{0.2,0.94,1.68},{0.2033,0.94,1.73},{0.2066,0.95,1.78},{0.2101,0.97,1.82},{0.2138,0.96,1.86},{0.2175,0.97,1.92},{0.2214,0.98,1.98},{0.2254,0.99,2.03},{0.2296,1.0,2.09},{0.2339,1.02,2.15},{0.2384,1.04,2.22},{0.2431,1.07,2.29},{0.248,1.1,2.35},{0.253,1.13,2.42},{0.2583,1.18,2.49},{0.2638,1.24,2.56},{0.2695,1.31,2.6},{0.2755,1.37,2.65},{0.2818,1.45,2.68},{0.2883,1.51,2.7},{0.2952,1.58,2.71},{0.3024,1.62,2.7},{0.31,1.64,2.68},{0.3179,1.64,2.67},{0.3263,1.61,2.69},{0.3351,1.57,2.72},{0.3444,1.52,2.81},{0.3542,1.5,2.93},{0.3647,1.53,3.05},{0.3757,1.57,3.15},{0.3875,1.62,3.26},{0.4,1.68,3.35},{0.4133,1.73,3.43},{0.4275,1.77,3.51},{0.4428,1.81,3.61},{0.4592,1.85,3.73},{0.4769,1.91,3.86},{0.4959,1.98,4.00},{0.5166,2.07,4.14},{0.5391,2.18,4.26},{0.5636,2.29,4.38},{0.5904,2.4,4.48},{0.6199,2.5,4.57},{0.6526,2.57,4.68},{0.6888,2.64,4.81},{0.7293,2.69,4.92},{0.7749,2.68,5.08},{0.8266,2.65,5.39},{0.8856,2.72,5.74},{0.9537,2.85,6.07},{1.033,2.96,6.41},{1.127,3.04,6.84},{1.24,3.15,7.31},{1.378,3.19,7.88},{1.55,3.14,8.61},{1.771,2.93,9.78},{2.066,2.79,11.6}};

RhodiumIndices={{0.2000,0.78,1.85},{0.2066,0.76,1.93},{0.2138,0.79,2.00},{0.2214,0.8,2.06},{0.2296,0.8,2.14},{0.2384,0.79,2.23},{0.248,0.79,2.34},{0.253,0.79,2.39},{0.2583,0.79,2.46},{0.2638,0.79,2.52},{0.2695,0.78,2.6},{0.2755,0.79,2.68},{0.2818,0.8,2.76},{0.2883,0.83,2.85},{0.2952,0.84,2.94},{0.3024,0.84,3.03},{0.31,0.86,3.12},{0.3179,0.88,3.23},{0.3263,0.91,3.34},{0.3351,0.95,3.45},{0.3444,0.99,3.58},{0.3542,1.04,3.71},{0.3647,1.11,3.84},{0.3757,1.2,3.97},{0.3875,1.3,4.09},{0.4000,1.41,4.2},{0.4133,1.53,4.29},{0.4275,1.63,4.36},{0.4592,1.8,4.49},{0.4769,1.85,4.55},{0.4959,1.88,4.65},{0.5166,1.9,4.78},{0.5391,1.94,4.94},{0.5636,2.0,5.11},{0.5904,2.05,5.3},{0.6199,2.12,5.51},{0.6526,2.2,5.76},{0.6888,2.3,6.02},{0.7293,2.42,6.33},{0.7749,2.6,6.64},{0.8266,2.78,6.97},{0.8856,3.01,7.31},{0.9537,3.26,7.63},{0.9919,3.38,7.8},{1.033,3.51,7.94},{1.078,3.6,8.09},{1.127,3.67,8.26},{1.181,3.72,8.44},{1.240,3.71,8.67},{1.305,3.69,8.94},{1.378,3.62,9.36},{1.459,3.64,9.81},{1.550,3.63,10.3},{1.653,3.63,11.0},{1.771,3.67,11.7},{1.907,3.77,12.6},{2.066,3.87,13.5}};

GoldIndices2 and SilverIndices2 are from Johnson and Christy, "Optical constants of the noble metals," Phys. Rev. B, 6 (1972) 4370-4379.  Note that at short wavelengths there is a substantial difference between the refractive indices of Ag in this list and that from the Handbook of Optical Constants of Solids.  When comparing to examples in the literature it is important to use the same source for the optical constants!

GoldIndices2 = {{0.18788, 1.28, 1.188}, {0.19165, 1.32, 1.203}, {0.19528, 1.34, 1.226}, {0.19936, 1.33, 1.251}, {0.20328, 1.33, 1.277}, {0.20736, 1.3, 1.304}, {0.21197, 1.3, 1.35}, {0.2164, 1.3, 1.387}, {0.22143, 1.3, 1.427}, {0.22628, 1.31, 1.46}, {0.23134, 1.3, 1.497}, {0.23709, 1.32, 1.536}, {0.24266, 1.32, 1.577}, {0.2490, 1.33, 1.631}, {0.25514, 1.33, 1.688}, {0.2616, 1.35, 1.749}, {0.26898, 1.38, 1.803}, {0.27617, 1.43, 1.847}, {0.2844, 1.47, 1.869}, {0.29245, 1.49, 1.878}, {0.30097, 1.53, 1.889}, {0.31078, 1.53, 1.893}, {0.32041, 1.54, 1.898}, {0.33155, 1.48, 1.883}, {0.34254, 1.48, 1.871}, {0.35429, 1.5, 1.866}, {0.36795, 1.48, 1.895}, {0.38154, 1.46, 1.933}, {0.39744, 1.47, 1.952}, {0.41333, 1.46, 1.958}, {0.43056, 1.45, 1.948}, {0.45091, 1.38, 1.914}, {0.47148, 1.31, 1.849}, {0.4960, 1.04, 1.833}, {0.52101, 0.62, 2.081}, {0.54867, 0.43, 2.455}, {0.58216, 0.29, 2.863}, {0.61692, 0.21, 3.272}, {0.65957, 0.14, 3.697}, {0.70455, 0.13, 4.103}, {0.7561, 0.14, 4.542}, {0.82119, 0.16, 5.083}, {0.89209, 0.17, 5.663}, {0.98413, 0.22, 6.35}, {1.08772, 0.27, 7.15}, {1.21569, 0.35, 8.145}, {1.39326, 0.43, 9.519}, {1.61039, 0.56, 11.21}, {1.9375, 0.92, 13.78}};

SilverIndices2={{0.18788,1.07,1.212},{0.191.65,1.1,1.232},{0.19528,1.12,1.255},{0.19936,1.14,1.277},{0.20328,1.15,1.296},{0.20736,1.18,1.312},{0.21197,1.2,1.325},{0.21640,1.22,1.336},{0.22143,1.25,1.342},{0.22628,1.26,1.344},{0.23134,1.28,1.357},{0.23709,1.28,1.367},{0.24266,1.3,1.378},{0.24900,1.31,1.389},{0.25514,1.33,1.393},{0.26160,1.35,1.387},{0.26898,1.38,1.372},{0.27617,1.41,1.331},{0.28440,1.41,1.264},{0.29245,1.39,1.161},{0.30097,1.34,0.964},{0.31078,1.13,0.616},{0.32041,0.81,0.392},{0.33155,0.17,0.829},{0.34254,0.14,1.142},{0.35429,0.1,1.419},{0.36795,0.07,1.657},{0.38154,0.05,1.864},{0.39744,0.05,2.07},{0.41333,0.05,2.275},{0.43056,0.04,2.462},{0.45091,0.04,2.657},{0.47148,0.05,2.869},{0.49600,0.05,3.093},{0.52101,0.05,3.324},{0.54867,0.06,3.586},{0.58216,0.05,3.858},{0.61692,0.06,4.152},{0.65957,0.05,4.483},{0.70455,0.04,4.838},{0.75610,0.03,5.242},{0.82119,0.04,5.727},{0.89209,0.04,6.312},{0.98413,0.04,6.992},{1.08772,0.04,7.795},{1.21569,0.09,8.828},{1.39326,0.13,10.1},{1.61039,0.15,11.85},{1.9375,0.24,14.08}};

IronIndices2 and CobaltIndices2 are from Johnson and Christy, "Optical constants of transition metals: Ti, V, Cr, Mn, Fe, Co, Ni, and Pd," Phys. Rev. B, 9 (1974) 5056 - 5070.

IronIndices2 = {{0.18788,1.29,1.35},{0.19165,1.35,1.37},{0.19528,1.42,1.39},{0.19936,1.45,1.4},{0.20328,1.47,1.4},{0.20736,1.49,1.41},{0.21197,1.47,1.43},{0.2164,1.47,1.44},{0.22143,1.47,1.47},{0.22628,1.47,1.49},{0.23134,1.48,1.53},{0.23709,1.48,1.57},{0.24266,1.5,1.61},{0.249,1.51,1.66},{0.25514,1.53,1.7},{0.2616,1.56,1.75},{0.26898,1.59,1.79},{0.27617,1.62,1.84},{0.2844,1.64,1.88},{0.29245,1.65,1.94},{0.30097,1.67,2},{0.31078,1.69,2.06},{0.32041,1.74,2.12},{0.33155,1.78,2.19},{0.34254,1.85,2.27},{0.35429,1.93,2.35},{0.36795,2.02,2.43},{0.38154,2.12,2.5},{0.39744,2.24,2.58},{0.41333,2.35,2.65},{0.43056,2.48,2.71},{0.45091,2.59,2.77},{0.47148,2.67,2.82},{0.496,2.74,2.88},{0.52101,2.86,2.91},{0.54867,2.95,2.93},{0.58216,2.94,2.99},{0.61692,2.88,3.05},{0.65957,2.92,3.1},{0.70455,2.86,3.19},{0.7561,2.87,3.28},{0.82119,2.94,3.39},{0.89209,2.96,3.56},{0.98413,2.92,3.79},{1.08772,2.97,4.06},{1.21569,3.03,4.39},{1.39326,3.09,4.83},{1.61039,3.11,5.39},{1.9375,3.17,6.12}};

CobaltIndices2={{0.18788,1.16,1.59},{0.19165,1.21,1.63},{0.19528,1.26,1.67},{0.19936,1.29,1.71},{0.20328,1.32,1.75},{0.20736,1.36,1.78},{0.21197,1.38,1.82},{0.2164,1.41,1.84},{0.22143,1.43,1.85},{0.22628,1.45,1.86},{0.23134,1.47,1.87},{0.23709,1.47,1.89},{0.24266,1.46,1.91},{0.249,1.45,1.93},{0.25514,1.45,1.97},{0.2616,1.44,2.01},{0.26898,1.44,2.04},{0.27617,1.44,2.09},{0.2844,1.44,2.14},{0.29245,1.44,2.19},{0.30097,1.44,2.25},{0.31078,1.44,2.31},{0.32041,1.45,2.37},{0.33155,1.46,2.44},{0.34254,1.48,2.52},{0.35429,1.49,2.61},{0.36795,1.5,2.71},{0.38154,1.53,2.82},{0.39744,1.57,2.93},{0.41333,1.61,3.05},{0.43056,1.67,3.17},{0.45091,1.74,3.28},{0.47148,1.81,3.41},{0.496,1.88,3.55},{0.52101,1.97,3.68},{0.54867,2.05,3.82},{0.58216,2.13,3.96},{0.61692,2.19,4.11},{0.65957,2.25,4.27},{0.70455,2.31,4.45},{0.7561,2.4,4.64},{0.82119,2.53,4.88},{0.89209,2.65,5.16},{0.98413,2.78,5.5},{1.08772,2.94,5.88},{1.21569,3.17,6.31},{1.39326,3.42,6.77},{1.61039,3.61,7.26},{1.9375,3.87,7.79}};

Routines to compute refractive indices of various metals

Drude[λ_]:=Module[{n,ω,ϵ,c,h,ωp,γ,ϵ∞},
c=2.99792456×(10^14);          (*  speed of light in μm/s  *)
h=4.1359×(10^(-15));                (*  Plank's constant in eV s  *)
ϵ∞=4;                                     (*  background dielectric constant  *)
ωp=15;                                       (*  ωp in eV  *)
γ = 0.1;                               (*  damping coefficient in eV  *)
ω = 7.7905/λ;                     (*  ang. frequency in eV when λ is in μm  *)
ϵ=ϵ∞-(ωp^2)/(ω^2 + I γ ω);
n = (Sqrt[ϵ]);
(*  Print["λ=",λ," n=",n];  *)
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271 for 0.1 to 5 eV  *)
LorentzDrudeAg[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.065,0.124,0.011,0.840,5.646};
γ={3.886,0.452,0.065,0.916,2.419};
ωω={0.816,4.481,8.185,9.083,20.29};
ωp=9.01;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.845(ωp^2)/(ω^2+I 0.048ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,5}];
n = (Sqrt[ϵ]);
Return[n]];


(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271 for 0.2 to 5 eV  *)
LorentzDrudeAu[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.024,0.010,0.071,0.601,4.384};
γ={0.241,0.345,0.870,2.494,2.214};
ωω={0.415,0.830,2.969,4.304,13.32};
ωp=9.03;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.760(ωp^2)/(ω^2+I 0.053ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,5}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271 for 0.1 to 6 eV  *)
LorentzDrudeCu[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.061,0.104,0.723,0.638};
γ={0.378,1.056,3.213,4.305};
ωω={0.291,2.957,5.300,11.18};
ωp=10.83;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.575(ωp^2)/(ω^2+I 0.030ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271 for 0.01 to 10 eV  *)
LorentzDrudeAl[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.227,0.050,0.166,0.030};
γ={0.333,0.312,1.351,3.382};
ωω={0.162,1.544,1.808,3.473};
ωp=14.98;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.523(ωp^2)/(ω^2+I 0.047ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271 for 0.02 to 5 eV  *)
LorentzDrudeBe[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.031,0.140,0.530,0.130};
γ={1.664,3.395,4.454,1.802};
ωω={0.100,1.032,3.183,4.604};
ωp=18.51;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.084(ωp^2)/(ω^2+I 0.035ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271 0.02 to 5 eV  *)
LorentzDrudeCr[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.151,0.150,1.149,0.825};
γ={3.175,1.305,2.676,1.335};
ωω={0.121,0.543,1.970,8.775};
ωp=10.75;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.168(ωp^2)/(ω^2+I 0.047ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271 for 0.2 to 5 eV  *)
LorentzDrudeNi[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.100,0.135,0.106,0.729};
γ={4.511,1.334,2.178,6.292};
ωω={0.174,0.582,1.597,6.089};
ωp=15.92;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.096(ωp^2)/(ω^2+I 0.048ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271 for 0.1 to 5 eV  *)
LorentzDrudePd[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.649,0.121,0.638,0.453};
γ={2.950,1.555,4.621,3.236};
ωω={0.336,0.501,1.659,5.715};
ωp=15.92;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.330(ωp^2)/(ω^2+I 0.008ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271 for 0.1 to 5 eV  *)
LorentzDrudePt[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.191,0.659,0.547,3.576};
γ={0.517,1.838,3.668,8.517};
ωω={0.780,1.314,3.141,9.249};
ωp=15.92;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.333(ωp^2)/(ω^2+I 0.080ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271 for 0.06 to 5 eV  *)
LorentzDrudeTi[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.899,0.393,0.187,0.001};
γ={2.276,2.518,1.663,1.762};
ωω={0.777,1.545,2.509,19.43};
ωp=15.92;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.148(ωp^2)/(ω^2+I 0.082ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271 for 0.1 to 5 eV  *)
LorentzDrudeW[λ_]:=Module[{n,ω,f,ωω,ϵ,c,ωp,γ,ϵ∞},
f={0.054,0.166,0.706,2.590};
γ={0.530,1.281,3.332,5.836};
ωω={1.004,1.917,3.580,7.498};
ωp=15.92;
ω = 1.24/λ;                 (*  ang. frequency in eV when λ is in μm  *)
Do[If[i==0,ϵ=1-(0.206(ωp^2)/(ω^2+I 0.064ω )),
ϵ=ϵ+f[[i]] ((ωp)^2)/((ωω[[i]])^2-ω^2-I ω γ[[i]])],{i,0,4}];
n = (Sqrt[ϵ]);
Return[n]];

(*  from D.W.Lynch and W.R.Hunter,"Gold" in Handbook of Optical Constants of Solids,ed.E.D.Palik (Academic,San Diego,1998) 293-295  *)
Gold[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<GoldIndices[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>GoldIndices[[Length[GoldIndices],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=GoldIndices[[i,1]];
n2=GoldIndices[[i,2]];
K2=GoldIndices[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=GoldIndices[[i-1,1]];
n1=GoldIndices[[i-1,2]];
K1=GoldIndices[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[GoldIndices]}]];
Return[n]];

(*  from D.W.Lynch and W.R.Hunter "Iridium" in Handbook of Optical Constants of Solids,ed.E.D.Palik (Academic,San Diego,1998) 302-303  *)
Iridium[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<IridiumIndices[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>IridiumIndices[[Length[IridiumIndices],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=IridiumIndices[[i,1]];
n2=IridiumIndices[[i,2]];
K2=IridiumIndices[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=IridiumIndices[[i-1,1]];
n1=IridiumIndices[[i-1,2]];
K1=IridiumIndices[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[IridiumIndices]}]];
Return[n]];

(*  from D.W.Lynch and W.R.Hunter "Rhodium" in Handbook of Optical Constants of Solids,ed.E.D.Palik (Academic,San Diego,1998) 348-349  *)
Rhodium[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<RhodiumIndices[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>RhodiumIndices[[Length[RhodiumIndices],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=RhodiumIndices[[i,1]];
n2=RhodiumIndices[[i,2]];
K2=RhodiumIndices[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=RhodiumIndices[[i-1,1]];
n1=RhodiumIndices[[i-1,2]];
K1=RhodiumIndices[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[RhodiumIndices]}]];
Return[n]];

(*  from Johnson and Christy,"Optical constants of the noble metals,"Phys. Rev. B,6(1972) 4370-4379  *)
Gold2[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<GoldIndices2[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>GoldIndices2[[Length[GoldIndices2],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=GoldIndices2[[i,1]];
n2=GoldIndices2[[i,2]];
K2=GoldIndices2[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=GoldIndices2[[i-1,1]];
n1=GoldIndices2[[i-1,2]];
K1=GoldIndices2[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[GoldIndices2]}]];
Return[n]];

(*  from D.W.Lynch and W.R.Hunter "Silver" in Handbook of Optical Constants of Solids ed.E.D.Palik (Academic,San Diego,1998) 355-357  *)
Silver[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<SilverIndices[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>SilverIndices[[Length[SilverIndices],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=SilverIndices[[i,1]];
n2=SilverIndices[[i,2]];
K2=SilverIndices[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=SilverIndices[[i-1,1]];
n1=SilverIndices[[i-1,2]];
K1=SilverIndices[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[SilverIndices]}]];
Return[n]];

(*  from Johnson and Christy,"Optical constants of the noble metals,"Phys. Rev. B,6(1972) 4370-4379  *)
Silver2[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<SilverIndices2[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>SilverIndices2[[Length[SilverIndices2],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=SilverIndices2[[i,1]];
n2=SilverIndices2[[i,2]];
K2=SilverIndices2[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=SilverIndices2[[i-1,1]];
n1=SilverIndices2[[i-1,2]];
K1=SilverIndices2[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[SilverIndices2]}]];
Return[n]];

(*  from Johnson and Christy, Phys. Rev. B,9(1974) 5056-5070  *)
Iron2[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<IronIndices2[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>IronIndices2[[Length[IronIndices2],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=IronIndices2[[i,1]];
n2=IronIndices2[[i,2]];
K2=IronIndices2[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=IronIndices2[[i-1,1]];
n1=IronIndices2[[i-1,2]];
K1=IronIndices2[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[IronIndices2]}]];
Return[n]];

(*  from Johnson and Christy, Phys. Rev. B,9(1974) 5056-5070  *)
Cobalt2[λ_]:=Module[{n,n1,n2,K1,K2,λ1,λ2},
If[λ<CobaltIndices2[[1,1]],Print[" Input wavelength is too short."];Abort[]];
If[λ>CobaltIndices2[[Length[CobaltIndices2],1]],Print[" Input wavelength is too long."];Abort[]];
Catch[Do[λ2=CobaltIndices2[[i,1]];
n2=CobaltIndices2[[i,2]];
K2=CobaltIndices2[[i,3]];
If[λ==λ2,n=n2+I K2;Throw[n],
If[λ<λ2,λ1=CobaltIndices2[[i-1,1]];
n1=CobaltIndices2[[i-1,2]];
K1=CobaltIndices2[[i-1,3]];
n=n1 + (n2-n1)(λ-λ1)/(λ2-λ1) + I(K1 + (K2-K1)(λ-λ1)/(λ2-λ1));
Throw[n]]],{i,Length[CobaltIndices2]}]];
Return[n]];

(*  from D.W.Lynch and W.R.Hunter "Copper" in Handbook of Optical Constants of Solids,ed.E.D.Palik (Academic,San Diego,1998) 284-285  *)
Copper[λ_] := Module[{n, n1, n2, K1, K2, λ1, λ2},
   If[λ < CopperIndices[[1, 1]], Print[" Input wavelength is too short."]; Abort[]];
   If[λ > CopperIndices[[Length[CopperIndices], 1]], Print[" Input wavelength is too long."]; Abort[]];
   Catch[Do[λ2 = CopperIndices[[i, 1]];
     n2 = CopperIndices[[i, 2]];
     K2 = CopperIndices[[i, 3]];
     If[λ == λ2, n = n2 + I K2; Throw[n],
      If[λ < λ2, λ1 = CopperIndices[[i - 1, 1]];
       n1 = CopperIndices[[i - 1, 2]];
       K1 = CopperIndices[[i - 1, 3]];
       n = n1 + (n2 - n1) (λ - λ1)/(λ2 - λ1) + I (K1 + (K2 - K1) (λ - λ1)/(λ2 - λ1));
       Throw[n]]], {i, Length[CopperIndices]}]];
   Return[n]];

(*  from D.Y.Smith E.Shiles and M.Inokuti "The optical properties of metallic aluminum" in Handbook of Optical Constants of Solids ed.E.D.Palik (Academic,San Diego,1998) 397-400  *)
Aluminum[λ_] := Module[{n, n1, n2, K1, K2, λ1, λ2},
   If[λ < AluminumIndices[[1, 1]], Print[" Input wavelength is too short."]; Abort[]];
   If[λ > AluminumIndices[[Length[AluminumIndices], 1]], Print[" Input wavelength is too long."]; Abort[]];
   Catch[Do[λ2 = AluminumIndices[[i, 1]];
     n2 = AluminumIndices[[i, 2]];
     K2 = AluminumIndices[[i, 3]];
     If[λ == λ2, n = n2 + I K2; Throw[n],
      If[λ < λ2, λ1 = AluminumIndices[[i - 1, 1]];
       n1 = AluminumIndices[[i - 1, 2]];
       K1 = AluminumIndices[[i - 1, 3]];
       n = n1 + (n2 - n1) (λ - λ1)/(λ2 - λ1) + I (K1 + (K2 - K1) (λ - λ1)/(λ2 - λ1));
       Throw[n]]], {i, Length[AluminumIndices]}]];
   Return[n]];

Various routines

The dielectric function is computed from the refractive index in this routine.

ϵ[λ_,n_]:=Module[{nmetal,ϵm},
nmetal=Which[n==-1,Gold[λ],n==-2,Silver[λ],n==-3,Copper[λ],n==-4,Aluminum[λ],n==-5,Drude[λ],n==-6,LorentzDrudeAu[λ],n==-7,LorentzDrudeAg[λ],n==-8,LorentzDrudeCu[λ],n==-9,LorentzDrudeAl[λ],n==-10,Gold2[λ],n==-11,Silver2[λ],n==-12,LorentzDrudeBe[λ],n==-13,LorentzDrudeCr[λ],n==-14,LorentzDrudeNi[λ],n==-15,LorentzDrudePd[λ],n==-16,LorentzDrudePt[λ],n==-17,LorentzDrudeTi[λ],n==-18,LorentzDrudeW[λ],n==-19,Iron2[λ],n==-20,Cobalt2[λ],n==-21,Rhodium[λ],n==-22,Iridium[λ],True,n];
ϵm=(nmetal^2);
Return[ϵm]];

The observability function defined in Eq. (5.4) is computed in this routine.

f[λ_,n_]:=Module[{nmetal,ϵm},
nmetal=Which[n==-1,Gold[λ],n==-2,Silver[λ],n==-3,Copper[λ],n==-4,Aluminum[λ],n==-5,Drude[λ],n==-6,LorentzDrudeAu[λ],n==-7,LorentzDrudeAg[λ],n==-8,LorentzDrudeCu[λ],n==-9,LorentzDrudeAl[λ],n==-10,Gold2[λ],n==-11,Silver2[λ],n==-12,LorentzDrudeBe[λ],n==-13,LorentzDrudeCr[λ],n==-14,LorentzDrudeNi[λ],n==-15,LorentzDrudePd[λ],n==-16,LorentzDrudePt[λ],n==-17,LorentzDrudeTi[λ],n==-18,LorentzDrudeW[λ],n==-19,Iron2[λ],n==-20,Cobalt2[λ],n==-21,Rhodium[λ],n==-22,Iridium[λ],True,n];
ϵm=(nmetal^2);
Return[Im[Sqrt[ϵm/(ϵm+1)]]]];

The figure of merit function defined in Eq. (5.5) is computed in this routine.

g[λ_,n_]:=Module[{nmetal,ϵm},
nmetal=Which[n==-1,Gold[λ],n==-2,Silver[λ],n==-3,Copper[λ],n==-4,Aluminum[λ],n==-5,Drude[λ],n==-6,LorentzDrudeAu[λ],n==-7,LorentzDrudeAg[λ],n==-8,LorentzDrudeCu[λ],n==-9,LorentzDrudeAl[λ],n==-10,Gold2[λ],n==-11,Silver2[λ],n==-12,LorentzDrudeBe[λ],n==-13,LorentzDrudeCr[λ],n==-14,LorentzDrudeNi[λ],n==-15,LorentzDrudePd[λ],n==-16,LorentzDrudePt[λ],n==-17,LorentzDrudeTi[λ],n==-18,LorentzDrudeW[λ],n==-19,Iron2[λ],n==-20,Cobalt2[λ],n==-21,Rhodium[λ],n==-22,Iridium[λ],True,n];
ϵm=(nmetal^2);
Return[3 Abs[Re[ϵm]/Im[ϵm]]]];

The reflectivity in the Kretschmann configuration is calculated in this routine.  The input parameters are the wavelength λ, the angle of incidence θ, the index of the incident medium (Subscript[n, 0]), the metal index (Subscript[n, 1]), the index of the transmitting medium (Subscript[n, 2]), and the thickness of the metal film d.

r[λ_,θ_,n0_,n1_,n2_,d_]:=Module[{nmetal,k0,cos0,cos1,cos2,sin1,sin2,a1,a2,b1,b2,exp1,refl,tran,R,T},
k0=2 π/λ;
nmetal=Which[n1==-1,Gold[λ],n1==-2,Silver[λ],n1==-3,Copper[λ],n1==-4,Aluminum[λ],n1==-5,Drude[λ],n1==-6,LorentzDrudeAu[λ],n1==-7,LorentzDrudeAg[λ],n1==-8,LorentzDrudeCu[λ],n1==-9,LorentzDrudeAl[λ],n1==-10,Gold2[λ],n1==-11,Silver2[λ],n1==-12,LorentzDrudeBe[λ],n1==-13,LorentzDrudeCr[λ],n1==-14,LorentzDrudeNi[λ],n1==-15,LorentzDrudePd[λ],n1==-16,LorentzDrudePt[λ],n1==-17,LorentzDrudeTi[λ],n1==-18,LorentzDrudeW[λ],n1==-19,Iron2[λ],n1==-20,Cobalt2[λ],n1==-21,Rhodium[λ],n1==-22,Iridium[λ],True,n1];
cos0=Cos[2 π θ/360];
sin1=n0 Sin[2 π θ/360]/nmetal;
cos1=Sqrt[1-(sin1^2)];
sin2=n0 Sin[2 π θ/360]/n2;
cos2=Sqrt[1-(sin2^2)];
a1=(n0 cos1-nmetal cos0)/(n2 cos1-nmetal cos2);
a2=(nmetal cos0+n0 cos1)/(n2 cos1+nmetal cos2);
b1=(n0 cos1-nmetal cos0)/(n2 cos1+nmetal cos2);
b2=(nmetal cos0+n0 cos1)/(n2 cos1-nmetal cos2);
exp1=Exp[2 I nmetal k0 d cos1];
refl=(a1-a2 exp1)/(b1 exp1-b2);
R=(Abs[refl])^2;
tran=(a2+b1 refl)Exp[ I nmetal k0 d cos1];
T= n2  Re[cos2] (Abs[tran])^2/(n0 cos0);
Return[R]];

The field in the medium outside of the prism in the Kretschmann configuration is calculated in this routine.  The input parameters are the wavelength λ, the angle of incidence θ, the index of the incident medium (Subscript[n, 0]), the metal index (Subscript[n, 1]), the index of the transmitting medium (Subscript[n, 2]), the thickness of the metal film d, and the distance below the surface z.

field[λ_,θ_,n0_,n1_,n2_,d_,z_]:=Module[{nmetal,k0,cos0,cos1,cos2,sin1,sin2,a1,a2,b1,b2,exp1,refl,tran,R,T,Ex,Ez,E},
k0=2 π/λ;
nmetal=Which[n1==-1,Gold[λ],n1==-2,Silver[λ],n1==-3,Copper[λ],n1==-4,Aluminum[λ],n1==-5,Drude[λ],n1==-6,LorentzDrudeAu[λ],n1==-7,LorentzDrudeAg[λ],n1==-8,LorentzDrudeCu[λ],n1==-9,LorentzDrudeAl[λ],n1==-10,Gold2[λ],n1==-11,Silver2[λ],n1==-12,LorentzDrudeBe[λ],n1==-13,LorentzDrudeCr[λ],n1==-14,LorentzDrudeNi[λ],n1==-15,LorentzDrudePd[λ],n1==-16,LorentzDrudePt[λ],n1==-17,LorentzDrudeTi[λ],n1==-18,LorentzDrudeW[λ],n1==-19,Iron2[λ],n1==-20,Cobalt2[λ],n1==-21,Rhodium[λ],n1==-22,Iridium[λ],True,n1];
cos0=Cos[2 π θ/360];
sin1=n0 Sin[2 π θ/360]/nmetal;
cos1=Sqrt[1-(sin1^2)];
sin2=n0 Sin[2 π θ/360]/n2;
cos2=Sqrt[1-(sin2^2)];
a1=(n0 cos1-nmetal cos0)/(n2 cos1-nmetal cos2);
a2=(nmetal cos0+n0 cos1)/(n2 cos1+nmetal cos2);
b1=(n0 cos1-nmetal cos0)/(n2 cos1+nmetal cos2);
b2=(nmetal cos0+n0 cos1)/(n2 cos1-nmetal cos2);
exp1=Exp[2 I nmetal k0 d cos1];
refl=(a1-a2 exp1)/(b1 exp1-b2);
R=(Abs[refl])^2;
tran=(a2+b1 refl)Exp[ I nmetal k0 d cos1];
T= n2  Re[cos2] (Abs[tran])^2/(n0 cos0);
Ex=cos2 tran Exp[I n2 k0 cos2 z];
Ez=sin2 tran Exp[I n2 k0 cos2 z];
E=Sqrt[Abs[Ex]^2 + Abs[Ez]^2];
Return[E]];

The SP propagation distance is computed in this function.  The parameters passed to the function are the wavelength, the index of the dielectric, and the index of the metal.

dist[λ_, nd_, nm_] := Module[{nmetal, ϵ1, ϵ2, sin1, k1},
      nmetal = Which[nm == -1, Gold[λ], nm == -2, Silver[λ], nm == -3, Copper[λ], nm == -4, Aluminum[λ], nm == -5, Drude[λ], nm == -6, LorentzDrudeAu[λ], nm == -7, LorentzDrudeAg[λ], nm == -8, LorentzDrudeCu[λ], nm == -9, LorentzDrudeAl[λ], nm == -10, Gold2[λ], nm == -11, Silver2[λ], nm == -12, LorentzDrudeBe[λ], nm == -13, LorentzDrudeCr[λ], nm == -14, LorentzDrudeNi[λ], nm == -15, LorentzDrudePd[λ], nm == -16, LorentzDrudePt[λ], nm == -17, LorentzDrudeTi[λ], nm == -18, LorentzDrudeW[λ], nm == -19, Iron2[λ], nm == -20, Cobalt2[λ], nm == -21, Rhodium[λ], nm == -22, Iridium[λ], True, nm];
      ϵ1 = nd^2; ϵ2 = nmetal^2;
      sin1 = Sqrt[ϵ2/(ϵ1 + ϵ2)];
      k1 = nd (2 π/λ);
      Return[1/Im[k1 sin1]]];

The SP lifetime is computed in this function.  The parameters passed to the function are the wavelength, the index of the dielectric, and the index of the metal.

τsp[λ_, nd_, nm_] := Module[{nmetal, ϵ1, ϵ2,x, k0,c},
      nmetal = Which[nm == -1, Gold[λ], nm == -2, Silver[λ], nm == -3, Copper[λ], nm == -4, Aluminum[λ], nm == -5, Drude[λ], nm == -6, LorentzDrudeAu[λ], nm == -7, LorentzDrudeAg[λ], nm == -8, LorentzDrudeCu[λ], nm == -9, LorentzDrudeAl[λ], nm == -10, Gold2[λ], nm == -11, Silver2[λ], nm == -12, LorentzDrudeBe[λ], nm == -13, LorentzDrudeCr[λ], nm == -14, LorentzDrudeNi[λ], nm == -15, LorentzDrudePd[λ], nm == -16, LorentzDrudePt[λ], nm == -17, LorentzDrudeTi[λ], nm == -18, LorentzDrudeW[λ], nm == -19, Iron2[λ], nm == -20, Cobalt2[λ], nm == -21, Rhodium[λ], nm == -22, Iridium[λ], True, nm];
      ϵ1 = nd^2; ϵ2 = nmetal^2;
      x = Sqrt[ϵ1 ϵ2/(ϵ1 + ϵ2)];
      k0 = (2 π/λ)×(10^6);   (*  free space wavevector in SI units  *)
c=2.99792456×(10^8);   (*  speed of light in SI units  *)
      Return[((10^15))Re[x/(c k0 Im[x])]]];    (* return lifetime in fs  *)

Calculation of the complex SP wavevector for a Drude metal with damping assuming the SP frequency is real.

Ksp:=Ω Sqrt[(Ω^2 -1+I Γ Ω)/(2 Ω^2 -1 + 2 I Γ Ω)];

Calculation of the complex SP frequency for a Drude metal with damping assuming the SP wavevector is real.

Ωsp=Module[{Ω,func,t},
func=Ω^4+I Γ Ω^3-(1+2K^2)Ω^2-I 2 Γ K^2 Ω+K^2;
t=NSolve[func==0,Ω];
Ω/.t];

This program computes the extinction and scattering coefficients for a sphere from Mie theory.  It has been adapted from the Fortran code in C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, New York, 1983) Appendix A.  The SphereScattering function can be used to compute extinction, scattering, backscattering, and near field efficiencies as a function of wavelength or sphere radius.  It can also be used to compute polar diagrams of the scattered light intensity.  For calculations as a function of wavelength or sphere radius, the function returns four lists of (x, y) coordinates for plotting.  For each list, x is the wavelength in microns. y is the scattering efficiency in the first list, the extinction efficiency in the second list, the backscattering efficiency in the third list, and the near field efficiency in the fourth list.  For polar scattering diagrams, the function also returns four lists.  Let (Subscript[i, ⊥]) be the scattered irradiance per unit incident irradiance when the incident light is polarized perpendicular to the scattering plane, and let (Subscript[i, ∥]) be the scattered irradiance per unit incident irradiance when the incident light is polarized parallel to the scattering plane.  The first list is the set of coordinates (θ, (Subscript[i, ⊥])) as θ goes from 0° to 180°.  The second list is the same set of data expressed in a Cartesian frame, ((Subscript[i, ⊥]) cos(θ), (Subscript[i, ⊥])sin(θ)).  The third list is the set of coordinates (θ, (Subscript[i, ∥])) as θ goes from 0° to 180°, and the fourth list is the same set of data expressed in a Cartesian frame, ((Subscript[i, ∥]) cos(θ), (Subscript[i, ∥])sin(θ)).

The parameters that are set by the user before a calculation are as follows:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[θ, inc]): increment of scattering angle in degrees
	(Subscript[r, 1]): initial radius of sphere in microns
	(Subscript[r, 2]): final radius of sphere in microns
	(Subscript[r, inc]): increment of radius of sphere in microns
	(Subscript[n, 1]): refractive index of outer medium
	(Subscript[n, 2]): refractive index of sphere

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes efficiencies vs. wavelength.  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes efficiencies vs. sphere radius.  Finally, if (Subscript[λ, 1])=(Subscript[λ, 2]) and (Subscript[r, 1])=(Subscript[r, 2]), the function computes the polar diagrams wiht an angle increment of (Subscript[θ, inc]).
The user can either set the refractive index of the sphere, (Subscript[n, 2]), to a specific (complex) value, or the user can set (Subscript[n, 2]) to a negative integer which tells the function to insert the wavelength-dependent refractive index for a particular metal as follows: 
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.

mie[λ_,r_,n1_,n2_,nang_ ]:=Module[{s1,s2,Qsct,Qext,Qbck,Qnf,nsph,x,y,nstop,nang1,nmx,amu,d,rn,π0,π1,ψ0,ψ1,ψ2,ψ3,χ0,χ1,χ2,χ3,ξ0,ξ1,ξ2,ξ3,h0,h1,h2,h3,fn,jj,pi,τ,p,t,an,bn,m},
If[nang==1,nang1=2,nang1=nang];
nsph=Which[n2==-1,Gold[λ],n2==-2,Silver[λ],n2==-3,Copper[λ],n2==-4,Aluminum[λ],n2==-5,Drude[λ],n2==-6,LorentzDrudeAu[λ],n2==-7,LorentzDrudeAg[λ],n2==-8,LorentzDrudeCu[λ],n2==-9,LorentzDrudeAl[λ],n2==-10,Gold2[λ],n2==-11,Silver2[λ],n2==-12,LorentzDrudeBe[λ],n2==-13,LorentzDrudeCr[λ],n2==-14,LorentzDrudeNi[λ],n2==-15,LorentzDrudePd[λ],n2==-16,LorentzDrudePt[λ],n2==-17,LorentzDrudeTi[λ],n2==-18,LorentzDrudeW[λ],n2==-19,Iron2[λ],n2==-20,Cobalt2[λ],n2==-21,Rhodium[λ],n2==-22,Iridium[λ],True,n2];
m=nsph/n1;
x=2π r n1/λ;
y=m x;

nstop=Round[x+4 (x^(1/3))+2];
nmx=Round[Max[nstop,Abs[y]]]+15;
amu=Table[Cos[π(j-1)/(2(nang1-1))],{j,nang1}];
d[nmx]=0;
Table[rn=nmx-n+1;d[nmx-n]=rn/y - 1/(d[nmx-n+1] + rn/y),{n,nmx-1}];
Do[π0[j]=0;π1[j]=1,{j,nang1}];
Do[s1[j]=0;s2[j]=0,{j,2 nang1-1}];
ψ0=Cos[x];ψ1=Sin[x];ψ2=ψ1/x - ψ0;
χ0=-Sin[x];χ1=Cos[x];χ2=χ1/x - χ0;
ξ0=ψ0-I χ0;ξ1=ψ1-I χ1;ξ2=ψ2-I χ2;          (*spherical Hankel fns of 1st kind*)
h0=(ψ0+I χ0)/x;h1=(ψ1+I χ1)/x;h2=(ψ2+I χ2)/x;  (*Hankel fns of 2nd kind*)
Qsct=0;
Qnf=0;
Do[fn=(2 n+1)/(n(n+1));
ψ3=(2 (n+1)-1) ψ2/x - ψ1;
χ3=(2 (n+1)-1) χ2/x - χ1;
ξ3=ψ3-I χ3;
h3=(ψ3+I χ3)/x;
an=((d[n]/m + n/x)ψ2 - ψ1)/((d[n]/m + n/x)ξ2 - ξ1);
bn=((m d[n] + n/x)ψ2 - ψ1)/((m d[n] + n/x)ξ2 - ξ1);
Qsct=Qsct+(2 n+1)(Abs[an]^2+Abs[bn]^2);
Qnf=Qnf+2(Abs[an]^2 ((n+1) Abs[h1]^2 + n Abs[h3]^2) + (2n+1)Abs[bn]^2 Abs[h2]^2);
Do[jj=Round[2 nang1-j];
pi[j]=π1[j];
τ[j]=n amu[[j]] pi[j] - (n+1)π0[j];
p=(-1)^(n-1);
s1[j]=s1[j] + fn (an pi[j] + bn τ[j]);
t=(-1)^n;
s2[j]=s2[j]+fn(an τ[j]+ bn pi[j]);
If[j!=jj,s1[jj]=s1[jj] + fn (an pi[j] p + bn τ[j] t);
s2[jj]=s2[jj] + fn (an τ[j] t + bn pi[j] p)],{j,nang1}];
ψ0=ψ1;ψ1=ψ2;ψ2=ψ3;
χ0=χ1;χ1=χ2;χ2=χ3;
ξ0=ξ1;ξ1=ξ2;ξ2=ξ3;
h0=h1;h1=h2;h2=h3;
Do[π1[j]=((2 n+1)/n) amu[[j]] pi[j]-(n+1) π0[j]/n;
π0[j]=pi[j],{j,nang1}],{n,nstop}];
Qsct=2Qsct/x^2;
Qext=4 Re[s1[1]]/x^2;
Qbck=4(Abs[s1[2 nang1-1]]^2/x^2);
Return[{Qsct,Qext,Qbck,Qnf,s1,s2}]];

SphereScattering[λ1_,λ2_,λinc_,r1_,r2_,rinc_,n1_,n2_,θinc_]:=Module[{Qs,Qe,Qb,Qn,Qsct,Qbck,Qext,Qnf,p1,pp1,p2,pp2,p3,pp3,rinc1,nr,r,nang,s1,s2,s11norm,θ,s11,s12,s33,s34,perp,perp1,par,par1,pol,result,λinc1,nλ,λ},
If[λ2≠λ1,            (*  find wavelength dependence  *)
If[λ2>λ1,λinc1=Abs[λinc],λinc1=-Abs[λinc] ];
nλ=Round[Abs[(λ2-λ1)/λinc1]+1];
result=Table[λ=λ1+(j-1)λinc1;
mie[λ,r1,n1,n2,1 ],{j,nλ}];
Qs=Table[{λ1+(j-1)λinc1,result[[j,1]]},{j,nλ}];
Qe=Table[{λ1+(j-1)λinc1,result[[j,2]]},{j,nλ}];
Qb=Table[{λ1+(j-1)λinc1,result[[j,3]]},{j,nλ}];
Qn=Table[{λ1+(j-1)λinc1,result[[j,4]]},{j,nλ}];
Return[{Qs,Qe,Qb,Qn}],

If[r1≠r2,            (*  find radius dependence  *)
If[r2>r1,rinc1=Abs[rinc],rinc1=-Abs[rinc] ];
nr = Round[Abs[(r2-r1)/rinc1]+1];
result=Table[r=r1+(j-1)rinc1;
mie[λ1,r,n1,n2,1 ] ,{j,nr}];
Qs=Table[{r1+(j-1)rinc1,result[[j,1]]},{j,nr}];
Qe=Table[{r1+(j-1)rinc1,result[[j,2]]},{j,nr}];
Qb=Table[{r1+(j-1)rinc1,result[[j,3]]},{j,nr}];
Qn=Table[{r1+(j-1)rinc1,result[[j,4]]},{j,nr}];
Return[{Qs,Qe,Qb,Qn}],

nang = Round[Abs[90/θinc]+1];          (*  find angle dependence  *)
{Qsct,Qext,Qbck,Qnf,s1,s2}=mie[λ1,r1,n1,n2,nang ] ;
Print["λ = ",λ1," Qsct=",Qsct, " Qext=",Qext, " Qbck=",Qbck];
s11norm=(((Abs[s1[1]])^2)+ ((Abs[s2[1]])^2) )/2;
Do[θ=θinc(j-1);
s11=(Abs[s1[j]]^2+Abs[s2[j]]^2)/2;
s12=(Abs[s2[j]]^2-Abs[s1[j]]^2)/2;
perp[j]={θ,s11-s12};
perp1[j]={(s11-s12)Cos[π θ/180.0],(s11-s12)Sin[π θ/180.0]};
par[j]={θ,s11+s12};
par1[j]={(s11+s12)Cos[π θ/180.0],(s11+s12)Sin[π θ/180.0]};
pol=-s12/s11;
s33=Re[s2[j] Conjugate[s1[j]]]/s11;
s34=Im[s2[j] Conjugate[s1[j]]]/s11;
s11=s11/s11norm,{j,2 nang-1}] ;
p1=Table[perp[j],{j,2 nang-1}];   (* perp scattering in polar coords *)
pp1=Table[perp1[j],{j,2 nang-1}]; (* perp scattering in x-y coords  *)
p2=Table[par[j],{j,2 nang-1}];      (* par scattering in polar coords *)
pp2=Table[par1[j],{j,2 nang-1}];  (* par scattering in x-y coords  *)
p3=Table[0.5(perp[j]+par[j]),{j,2 nang-1}];      (* perp + par scattering in polar coords *)
pp3=Table[0.5(perp1[j]+par1[j]),{j,2 nang-1}];  (* perp + par scattering in x-y coords  *)
Return[{p1,pp1,p2,pp2,p3,pp3}]]]];

Code for generating figures

The routine for computinig the dielectric function is "ϵ."  It requires two input parameters, the wavelength λ and the metal index.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

λ1 = 0.4;      (* initial wavelength in microns *)
λ2 = 1.2;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavleneght  *)
num = Round[Abs[λ2 - λ1]/λinc] + 1;
ReAu=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-6]]},{i,num}];
ReAg=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-7]]},{i,num}];
ReAl=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-9]]},{i,num}];
ReCu=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-8]]},{i,num}];
Fig1a=ListLinePlot[{ReAu,ReAg,ReAl,ReCu}, LegendPosition -> {-0.45, -0.25}, PlotLegend->{"Au","Ag","Al","Cu"},LegendSize -> {0.4,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^′)"},PlotRange->All,Evaluate[option]];
Fig1aBW=ListLinePlot[{ReAu,ReAg,ReAl,ReCu}, LegendPosition -> {-0.45, -0.25}, PlotLegend->{"Au","Ag","Al","Cu"},LegendSize -> {0.4,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^′)"},PlotRange->All,Evaluate[optB]];
Clear[λ,λ1,λ2,λinc,num,ReAu,ReAg,ReAl,ReCu];

λ1 = 0.4;      (* initial wavelength in microns *)
λ2 = 1.2;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavleneght  *)
num = Round[Abs[λ2 - λ1]/λinc] + 1;
ImAu=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-6]]},{i,num}];
ImAg=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-7]]},{i,num}];
ImAl=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-9]]},{i,num}];
ImCu=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-8]]},{i,num}];
Fig1b=ListLinePlot[{ImAu,ImAg,ImAl,ImCu},  FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^″)"},Evaluate[option]] ;
Fig1bBW=ListLinePlot[{ImAu,ImAg,ImAl,ImCu},  FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^″)"},Evaluate[optB]] ;

Fig1=Framed[GraphicsGrid[{{Fig1a, Fig1b}},ImageSize->gr2,
Frame->False,Spacings -> {-30,0}],FrameMargins->1]
Fig1BW=Framed[GraphicsGrid[{{Fig1aBW, Fig1bBW}},ImageSize->gr2,
Frame->False,Spacings -> {-30,0}],FrameMargins->1];
Clear[λ,λ1,λ2,λinc,num,ImAu,ImAg,ImAl,ImCu];

The routine for computinig the dielectric function is "ϵ."  It requires two input parameters, the wavelength λ and the metal index.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for Au interpolated from values in Palik 
	-2: complex index for Ag interpolated from values in Palik  
	-3: complex index for Cu interpolated from values in Palik 
	-4: complex index for Al interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for Au from Drude-Lorentz fit from Rakik et al.
	-7: complex index for Ag from Drude-Lorentz fit from Rakik et al.
	-8: complex index for Cu from Drude-Lorentz fit from Rakik et al.
	-9: complex index for Al from Drude-Lorentz fit from Rakik et al.
	-10: complex index for Au from Johnson and Christy
	-11: complex index for Ag from Johnson and Christy
	-12: complex index for Be from Drude-Lorentz fit from Rakik et al.
	-13: complex index for Cr from Drude-Lorentz fit from Rakik et al.
	-14: complex index for Ni from Drude-Lorentz fit from Rakik et al.
	-15: complex index for Pd from Drude-Lorentz fit from Rakik et al.
	-16: complex index for Pt from Drude-Lorentz fit from Rakik et al.
	-17: complex index for Ti from Drude-Lorentz fit from Rakik et al.
	-18: complex index for W from Drude-Lorentz fit from Rakik et al.
	-19: complex index for Fe from Johnson and Christy
	-20: complex index for Co from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

λ1 = 0.4;
      
(* initial wavelength in microns *)
λ2 = 1.2;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavleneght  *)

num = Round[Abs[λ2 - λ1]/λinc] + 1;
RePd=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-15]]},{i,num}];
RePt=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-16]]},{i,num}];
ReRh=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-21]]},{i,num}];
ReIr=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-22]]},{i,num}];

Fig2a=ListLinePlot[{RePd,RePt,ReRh,ReIr}, LegendPosition -> {-.55, -.2}, PlotLegend->{"Pd","Pt","Rh","Ir"},LegendSize -> {0.3,0.5}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^′)"}, PlotRange->All,Evaluate[option]] ;
Fig2aBW=ListLinePlot[{RePd,RePt,ReRh,ReIr}, LegendPosition -> {-.55, -.2}, PlotLegend->{"Pd","Pt","Rh","Ir"},LegendSize -> {0.3,0.5}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^′)"}, PlotRange->All,Evaluate[optB]] ;
Clear[λ,λ1,λ2,λinc,num,RePd,RePt,ReRh,ReIr];

λ1 = 0.4;    (* initial wavelength in microns *)
λ2 = 1.2;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavelength  *)
num = Round[Abs[λ2 - λ1]/λinc] + 1;
ImPd=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-15]]},{i,num}];
ImPt=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-16]]},{i,num}];
ImRh=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-21]]},{i,num}];
ImIr=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-22]]},{i,num}];
Fig2b=ListLinePlot[{ImPd,ImPt,ImRh,ImIr},  FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^″)"},PlotRange->All,Evaluate[option]] ;
Fig2bBW=ListLinePlot[{ImPd,ImPt,ImRh,ImIr},  FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^″)"},PlotRange->All,Evaluate[optB]] ;

Fig2=Framed[GraphicsGrid[{{Fig2a, Fig2b}}, ImageSize -> gr2,
 Frame -> False,Spacings -> {-10,0}],FrameMargins->1]
Fig2BW=Framed[GraphicsGrid[{{Fig2aBW, Fig2bBW}}, ImageSize -> gr2,
 Frame -> False,Spacings -> {-10,0}],FrameMargins->1];
Clear[λ,λ1,λ2,λinc,num,ImPd,ImPt,ImRh,ImIr];

The routine for computing the dielectric function is "ϵ."  It requires two input parameters, the wavelength λ and the metal index.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy

λ1 = 0.4;      (* initial wavelength in microns *)
λ2 = 1.2;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavleneght  *)
num = Round[Abs[λ2 - λ1]/λinc] + 1;
ReFe=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-19]]},{i,num}];
ReCo=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-20]]},{i,num}];
ReCr=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-13]]},{i,num}];
ReW=Table[{λ=λ1+(i-1)λinc,Re[ϵ[λ,-18]]},{i,num}];
Fig3a=ListLinePlot[{ReFe,ReCo,ReCr,ReW}, LegendPosition -> {-0.5, -0.3}, PlotLegend->{"Fe","Co","Cr","W"},LegendSize -> {0.3,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^′)"}, Evaluate[option]] ;
Fig3aBW=ListLinePlot[{ReFe,ReCo,ReCr,ReW}, LegendPosition -> {-0.5, -0.3}, PlotLegend->{"Fe","Co","Cr","W"},LegendSize -> {0.3,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^′)"}, Evaluate[optB]] ;
Clear[λ,λ1,λ2,λinc,num,ReFe,ReCo,ReCr,ReW];

λ1 = 0.4;      (* initial wavelength in microns *)
λ2 = 1.2;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavleneght  *)
num = Round[Abs[λ2 - λ1]/λinc] + 1;
ImFe=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-19]]},{i,num}];
ImCo=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-20]]},{i,num}];
ImCr=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-13]]},{i,num}];
ImW=Table[{λ=λ1+(i-1)λinc,Im[ϵ[λ,-18]]},{i,num}];
Fig3b=ListLinePlot[{ImFe,ImCo,ImCr,ImW}, 
FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^″)"}, Evaluate[option]] ;
Fig3bBW=ListLinePlot[{ImFe,ImCo,ImCr,ImW}, 
FrameLabel -> {"λ (μm)", "(Subscript[ϵ, m]^″)"}, Evaluate[optB]] ;

Fig3=Framed[GraphicsGrid[{{Fig3a, Fig3b}}, ImageSize ->gr2,
 Frame -> False,Spacings -> {-40,0}],FrameMargins->1]
Fig3BW=Framed[GraphicsGrid[{{Fig3aBW, Fig3bBW}}, ImageSize ->gr2,
 Frame -> False,Spacings -> {-40,0}],FrameMargins->1];
Clear[λ,λ1,λ2,λinc,num,ImFe,ImCo,ImCr,ImW];

The observability function is "f."  It requires two input parameters, the wavelength λ and the metal index.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy

λ1 = 0.4;      (* initial wavelength in microns *)
λ2 = 1.2;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavelength  *)
num = Round[Abs[λ2 - λ1]/λinc] + 1;
Au=Table[{λ=λ1+(i-1)λinc,f[λ,-6]},{i,num}];
Ag=Table[{λ=λ1+(i-1)λinc,f[λ,-7]},{i,num}];
Ir=Table[{λ=λ1+(i-1)λinc,f[λ,-22]},{i,num}];
Pt=Table[{λ=λ1+(i-1)λinc,f[λ,-16]},{i,num}];
Fe=Table[{λ=λ1+(i-1)λinc,f[λ,-19]},{i,num}];
Cr=Table[{λ=λ1+(i-1)λinc,f[λ,-13]},{i,num}];
Fig4=ListLinePlot[{Au,Ag,Ir,Pt,Fe,Cr}, PlotRange->{0,0.03}, FrameLabel -> {"λ (μm)", "\!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)"}, Evaluate[option]] 
Fig4BW=ListLinePlot[{Au,Ag,Ir,Pt,Fe,Cr}, PlotRange->{0,0.03}, FrameLabel -> {"λ (μm)", "\!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)"}, Evaluate[optB]] ;
Clear[λ,λ1,λ2,λinc,num,Au,Ag,Ir,Pt,Fe,Cr];

The routine for computing reflectivity in the Kretschmann configuration requires six input parameters, the wavelength λ, the angle of incidence θ, the index of the incident medium (Subscript[n, 0]), the metal index (Subscript[n, 1]), the index of the transmitting medium (Subscript[n, 2]), and the thickness of the metal film d.  The metal index may be entered as a complex number, or the metal index can be directly calculated for any wavelength for several different metals according to which option is selected below for (Subscript[n, 1]).
	-1: complex index for Au interpolated from values in Palik 
	-2: complex index for Ag interpolated from values in Palik  
	-3: complex index for Cu interpolated from values in Palik 
	-4: complex index for Al interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for Au from Drude-Lorentz fit from Rakik et al.
	-7: complex index for Ag from Drude-Lorentz fit from Rakik et al.
	-8: complex index for Cu from Drude-Lorentz fit from Rakik et al.
	-9: complex index for Al from Drude-Lorentz fit from Rakik et al.
	-10: complex index for Au from Johnson and Christy
	-11: complex index for Ag from Johnson and Christy
	-12: complex index for Be from Drude-Lorentz fit from Rakik et al.
	-13: complex index for Cr from Drude-Lorentz fit from Rakik et al.
	-14: complex index for Ni from Drude-Lorentz fit from Rakik et al.
	-15: complex index for Pd from Drude-Lorentz fit from Rakik et al.
	-16: complex index for Pt from Drude-Lorentz fit from Rakik et al.
	-17: complex index for Ti from Drude-Lorentz fit from Rakik et al.
	-18: complex index for W from Drude-Lorentz fit from Rakik et al.
	-19: complex index for Fe from Johnson and Christy
	-20: complex index for Co from Johnson and Christy

θ1 = 40;      (* initial angle of incidence in microns *)
θ2 = 60;      (* final angle of incidence in microns *)
θinc=0.1;    (* increment of angle  *)
n0=1.5;         (*  index of incident medium - glass prism  *)
n1=-7;            (*  index of silver  *)
n2=1.0;          (*  index of transmitting medium  *)
d=0.050;        (*  thickness of metal film in microns  *)

num = Round[Abs[θ2 - θ1]/θinc] + 1;
λ=0.4;            (*  wavelength in microns  *)

Ag5=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.5;            (*  wavelength in microns  *)
Ag4=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.6;            (*  wavelength in microns  *)
Ag3=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.7;            (*  wavelength in microns  *)
Ag2=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.8;            (*  wavelength in microns  *)
Ag1=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}]; 

Fig5=ListLinePlot[{Ag1,Ag2,Ag3,Ag4,Ag5},FrameLabel -> {"θ (°)", "R"}, Evaluate[option]]
Fig5BW=ListLinePlot[{Ag1,Ag2,Ag3,Ag4,Ag5}, FrameLabel -> {"θ (°)", "R"}, Evaluate[optB]];
Clear[λ,θ,θ1,θ2,θinc,n0,n1,n2,d,num,Ag1,Ag2,Ag3,Ag4,Ag5];

θ1 = 40;      (* initial angle of incidence in microns *)
θ2 = 60;      (* final angle of incidence in microns *)
θinc=0.1;    (* increment of angle  *)
n0=1.5;         (*  index of incident medium - glass prism  *)
n1=-6;            (*  index of gold  *)
n2=1.0;          (*  index of transmitting medium  *)
d=0.050;        (*  thickness of metal film in microns  *)

num = Round[Abs[θ2 - θ1]/θinc] + 1;
λ=0.4;            (*  wavelength in microns  *)

Au5=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.5;            (*  wavelength in microns  *)
Au4=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.6;            (*  wavelength in microns  *)
Au3=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.7;            (*  wavelength in microns  *)
Au2=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.8;            (*  wavelength in microns  *)
Au1=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}]; 

Fig6=ListLinePlot[{Au1,Au2,Au3,Au4,Au5},  FrameLabel -> {"θ (°)", "R"}, Evaluate[option]]
Fig6BW=ListLinePlot[{Au1,Au2,Au3,Au4,Au5},  FrameLabel -> {"θ (°)", "R"}, Evaluate[optB]];
Clear[λ,θ,θ1,θ2,θinc,n0,n1,n2,d,num,Au1,Au2,Au3,Au4,Au5];

θ1 = 40;      (* initial angle of incidence in microns *)
θ2 = 60;      (* final angle of incidence in microns *)
θinc=0.1;    (* increment of angle  *)
n0=1.5;         (*  index of incident medium - glass prism  *)
n1=-9;            (*  index of aluminum  *)
n2=1.0;          (*  index of transmitting medium  *) 

num = Round[Abs[θ2 - θ1]/θinc] + 1;
λ=0.4;            (*  wavelength in microns  *)
d=0.015;        (*  thickness of metal film in microns  *)
Al5=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.5;            (*  wavelength in microns  *)
d=0.015;        (*  thickness of metal film in microns  *)
Al4=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.6;            (*  wavelength in microns  *)
d=0.011;        (*  thickness of metal film in microns  *)
Al3=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.7;            (*  wavelength in microns  *)
d=0.009;        (*  thickness of metal film in microns  *)
Al2=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
λ=0.8;            (*  wavelength in microns  *)
d=0.007;        (*  thickness of metal film in microns  *)
Al1=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}]; 

Fig7=ListLinePlot[{Al1,Al2,Al3,Al4,Al5}, FrameLabel -> {"θ (°)", "R"}, Evaluate[option]]
Fig7BW=ListLinePlot[{Al1,Al2,Al3,Al4,Al5}, FrameLabel -> {"θ (°)", "R"}, Evaluate[optB]];
Clear[λ,θ,θ1,θ2,θinc,n0,n1,n2,d,num,Al1,Al2,Al3,Al4,Al5];

λ=0.8;            (*  wavelength in microns  *)
θ1 = 40;      (* initial angle of incidence in microns *)
θ2 = 60;      (* final angle of incidence in microns *)
θinc=0.1;    (* increment of angle  *)
n0=1.5;         (*  index of incident medium - glass prism  *)
n2=1.0;          (*  index of transmitting medium  *)

num = Round[Abs[θ2 - θ1]/θinc] + 1;
n1=-19;            (*  index of iron  *)
d=0.010;        (*  thickness of metal film in microns  *)
Fe=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
n1=-20;            (*  index of cobalt  *)
d=0.010;        (*  thickness of metal film in microns  *)
Co=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
n1=-15;            (*  index of palladium  *)
d=0.010;        (*  thickness of metal film in microns  *)
Pd=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];
n1=-16;            (*  index of platinum  *)
d=0.010;        (*  thickness of metal film in microns  *)
Pt=Table[{θ=θ1+(i-1)θinc,r[λ,θ,n0,n1,n2,d]},{i,num}];

Fig8=ListLinePlot[{Fe,Co,Pd,Pt}, FrameLabel -> {"θ (°)", "R"}, Evaluate[option]]
Fig8BW=ListLinePlot[{Fe,Co,Pd,Pt}, FrameLabel -> {"θ (°)", "R"}, Evaluate[optB]];
Clear[λ,θ,θ1,θ2,θinc,n0,n1,n2,d,num,Fe,Co,Pd,Pt];

The routine for computing the field amplitude in the medium below the prism in the Kretschmann configuration requires seven input parameters, the wavelength λ, the angle of incidence θ, the index of the incident medium (Subscript[n, 0]), the metal index (Subscript[n, 1]), the index of the transmitting medium (Subscript[n, 2]), the thickness of the metal film d, and the distance below the surface of the metal film z.  The metal index may be entered as a complex number, or the metal index can be directly calculated for any wavelength for several different metals according to which option is selected below for (Subscript[n, 1]).
	-1: complex index for Au interpolated from values in Palik 
	-2: complex index for Ag interpolated from values in Palik  
	-3: complex index for Cu interpolated from values in Palik 
	-4: complex index for Al interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for Au from Drude-Lorentz fit from Rakik et al.
	-7: complex index for Ag from Drude-Lorentz fit from Rakik et al.
	-8: complex index for Cu from Drude-Lorentz fit from Rakik et al.
	-9: complex index for Al from Drude-Lorentz fit from Rakik et al.
	-10: complex index for Au from Johnson and Christy
	-11: complex index for Ag from Johnson and Christy
	-12: complex index for Be from Drude-Lorentz fit from Rakik et al.
	-13: complex index for Cr from Drude-Lorentz fit from Rakik et al.
	-14: complex index for Ni from Drude-Lorentz fit from Rakik et al.
	-15: complex index for Pd from Drude-Lorentz fit from Rakik et al.
	-16: complex index for Pt from Drude-Lorentz fit from Rakik et al.
	-17: complex index for Ti from Drude-Lorentz fit from Rakik et al.
	-18: complex index for W from Drude-Lorentz fit from Rakik et al.
	-19: complex index for Fe from Johnson and Christy
	-20: complex index for Co from Johnson and Christy

z = 0;      (* initial distance below metal film in microns *)
n0=1.5;         (*  index of incident medium - glass prism  *)
n2=1.0;          (*  index of transmitting medium - air  *)
λ1=0.3;            (*  initial wavelength in microns  *)
λ2=1.2;            (*  final wavelength in microns  *)
λinc = 0.01;    (*  increment of wavelength  *)

num = Round[Abs[λ2 -λ1]/λinc] + 1;
n1=-7;             (*  silver  *)
d=0.045;        (*  thickness of metal film in microns  *)
θ=42.7;         (* angle of incidence  *)
Ag=Table[{λ=λ1+(i-1)λinc,field[λ,θ,n0,n1,n2,d,z]},{i,num}];

n1=-6;             (*  gold  *)
d=0.045;        (*  thickness of metal film in microns  *)
θ=42.7;         (*  angle of incidence  *)
Au=Table[{λ=λ1+(i-1)λinc,field[λ,θ,n0,n1,n2,d,z]},{i,num}];

n1=-8;             (*  copper  *)
d=0.045;        (*  thickness of metal film in microns  *)
θ=42.7;         (* angle of incidence  *)
Cu=Table[{λ=λ1+(i-1)λinc,field[λ,θ,n0,n1,n2,d,z]},{i,num}];

n1=-9;             (*  aluminum  *)
d=0.018;        (*  thickness of metal film in microns  *)
θ=43.1;          (*  angle of incidence  *)
Al=Table[{λ=λ1+(i-1)λinc,field[λ,θ,n0,n1,n2,d,z]},{i,num}];

Fig9a=ListLinePlot[{Ag,Au,Cu,Al}, LegendPosition -> {0.8, -.05}, PlotLegend->{"Ag","Au","Cu","Al"},LegendSize -> {0.3,0.5}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "|E|"}, Evaluate[option]] ;
Fig9aBW=ListLinePlot[{Ag,Au,Cu,Al}, LegendPosition -> {0.8, -.05}, PlotLegend->{"Ag","Au","Cu","Al"},LegendSize -> {0.3,0.5}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "|E|"}, Evaluate[optB]] ;

n1=-20;           (*  cobalt  *)
d=0.015;        (*  thickness of metal film in microns  *)
θ=42.7;         (* angle of incidence  *)
Co=Table[{λ=λ1+(i-1)λinc,field[λ,θ,n0,n1,n2,d,z]},{i,num}];

n1=-13;           (*  chromium  *)
d=0.015;        (*  thickness of metal film in microns  *)
θ=42.7;         (* angle of incidence  *)
Cr=Table[{λ=λ1+(i-1)λinc,field[λ,θ,n0,n1,n2,d,z]},{i,num}];

n1=-15;            (*  palladium  *)
d=0.006;        (*  thickness of metal film in microns  *)
θ=42.7;         (* angle of incidence  *)
Pd=Table[{λ=λ1+(i-1)λinc,field[λ,θ,n0,n1,n2,d,z]},{i,num}];

n1=-16;            (*  platinum  *)
d=0.004;        (*  thickness of metal film in microns  *)
θ=42.7;         (* angle of incidence  *)
Pt=Table[{λ=λ1+(i-1)λinc,field[λ,θ,n0,n1,n2,d,z]},{i,num}];

Fig9b=ListLinePlot[{Co,Cr,Pd,Pt}, LegendPosition -> {0.85, -.05}, PlotLegend->{"Co","Cr","Pd","Pt"},LegendSize -> {0.3,0.5}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "|E|"}, Evaluate[option]] ;
Fig9bBW=ListLinePlot[{Co,Cr,Pd,Pt}, LegendPosition -> {0.85, -.05}, PlotLegend->{"Co","Cr","Pd","Pt"},LegendSize -> {0.3,0.5}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "|E|"}, Evaluate[optB]] ;

Fig9=Framed[GraphicsGrid[{{Fig9a, Fig9b}}, ImageSize ->gr2,
 Frame -> False,Spacings -> {10,0}],FrameMargins->1]
Fig9BW=Framed[GraphicsGrid[{{Fig9aBW, Fig9bBW}}, ImageSize ->gr2,
 Frame -> False,Spacings -> {10,0}],FrameMargins->1];
Clear[λ1,λ2,λinc,θ,θ1,θ2,θinc,n0,n1,n2,d,num,Ag,Au,Cu,Al,Co,Cr,Pd,Pt];

The routine for computing the propagation distance is "dist."  It requires three input parameters, the wavelength λ, the dielectric index, and the metal index.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

λ1 = 0.4;      (* initial wavelength in microns *)
λ2 = 1.0;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavleneght  *)

num = Round[Abs[λ2 - λ1]/λinc] + 1;
Au=Table[{λ=λ1+(i-1)λinc,Re[dist[λ,1.0,-6]]},{i,num}];
Ag=Table[{λ=λ1+(i-1)λinc,Re[dist[λ,1.0,-7]]},{i,num}];
Al=Table[{λ=λ1+(i-1)λinc,Re[dist[λ,1.0,-9]]},{i,num}];
Pt=Table[{λ=λ1+(i-1)λinc,Re[dist[λ,1.0,-16]]},{i,num}];
Fe=Table[{λ=λ1+(i-1)λinc,Re[dist[λ,1.0,-19]]},{i,num}];

Fig10=ListLinePlot[{Au,Ag,Al,Pt,Fe}, LegendPosition -> {-.6, 0.15}, PlotLegend->{"Au","Ag","Al","Pt","Fe"},LegendSize ->{0.3,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "d (μm)"}, Evaluate[option]] 
Fig10BW=ListLinePlot[{Au,Ag,Al,Pt,Fe}, LegendPosition -> {-.6, 0.15}, PlotLegend->{"Au","Ag","Al","Pt","Fe"},LegendSize -> {0.3,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "d (μm)"}, Evaluate[optB]] ;

Clear[λ,λ1,λ2,λinc,num,Au,Ag,Al,Cu];

The routine for computing the SP lifetime is (Subscript[τ, SP]).  It requires three input parameters, the wavelength λ, the dielectric index, and the metal index.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy
	-21: complex index for Rh interpolated from values in Palik 
	-22: complex index for Ir interpolated from values in Palik

λ1 = 0.4;      (* initial wavelength in microns *)
λ2 = 1.0;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavlength  *)

num = Round[Abs[λ2 - λ1]/λinc] + 1;
Au=Table[{λ=λ1+(i-1)λinc,Re[τsp[λ,1.0,-6]]},{i,num}];
Ag=Table[{λ=λ1+(i-1)λinc,Re[τsp[λ,1.0,-7]]},{i,num}];
Al=Table[{λ=λ1+(i-1)λinc,Re[τsp[λ,1.0,-9]]},{i,num}];
Pt=Table[{λ=λ1+(i-1)λinc,Re[τsp[λ,1.0,-16]]},{i,num}];
Fe=Table[{λ=λ1+(i-1)λinc,Re[τsp[λ,1.0,-19]]},{i,num}];

Fig11=ListLinePlot[{Au,Ag,Al,Pt,Fe}, LegendPosition -> {-.6, 0.15}, PlotLegend->{"Au","Ag","Al","Pt","Fe"},LegendSize -> {0.3,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "τ (fs)"}, Evaluate[option]] 
Fig11BW=ListLinePlot[{Au,Ag,Al,Pt,Fe}, LegendPosition -> {-.6, 0.15}, PlotLegend->{"Au","Ag","Al","Pt","Fe"},LegendSize -> {0.3,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "τ (fs)"}, Evaluate[optB]] ;
Clear[λ,λ1,λ2,λinc,num,Au,Ag,Al,Cu];

λ1=0.25;                  (*  initial wavelength  *)
λ2=0.7;                  (*  final wavelength  *)
λinc=0.001;                  (*  increment of wavelength  *)
θinc = 1;                   (*  increment of angle  *)   
r1=0.02;           (*  initial radius of sphere in microns  *)
r2=0.02;                  (*  final radius of sphere in microns  *)
rinc = 0.0;              (*  increment of radius  *)
n1 = 1.5;                   (* refractive index of glass medium  *)  

n2=-7;                   (*  refractive index of silver sphere  *){QsAg,QeAg,QbAg,QnAg}=SphereScattering[λ1,λ2,λinc,r1,r2,rinc,n1,n2,θinc];n2=-6;                   (*  refractive index of gold sphere  *)
{QsAu,QeAu,QbAu,QnAu}=SphereScattering[λ1,λ2,λinc,r1,r2,rinc,n1,n2,θinc];
n2=-8;                   (*  refractive index of copper sphere  *)
{QsCu,QeCu,QbCu,QnCu}=SphereScattering[λ1,λ2,λinc,r1,r2,rinc,n1,n2,θinc];
n2=-9;                   (*  refractive index of aluminum sphere  *)
{QsAl,QeAl,QbAl,QnAl}=SphereScattering[λ1,λ2,λinc,r1,r2,rinc,n1,n2,θinc];

Fig12a=ListLinePlot[{QeAg,QeAu,QeCu,QeAl},LegendPosition -> {.35, -.05}, PlotLegend->{"Ag","Au","Cu","Al"},LegendSize -> {0.3,0.6}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[Q, ext])"}, Evaluate[option]] ;
Fig12aBW=ListLinePlot[{QeAg,QeAu,QeCu,QeAl},LegendPosition -> {.35, -.05}, PlotLegend->{"Ag","Au","Cu","Al"},LegendSize -> {0.3,0.6}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[Q, ext])"}, Evaluate[optB]] ;

n2=-20;                   (*  refractive index of cobalt sphere  *){QsCo,QeCo,QbCo,QnCo}=SphereScattering[λ1,λ2,λinc,r1,r2,rinc,n1,n2,θinc];
n2=-13;                   (*  refractive index of chromium sphere  *){QsCr,QeCr,QbCr,QnCr}=SphereScattering[λ1,λ2,λinc,r1,r2,rinc,n1,n2,θinc];n2=-15;                   (*  refractive index of palladium sphere  *)
{QsPd,QePd,QbPd,QnPd}=SphereScattering[λ1,λ2,λinc,r1,r2,rinc,n1,n2,θinc];
n2=-16;                   (*  refractive index of platinum sphere  *)
{QsPt,QePt,QbPt,QnPt}=SphereScattering[λ1,λ2,λinc,r1,r2,rinc,n1,n2,θinc];
Fig12b=ListLinePlot[{QeCo,QeCr,QePd,QePt},LegendPosition -> {.3,-0.05}, PlotLegend->{"Co","Cr","Pd","Pt"},LegendSize -> {0.3,0.6}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[Q, ext])"}, Evaluate[option]] ;
Fig12bBW=ListLinePlot[{QeCo,QeCr,QePd,QePt},LegendPosition -> {.3,-0.05}, PlotLegend->{"Co","Cr","Pd","Pt"},LegendSize -> {0.3,0.6}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[Q, ext])"}, Evaluate[optB]] ;

Fig12=Framed[GraphicsGrid[{{Fig12a, Fig12b}},ImageSize -> gr2,Frame -> False, Spacings -> {-40,0}],FrameMargins->1]
Fig12BW=Framed[GraphicsGrid[{{Fig12aBW, Fig12bBW}},ImageSize -> gr2,Frame -> False, Spacings -> {-40,0}],FrameMargins->1];
Clear[λ1,λ2,λinc,r1,r2,rinc,n1,n2,θinc];

The observability function is "f."  It requires two input parameters, the wavelength λ and the metal index.  The metal index may be entered as a complex number, or the metal index can be directly calculated vs. wavelength for several different metals according to which option is selected below.
	-1: complex index for gold interpolated from values in Palik 
	-2: complex index for silver interpolated from values in Palik  
	-3: complex index for copper interpolated from values in Palik 
	-4: complex index for aluminum interpolated from values in Palik
	-5: complex index of a Drude metal
	-6: complex index for gold from Drude-Lorentz fit from Rakik et al.
	-7: complex index for silver from Drude-Lorentz fit from Rakik et al.
	-8: complex index for copper from Drude-Lorentz fit from Rakik et al.
	-9: complex index for aluminum from Drude-Lorentz fit from Rakik et al.
	-10: complex index for gold from Johnson and Christy
	-11: complex index for silver from Johnson and Christy
	-12: complex index for beryllium from Drude-Lorentz fit from Rakik et al.
	-13: complex index for chromium from Drude-Lorentz fit from Rakik et al.
	-14: complex index for nickel from Drude-Lorentz fit from Rakik et al.
	-15: complex index for palladium from Drude-Lorentz fit from Rakik et al.
	-16: complex index for platinum from Drude-Lorentz fit from Rakik et al.
	-17: complex index for titanium from Drude-Lorentz fit from Rakik et al.
	-18: complex index for tungsten from Drude-Lorentz fit from Rakik et al.
	-19: complex index for iron from Johnson and Christy
	-20: complex index for cobalt from Johnson and Christy

λ1 = 0.4;        (* initial wavelength in microns *)
λ2 = 1.2;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavelength  *)
num = Round[Abs[λ2 - λ1]/λinc] + 1;
Au=Table[{λ=λ1+(i-1)λinc,g[λ,-6]},{i,num}];
Ag=Table[{λ=λ1+(i-1)λinc,g[λ,-7]},{i,num}];
Al=Table[{λ=λ1+(i-1)λinc,g[λ,-9]},{i,num}];
Cu=Table[{λ=λ1+(i-1)λinc,g[λ,-8]},{i,num}];
Fig13a=ListLinePlot[{Au,Ag,Al,Cu}, LegendPosition -> {.8, .05}, PlotLegend->{"Au","Ag","Al","Cu"},LegendSize -> {0.3,0.5}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "\!\(\*
StyleBox[\"g\",\nFontSlant->\"Italic\"]\)"}, Evaluate[option]] ;
Fig13aBW=ListLinePlot[{Au,Ag,Al,Cu}, LegendPosition -> {.8, .05}, PlotLegend->{"Au","Ag","Al","Cu"},LegendSize -> {0.3,0.5}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "\!\(\*
StyleBox[\"g\",\nFontSlant->\"Italic\"]\)"}, Evaluate[optB]] ;
Clear[λ,λ1,λ2,λinc,num,Au,Ag,Al,Cu];

λ1 = 0.4;      (* initial wavelength in microns *)
λ2 = 1.2;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavelength  *)
num = Round[Abs[λ2 - λ1]/λinc] + 1;
Fe=Table[{λ=λ1+(i-1)λinc,g[λ,-19]},{i,num}];
Co=Table[{λ=λ1+(i-1)λinc,g[λ,-20]},{i,num}];
Cr=Table[{λ=λ1+(i-1)λinc,g[λ,-13]},{i,num}];
W=Table[{λ=λ1+(i-1)λinc,g[λ,-18]},{i,num}];
Fig13b=ListLinePlot[{Fe,Co,Cr,W}, LegendPosition -> {.85, .05}, PlotLegend->{"Fe","Co","Cr","W"},LegendSize -> {0.3,0.5}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "\!\(\*
StyleBox[\"g\",\nFontSlant->\"Italic\"]\)"}, Evaluate[option]] ;
Fig13bBW=ListLinePlot[{Fe,Co,Cr,W}, LegendPosition -> {.85, .05}, PlotLegend->{"Fe","Co","Cr","W"},LegendSize -> {0.3,0.5}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "\!\(\*
StyleBox[\"g\",\nFontSlant->\"Italic\"]\)"}, Evaluate[optB]] ;

Fig13=Framed[GraphicsGrid[{{Fig13a, Fig13b}}, ImageSize ->gr2,
 Frame -> False, Spacings -> {-10,0}],FrameMargins->1]
Fig13BW=Framed[GraphicsGrid[{{Fig13aBW, Fig13bBW}}, ImageSize ->gr2,
 Frame -> False, Spacings -> {-10,0}],FrameMargins->1];
Clear[λ,λ1,λ2,λinc,num,Fe,Co,Cr,W];

Γ=0;

Ω1 = 0.0;      (* initial normalized wavevector for lower branch *)
Ω2 = 1/Sqrt[2]-0.01;      (* final normalized wavelength *)
Ωinc=0.01;    (* increment of wavevector  *)
num = Round[Abs[Ω2 - Ω1]/Ωinc] + 1;
ReΩ1=Table[Ω=Ω1+(i-1)Ωinc;{Re[Ksp],Ω},{i,num}];
Ω1 = 1.0;      (* initial normalized wavevector for upper branch *)
Ω2 = 2.0;      (* final normalized wavelength *)
Ωinc=0.01;    (* increment of wavevector  *)
num = Round[Abs[Ω2 - Ω1]/Ωinc] + 1;
ReΩ2=Table[Ω=Ω1+(i-1)Ωinc;{Re[Ksp],Ω},{i,num}];

Fig15=ListLinePlot[{ReΩ1,ReΩ2}, LegendPosition -> {0, 0.1}, PlotLegend->{"Nonradiative","Radiative"},LegendSize -> {0.5,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"K = c (Subscript[k, sp])/(Subscript[ω, p])", "Ω = ω/(Subscript[ω, p])"}, Evaluate[option]] 
Fig15BW=ListLinePlot[{ReΩ1,ReΩ2}, LegendPosition -> {0, 0.1}, PlotLegend->{"Nonradiative","Radiative"},LegendSize -> {0.6,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"K = c (Subscript[k, sp])/(Subscript[ω, p])", "Ω = ω/(Subscript[ω, p])"}, Evaluate[optB]] ;
Clear[Ω,Ω1,Ω2,Ωinc,Γ,num,ReΩ];

The function for computing the dispersion of a lossy Drude metal with a real frequency is Ksp.  The parameters that get pased to the routine are the normalized frequency Ω and the normalized damping coefficient Γ.

Γ=0;
Ω1 = 0.0;      (* initial normalized wavevector *)
Ω2 = 1/Sqrt[2]-0.01;      (* final normalized wavelength *)
Ωinc=0.01;    (* increment of wavevector  *)
num = Round[Abs[Ω2 - Ω1]/Ωinc] + 1;
ReΩ1=Table[Ω=Ω1+(i-1)Ωinc;{Re[Ksp],Ω},{i,num}];
Frohlich={{0.0,1/Sqrt[2]},{10,1/Sqrt[2]}};
LightLine={{0.0,0.0},{10,10}};

Fig16=ListLinePlot[{ReΩ1,Frohlich,LightLine}, LegendPosition -> {-0.75, 0.25}, PlotLegend->{"Nonradiative","Fröhlich mode","Light Line"},LegendSize -> {0.6,0.3}, LegendShadow -> {.01, -.01}, PlotRange->{{0,3},{0,2}}, FrameLabel -> {"K = c (Subscript[k, sp])/(Subscript[ω, p])", "Ω = ω/(Subscript[ω, p])"}, Evaluate[option]]
Fig16BW=ListLinePlot[{ReΩ1,Frohlich,LightLine}, LegendPosition -> {-0.75, 0.25}, PlotLegend->{"Nonradiative","Fröhlich mode","Light Line"},LegendSize -> {0.65,0.3}, LegendShadow -> {.01, -.01}, PlotRange->{{0,3},{0,2}}, FrameLabel -> {"K = c (Subscript[k, sp])/(Subscript[ω, p])", "Ω = ω/(Subscript[ω, p])"}, Evaluate[optB]];
Clear[Ω,Ω1,Ω2,Ωinc,Γ,num,ReΩ];

Γ=0;
Ω1 = 1.0;      (* initial normalized wavevector *)
Ω2 = 2.0;      (* final normalized wavelength *)
Ωinc=0.01;    (* increment of wavevector  *)
num = Round[Abs[Ω2 - Ω1]/Ωinc] + 1;
ReΩ2=Table[Ω=Ω1+(i-1)Ωinc;{Re[Ksp],Ω},{i,num}];
Asymptote={{0.0,0.0},{10,10Sqrt[2]}};

Fig17=ListLinePlot[{ReΩ2,Asymptote}, LegendPosition -> {0.1, -0.2}, PlotLegend->{"Radiative","Asymptote"},LegendSize -> {0.5,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"K = c (Subscript[k, sp])/(Subscript[ω, p])", "Ω = ω/(Subscript[ω, p])"},  PlotRange->{{0,3},{0,2}},Evaluate[option]] 
Fig17BW=ListLinePlot[{ReΩ2,Asymptote}, LegendPosition -> {0.1, -0.2}, PlotLegend->{"Radiative","Asymptote"},LegendSize -> {0.5,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"K = c (Subscript[k, sp])/(Subscript[ω, p])", "Ω = ω/(Subscript[ω, p])"},  PlotRange->{{0,3},{0,2}},Evaluate[optB]] ;
Clear[Ω,Ω1,Ω2,Ωinc,Γ,num,ReΩ];

The function for computing the dispersion of a lossy Drude metal with a real frequency is Ksp.  The parameters that get passed to the routine are the normalized frequency Ω and the normalized damping coefficient Γ.

Ω1 = 0.0;      (* initial normalized wavevector *)
Ω2 = 2.0;      (* final normalized wavelength *)
Ωinc=0.01;    (* increment of wavevector  *)
num = Round[Abs[Ω2 - Ω1]/Ωinc] + 1;

Γ=0.01;
ReΩ1=Table[Ω=Ω1+(i-1)Ωinc;{Re[Ksp],Ω},{i,num}];
ImΩ1=Table[Ω=Ω1+(i-1)Ωinc;{Im[Ksp],Ω},{i,num}];
Γ=0.4;
ReΩ2=Table[Ω=Ω1+(i-1)Ωinc;{Re[Ksp],Ω},{i,num}];
ImΩ2=Table[Ω=Ω1+(i-1)Ωinc;{Im[Ksp],Ω},{i,num}];
Γ=0.8;
ReΩ3=Table[Ω=Ω1+(i-1)Ωinc;{Re[Ksp],Ω},{i,num}];
ImΩ3=Table[Ω=Ω1+(i-1)Ωinc;{Im[Ksp],Ω},{i,num}];
Γ=1.2;
ReΩ4=Table[Ω=Ω1+(i-1)Ωinc;{Re[Ksp],Ω},{i,num}];
ImΩ4=Table[Ω=Ω1+(i-1)Ωinc;{Im[Ksp],Ω},{i,num}];

Fig19a=ListLinePlot[{ReΩ1,ReΩ2,ReΩ3,ReΩ4}, LegendPosition -> {0.3, 0.05}, PlotLegend->{"0.01","0.4","0.8","1.2",},LegendSize -> 0.5, LegendShadow -> {.01, -.01}, FrameLabel -> {"K'", "Ω"}, PlotRange->{0,2},Evaluate[option]] ;
Fig19aBW=ListLinePlot[{ReΩ1,ReΩ2,ReΩ3,ReΩ4}, LegendPosition -> {0.3, 0.05}, PlotLegend->{"0.01","0.4","0.8","1.2",},LegendSize -> 0.5, LegendShadow -> {.01, -.01}, FrameLabel -> {"K'", "Ω"}, PlotRange->{0,2},Evaluate[optB]] ;

Fig19b=ListLinePlot[{ImΩ1,ImΩ2,ImΩ3,ImΩ4},  FrameLabel -> {"K''", "Ω"}, PlotRange->{0,2},Evaluate[option]] ;
Fig19bBW=ListLinePlot[{ImΩ1,ImΩ2,ImΩ3,ImΩ4},  FrameLabel -> {"K''", "Ω"}, PlotRange->{0,2},Evaluate[optB]] ;

Fig19=Framed[GraphicsGrid[{{Fig19a, Fig19b}}, ImageSize ->gr2,
 Frame -> False, Spacings -> {-10,0}],FrameMargins->1]
Fig19BW=Framed[GraphicsGrid[{{Fig19aBW, Fig19bBW}}, ImageSize ->gr2,
 Frame -> False, Spacings -> {-10,0}],FrameMargins->1];
Clear[Ω,Ω1,Ω2,Ωinc,Γ,num,ReΩ1,ReΩ2,ReΩ3,ReΩ4,ImΩ1,ImΩ2,ImΩ3,ImΩ4];

The function for computing the SP wavevector in a lossy Drude metal is Ksp.  The parameters that get passed to the routine are the normalized frequency Ω and the normalized damping coefficient Γ.

Ω1 = 0.01;      (* initial normalized wavevector *)
Ω2 = 2.0;      (* final normalized wavelength *)
Ωinc = 0.01;    (* increment of wavevector  *)
num = Round[Abs[Ω2 - Ω1]/Ωinc] + 1;

Γ = 0.01;
d1 = Table[Ω = Ω1 + (i - 1) Ωinc; {Ω,1/(2 π Im[Ksp]) }, {i, num}];
Γ = 0.4;
d2 = Table[Ω = Ω1 + (i - 1) Ωinc; {Ω,1/(2 π Im[Ksp])}, {i, num}];
Γ = 0.8;
d3 = Table[Ω = Ω1 + (i - 1) Ωinc; {Ω,1/(2 π Im[Ksp])}, {i, num}];
Γ = 1.2;
d4 = Table[Ω = Ω1 + (i - 1) Ωinc; {Ω,1/(2 π Im[Ksp])}, {i, num}];

Fig20=ListLinePlot[{d1, d2, d3,d4}, LegendPosition -> {0.9, -0.2}, LegendLabel -> "Γ", PlotLegend -> {"0.01", "0.4", "0.8", "1.2",}, LegendSize -> 0.5, LegendShadow -> {.01, -.01}, FrameLabel -> { "Ω","d/(Subscript[λ, p])"},PlotRange -> {0, 10}, Evaluate[option]] 
Fig20BW=ListLinePlot[{d1, d2, d3,d4}, LegendPosition -> {0.9, -0.2}, LegendLabel -> "Γ", PlotLegend -> {"0.01", "0.4", "0.8", "1.2",}, LegendSize -> 0.5, LegendShadow -> {.01, -.01}, FrameLabel -> { "Ω","d/(Subscript[λ, p])"},PlotRange -> {0, 10}, Evaluate[optB]] ;
Clear[Ω, Ω1, Ω2, Ωinc, Γ, num, d1, d2, d3, d4];

The routine for computing the dispersion of a lossy Drude metal with a real wavevector is Ω.  The parameters that get pased to the routine are the normalized wavevector K and the normalized damping coefficient Γ.

K1 = 0.0;      (* initial normalized wavevector *)
K2 = 3.0;      (* final normalized wavelength *)
Kinc=0.01;    (* increment of wavevector  *)
num = Round[Abs[K2 - K1]/Kinc] + 1;

Γ=0.0;
ReΩup1=Table[{K=K1+(i-1)Kinc,Re[Ωsp[[4]]]},{i,num}];
ReΩdn1=Table[{K=K1+(i-1)Kinc,Re[Ωsp[[3]]]},{i,num}];
ImΩup1=Table[{K=K1+(i-1)Kinc,Im[Ωsp[[4]]]},{i,num}];
ImΩdn1=Table[{K=K1+(i-1)Kinc,Im[Ωsp[[3]]]},{i,num}];
Γ=0.4;
ReΩup2=Table[{K=K1+(i-1)Kinc,Re[Ωsp[[4]]]},{i,num}];
ReΩdn2=Table[{K=K1+(i-1)Kinc,Re[Ωsp[[3]]]},{i,num}];
ImΩup2=Table[{K=K1+(i-1)Kinc,Im[Ωsp[[4]]]},{i,num}];
ImΩdn2=Table[{K=K1+(i-1)Kinc,Im[Ωsp[[3]]]},{i,num}];
Γ=0.8;
ReΩup3=Table[{K=K1+(i-1)Kinc,Re[Ωsp[[4]]]},{i,num}];
ReΩdn3=Table[{K=K1+(i-1)Kinc,Re[Ωsp[[3]]]},{i,num}];
ImΩup3=Table[{K=K1+(i-1)Kinc,Im[Ωsp[[4]]]},{i,num}];
ImΩdn3=Table[{K=K1+(i-1)Kinc,Im[Ωsp[[3]]]},{i,num}];
Γ=1.2;
ReΩup4=Table[{K=K1+(i-1)Kinc,Re[Ωsp[[4]]]},{i,num}];
ReΩdn4=Table[{K=K1+(i-1)Kinc,Re[Ωsp[[3]]]},{i,num}];
ImΩup4=Table[{K=K1+(i-1)Kinc,Im[Ωsp[[4]]]},{i,num}];
ImΩdn4=Table[{K=K1+(i-1)Kinc,Im[Ωsp[[3]]]},{i,num}];

Fig22a=ListLinePlot[{ReΩdn1,ReΩdn2,ReΩdn3,ReΩdn4},  FrameLabel -> {"K", "Ω'"},PlotRange->{0,2},Evaluate[option]] ;
Fig22aBW=ListLinePlot[{ReΩdn1,ReΩdn2,ReΩdn3,ReΩdn4},  FrameLabel -> {"K", "Ω'"},PlotRange->{0,2},Evaluate[optB]] ;

Fig22b=ListLinePlot[{ReΩup1,ReΩup2,ReΩup3,ReΩup4},  FrameLabel -> {"K", "Ω'"},PlotRange->{0,2},Evaluate[option]]; 
Fig22bBW=ListLinePlot[{ReΩup1,ReΩup2,ReΩup3,ReΩup4},  FrameLabel -> {"K", "Ω'"},PlotRange->{0,2},Evaluate[optB]] ;

Fig22=Framed[GraphicsGrid[{{Fig22a, Fig22b}}, ImageSize ->gr2,
 Frame -> False, Spacings -> {-10,0}],FrameMargins->1]
Fig22BW=Framed[GraphicsGrid[{{Fig22aBW, Fig22bBW}}, ImageSize ->gr2,
 Frame -> False, Spacings -> {-10,0}],FrameMargins->1];

Fig23a=ListLinePlot[{ImΩdn1,ImΩdn2,ImΩdn3,ImΩdn4},  FrameLabel -> {"K", "Ω''"}, Evaluate[option]] ;
Fig23aBW=ListLinePlot[{ImΩdn1,ImΩdn2,ImΩdn3,ImΩdn4},  FrameLabel -> {"K", "Ω''"}, Evaluate[optB]] ;

Fig23b=ListLinePlot[{ImΩup1,ImΩup2,ImΩup3,ImΩup4},FrameLabel -> {"K", "Ω''"}, Evaluate[option]];
Fig23bBW=ListLinePlot[{ImΩup1,ImΩup2,ImΩup3,ImΩup4},FrameLabel -> {"K", "Ω''"}, Evaluate[optB]];

Fig23=Framed[GraphicsGrid[{{Fig23a, Fig23b}}, ImageSize ->gr2,
 Frame -> False, Spacings -> {-10,0}],FrameMargins->1]
Fig23BW=Framed[GraphicsGrid[{{Fig23aBW, Fig23bBW}}, ImageSize ->gr2,
 Frame -> False, Spacings -> {-10,0}],FrameMargins->1];
Clear[K,K1,K2,Kinc,Γ,num,ReΩup1,ReΩdn1,ReΩup2,ReΩdn2,ReΩup3,ReΩdn3,ReΩup4,ReΩdn4,ImΩup1,ImΩup2,ImΩup3,ImΩup4,ImΩdn1,ImΩdn2,ImΩdn3,ImΩdn4];

λ1 = 0.22;      (* initial wavelength in microns *)
λ2 =0.8;       (* final wavelength in microns *)
λinc=0.01;    (* increment of wavelength  *)

num = Round[Abs[λ2 - λ1]/λinc] + 1;
Cr=Table[{λ=λ1+(i-1)λinc,g[λ,-13]},{i,num}];

ex1=ListLinePlot[{Cr}, FrameLabel -> {"λ (μm)", "g"}, Evaluate[option]]
ex1BW=ListLinePlot[{Cr}, FrameLabel -> {"λ (μm)", "g"}, Evaluate[optB]];
Clear[λ,λ1,λ2,λinc,num,Cr];

λ1 = 1;      (* initial wavelength in microns *)
λ2 = 10;      (* final wavelength in microns *)
λinc=0.01;    (* increment of wavelength  *)
num = Round[Abs[λ2 - λ1]/λinc] + 1;
W=Table[{λ=λ1+(i-1)λinc,g[λ,-18]},{i,num}];

ex2=ListLinePlot[{W}, FrameLabel -> {"λ (μm)", "g"}, Evaluate[option]]
ex2BW=ListLinePlot[{W}, FrameLabel -> {"λ (μm)", "g"}, Evaluate[optB]];
Clear[λ,λ1,λ2,λinc,num,W];

Figures

11.2  Real metals

11.2.1  Observability function

11.2.2  Kretschmann resonance

11.2.3  Propagation distance

11.2.4  Nanoparticle field enhancement

11.3  Drude metals

11.3.2  Lossless Drude model

Nonradiative SPs

Radiative SPs

11.3.3  Lossy Drude model

Exercises

1. In Fig. 11.11(b) the graph for chromium indicates an increasing field enhancement function at short wavelengths.  Plot g for chromium from 250 nm to 800 nm.  Does chromium become a plasmonic metal in the UV?

ex1

2. Tungsten is a highly resistive metal with a high melting point.  These are properties that make it useful for heating elements and lamp filaments, but not useful for supporting SPs at optical frequencies as seen in Fig. 11(b).  By plotting the field enhancement function from a wavelength of 1 μm to 10 μm, decide whether or not this metal could support SPs at these longer IR wavelengths.

ex2

References

[1] J. R. Sambles, G. W. Bradbery and F. Yang.  Optical excitation of surface plasmons: an introduction.  Contemp. Phys. 32 (1991) 173.
[2] O. Stephan, D. Taverna, M. Kociak, L. Henrard, K. Suenaga and C. Colliex.  Surface plasmon coupling in nanotubes.  In Structural and Electronic Properties of Molecular Nanostructures, ed. H. Kuzmany, J. Fink, M. Mehring and S. Roth.  (New York:  AIP, 2002), p. 477.
[3] D. J. Nash and J. R. Sambles.  Surface plasmon-polariton study of the optical dielectric function of silver.  J. Mod. Opt. 43 (1996) 81.
[4] R. A. Innes and J. R. Sambles.  Optical characterisation of gold using surface plasmon-polaritons.  J. Phys. F: Met. Phys. 17 (1987) 277.
[5] M. D. Tillin and J. R. Sambles.  A surface plasmon-polariton study of the dielectric constants of reactive metals: aluminum.  Thin Solid Films 167 (1988) 73.
[6] D. Y. Smith, E. Shiles and M. Inokuti.  The optical properties of metallic aluminum.  In Handbook of Optical Constants of Solids, ed. E. D. Palik.  (San Diego: Academic Press, 1998), pp. 369.
[7] M. van Exeter and A. Lagendijk.  Ultrashort surface-plasmon and phonon dynamics.  Phys. Rev. Lett. 60 (1988) 49.
[8] N. Marschall, B. Fischer and H. J. Queisser.  Dispersion of surface plasmons in InSb.  Phys. Rev. Lett. 27 (1971) 95.
[9] J. G. Rivas, M. Kuttge, P. H. Bolivar, H. Kurz and J. A. Sánchez-Gil.  Propagation of surface plasmon polaritons on semiconductor gratings.  Phys. Rev. Lett. 93 (2004) 256804.
[10] A. S. Barker, Jr.  Direct optical coupling to surface excitations.  Phys. Rev. Lett. 28 (1972) 892.
[11] A. Dereux, J-P. Vigneron, P. Lambin and A. A. Lucas.  Polaritons in semiconductor multilayered materials.  Phys. Rev. B 38 (1988) 5438.
[12] D. J. Nash and J. R. Sambles.  Surface plasmon-polariton study of the optical dielectric function of copper.  J. Mod. Opt. 42 (1995) 1639.
[13] T. Kloos and H. Raether.  The dispersion of surface plasmons of Al and Mg.  Phys. Lett. 44A (1973) 157.
[14] M. D. Tillin and J. R. Sambles.  Surface plasmon-polariton study of the dielectric function of magnesium.  Thin Solid Films 172 (1989) 27.
[15] G. J. Kovacs.  Optical excitation of surface plasma waves in an indium film bounded by dielectric layers.  Thin Solid Films 60 (1979) 33.
[16] M. G. Blaber, M. D. Arnold, N. Harris, M. J. Ford and M. B. Cortie.  Plasmon absorption in nanospheres:  a comparison of sodium, potassium, aluminium, silver and gold.  Physica B: Cond. Mat. 394 (2007) 184.
[17] H. J. Simon, D. E. Mitchell and J. G. Watson.  Second harmonic generation with surface plasmons in alkali metals.  Opt. Commun. 13 (1975) 294.
[18] Y. Xiong, J. Chen, B. Wiley, Y. Xia, Y. Yin and Z-Y. Li.  Size-dependence of surface plasmon resonance and oxidation for Pd nanocubes synthesized via a seed etching process.  Nano Letters 5 (2005) 1237.
[19] N. Zettsu, J. M. McLellan, B. Wiley, Y. Yin, Z-Y. Li and Y. Xia.  Synthesis, stability, and surface plasmonic properties of rhodium multipods, and their use as substrates for surface-enhanced Raman scattering.  Ang. Chem. - Intern. Ed. 45 (2006) 1288.
[20] F. Yang, G. W. Bradberry and J. R. Sambles.  Study of the optical properties of obliquely evaporated nickel films using IR surface plasmons.  Thin Solid Films 196 (1991) 35.
[21] R. W. Brown, P. Wessel and E. P. Trounson.  Plasmon radiation from silver films.  Phys Rev. Lett. 5 (1960) 472.
[22] R. A. Ferrell.  Predicted radiation of plasma oscillations in metal films.  Phys. Rev. 111 (1958) 1214.
[23] E. T. Arakawa, M. W. Williams, R. N. Hamm and R. H. Ritchie.  Effect of damping on surface plasmon dispersion.  Phys. Rev. Lett. 31 (1973) 1127.