8 
Quasi-one-dimensional surface plasmons

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;
base="C:\Users\Bill\Desktop\Surface Plasmons\Final figures\Ch 8 nanowires\jpg";

Packages

ϵ0=8.8541878176×10^-12;
μ0=4 π ×10^-7;
c=1/Sqrt[ϵ0 μ0];
η=Sqrt[μ0/ϵ0];
ϵ∞=4;                                     (*  background dielectric constant for Drude metal  *)
ωp=15.565;                         (*  ωp is plasma freq. in eV for Drude function  *)      
kp=2 π ωp/7.7905;            (* kp is plasma wavevector in ((μm)^(-1)) for Drude fn. *) 
γp = 0.608;                               (*  Drude damping coefficient in eV  *)

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={
LabelStyle->{FontSize -> fs},
GridLines -> Automatic,
Frame -> True,
ImageSize->is,
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[optC]]

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[optC]]

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}]

8.0.2  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 over certain wavelength ranges have been included for use in the 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 from the Handbook of Optical Constants of Solids are included for gold from the article by D. W. Lynch and W. R. Hunter, "Gold" pp. 293-295, for silver from the article by D. W. Lynch and W. R. Hunter, "Silver" pp. 355-357, for copper from the article by D. W. Lynch and W. R. Hunter, "Copper" pp. 284-285, for aluminum from the article by D. Y. Smith, E. Shiles and M. Inokuti, "The optical properties of metallic aluminum" pp. 397-400, for iridium form the article by D. W. Lynch and W. R. Hunter, "Iridium" pp. 302-303., and for rhodium from the article by D. W. Lynch and W. R. Hunter, "Rhodium" pp. 348-349.

GoldIndices = {{0.195,1.424,1.17},{0.2,1.427,1.215},{0.2066,1.422,1.306},{0.2101,1.43,1.334},{0.2138,1.432,1.364},{0.2175,1.438,1.388},{0.2254,1.452,1.442},{0.2296,1.454,1.478},{0.2339,1.462,1.51},{0.2384,1.47,1.55},{0.2431,1.478,1.59},{0.2480,1.484,1.636},{0.2530,1.49,1.698},{0.2583,1.504,1.748},{0.2638,1.546,1.784},{0.2695,1.598,1.822},{0.2755,1.648,1.852},{0.2818,1.69,1.882},{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}};

8.0.3  Routines to compute refractive indices of various metals

Drude[λ_]:=Module[{n,ω,ϵ},
ω = 7.7905/λ;                     (*  ang. frequency in eV when λ is in μm  *)
ϵ=ϵ∞-(ωp^2)/(ω^2 + I γp ω);
n = (Sqrt[ϵ]);
Return[n]];

(* from Rakic et al., Appl. Opt. vol. 37 (1998) 5271   *)
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   *)
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   *)
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   *)
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   *)
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   *)
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   *)
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   *)
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   *)
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   *)
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   *)
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,"Gold" in Handbook of Optical Constants of Solids,ed.E.D.Palik (Academic,San Diego,1998) 293-295.  However,the refractive index is modified so that Im(ϵ) is zero.  *)
Gold3[λ_]:=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[Sqrt[Re[n^2]]]];

(*  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]];

k0[λ_] := 2 π/λ;
nn[n_,λ_]=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[λ],n==-23,Gold3[λ],True,n];
ϵ[n_] := n^2;
γ[ϵ_, k0_, kz_] := Sqrt[kz^2-ϵ k0^2 ];

8.0.4  Normal modes of a planar film

err1[kz_,ϵ1_,ϵ2_,p_,h_,t_]:=Tan[h t]-2 p h/(ϵ2 h^2/ϵ1-ϵ1 p^2/ϵ2);

wv1[λ_,t_,nmid_,nout_,k1_]:=Module[{k,ϵ1,ϵ2,h,p,Δ},
ϵ1=ϵ[nn[nmid,λ]];
ϵ2=ϵ[nn[nout,λ]];
h[kz_]:=I γ[ϵ1, k0[λ], kz];
p[kz_]:=γ[ϵ2, k0[λ], kz];
Return[x/.FindRoot[err1[x,ϵ1,ϵ2,p[x],h[x],t],{x,k1}]]];

planarFilm[λ1_,λ2_,λinc_,x1_,x2_,xinc_,t_,nmid_,nout_]:=Module[{nλ,λ,λinc1,kz,wv,xinc1,x,nx,fields,Ex,Ez,sin0,cos0,sin1,cos1,h,p,c},
If[λ2≠λ1,            (*  find wavelength dependence of SP wavevector *)
If[λ2>λ1,λinc1=Abs[λinc],λinc1=-Abs[λinc] ];
nλ=Round[Abs[(λ2-λ1)/λinc1]+1];
kz=Table[λ=λ1+(j-1)λinc1;
wv=wv1[λ, t,  nn[nmid,λ], nn[nout,λ],1.1 2 π/λ+0.0 I ];
{{λ,wv},{λ,Re[wv]},{λ,Im[wv]},{λ,1/Im[wv]}},{j,nλ}];
Return[kz],

If[x2>x1,xinc1=Abs[xinc],xinc1=-Abs[xinc] ];
nx=Round[Abs[(x2-x1)/xinc1]+1];
kz=wv1[λ1, t, nn[nmid,λ1], nn[nout,λ1],1.1 2 π/λ1 ];
Print["kz=",kz,"  neff=",kz λ1/(2 π)];
sin0=kz/(nout k0[λ1]);
cos0=Sqrt[1-sin0^2];
sin1=kz/(nmid k0[λ1]);
cos1=Sqrt[1-sin1^2];
h=I γ[nmid^2,k0[λ1],kz];
p=-γ[nout^2, k0[λ1], kz];
c=-Cos[0.5 h t]/cos0;

fields=Table[x=x1+(j-1)xinc1;
If[x≤ 0.5 t,Ex=-I sin1 Sin[h x]/cos1 ;
Ez=Cos[h x],

Ex=+c sin0 Exp[p(x-0.5 t)];
Ez=-c cos0 Exp[p(x-0.5 t)]];
{{x,Abs[Ez]},{x,Abs[Ex]},{x,Sqrt[Abs[Ex]^2+Abs[Ez]^2]}},{j,nx}]];
Return[fields]];

8.0.5  Routines to compute refractive indices of various metals

err1[kz_,ϵ1_,ϵ2_,p_,h_,t_]:=Tan[h t]-2 p h/(ϵ2 h^2/ϵ1-ϵ1 p^2/ϵ2);

wv1[λ_,t_,nmid_,nout_,k1_]:=Module[{k,ϵ1,ϵ2,h,p,Δ},
ϵ1=ϵ[nn[nmid,λ]];
ϵ2=ϵ[nn[nout,λ]];
h[kz_]:=I γ[ϵ1, k0[λ], kz];
p[kz_]:=γ[ϵ2, k0[λ], kz];
Return[x/.FindRoot[err1[x,ϵ1,ϵ2,p[x],h[x],t],{x,k1}]]];

planarFilm[λ1_,λ2_,λinc_,x1_,x2_,xinc_,t_,nmid_,nout_]:=Module[{nλ,λ,λinc1,kz,wv,xinc1,x,nx,fields,Ex,Ez,sin0,cos0,sin1,cos1,h,p,c},
If[λ2≠λ1,            (*  find wavelength dependence of SP wavevector *)
If[λ2>λ1,λinc1=Abs[λinc],λinc1=-Abs[λinc] ];
nλ=Round[Abs[(λ2-λ1)/λinc1]+1];
kz=Table[λ=λ1+(j-1)λinc1;
wv=wv1[λ, t,  nn[nmid,λ], nn[nout,λ],1.1 2 π/λ+0.0 I ];
{{λ,wv},{λ,Re[wv]},{λ,Im[wv]},{λ,1/Im[wv]}},{j,nλ}];
Return[kz],

If[x2>x1,xinc1=Abs[xinc],xinc1=-Abs[xinc] ];
nx=Round[Abs[(x2-x1)/xinc1]+1];
kz=wv1[λ1, t, nn[nmid,λ1], nn[nout,λ1],1.1 2 π/λ1 ];
Print["kz=",kz,"  neff=",kz λ1/(2 π)];
sin0=kz/(nout k0[λ1]);
cos0=Sqrt[1-sin0^2];
sin1=kz/(nmid k0[λ1]);
cos1=Sqrt[1-sin1^2];
h=I γ[nmid^2,k0[λ1],kz];
p=-γ[nout^2, k0[λ1], kz];
c=-Cos[0.5 h t]/cos0;

fields=Table[x=x1+(j-1)xinc1;
If[x≤ 0.5 t,Ex=-I sin1 Sin[h x]/cos1 ;
Ez=Cos[h x],

Ex=+c sin0 Exp[p(x-0.5 t)];
Ez=-c cos0 Exp[p(x-0.5 t)]];
{{x,Abs[Ez]},{x,Abs[Ex]},{x,Sqrt[Abs[Ex]^2+Abs[Ez]^2]}},{j,nx}]];
Return[fields]];

8.0.6  Normal modes of a cylindrical wire or hollow waveguide (2 regions)

The value of (Subscript[k, z]) for which the function "err2" is zero is the wavevector for the SP propagating along the nanowire.  The function "wv2" uses err2 to solve for the (Subscript[k, z]) component of the wavevector.  The parameters passed to the function are the wavelength λ, the radius of the cylinder R, the refractive index inside the cylinder (Subscript[n, in]), the refractive index outside the cylinder (Subscript[n, out]), and a starting guess for the wavevector k1.  The value for (Subscript[k, z]) is returned.

The function "cylinder" tabulates the normal modes for a cylindrical wire or metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): initial radius of cylinder in microns
	(Subscript[R, 2]): final radius of cylinder in microns
	(Subscript[R, inc]): increment of cylinder radius
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for SP effective index

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes effective index of the SP vs. wavelength for a cylinder of radius (Subscript[R, 1]) using wavelength increment (Subscript[λ, inc]).  The function outputs six tables of values for plotting, {λ, (Subscript[k, z])/k}, {λ, Re((Subscript[k, z]))/(Subscript[k, 0])}, {λ, Im((Subscript[k, z]))/(Subscript[k, 0])}, {λ, 1/Im ((Subscript[k, z]))}, {Re((Subscript[k, z]))/(Subscript[k, p]), ω/(Subscript[ω, p])}, and {Im((Subscript[k, z]))/(Subscript[k, p]), ω/(Subscript[ω, p])}.  The last two tables are for use with a Drude metal.  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[R, 1]) ≠ (Subscript[R, 2]) then the function computes effective index of the SP  vs. radius of the cylinder at wavelength (Subscript[λ, 1]).  The function outputs four tables of values for plotting, {R, (Subscript[k, z])/(Subscript[k, 0])}, {R, Re((Subscript[k, z]))/(Subscript[k, 0])}, {R, Im((Subscript[k, z]))/(Subscript[k, 0])}, and {R, 1/Im((Subscript[k, z]))}.  If (Subscript[λ, 1])=(Subscript[λ, 2]) and (Subscript[R, 1]) = (Subscript[R, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of radius (Subscript[R, 1]) at wavelength (Subscript[λ, 1]) from (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The function outputs three tables of values for plotting, {r, | Ez |}, {r, |Eρ|}, and {r, Sqrt[|Eρ((|)^2) + |Ez((|)^2)]}.  The user can either set the refractive index to a specific (complex) value, or the user can set it 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 (approx. aluminum)
	-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

err2[kz_,ϵ1_,ϵ2_,k0_,R_]:=Module[{γ1,γ2},
γ1=γ[ϵ1,k0,kz];
γ2=γ[ϵ2,k0,kz];
If[Im[γ2]>0.0,γ2=-γ2];Return[(γ2/γ1)(BesselI[1,R γ1]/BesselI[0,R γ1])(BesselK[0,R γ2]/BesselK[1,R γ2])+(ϵ2/ϵ1)]];

wv2[λ_,R_,n1_,n2_,k1_]:=Module[{ϵ1,ϵ2,k,x},
k=k0[λ];
ϵ1=ϵ[n1];
ϵ2=ϵ[n2];
Return[x/.FindRoot[err2[x,ϵ1,ϵ2,k,R],{x,k1}]]];

cylinder[λ1_,λ2_,λinc_,r1_,r2_,rinc_,R1_,R2_,Rinc_,nin_,nout_,ng_]:=Module[{nλ,λ,λinc1,kz,wv,wv1,Rinc1,R,nR,rinc1,r,nr,fields,Ez,Eρ,γ1,γ2,c,k,ν,nn1,nn2},
If[λ2≠λ1,            (*  find wavelength dependence of SP wavevector *)
If[λ2>λ1,λinc1=Abs[λinc],λinc1=-Abs[λinc] ];
nλ=Round[Abs[(λ2-λ1)/λinc1]+1];
wv1=ng k0[λ1];
kz=Table[λ=λ1+(j-1)λinc1;
k=k0[λ];
If[Re[nin]<0,nn1=nn[nin,λ],nn1=nin];
If[Re[nout]<0,nn2=nn[nout,λ],nn2=nout];
wv=wv2[λ,R1,nn1 ,nn2,wv1 ];
wv1=wv;               (*  use this kz as best guess for next calculation  *)
ν = 1.24/λ;    (*  angular freq. in eV when λ is in microns  *)
{{λ,wv/k},{λ,Re[wv]/k},{λ,Im[wv]/k},{λ,1/Im[wv]},{Re[wv]/k,ν},{Im[wv]/k,ν}},{j,nλ}];
Return[kz],

If[R2≠R1,            (*  find cylinder radius dependence of SP wavevector *)
If[R2>R1,Rinc1=Abs[Rinc],Rinc1=-Abs[Rinc] ];
nR=Round[Abs[(R2-R1)/Rinc1]+1];
k=k0[λ1];
wv1=ng k;    (*  starting guess for SP refractive index  *)
kz=Table[R=R1+(j-1)Rinc1;
If[Re[nin]<0,nn1=nn[nin,λ1],nn1=nin];
If[Re[nout]<0,nn2=nn[nout,λ1],nn2=nout];
wv=wv2[λ1,R,nn1,nn2,wv1];
If[(j==1)||(j==nR),Print["kz(r=",R,") = ",wv,"   neff = ",wv/k]];
wv1=wv;   (*  use this kz as best guess for next calculation  *)
{{R,wv/k},{R,Re[wv]/k},{R,Im[wv]/k},{R,1/Im[wv]}},{j,nR}];
Return[kz],

If[r2>r1,rinc1=Abs[rinc],rinc1=-Abs[rinc] ];  (*  calculate fields vs. radial distance  *)
nr=Round[Abs[(r2-r1)/rinc1]+1];
k=k0[λ1];
If[Re[nin]<0,nn1=nn[nin,λ1],nn1=nin];
If[Re[nout]<0,nn2=nn[nout,λ1],nn2=nout];
kz=wv2[λ1,R1,nn1,nn2,ng k];
Print["kz=",kz,"  neff = ",kz/k];
γ1=γ[ϵ[nn1],k,kz];
γ2=γ[ϵ[nn2],k,kz];
(*If[Im[γ2]>0.0,γ2=-γ2];*)
c=BesselI[0,R1 γ1]/BesselK[0,R1 γ2];
fields=Table[r=r1+(j-1)rinc1;
If[r≤ R1,Ez=BesselI[0,r γ1];
Eρ=(I kz/γ1)BesselI[1,r γ1],

Ez=c BesselK[0,r γ2];
Eρ=(-I kz/γ2)c BesselK[1,r γ2]];
{{r,Abs[Ez]},{r,Abs[Eρ]},{r,Sqrt[Abs[Eρ]^2+Abs[Ez]^2]}},{j,nr}];
Return[fields]]]];

8.0.7  Normal modes of cylindrical shell (3 regions)

The value of (Subscript[k, z]) for which the function "err3" is zero is the wavevector for the SP propagating along the hollow nanotube.
	The function "wv3" solves for the (Subscript[k, z]) component of the wavevector.  The parameters passed to the function are the wavelength λ, the inner radius of the cylinder (Subscript[R, 1]), the outer radius of the cylinder (Subscript[R, 2]), the refractive index inside the hollow region of the cylinder (Subscript[n, 1]), the refractive index of the cylinder itself (Subscript[n, 2]), the refractive index outside the cylinder (Subscript[n, 3]), and a starting guess for the wavevector k1.  The value for (Subscript[k, z]) is returned.
	The function "tube" tabulates the normal modes for a hollow metallic cylinder.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): inner radius of cylinder in microns
	(Subscript[R, 2]): outer radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, 1]): refractive index inside hollow region of cylinder
	(Subscript[n, 2]): refractive index of cylinder
	(Subscript[n, 3]): refractive index outside cylinder
	(Subscript[n, g]): inritial guess for refractive index of SP

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes effective index of the SP vs. wavelength for a cylinder of inner radius (Subscript[R, 1]) and outer radius (Subscript[R, 2]) using wavelength increment (Subscript[λ, inc]).  The function outputs six tables of values for plotting, {λ, (Subscript[k, z])/k}, {λ, Re((Subscript[k, z]))/(Subscript[k, 0])}, {λ, Im((Subscript[k, z]))/(Subscript[k, 0])}, {λ, 1/Im((Subscript[k, z]))}, {Re((Subscript[k, z]))/(Subscript[k, p]), ω/(Subscript[ω, p])}, and {Im((Subscript[k, z]))/(Subscript[k, p]), ω/(Subscript[ω, p])}.  The last two tables are for use with a Drude metal.  If (Subscript[λ, 1])=(Subscript[λ, 2]) but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of inner radius (Subscript[R, 1]) and outer radius (Subscript[R, 2]) at wavelength (Subscript[λ, 1]) from (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The function outputs three tables of values for plotting, {r, | Ez |}, {r, |Eρ|}, and {r, Sqrt[|Eρ((|)^2) + |Ez((|)^2)]}.  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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
	-23: complex index for gold interpolated from values in Palik, but modified so that Im(ϵ) = 0

err3[kz_,ϵ1_,ϵ2_,ϵ3_,k0_,R1_,R2_,sign_]:=Module[{γ1,γ2,γ3,α1,α2,I011,I012,I022,I111,I112,I122,K012,K022,K023,K112,K122,K123},
γ1:=γ[ϵ1,k0,kz];
γ3:=γ[ϵ3,k0,kz];
If[sign≥0,
γ2:=+γ[ϵ2,k0,kz];
If[Im[γ3]<0.0,γ3=-γ3],
γ2:=-γ[ϵ2,k0,kz];
If[Im[γ3]>0.0,γ3=-γ3]];
α1:=(ϵ1 γ2)/(ϵ2 γ1);
α2:=(ϵ3 γ2)/(ϵ2 γ3);
I011:=BesselI[0,R1 γ1];
I012:=BesselI[0,R1 γ2];
I022:=BesselI[0,R2 γ2];
I111:=BesselI[1,R1 γ1];
I112:=BesselI[1,R1 γ2];
I122:=BesselI[1,R2 γ2];
K012:=BesselK[0,R1 γ2];
K022:=BesselK[0,R2 γ2];
K023:=BesselK[0,R2 γ3];
K112:=BesselK[1,R1 γ2];
K122:=BesselK[1,R2 γ2];
K123:=BesselK[1,R2 γ3];
Return[I011 K023(I122 K112-I112 K122)+α1 I111 K023(I122 K012+I012 K122)+α2 I011 K123 (I022 K112+I112 K022)+α1 α2 I111 K123(I022 K012-I012 K022)]];

wv3[k_,R1_,R2_,n1_,n2_,n3_,k1_,sign_]:=Module[{x},
Return[x/.FindRoot[err3[x,n1^2,n2^2,n3^2,k,R1,R2,sign],{x,k1}]]];

tube[λ1_,λ2_,λinc_,r1_,r2_,rinc_,R1_,R2_,n1_,n2_,n3_,ng_,sign_]:=Module[{nλ,λ,λinc1,ν,kz,wv,rinc1,r,nr,fields,Ez,Eρ,Hϕ,kz1,nn1,nn2,nn3,ϵ1,ϵ2,ϵ3,γ1,γ2,γ3,α1,α2,b,c,d,k,wv1,I011,I012,I022,I111,I112,I122,K011,K012,K022,K023,K112,K122,K123,error},
If[λ2≠λ1,            (*  find wavelength dependence of SP wavevector *)
If[λ2>λ1,λinc1=Abs[λinc],λinc1=-Abs[λinc] ];
nλ=Round[Abs[(λ2-λ1)/λinc1]+1];
wv1=ng k0[λ1];
kz=Table[λ=λ1+(j-1)λinc1;
If[Re[n1]<0,nn1=nn[n1,λ],nn1=n1];
If[Re[n2]<0,nn2=nn[n2,λ],nn2=n2];
If[Re[n3]<0,nn3=nn[n3,λ],nn3=n3];
k=k0[λ];
ν = 1.24/λ;    (*  freq. in eV when λ is in microns  *)
wv=wv3[k, R1,R2,nn1, nn2,nn3,wv1,sign];
error=err3[wv,nn1^2,nn2^2,nn3^2,k,R1,R2,sign];
If[error>0.0001,Print["λ=",λ,"  err=",error];wv=0];
wv1=wv;
{{λ,wv/k},{λ,Re[wv]/k},{λ,Im[wv]/k},{λ,1/Im[wv]},{Re[wv]/k,ν},{Im[wv]/k,ν}},{j,nλ}];
Return[kz],

If[r2>r1,rinc1=Abs[rinc],rinc1=-Abs[rinc] ];  (*  calculate fields vs. radial distance  *)
nr=Round[Abs[(r2-r1)/rinc1]+1];
If[Re[n1]<0,nn1=nn[n1,λ1],nn1=n1];
If[Re[n2]<0,nn2=nn[n2,λ1],nn2=n2];
If[Re[n3]<0,nn3=nn[n3,λ1],nn3=n3];
k=k0[λ1];
ϵ1=ϵ[nn1];
ϵ2=ϵ[nn2];
ϵ3=ϵ[nn3];
kz=wv3[k,R1,R2,nn1,nn2,nn3,ng k,sign]; 
Print["kz=",kz,"  neff=",kz/k];
γ1:=γ[ϵ1,k,kz];
γ3:=γ[ϵ3,k,kz]; 
If[sign≥0,
γ2:=+γ[ϵ2,k,kz];
If[Im[γ3]<0.0,γ3=-γ3],
γ2:=-γ[ϵ2,k,kz];
If[Im[γ3]>0.0,γ3=-γ3]];
α1:=(ϵ1 γ2)/(ϵ2 γ1);
α2:=(ϵ3 γ2)/(ϵ2 γ3);
I011:=BesselI[0,R1 γ1];
I012:=BesselI[0,R1 γ2];
I022:=BesselI[0,R2 γ2];
I111:=BesselI[1,R1 γ1];
I112:=BesselI[1,R1 γ2];
I122:=BesselI[1,R2 γ2];
K012:=BesselK[0,R1 γ2];
K022:=BesselK[0,R2 γ2];
K023:=BesselK[0,R2 γ3];
K112:=BesselK[1,R1 γ2];
K122:=BesselK[1,R2 γ2];
K123:=BesselK[1,R2 γ3];

b=R1 γ2(α1 I111 K012+I011 K112);
c=R1 γ2(I011 I112-α1 I012 I111);
d=(I022 b+K022 c)/K023;
fields=Table[r=r1+(j-1)rinc1;
If[r≤ R1,           (*  fields inside hollow of cylinder  *)
Ez=BesselI[0,r γ1];
Eρ=(I kz/γ1)BesselI[1,r γ1];
Hϕ=(k ϵ1/(kz η)) Eρ,

If[r≤ R2,          (*  fields inside cylinder itself  *)
Ez=b BesselI[0,r γ2]+c BesselK[0,r γ2];
Eρ=(I kz/γ2)(b BesselI[1,r γ2]-c BesselK[1,r γ2]);
Hϕ=(k ϵ2/(kz η)) Eρ,

Ez=d BesselK[0,r γ3];    (*  fields outside cylinder  *)
Eρ=(-I kz/γ3)d BesselK[1,r γ3]]];
Hϕ=(k ϵ3/(kz η)) Eρ;
{{r,Re[Ez]},{r,Im[Eρ]},{r,Sqrt[Abs[Eρ]^2+Abs[Ez]^2]}},{j,nr}];
Return[fields]]];

8.0.8  Plane wave incident on metallic cylinder

The function "local" computes the local fields around a metallic cylinder when illuminated by a plane wave such that the cylinder axis is normal to the plane of incidence and the incident field is polarized in the plane of incidence at right angles to the cylinder axis.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): initial radius of cylinder in microns
	(Subscript[R, 2]): final radius of cylinder in microns
	(Subscript[R, inc]): increment of cyinder radius
	(Subscript[r, 1]): initial radius in μm for calculation of local fields
	(Subscript[r, 2]): final radius in μm for calculation of local fields
	ϕ: angle in degrees for field calculation in xy plane measured from the x axis.  0 degrees is the backward scattered direction.
	(Subscript[x, 1]): initial x axis point in μm for calculation of local fields
	(Subscript[x, 2]): final x axis point in μm for calculation of local fields
	(Subscript[y, 1]): initial x axis point in μm for calculation of local fields
	(Subscript[y, 2]): final x axis point in μm for calculation of local fields
	inc: increment in μm of radius, x position, or y position for calculation of local fields
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes peak field intensity at the surface of the cylinder and scattering efficiency vs. wavelength for a cylinder of radius (Subscript[R, 1]) using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]) or (Subscript[x, 1])≠(Subscript[x, 2]) or (Subscript[y, 1])≠(Subscript[y, 2]), then the function computes fields at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) at angle ϕ or along the x axis from (Subscript[x, 1]) to (Subscript[x, 2]) or along the y axis from (Subscript[y, 1]) to (Subscript[y, 2]) by increment inc.  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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

local[n1_,n2_,R_,nmax_,λ1_,λ2_,λinc_,r1_,r2_,ϕ_,x1_,x2_,y1_,y2_,inc_]:=Module[{ξ1,ξ2,J1,J1p,J2,J2p,H2,H2p,an,gn,num,λ,r,ϕ1,x,y,t,inc1,En,Eρ,Eϕ,Eiρ,Eiϕ,Esρ,Esϕ,J1neg,J1pneg,J2neg,J2pneg,H2neg,H2pneg,nn1,nn2,Imax,intensity,ϕmax,Qsct},
If[λ1≠λ2,num = Round[Abs[(λ2-λ1)/λinc]] + 1;
Table[λ=λ1+(i-1)λinc;
If[Re[n1]<0,nn1=nn[n1,λ],nn1=n1];
If[Re[n2]<0,nn2=nn[n2,λ],nn2=n2];
ξ1=nn1 k0[λ]R ;
ξ2=nn2 k0[λ] R ;
J1=Table[BesselJ[n,ξ1],{n,0,nmax}];     (*  Note that (Subscript[J, n-1]) is J1[[n]].  *)
J1p=Table[If[n==0,-J1[[2]],J1[[n]]-n J1[[n+1]]/ξ1],{n,0,nmax}];   (*  derivative  *)
J1neg=Table[((-1)^n)J1[[n+1]],{n,0,nmax}];     (*  (Subscript[J, -n])(ξ) = (((-1))^n) (Subscript[J, +n])(ξ)  *)
J1pneg=Table[((-1)^n)J1p[[n+1]],{n,0,nmax}];     J2=Table[BesselJ[n,ξ2],{n,0,nmax}];
J2p=Table[If[n==0,-J2[[2]],J2[[n]]-n J2[[n+1]]/ξ2],{n,0,nmax}];
J2neg=Table[((-1)^n)J2[[n+1]],{n,0,nmax}];     (*  (Subscript[J, -n])(ξ) = (((-1))^n) (Subscript[J, +n])(ξ)  *)
J2pneg=Table[((-1)^n)J2p[[n+1]],{n,0,nmax}];     H2=Table[HankelH1[n,ξ2],{n,0,nmax}];  (*  Note that (Subscript[H, n-1]^((1))) is H1[[n]].  *)
H2p=Table[If[n==0,-H2[[2]],H2[[n]]-n H2[[n+1]]/ξ2],{n,0,nmax}];
H2neg=Table[((-1)^n)H2[[n+1]],{n,0,nmax}];     (*  (Subscript[H, -n])(ξ) = (((-1))^n) (Subscript[H, +n])(ξ)  *)
H2pneg=Table[((-1)^n)H2p[[n+1]],{n,0,nmax}];     an=Table[I(ξ1 J1[[n+1]] J2p[[n+1]]-ξ2 J2[[n+1]] J1p[[n+1]])/(ξ1 J1[[n+1]] H2p[[n+1]]-ξ2 J1p[[n+1]] H2[[n+1]]),{n,0,nmax}];
gn=Table[(I ξ2^2/ξ1)(J2p[[n+1]] H2[[n+1]]-J2[[n+1]] H2p[[n+1]])/(ξ1 J1[[n+1]] H2p[[n+1]]-ξ2 J1p[[n+1]] H2[[n+1]]),{n,0,nmax}];

Qsct=Abs[an[[1]]]^2;                (*  compute scattering efficiency  *)
Do[Qsct=Qsct+2 Abs[an[[i+1]]]^2,{i,1,nmax}];
Qsct=Qsct/(π R nn2 k0[λ]);

Imax=0;
ϕmax=-10;
r=R;
Do[ϕ1=j;   (*  compute |E((|)^2) along entire top surface of cylinder  *)
Eρ=0;  (*  field components on outside surface of the cylinder  *)
Eϕ=0;
Eiρ=0;
Eiϕ=0;
Esρ=0;
Esϕ=0;
Do[En=((-I)^n) r Exp[I n 2 π ϕ1/360]/ξ2;
Eiρ=Eiρ+En n J2[[n+1]]/r;            
Eiϕ=Eiϕ+En I ξ2 J2p[[n+1]]/r;
Esρ=Esρ+En (an[[n+1]] I n/r)H2[[n+1]];
Esϕ=Esϕ-En(an[[n+1]] ξ2/r)H2p[[n+1]];
If[n≠0,En=((+I)^n) r Exp[-I n 2 π ϕ1/360]/ξ2;
Eiρ=Eiρ-En n J2neg[[n+1]]/r;
Eiϕ=Eiϕ+En I ξ2 J2pneg[[n+1]]/r;
Esρ=Esρ-En (an[[n+1]] I n/r)H2neg[[n+1]];
Esϕ=Esϕ-En(an[[n+1]] ξ2/r)H2pneg[[n+1]];
Eρ=Eiρ+Esρ;
Eϕ=Eiϕ+Esϕ],{n,0,nmax}];
intensity=Abs[Eρ]^2+Abs[Eϕ]^2;
If[intensity>Imax,Imax=intensity;ϕmax=ϕ1],{j,1,180}];
{{λ,Imax},{λ,Qsct}},{i,num}],

If[Re[n1]<0,nn1=nn[n1,λ1],nn1=n1];    (*  compute fields along a line  *)
If[Re[n2]<0,nn2=nn[n2,λ1],nn2=n2];
ξ1=nn1 k0[λ1]R ; 
ξ2=nn2 k0[λ1] R ;
J1=Table[BesselJ[n,ξ1],{n,0,nmax}];     (*  Note that (Subscript[J, n-1]) is J1[[n]].  *)
J1p=Table[If[n==0,-J1[[2]],J1[[n]]-n J1[[n+1]]/ξ1],{n,0,nmax}];   (*  derivative  *)
J2=Table[BesselJ[n,ξ2],{n,0,nmax}];
J2p=Table[If[n==0,-J2[[2]],J2[[n]]-n J2[[n+1]]/ξ2],{n,0,nmax}];
H2=Table[HankelH1[n,ξ2],{n,0,nmax}];  (*  Note that (Subscript[H, n-1]^((1))) is H1[[n]].  *)
H2p=Table[If[n==0,-H2[[2]],H2[[n]]-n H2[[n+1]]/ξ2],{n,0,nmax}];
an=Table[I(ξ1 J1[[n+1]] J2p[[n+1]]-ξ2 J2[[n+1]] J1p[[n+1]])/(ξ1 J1[[n+1]] H2p[[n+1]]-ξ2 J1p[[n+1]] H2[[n+1]]),{n,0,nmax}];
gn=Table[(I ξ2^2/ξ1)(J2p[[n+1]] H2[[n+1]]-J2[[n+1]] H2p[[n+1]])/(ξ1 J1[[n+1]] H2p[[n+1]]-ξ2 J1p[[n+1]] H2[[n+1]]),{n,0,nmax}];

(*Print["J1:  ",J1];
Print["J1p: ",J1p];
Print["J2:  ",J2];
Print["J2p: ",J2p];
Print["H2:  ",H2];
Print["H2p: ",H2p];
Print["an:  ",an];
Print["gn:  ",gn];*)

If[r1≠r2,num = Round[Abs[(r2-r1)/inc]] + 1;
If[r2>r1,inc1=Abs[inc],inc1=-Abs[inc]],
If[x1≠x2,num = Round[Abs[(x2-x1)/inc]] + 1;
If[x2>x1,inc1=Abs[inc],inc1=-Abs[inc]],
num = Round[Abs[(y2-y1)/inc]] + 1;
If[y2>y1,inc1=Abs[inc],inc1=-Abs[inc]]]];

Table[If[r2≠r1,r=r1+(i-1) inc1;If[r==0,r=0.000001];t=r;ϕ1=ϕ,
If[x2≠x1,x=x1+(i-1) inc1;If[x==0,x=0.000001];t=x;
r=Sqrt[x^2+y1^2];ϕ1=(180/π)ArcTan[x,y1],
y=y1+(i-1) inc1;If[y==0,y=0.000001];t=y;
r=Sqrt[x1^2+y^2];ϕ1=(180/π)ArcTan[x1,y]]];
ξ1=nn1 k0[λ1] r ;
ξ2=nn2 k0[λ1] r ;
J1=Table[BesselJ[n,ξ1],{n,0,nmax}];     (*  Note that (Subscript[J, n-1]) is J1[[n]].  *)
J1p=Table[If[n==0,-J1[[2]],J1[[n]]-n J1[[n+1]]/ξ1],{n,0,nmax}];   (*  derivative  *)
J1neg=Table[((-1)^n)J1[[n+1]],{n,0,nmax}];     (*  (Subscript[J, -n])(ξ) = (((-1))^n) (Subscript[J, +n])(ξ)  *)
J1pneg=Table[((-1)^n)J1p[[n+1]],{n,0,nmax}];     (*  (Subscript[J, -n]^′)(ξ) = (((-1))^n) (Subscript[J, +n]^′)(ξ)  *)
J2=Table[BesselJ[n,ξ2],{n,0,nmax}];
J2p=Table[If[n==0,-J2[[2]],J2[[n]]-n J2[[n+1]]/ξ2],{n,0,nmax}];
J2neg=Table[((-1)^n)J2[[n+1]],{n,0,nmax}];     (*  (Subscript[J, -n])(ξ) = (((-1))^n) (Subscript[J, +n])(ξ)  *)
J2pneg=Table[((-1)^n)J2p[[n+1]],{n,0,nmax}];     (*  (Subscript[J, -n]^′)(ξ) = (((-1))^n) (Subscript[J, +n]^′)(ξ)  *)
H2=Table[HankelH1[n,ξ2],{n,0,nmax}];  (*  Note that (Subscript[H, n-1]^((1))) is H1[[n]].  *)
H2p=Table[If[n==0,-H2[[2]],H2[[n]]-n H2[[n+1]]/ξ2],{n,0,nmax}];
H2neg=Table[((-1)^n)H2[[n+1]],{n,0,nmax}];     (*  (Subscript[H, -n])(ξ) = (((-1))^n) (Subscript[H, +n])(ξ)  *)
H2pneg=Table[((-1)^n)H2p[[n+1]],{n,0,nmax}];     (*  (Subscript[H, -n]^′)(ξ) = (((-1))^n) (Subscript[H, +n]^′)(ξ)  *)

Eρ=0;
Eϕ=0;
Eiρ=0;
Eiϕ=0;
Esρ=0;
Esϕ=0;
Do[If[r<R,
En=((-I)^n) r Exp[I n 2 π ϕ1/360]/ξ2;
Eρ=Eρ+En(gn[[n+1]] I n/r)J1[[n+1]]; (*  field components inside the cylinder  *)
Eϕ=Eϕ-En(gn[[n+1]] ξ1/r)J1p[[n+1]];
If[n≠0,En=((+I)^n) r Exp[-I n 2 π ϕ1/360]/ξ2;
 Eρ=Eρ-En(gn[[n+1]] I n/r)J1neg[[n+1]];
Eϕ=Eϕ-En(gn[[n+1]] ξ1/r)J1pneg[[n+1]]],

En=((-I)^n) r Exp[I n 2 π ϕ1/360]/ξ2;
Eiρ=Eiρ+En n J2[[n+1]]/r;            (*  field components outside the cylinder  *)
Eiϕ=Eiϕ+En I ξ2 J2p[[n+1]]/r;
Esρ=Esρ+En (an[[n+1]] I n/r)H2[[n+1]];
Esϕ=Esϕ-En(an[[n+1]] ξ2/r)H2p[[n+1]];
If[n≠0,En=((+I)^n) r Exp[-I n 2 π ϕ1/360]/ξ2;
Eiρ=Eiρ-En n J2neg[[n+1]]/r;
Eiϕ=Eiϕ+En I ξ2 J2pneg[[n+1]]/r;
Esρ=Esρ-En (an[[n+1]] I n/r)H2neg[[n+1]];
Esϕ=Esϕ-En(an[[n+1]] ξ2/r)H2pneg[[n+1]]];
Eρ=Eiρ+Esρ;
Eϕ=Eiϕ+Esϕ],{n,0,nmax}];
{{t,Sqrt[Abs[Eρ]^2+Abs[Eϕ]^2]},{t,Abs[Eρ]},{t,Abs[Eϕ]},{t,Re[Eρ]},{t,Re[Eϕ]},{t,Im[Eρ]},{t,Im[Eϕ]},{t,Sqrt[Abs[Esρ]^2+Abs[Esϕ]^2]},{t,Abs[Esρ]},{t,Abs[Esϕ]}},{i,num}]]];

Code for generating figures

8.2 Propagating modes on wires of circular cross section

The function cylinder tabulates the normal modes for a cylindrical wire or metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): initial radius of cylinder in microns
	(Subscript[R, 2]): final radius of cylinder in microns
	(Subscript[R, inc]): increment of cyinder radius
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder
	(Subscript[n, g]): guess for the initial SP refractive index, generally ~1

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of radius (Subscript[R, 1]) using wavelength increment (Subscript[λ, inc]).    If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[R, 1]) ≠ (Subscript[R, 2]) then the function computes the wavevector vs. radius of the cylinder at wavelength (Subscript[λ, 1]).  If (Subscript[λ, 1])=(Subscript[λ, 2]) and (Subscript[R, 1])=(Subscript[R, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of radius (Subscript[R, 1]) at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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

The function planarFilm computes the SP fields for propagation on a planar surface.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[x, 1]): initial point from center of film for calculation in microns
	(Subscript[x, 2]): final point from center of film for calculation in microns
	(Subscript[x, inc]): increment of radius for calculation in microns
	t: thickness of metallic film
	(Subscript[n, mid]): refractive index of metallic film
	(Subscript[n, out]): refractive index of dielectric outside metallic film

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a planar metallic film of thickness t using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[x, 1])≠(Subscript[x, 2]), then the function computes fields for a planar metallic film of thickness t from (Subscript[x, 1]) to (Subscript[x, 2]) by increment (Subscript[x, inc]) measured from the center of the film.  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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.633;         (*  initial wavelength  *)
λ2 = 0.633;         (*  final wavelength  *)
λinc = 0;             (*  increment of wavelength  *)
R1 =0.025;          (*  radius of cylinder in microns  *)
R2 =0.025;         (*  radius of cylinder in microns  *)
Rinc = 0.0;
r1 = 0.001;         (*  initial radius in microns for calculation  *)
r2 = 0.300;                  (*  final radius in microns for calculation  *)
rinc = 0.001;              (*  increment of radius for calculation  *)
nin = 0.135+I 3.99;     (* refractive index of silver wire  *)  
nout = 1;                   (*  refractive index of air around wire  *)
ng=1.1;                         (*  initial guess for SP refractive index  *)
fields= cylinder[λ1,λ2,λinc,r1,r2,rinc,R1,R2,Rinc,nin,nout,ng];

Fig3a=ListLinePlot[{fields[[All,1]],fields[[All,2]],fields[[All,3]]}, PlotLegend -> {"|(Subscript[E, z])|", "|(Subscript[E, ρ])|", "|E|"}, LegendPosition -> {0.35, 0.2}, LegendSize ->{0.4,0.35}, LegendShadow -> {.01, -.01},FrameLabel -> {"ρ (μm)", "Field Amplitudes"},GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[option]];
Fig3aBW=ListLinePlot[{fields[[All,1]],fields[[All,2]],fields[[All,3]]}, PlotLegend -> {"|(Subscript[E, z])|", "|(Subscript[E, ρ])|", "|E|"}, LegendPosition -> {0.35, 0.2}, LegendSize ->{0.4,0.35}, LegendShadow -> {.01, -.01},FrameLabel -> {"ρ (μm)", "Field Amplitudes"},GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[optB]];
Clear[λ1,λ2,λinc,R1,R2,Rinc,r1,r2,rinc,nin,nout,ng,fields];

λ1 = 0.633;               (*  initial wavelength  *)
λ2 = 0.633;                  (*  final wavelength  *)
λinc = 0;             (*  increment of wavelength  *)
t =0.05;                       (*  thickness of film in microns  *)
x1 = 0.0;                  (*  initial point on x axis for calculation  *)
x2 = 0.3;                  (*  final point on x axis for calculation  *)
xinc = 0.001;              (*  increment of distance on x axis for calculation  *)
nmid = 0.135+I 3.99;     (* refractive index of silver wire  *)  
nout = 1;                   (*  refractive index of air around wire  *)
fields= planarFilm[λ1, λ2, λinc, x1, x2, xinc,t, nmid, nout];

Fig3b=ListLinePlot[{fields[[All,1]],fields[[All,2]],fields[[All,3]]}, PlotLegend -> {"|(Subscript[E, z])|","|(Subscript[E, x])|","|E|"}, LegendPosition -> {0.3, 0.25}, LegendSize -> {0.4,0.35}, LegendShadow -> {.01, -.01}, FrameLabel -> {"x (μm)", "Field Amplitudes"},GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[option]];
Fig3bBW=ListLinePlot[{fields[[All,1]],fields[[All,2]],fields[[All,3]]}, PlotLegend -> {"|(Subscript[E, z])|","|(Subscript[E, x])|","|E|"}, LegendPosition -> {0.3, 0.25}, LegendSize -> {0.4,0.35}, LegendShadow -> {.01, -.01}, FrameLabel -> {"x (μm)", "Field Amplitudes"},GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[optB]];
Fig3=GraphicsRow[{Fig3a,Fig3b},ImageSize->gr2,Spacings->{-30,30},Frame->True]
Fig3BW=GraphicsRow[{Fig3aBW,Fig3bBW},ImageSize->gr2,Spacings->{-30,30},Frame->True];

Export[FileNameJoin[{base,"Fig8.3.JPG"}],Fig3BW,ImageResolution->600];
Clear[λ1,λ2,λinc,t,x1,x2,xinc,nmid,nout,fields];

The function cylinder tabulates the normal modes for a cylindrical wire or metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	R: radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for SP effective refractive index

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of radius R using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of radius R from (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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

λinc = 0.01;         (*  increment of wavelength  *)
R1 =0.1;               (*  initial radius of cylinder in microns  *)
Rinc=0.0005;            (*  increment of radius  *)
r1 = 0.2;                  (*  initial radius in microns for calculation  *)
r2 = 0.2;                  (*  final radius in microns for calculation  *)
rinc = 0.0;              (*  increment of radius for calculation  *)
nin = -7;                   (* refractive index of silver wire  *)  
nout = 1;                   (*  refractive index of air around wire  *)
ng=1.1;

λ1 = 0.5;                  (*  initial wavelength  *)
λ2 = 0.5;                  (*  final wavelength  *)
R2=0.012;                    (*  final radius of cylinder in microns  *)
kz1= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];
λ1 = 0.6;                  (*  initial wavelength  *)
λ2 = 0.6;                  (*  final wavelength  *)
R2=0.008;                    (*  final radius of cylinder in microns  *)
kz2= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];
λ1 = 0.8;                  (*  initial wavelength  *)
λ2 = 0.8;                  (*  final wavelength  *)
R2=0.005;                    (*  final radius of cylinder in microns  *)
kz3= cylinder[λ1,λ2,λinc,r1,r2,rinc,R1,R2,Rinc,nin,nout,ng];

Fig4a=ListLinePlot[{kz1[[All,2]],kz2[[All,2]],kz3[[All,2]]}, PlotLegend -> {"λ=500 nm", "λ=600 nm", "λ=800 nm"}, LegendPosition -> {-0.3, 0.05}, LegendSize -> {0.7,0.4}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "(Subscript[n, eff]^′)"}, GridLines -> Automatic, Frame -> True,PlotRange -> All, Evaluate[option]];
Fig4aBW=ListLinePlot[{kz1[[All,2]],kz2[[All,2]],kz3[[All,2]]},  FrameLabel -> {"R (μm)", "(Subscript[n, eff]^′)"}, GridLines -> Automatic, Frame -> True,PlotRange -> All, Evaluate[optB]];

Fig4b=ListLinePlot[{kz1[[All,3]],kz2[[All,3]],kz3[[All,3]]},  FrameLabel -> {"R (μm)", "(Subscript[n, eff]^″)"}, GridLines -> Automatic, Frame -> True,PlotRange -> All,Evaluate[option]];
Fig4bBW=ListLinePlot[{kz1[[All,3]],kz2[[All,3]],kz3[[All,3]]},  FrameLabel -> {"R (μm)", "(Subscript[n, eff]^″)"}, GridLines -> Automatic, Frame -> True,PlotRange -> All,Evaluate[optB]];

Fig4=GraphicsRow[{Fig4a,Fig4b},ImageSize->gr2,Spacings->{-30,10},Frame->True]
Fig4BW=GraphicsRow[{Fig4aBW,Fig4bBW},ImageSize->gr2,Spacings->{-30,10},Frame->True];

Export[FileNameJoin[{base,"Fig8.4.JPG"}],Fig4BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,Rinc,r1,r2,rinc,nin,nout,ng,kz1,kz2,kz3];

The function cylinder tabulates the normal modes for a cylindrical wire or metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	R: radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of radius R using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of radius R from (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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

λinc = 0.01;         (*  increment of wavelength  *)
R1 =0.1;               (*  initial radius of cylinder in microns  *)
Rinc=0.0005;            (*  increment of radius  *)
r1 = 0.2;                  (*  initial radius in microns for calculation  *)
r2 = 0.2;                  (*  final radius in microns for calculation  *)
rinc = 0.0;              (*  increment of radius for calculation  *)
nin = -7;                   (* refractive index of silver wire  *)  
nout = 1;                   (*  refractive index of air around wire  *)
ng=1.1;                    (*  initial guess for SP refractive index  *)

λ1 = 0.5;                  (*  initial wavelength  *)
λ2 = 0.5;                  (*  final wavelength  *)
R2=0.012;                    (*  final radius of cylinder in microns  *)
kz1= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];
λ1 = 0.6;                  (*  initial wavelength  *)
λ2 = 0.6;                  (*  final wavelength  *)
R2=0.008;                    (*  final radius of cylinder in microns  *)
kz2= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];
λ1 = 0.8;                  (*  initial wavelength  *)
λ2 = 0.8;                  (*  final wavelength  *)
R2=0.005;                    (*  final radius of cylinder in microns  *)
kz3= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];

Fig5=ListLinePlot[{kz1[[All,4]],kz2[[All,4]],kz3[[All,4]]}, PlotLegend -> {"λ=500 nm", "λ=600 nm", "λ=800 nm"}, LegendPosition -> {-0.7, 0.2}, LegendSize -> {0.7,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Propagation Distance (μm)"},  GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[option]]
Fig5BW=ListLinePlot[{kz1[[All,4]],kz2[[All,4]],kz3[[All,4]]}, FrameLabel -> {"R (μm)", "Propagation Distance (μm)"},  GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[optB]];

Export[FileNameJoin[{base,"Fig8.5.JPG"}],Fig5BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,Rinc,r1,r2,rinc,nin,nout,ng,kz1,kz2,kz3];

The function cylinder tabulates the normal modes for a cylindrical wire or metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	R: radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for SP refractive index

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of radius R using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of radius R from (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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

λinc = 0.1;         (*  increment of wavelength  *)
R1 =0.1;               (*  initial radius of cylinder in microns  *)
R2=10;                    (*  final radius of cylinder in microns  *)
Rinc=0.1;            (*  increment of radius  *)
r1 = 0.2;                  (*  initial radius in microns for calculation  *)
r2 = 0.2;                  (*  final radius in microns for calculation  *)
rinc = 0.0;              (*  increment of radius for calculation  *)
nin = -7;                   (* refractive index of silver wire  *)  
nout = 1;                   (*  refractive index of air around wire  *)
ng=1.1;                  (*  initial guess for SP refractive index  *)

λ1 = 0.8;                  (*  initial wavelength  *)
λ2 = 0.8;                  (*  final wavelength  *)
kz1= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];
n0=nout;
n1=nn[nin,λ1];
kzplane=n0 n1 k0[λ1]/Sqrt[n0^2+n1^2];
d=1/Im[kzplane];
kzp={{R1,d},{R2,d}};

Fig6=ListLinePlot[{kz1[[All,4]],kzp}, PlotLegend -> {"nanowire","planar limit"}, LegendPosition -> {0.1, -0.2}, LegendSize -> {0.5,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Propagation Distance (μm)"}, GridLines -> Automatic, Frame -> True, PlotRange -> All,  Evaluate[option]]
Fig6BW=ListLinePlot[{kz1[[All,4]],kzp}, PlotLegend -> {"nanowire","planar limit"}, LegendPosition -> {0.1, 0}, LegendSize -> {0.5,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Propagation Distance (μm)"}, GridLines -> Automatic, Frame -> True, PlotRange -> All,  Evaluate[optB]];

Export[FileNameJoin[{base,"Fig8.6.JPG"}],Fig6BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,Rinc,r1,r2,rinc,nin,nout,ng,kz1,kzplane,kzp,n0,n1,d];

The function cylinder tabulates the normal modes for a cylindrical wire or metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	R: radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for refractive index of SP

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of radius R using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of radius R from (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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;                                     (*  background dielectric constant for Drude metal  *)
ωp=15;                         (*  ωp is plasma freq. in eV for Drude function  *)      
kp=2 π ωp/7.7905;            (* kp is plasma wavevector in ((μm)^(-1)) for Drude fn. *) 
γp = 0.1;                               (*  Drude damping coefficient in eV  *)
Fmode={{0,ωp/(2 π Sqrt[2])},{60,ωp/(2 π Sqrt[2])}};

λ1 = 2.0;                  (*  initial wavelength  *)
λ2 = 1.24/1.68;     (*  final wavelength  *)
λinc = 0.01;         (*  increment of wavelength  *)
R1 =0.05;               (*  initial radius of cylinder in microns  *)
R2=0.05;                    (*  final radius of cylinder in microns  *)
Rinc=0.0;            (*  increment of radius  *)
r1 = 0.2;                  (*  initial radius in microns for calculation  *)
r2 = 0.2;                  (*  final radius in microns for calculation  *)
rinc = 0.0;              (*  increment of radius for calculation  *)
nin = -5;                   (* refractive index of Drude metal  *)  
nout = 1;                   (*  refractive index of air around wire  *)
ng=1.1;                    (*  initial guess for SP refractive index  *)

Off[FindRoot::lstol]
kz1= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];
Fig7=ListLinePlot[{kz1[[All,5]],kz1[[All,6]],Fmode},PlotLegend -> {"(Subscript[n, eff]^′)","(Subscript[n, eff]^″)","Fröhlich mode"}, LegendPosition -> {-0.2, -0.2}, LegendSize -> {0.6,0.3}, LegendShadow -> {.01, -.01},FrameLabel -> {"Effective Index", "Frequency (eV)"},  GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[option]]
Fig7BW=ListLinePlot[{kz1[[All,5]],kz1[[All,6]],Fmode},PlotLegend -> {"(Subscript[n, eff]^′)","(Subscript[n, eff]^″)","Fröhlich mode"}, LegendPosition -> {-0.2, 0}, LegendSize -> {0.6,0.3}, LegendShadow -> {.01, -.01},FrameLabel -> {"Effective Index", "Frequency (eV)"},  GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[optB]];

Export[FileNameJoin[{base,"Fig8.7.JPG"}],Fig7BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,Rinc,r1,r2,rinc,nin,nout,ng,Fmode,kz1];
On[FindRoot::lstol]

8.4  Propagating modes on hollow cylindrical waveguides

The function cylinder tabulates the normal modes for a cylindrical wire or metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): initial radius of cylinder in microns
	(Subscript[R, 2]): final radius of cylinder in microns
	(Subscript[R, inc]): increment of cyinder radius
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for the SP refractive index

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of radius (Subscript[R, 1]) using wavelength increment (Subscript[λ, inc]).    If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[R, 1]) ≠ (Subscript[R, 2]) then the function computes the wavevector vs. radius of the cylinder at wavelength (Subscript[λ, 1]).  If (Subscript[λ, 1])=(Subscript[λ, 2]) and (Subscript[R, 1])=(Subscript[R, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of radius (Subscript[R, 1]) at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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.633;                 (*  initial wavelength  *)
λ2 = 0.633;                  (*  final wavelength  *)
λinc = 0;                     (*  increment of wavelength  *)
R1 =1;                       (*  radius of cylinder in microns  *)
R2 =0.1;                       (*  radius of cylinder in microns  *)
Rinc = 0.001;
r1 = 0;                  (*  initial radius in microns for calculation  *)
r2 = 0;                  (*  final radius in microns for calculation  *)
rinc = 0;              (*  increment of radius for calculation  *)
nin = 1;                   (*  refractive index of air around wire  *)
nout = 0.135+I 3.99;     (* refractive index of silver wire  *)  

Off[FindRoot::lstol]
ng=1.0117084734882436`+0.0021518802114859365` I;
kz1= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];
R2 =0.4;                       (*  radius of cylinder in microns  *)
ng=0.8616609165562797`+0.0009855660003827405` I;
kz2= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];
On[FindRoot::lstol]

Fig18=ListLinePlot[{kz1[[All,2]],kz2[[All,2]]},PlotLegend -> {"(Subscript[TM, 01])", "(Subscript[TM, 02])"}, LegendPosition -> {0.35, -0.2}, LegendSize ->{0.3,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "(Subscript[n, SP]^′)"}, GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[option]]
Fig18BW=ListLinePlot[{kz1[[All,2]],kz2[[All,2]]},PlotLegend -> {"(Subscript[TM, 01])", "(Subscript[TM, 02])"}, LegendPosition -> {0.25, -0.2}, LegendSize ->{0.4,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "(Subscript[n, SP]^′)"}, GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[optB]];

Export[FileNameJoin[{base,"Fig8.18.JPG"}],Fig18BW,ImageResolution->600];

Fig19=ListLinePlot[{kz1[[All,4]],kz2[[All,4]]},PlotLegend -> {"(Subscript[TM, 01])", "(Subscript[TM, 02])"}, LegendPosition -> {-0.6, 0.2}, LegendSize -> {0.4,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "d (μm)"}, GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[option]]
Fig19BW=ListLinePlot[{kz1[[All,4]],kz2[[All,4]]},PlotLegend -> {"(Subscript[TM, 01])", "(Subscript[TM, 02])"}, LegendPosition -> {-0.6, 0.2}, LegendSize -> {0.4,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "d (μm)"}, GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[optB]];

Export[FileNameJoin[{base,"Fig8.19.JPG"}],Fig19BW,ImageResolution->600];

Clear[λ1,λ2,λinc,R1,R2,Rinc,r1,r2,rinc,nin,nout,ng];

The function cylinder tabulates the normal modes for a cylindrical wire or metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): initial radius of cylinder in microns
	(Subscript[R, 2]): final radius of cylinder in microns
	(Subscript[R, inc]): increment of cyinder radius
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for SP refractive index

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of radius (Subscript[R, 1]) using wavelength increment (Subscript[λ, inc]).    If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[R, 1]) ≠ (Subscript[R, 2]) then the function computes the wavevector vs. radius of the cylinder at wavelength (Subscript[λ, 1]).  If (Subscript[λ, 1])=(Subscript[λ, 2]) and (Subscript[R, 1])=(Subscript[R, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of radius (Subscript[R, 1]) at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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.633;                 (*  initial wavelength  *)
λ2 = 0.633;                  (*  final wavelength  *)
λinc = 0;             (*  increment of wavelength  *)
R1 =0.1;                       (*  radius of cylinder in microns  *)
R2 =0.1;                       (*  radius of cylinder in microns  *)
Rinc = 0.0;
r1 = 0.001;                  (*  initial radius in microns for calculation  *)
r2 = 0.2;                  (*  final radius in microns for calculation  *)
rinc = 0.001;              (*  increment of radius for calculation  *)
nin = 1;                   (*  refractive index of air around wire  *)
nout = 0.135+I 3.99;     (* refractive index of silver wire  *)  
ng=0.0014597920175993726`+2.1170308347596007` I;
fields= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];

Fig20=ListLinePlot[{fields[[All,1]],fields[[All,2]],fields[[All,3]]}, PlotLegend -> {"|(Subscript[E, z])|", "|(Subscript[E, ρ])|", "|E|"}, LegendPosition -> {0.25, 0.2}, LegendSize -> {0.4,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"ρ (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[option]]
Fig20BW=ListLinePlot[{fields[[All,1]],fields[[All,2]],fields[[All,3]]}, PlotLegend -> {"|(Subscript[E, z])|", "|(Subscript[E, ρ])|", "|E|"}, LegendPosition -> {0.25, 0.2}, LegendSize -> {0.4,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"ρ (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[optB]];


Export[FileNameJoin[{base,"Fig8.20.JPG"}],Fig20BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,Rinc,r1,r2,rinc,nin,nout,ng];

The function cylinder tabulates the normal modes for a cylindrical wire or metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): initial radius of cylinder in microns
	(Subscript[R, 2]): final radius of cylinder in microns
	(Subscript[R, inc]): increment of cyinder radius
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for the SP refractive index

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of radius (Subscript[R, 1]) using wavelength increment (Subscript[λ, inc]).    If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[R, 1]) ≠ (Subscript[R, 2]) then the function computes the wavevector vs. radius of the cylinder at wavelength (Subscript[λ, 1]).  If (Subscript[λ, 1])=(Subscript[λ, 2]) and (Subscript[R, 1])=(Subscript[R, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of radius (Subscript[R, 1]) at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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.633;                 (*  initial wavelength  *)
λ2 = 0.633;                  (*  final wavelength  *)
λinc = 0;                     (*  increment of wavelength  *)
R1 =1;                       (*  radius of cylinder in microns  *)
R2 =1;                       (*  radius of cylinder in microns  *)
Rinc = 0;
r1 = 0.001;                  (*  initial radius in microns for calculation  *)
r2 = 2;                  (*  final radius in microns for calculation  *)
rinc = 0.001;              (*  increment of radius for calculation  *)
nin = 1;                   (*  refractive index of air around wire  *)
nout = 0.135+I 3.99;     (* refractive index of silver wire  *)  
ng=1.0117084734882436`+0.0021518802114859365` I;  (*  first n = 0 mode  *)
fields1= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];
Fig21a=ListLinePlot[{fields1[[All,1]],fields1[[All,2]],fields1[[All,3]]}, PlotLegend -> {"|(Subscript[E, z])|", "|(Subscript[E, ρ])|", "|E|"}, LegendPosition -> {0.1, 0.1}, LegendSize -> {0.4,0.4}, LegendShadow -> {.01, -.01},FrameLabel -> {"ρ (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[option]];
Fig21aBW=ListLinePlot[{fields1[[All,1]],fields1[[All,2]],fields1[[All,3]]}, PlotLegend -> {"|(Subscript[E, z])|", "|(Subscript[E, ρ])|", "|E|"}, LegendPosition -> {0.1, 0.1}, LegendSize -> {0.4,0.4}, LegendShadow -> {.01, -.01},FrameLabel -> {"ρ (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[optB]];

ng=0.8616609165562797`+0.0009855660003827405` I;  (*  second n = 0 mode  *)
fields2= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];
Fig21b=ListLinePlot[{fields2[[All,1]],fields2[[All,2]],fields2[[All,3]]}, FrameLabel -> {"ρ (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[option]];
Fig21bBW=ListLinePlot[{fields2[[All,1]],fields2[[All,2]],fields2[[All,3]]}, FrameLabel -> {"ρ (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange -> All, Evaluate[optB]];
Fig21=GraphicsRow[{Fig21a,Fig21b},ImageSize->gr2,Spacings->{-30,30},Frame->True]
Fig21BW=GraphicsRow[{Fig21aBW,Fig21bBW},ImageSize->gr2,Spacings->{-30,30},Frame->True];

Export[FileNameJoin[{base,"Fig8.21.JPG"}],Fig21BW,ImageResolution->600];
Clear[fields1,fields2,nin,nout,Rinc];

8.5  Propagating surface plasmons on cylindrical shells

The function "tube" tabulates the normal modes for a hollow cylindrical metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): inner radius of cylinder in microns
	(Subscript[R, 2]): outer radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, 1]): refractive index inside hollow region of cylinder
	(Subscript[n, 2]): refractive index cylinder itself
	(Subscript[n, 3]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for refractive index of SP

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of inner radius (Subscript[R, 1]) and outer radius (Subscript[R, 2]) using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]) but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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

R1 =0.08;                       (*  inner radius of cylinder in microns  *)
R2 =0.112;                       (*  outer radius of cylinder in microns  *)
r1 = 0.001;                  (*  initial radius in microns for calculation  *)
r2 = 0.2;                  (*  final radius in microns for calculation  *)
rinc = 0.001;              (*  increment of radius for calculation  *)
n1 = 1.414;               (* refractive index of region inside tube  *)  
n2=-23;                    (*  refractive index of gold from Palik, but with Im(ϵ)=0   *)
n3 = 1;                   (*  refractive index of air outside tube  *)

Off[General::stop]
Off[Greater::nord]
Off[FindRoot::lstol]
Off[Power::infy]
λ1 = 0.405;                 (*  initial wavelength  *)
λ2 = .384;                  (*  final wavelength  *)
λinc = 0.001;                     (*  increment of wavelength  *)
ng=1.01;                   (*  initial guess for SP refractive index  *)
kz1= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,-1]; 
λ1 = 0.385;                 (*  initial wavelength  *)
λ2 = .38364;                  (*  final wavelength  *)
λinc = 0.00001;                     (*  increment of wavelength  *)
ng=1.4;                   (*  initial guess for SP refractive index  *)
kz2= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,-1]; 
λ1 = 0.385;                 (*  initial wavelength  *)
λ2 = .38364;                  (*  final wavelength  *)
ng=2;                   (*  initial guess for SP refractive index  *)
kz3= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,-1];λ1 = 0.405;                 (*  initial wavelength  *)
λ2 = .384;                  (*  final wavelength  *)
λinc = 0.001;                     (*  increment of wavelength  *)
ng=1.6;                   (*  initial guess for SP refractive index  *)
kz4= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,-1];
λ1 = 0.8;                 (*  initial wavelength  *)
λ2 = .485;                  (*  final wavelength  *)
ng=1.1;                   (*  initial guess for SP refractive index  *)
kz5= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,+1]; 
On[General::stop]
On[Greater::nord]
On[FindRoot::lstol]
On[Power::infy]

Fig24=ListLinePlot[{kz1[[All,5]],kz2[[All,5]],kz3[[All,5]],kz4[[All,5]],kz5[[All,5]]}, FrameLabel -> {"(Subscript[n, eff]^′)", "Frequency (eV)"},LabelStyle -> {FontSize -> 12}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]]
Fig24BW=ListLinePlot[{kz1[[All,5]],kz2[[All,5]],kz3[[All,5]],kz4[[All,5]],kz5[[All,5]]}, FrameLabel -> {"(Subscript[n, eff]^′)", "Frequency (eV)"},LabelStyle -> {FontSize -> 12}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]];

Export[FileNameJoin[{base,"Fig8.24.JPG"}],Fig24BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,r1,r2,rinc,n1,n2,n3,ng,kz1,kz2,kz3,kz4,kz5];

The function "tube" tabulates the normal modes for a hollow cylindrical metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): inner radius of cylinder in microns
	(Subscript[R, 2]): outer radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, 1]): refractive index inside hollow region of cylinder
	(Subscript[n, 2]): refractive index cylinder itself
	(Subscript[n, 3]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for refractive index of SP

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of inner radius (Subscript[R, 1]) and outer radius (Subscript[R, 2]) using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]) but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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
	-23: complex index for gold interpolated from values in Palik, but with Im(ϵ) set = 0

λ1 = 0.496;                 (*  initial wavelength  *)
λ2 = 0.496;                  (*  final wavelength  *)
λinc = 0;                     (*  increment of wavelength  *)
R1 =0.08;                       (*  inner radius of cylinder in microns  *)
R2 =0.112;                       (*  outer radius of cylinder in microns  *)
r1 = 0.001;                  (*  initial radius in microns for calculation  *)
r2 = 0.2;                  (*  final radius in microns for calculation  *)
rinc = 0.001;              (*  increment of radius for calculation  *)
n1 = 1.414;               (* refractive index of dielectric inside tube  *)  
n2=-23;           (*  refractive index of gold with Im(ϵ) = 0   *)
n3 = 1;                   (*  refractive index of air outside tube  *)
ng=3;                   (*  initial guess for SP refractive index  *)

fields1= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,+1]; 
Fig25a=ListLinePlot[{fields1[[All,1]],fields1[[All,2]],fields1[[All,3]]}, PlotLegend -> {"(Subscript[E, z]^′)","(Subscript[E, ρ]^″)","|E|"}, LegendPosition -> {0.4, 0.25}, LegendSize -> {0.65,0.4}, LegendShadow -> {.01, -.01},FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]];
Fig25aBW=ListLinePlot[{fields1[[All,1]],fields1[[All,2]],fields1[[All,3]]}, PlotLegend -> {"(Subscript[E, z]^′)","(Subscript[E, ρ]^″)","|E|"}, LegendPosition -> {0.4, 0.25}, LegendSize -> {0.55,0.35}, LegendShadow -> {.01, -.01},FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]];

λ1 = 0.396;                 (*  initial wavelength  *)
λ2 = 0.396;                  (*  final wavelength  *)
fields2= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,+1];
Fig25b=ListLinePlot[{fields2[[All,1]],fields2[[All,2]],fields2[[All,3]]}, FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]] ;
Fig25bBW=ListLinePlot[{fields2[[All,1]],fields2[[All,2]],fields2[[All,3]]}, FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]] ;
Fig25=GraphicsRow[{Fig25a,Fig25b},ImageSize->gr2,Spacings->{-30,30},Frame->True]
Fig25BW=GraphicsRow[{Fig25aBW,Fig25bBW},ImageSize->gr2,Spacings->{-30,30},Frame->True];

Export[FileNameJoin[{base,"Fig8.25.JPG"}],Fig25BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,r1,r2,rinc,n1,n2,n3,ng,fields1,fields2];

The function "tube" tabulates the normal modes for a hollow cylindrical metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): inner radius of cylinder in microns
	(Subscript[R, 2]): outer radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, 1]): refractive index inside hollow region of cylinder
	(Subscript[n, 2]): refractive index cylinder itself
	(Subscript[n, 3]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for refractive index of SP

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of inner radius (Subscript[R, 1]) and outer radius (Subscript[R, 2]) using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]) but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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
	-23: complex index for gold interpolated from values in Palik, but with Im(ϵ) set = 0

λ1 = 1;                 (*  initial wavelength  *)
λ2 = 0.26;                  (*  final wavelength  *)
λinc = 0.001;                     (*  increment of wavelength  *)
R1 =0.08;                       (*  inner radius of cylinder in microns  *)
R2 =0.112;                       (*  outer radius of cylinder in microns  *)
r1 = 0.001;                  (*  initial radius in microns for calculation  *)
r2 = 0.2;                  (*  final radius in microns for calculation  *)
rinc = 0.001;              (*  increment of radius for calculation  *)
n1 = 1.414;               (* refractive index of region inside tube  *)  
n2=-6;                    (*  refractive index of gold from Rakic   *)
n3 = 1;                   (*  refractive index of air outside tube  *)
(*ng=1.1;                   (*  initial guess for SP refractive index  *)
kz1= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,-1]; *)
ng=1.1;                   (*  initial guess for SP refractive index  *)
Off[General::stop]
Off[Greater::nord]
Off[FindRoot::lstol]
kz2= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,+1];
On[General::stop]
On[Greater::nord]
On[FindRoot::lstol]

Fig27a=ListLinePlot[{kz2[[All,5]]},FrameLabel -> {"(Subscript[n, eff]^′)", "Frequency (eV)"},LabelStyle -> {FontSize -> 12}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]];
Fig27b=ListLinePlot[{kz2[[All,6]]}, FrameLabel -> {"(Subscript[n, eff]^″)", "Frequency (eV)"},LabelStyle -> {FontSize -> 12}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]];
Fig27aBW=ListLinePlot[{kz2[[All,5]]},FrameLabel -> {"(Subscript[n, eff]^′)", "Frequency (eV)"},LabelStyle -> {FontSize -> 12}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]];
Fig27bBW=ListLinePlot[{kz2[[All,6]]}, FrameLabel -> {"(Subscript[n, eff]^″)", "Frequency (eV)"},LabelStyle -> {FontSize -> 12}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]];
Fig27=GraphicsRow[{Fig27a,Fig27b},ImageSize->gr2,Spacings->{-30,30},Frame->True]
Fig27BW=GraphicsRow[{Fig27aBW,Fig27bBW},ImageSize->gr2,Spacings->{-30,30},Frame->True];

Export[FileNameJoin[{base,"Fig8.27.JPG"}],Fig27BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,r1,r2,rinc,n1,n2,n3,ng,kz1,kz2];

The function "tube" tabulates the normal modes for a hollow cylindrical metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): inner radius of cylinder in microns
	(Subscript[R, 2]): outer radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, 1]): refractive index inside hollow region of cylinder
	(Subscript[n, 2]): refractive index cylinder itself
	(Subscript[n, 3]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for refractive index of SP

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of inner radius (Subscript[R, 1]) and outer radius (Subscript[R, 2]) using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]) but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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
	-23: complex index for gold interpolated from values in Palik, but with Im(ϵ) set = 0

λinc = 0;                     (*  increment of wavelength  *)
R1 =0.08;                       (*  inner radius of cylinder in microns  *)
R2 =0.112;                       (*  outer radius of cylinder in microns  *)
r1 = 0.001;                  (*  initial radius in microns for calculation  *)
r2 = 0.2;                  (*  final radius in microns for calculation  *)
rinc = 0.001;              (*  increment of radius for calculation  *)
n1 = 1.414;               (* refractive index of region inside tube  *)  
n2=-6;                       (*  refractive index of gold from Rakic  *)
n3 = 1;                   (*  refractive index of air outside tube  *)
ng=1.1;                   (*  initial guess for SP refractive index  *)

λ1 = 0.875;                 (*  initial wavelength  *)
λ2 = 0.875;                  (*  final wavelength  *)
fields1= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,+1]; 
Fig28a=ListLinePlot[{fields1[[All,1]],fields1[[All,2]],fields1[[All,3]]}, PlotLegend -> {"Re((Subscript[E, z]))","Im((Subscript[E, ρ]))","|E|"}, LegendPosition -> {-0.5, -0.2}, LegendSize -> {0.55,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]];
Fig28aBW=ListLinePlot[{fields1[[All,1]],fields1[[All,2]],fields1[[All,3]]}, PlotLegend -> {"Re((Subscript[E, z]))","Im((Subscript[E, ρ]))","|E|"}, LegendPosition ->{-0.5, -0.12}, LegendSize -> {0.55,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]];

λ1 = 0.506;                 (*  initial wavelength  *)
λ2 = 0.506;                  (*  final wavelength  *)
fields2= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,+1];
Fig28b=ListLinePlot[{fields2[[All,1]],fields2[[All,2]],fields2[[All,3]]},  FrameLabel -> {"R (μm)", "Field Amplitudes"},GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]] ;
Fig28bBW=ListLinePlot[{fields2[[All,1]],fields2[[All,2]],fields2[[All,3]]},  FrameLabel -> {"R (μm)", "Field Amplitudes"},GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]] ;
Fig28=GraphicsRow[{Fig28a,Fig28b},ImageSize->gr2,Spacings->{0,30},Frame->True]
Fig28BW=GraphicsRow[{Fig28aBW,Fig28bBW},ImageSize->gr2,Spacings->{0,30},Frame->True];

Export[FileNameJoin[{base,"Fig8.28.JPG"}],Fig28BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,r1,r2,rinc,n1,n2,n3,ng,fields1,fields2];

The function "tube" tabulates the normal modes for a hollow cylindrical metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): inner radius of cylinder in microns
	(Subscript[R, 2]): outer radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, 1]): refractive index inside hollow region of cylinder
	(Subscript[n, 2]): refractive index cylinder itself
	(Subscript[n, 3]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for refractive index of SP

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of inner radius (Subscript[R, 1]) and outer radius (Subscript[R, 2]) using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]) but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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
	-23: complex index for gold interpolated from values in Palik, but with Im(ϵ) set = 0

λ1 = 0.83;                 (*  initial wavelength  *)
λ2 = 0.83;                  (*  final wavelength  *)
λinc = 0;                     (*  increment of wavelength  *)
R1 =0.38;                       (*  inner radius of cylinder in microns  *)
R2 =0.42;                       (*  outer radius of cylinder in microns  *)
r1 = 0.001;                  (*  initial radius in microns for calculation  *)
r2 = 0.8;                  (*  final radius in microns for calculation  *)
rinc = 0.001;              (*  increment of radius for calculation  *)
n1 = 1.414;               (* refractive index of region inside tube  *)          
n2=-1;                   (* refractive index of gold tube from Palik  *)  
n3 = 1;                   (*  refractive index of air outside tube  *)

ng=1.29;                   (*  initial guess for SP refractive index  *)
fields1= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,+1]; 
Fig29a=ListLinePlot[{fields1[[All,1]],fields1[[All,2]],fields1[[All,3]]}, PlotLegend -> {"(Subscript[E, z]^′)","(Subscript[E, ρ]^″)","|E|"}, LegendPosition -> {0.25, 0.2}, LegendSize -> {0.3,0.35}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]];
Fig29aBW=ListLinePlot[{fields1[[All,1]],fields1[[All,2]],fields1[[All,3]]}, PlotLegend -> {"(Subscript[E, z]^′)","(Subscript[E, ρ]^″)","|E|"}, LegendPosition -> {0.25, 0.2}, LegendSize -> {0.3,0.35}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]];

ng=1.04;                   (*  initial guess for SP refractive index  *)
fields2= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,+1];
Fig29b=ListLinePlot[{fields2[[All,1]],fields2[[All,2]],fields2[[All,3]]}, FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]];
Fig29bBW=ListLinePlot[{fields2[[All,1]],fields2[[All,2]],fields2[[All,3]]}, FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]];
Fig29=GraphicsRow[{Fig29a,Fig29b},ImageSize->gr2,Spacings->{-30,30},Frame->True]
Fig29BW=GraphicsRow[{Fig29aBW,Fig29bBW},ImageSize->gr2,Spacings->{-30,30},Frame->True];

Export[FileNameJoin[{base,"Fig8.29.JPG"}],Fig29BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,r1,r2,rinc,n1,n2,n3,ng,fields1,fields2];

The function "tube" tabulates the normal modes for a hollow cylindrical metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[R, 1]): inner radius of cylinder in microns
	(Subscript[R, 2]): outer radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, 1]): refractive index inside hollow region of cylinder
	(Subscript[n, 2]): refractive index cylinder itself
	(Subscript[n, 3]): refractive index outside cylinder
	(Subscript[n, g]): initial guess for the refractive index of the SP

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes effective index of the SP vs. wavelength for a cylinder of inner radius (Subscript[R, 1]) and outer radius (Subscript[R, 2]) using wavelength increment (Subscript[λ, inc]).  The function outputs six tables of values for plotting, {λ, (Subscript[k, z])/k}, {λ, Re((Subscript[k, z]))/(Subscript[k, 0])}, {λ, Im((Subscript[k, z]))/(Subscript[k, 0])}, {λ, 1/Im((Subscript[k, z]))}, {Re((Subscript[k, z]))/(Subscript[k, p]), ω/(Subscript[ω, p])}, and {Im((Subscript[k, z]))/(Subscript[k, p]), ω/(Subscript[ω, p])}.  The last two tables are for use with a Drude metal.  If (Subscript[λ, 1])=(Subscript[λ, 2]) but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of inner radius (Subscript[R, 1]) and outer radius (Subscript[R, 2]) at wavelength (Subscript[λ, 1]) from (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The function outputs three tables of values for plotting, {r, | Ez |}, {r, |Eρ|}, and {r, Sqrt[|Eρ((|)^2) + |Ez((|)^2)]}.  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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 = 2.0;                  (*  initial wavelength  *)
λ2 = 1.16;                  (*  final wavelength  *)
λinc = 0.01;         (*  increment of wavelength  *)
R1 =0.02;               (*  inner radius of tube in microns  *)
R2=0.025;                (*  outer radius of cylinder/tube in microns  *)
Rinc=0;
r1 = 0.2;                  (*  initial radius in microns for calculation  *)
r2 = 0.2;                  (*  final radius in microns for calculation  *)
rinc = 0.0;              (*  increment of radius for calculation  *)
n1 = 1.5;                   (* refractive index of Drude metal  *) 
n2=-7;                    (*  refractive index of silver  *) 
n3 = 1;                   (*  refractive index of air around wire  *)
ng=1.1;                   (*  initial guess for SP refractive index  *)

Off[General::stop]
Off[Greater::nord]
kzwire= cylinder[λ1,λ2,λinc,r1,r2,rinc,R2,R2,Rinc,n2,n3,ng];
kztube= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,+1];
Fig30=ListLinePlot[{kztube[[All,4]],kzwire[[All,4]]}, PlotLegend -> {"shell","wire"}, LegendPosition -> {-0.7, 0.3}, LegendSize -> {0.3,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "ξ (μm)"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]]
Fig30BW=ListLinePlot[{kztube[[All,4]],kzwire[[All,4]]}, PlotLegend -> {"shell","wire"}, LegendPosition -> {-0.5, 0.3}, LegendSize -> {0.3,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "ξ (μm)"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]];

Export[FileNameJoin[{base,"Fig8.30.JPG"}],Fig30BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,Rinc,r1,r2,rinc,n1,n2,n3,ng,kzwire,kztube];
On[General::stop]
On[Greater::nord]

8.6  Interaction with plane waves; nonpropagating surface plasmons on nanowires

The function "local" computes the local fields around a metallic cylinder when illuminated by a plane wave such that the cylinder axis is normal to the plane of incidence and the incident field is polarized in the plane of incidence at right angles to the cylinder axis.  The parameters that are passed by the user to the function are:
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder
	R: radius of cylinder in microns
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[r, 1]): initial radius in μm for calculation of local fields
	(Subscript[r, 2]): final radius in μm for calculation of local fields
	ϕ: angle in degrees for field calculation in xy plane measured from the x axis.  0 degrees is the backward scattered direction.
	(Subscript[x, 1]): initial x axis point in μm for calculation of local fields
	(Subscript[x, 2]): final x axis point in μm for calculation of local fields
	(Subscript[y, 1]): initial x axis point in μm for calculation of local fields
	(Subscript[y, 2]): final x axis point in μm for calculation of local fields
	inc: increment in μm of radius, x position, or y position for calculation of local fields

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes the peak |E((|)^2) field intensity at the surface of the cylinder and the scattering efficiency vs. wavelength for a cylinder of radius (Subscript[R, 1]) using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]) or (Subscript[x, 1])≠(Subscript[x, 2]) or (Subscript[y, 1])≠(Subscript[y, 2]), then the function computes fields at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) at angle ϕ or along the x axis from (Subscript[x, 1]) to (Subscript[x, 2]) or along the y axis from (Subscript[y, 1]) to (Subscript[y, 2]) by increment inc.  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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 μm  *)
λ2= 2.0;                      (*  final wavelength in μm  *)
λinc= 0.01;               (*  increment of wavelength in μm  *)
n1 = -6;                      (*  refractive index of gold nanowire  *)
n2 = 1.0;                     (* refractive index of surrounding medium  *)
nmax = 20;                       (*  # modes to keep in the expansion  *)
r1 = 0;                       (*  initial radius in μm for field calculation  *)
r2 = 0;                       (*  final radius in μm for field calculation  *)
ϕ = 90;         (*  angle in xy plane measured from x axis for field calculation  *)
x1 = 0;                       (*  initial point in μm on x axis for calculation  *)
x2 = 0;                         (*  final point in μm on x axis for calculation  *)
y1 = 0;                       (*  initial point in μm on y axis for calculation  *)
y2 = 0;                         (*  final point in μm on y axis for calculation  *)
inc = 0;                 (*  increment of position (r, x, or y) in μm  *)

R = 0.025;                          (*  radius of cylinder in μm  *)
i1=local[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc];
Print["Completed 50 nm wire."]
R = 0.05;                          (*  radius of cylinder in μm  *)
i2=local[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc];
Print["Completed 100 nm wire."]
R = 0.075;                          (*  radius of cylinder in μm  *)
i3=local[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc];
Print["Completed 150 nm wire."]
R = 0.1;                          (*  radius of cylinder in μm  *)
i4=local[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc];

Fig32=ListLinePlot[{i1[[All,1]],i2[[All,1]],i3[[All,1]],i4[[All,1]]},PlotLegend -> { "50 nm","100 nm","150 nm", "200 nm"}, LegendPosition -> {0.35, 0.25}, LegendSize -> {0.45,0.35}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "|E((|)^2)"}, Evaluate[option]]
Fig32BW=ListLinePlot[{i1[[All,1]],i2[[All,1]],i3[[All,1]],i4[[All,1]]},PlotLegend -> { "50 nm","100 nm","150 nm", "200 nm"}, LegendPosition -> {0.35, 0.25}, LegendSize -> {0.45,0.35}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "|E((|)^2)"}, Evaluate[optB]];

Export[FileNameJoin[{base,"Fig8.32.JPG"}],Fig32BW,ImageResolution->600];
Clear[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc,i1,i2,i3];

The function "local" computes the local fields around a metallic cylinder when illuminated by a plane wave such that the cylinder axis is normal to the plane of incidence and the incident field is polarized in the plane of incidence at right angles to the cylinder axis.  The parameters that are passed by the user to the function are:
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder
	R: radius of cylinder in microns
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[r, 1]): initial radius in μm for calculation of local fields
	(Subscript[r, 2]): final radius in μm for calculation of local fields
	ϕ: angle in degrees for field calculation in xy plane measured from the x axis.  0 degrees is the backward scattered direction.
	(Subscript[x, 1]): initial x axis point in μm for calculation of local fields
	(Subscript[x, 2]): final x axis point in μm for calculation of local fields
	(Subscript[y, 1]): initial x axis point in μm for calculation of local fields
	(Subscript[y, 2]): final x axis point in μm for calculation of local fields
	inc: increment in μm of radius, x position, or y position for calculation of local fields

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes the peak |E((|)^2) field intensity at the surface of the cylinder vs. wavelength for a cylinder of radius (Subscript[R, 1]) using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]) or (Subscript[x, 1])≠(Subscript[x, 2]) or (Subscript[y, 1])≠(Subscript[y, 2]), then the function computes fields at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) at angle ϕ or along the x axis from (Subscript[x, 1]) to (Subscript[x, 2]) or along the y axis from (Subscript[y, 1]) to (Subscript[y, 2]) by increment inc.  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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.53;                      (*  initial wavelength in μm  *)
λ2= 0.53;                      (*  final wavelength in μm  *)
λinc= 0;                      (*  increment of wavelength in μm  *)
n1 = -6;                     (*  refractive index of gold nanowire  *)
n2 = 1.0;                     (* refractive index of surrounding medium  *)
R = 0.025;                          (*  radius of cylinder in μm  *)
nmax = 20;                       (*  # modes to keep in the expansion  *)
r1 = 0;                       (*  initial radius in μm for field calculation  *)
r2 = 0;                       (*  final radius in μm for field calculation  *)
ϕ = 90;         (*  angle in xy plane measured from x axis for field calculation  *)
x1 = 0;                       (*  initial point in μm on x axis for calculation  *)
x2 = 0;                         (*  final point in μm on x axis for calculation  *)
y1 = -0.1;                       (*  initial point in μm on y axis for calculation  *)
y2 = +0.1;                         (*  final point in μm on y axis for calculation  *)
inc = 0.001;                 (*  increment of position (r, x, or y) in μm  *)

fields=local[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc];

Fig34=ListLinePlot[{fields[[All,1]],fields[[All,2]],fields[[All,3]]}, PlotLegend -> { "|E|","|(Subscript[E, ρ])|", "|(Subscript[E, ϕ])|"}, LegendPosition -> {0.45, 0.3}, LegendSize -> {0.3,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Field Amplitudes"},Evaluate[option]]
Fig34BW=ListLinePlot[{fields[[All,1]],fields[[All,2]],fields[[All,3]]}, PlotLegend -> { "|E|","|(Subscript[E, ρ])|", "|(Subscript[E, ϕ])|"}, LegendPosition -> {0.45, 0.3}, LegendSize -> {0.3,0.3}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Field Amplitudes"},Evaluate[optB]];

Export[FileNameJoin[{base,"Fig8.34.JPG"}],Fig34BW,ImageResolution->600];
Clear[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc,fields];

The function "local" computes the local fields around a metallic cylinder when illuminated by a plane wave such that the cylinder axis is normal to the plane of incidence and the incident field is polarized in the plane of incidence at right angles to the cylinder axis.  The parameters that are passed by the user to the function are:
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder
	R: radius of cylinder in microns
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	(Subscript[r, 1]): initial radius in μm for calculation of local fields
	(Subscript[r, 2]): final radius in μm for calculation of local fields
	ϕ: angle in degrees for field calculation in xy plane measured from the x axis.  0 degrees is the backward scattered direction.
	(Subscript[x, 1]): initial x axis point in μm for calculation of local fields
	(Subscript[x, 2]): final x axis point in μm for calculation of local fields
	(Subscript[y, 1]): initial x axis point in μm for calculation of local fields
	(Subscript[y, 2]): final x axis point in μm for calculation of local fields
	inc: increment in μm of radius, x position, or y position for calculation of local fields

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes the peak |E((|)^2) field intensity at the surface of the cylinder and the scattering efficiency vs. wavelength for a cylinder of radius (Subscript[R, 1]) using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]) or (Subscript[x, 1])≠(Subscript[x, 2]) or (Subscript[y, 1])≠(Subscript[y, 2]), then the function computes fields at wavelength (Subscript[λ, 1]) from radial position (Subscript[r, 1]) to (Subscript[r, 2]) at angle ϕ or along the x axis from (Subscript[x, 1]) to (Subscript[x, 2]) or along the y axis from (Subscript[y, 1]) to (Subscript[y, 2]) by increment inc.  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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 μm  *)
λ2= 2.0;                      (*  final wavelength in μm  *)
λinc= 0.01;               (*  increment of wavelength in μm  *)
n1 = -6;                      (*  refractive index of gold nanowire  *)
n2 = 1.0;                     (* refractive index of surrounding medium  *)
nmax = 20;                       (*  # modes to keep in the expansion  *)
r1 = 0;                       (*  initial radius in μm for field calculation  *)
r2 = 0;                       (*  final radius in μm for field calculation  *)
ϕ = 90;         (*  angle in xy plane measured from x axis for field calculation  *)
x1 = 0;                       (*  initial point in μm on x axis for calculation  *)
x2 = 0;                         (*  final point in μm on x axis for calculation  *)
y1 = 0;                       (*  initial point in μm on y axis for calculation  *)
y2 = 0;                         (*  final point in μm on y axis for calculation  *)
inc = 0;                 (*  increment of position (r, x, or y) in μm  *)

R = 0.025;                          (*  radius of cylinder in μm  *)
i1=local[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc];
Print["Completed 50 nm wire."]
R = 0.05;                          (*  radius of cylinder in μm  *)
i2=local[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc];
Print["Completed 100 nm wire."]
R = 0.075;                          (*  radius of cylinder in μm  *)
i3=local[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc];
Print["Completed 150 nm wire."]
R = 0.1;                          (*  radius of cylinder in μm  *)
i4=local[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc];

Fig35=ListLinePlot[{i1[[All,2]],i2[[All,2]],i3[[All,2]],i4[[All,2]]},PlotLegend -> { "50 nm","100 nm","150 nm", "200 nm"}, LegendPosition -> {0.33, 0.1}, LegendSize -> {0.5,0.35}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[Q, sct])"}, Evaluate[option]]
Fig35BW=ListLinePlot[{i1[[All,2]],i2[[All,2]],i3[[All,2]],i4[[All,2]]},PlotLegend -> { "50 nm","100 nm","150 nm", "200 nm"}, LegendPosition -> {0.33, 0.1}, LegendSize -> {0.5,0.35}, LegendShadow -> {.01, -.01}, FrameLabel -> {"λ (μm)", "(Subscript[Q, sct])"}, Evaluate[optB]];

Export[FileNameJoin[{base,"Fig8.35.JPG"}],Fig35BW,ImageResolution->600];
Clear[n1,n2,R,nmax,λ1,λ2,λinc,r1,r2,ϕ,x1,x2,y1,y2,inc,i1,i2,i3];

Exercises

The function cylinder tabulates the normal modes for a cylindrical wire or metallic tube.  The parameters that are passed by the user to the function are:
	(Subscript[λ, 1]):  initial wavelength in microns	
	(Subscript[λ, 2]):  final wavelength in microns
	(Subscript[λ, inc]): increment of wavelength in microns
	R: radius of cylinder in microns
	(Subscript[r, 1]): initial radius for calculation in microns
	(Subscript[r, 2]): final radius for calculation in microns
	(Subscript[r, inc]): increment of radius for calculation in microns
	(Subscript[n, in]): refractive index inside cylinder
	(Subscript[n, out]): refractive index outside cylinder

If (Subscript[λ, 1])≠(Subscript[λ, 2]) then the function computes wavevector vs. wavelength for a cylinder of radius R using wavelength increment (Subscript[λ, inc]).  If (Subscript[λ, 1])=(Subscript[λ, 2]), but (Subscript[r, 1])≠(Subscript[r, 2]), then the function computes fields for a cylinder of radius R from (Subscript[r, 1]) to (Subscript[r, 2]) by increment (Subscript[r, inc]).  The user can either set the refractive index to a specific (complex) value, or the user can set it 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.
	-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

λinc = 0.01;         (*  increment of wavelength  *)
R1 =0.1;               (*  initial radius of cylinder in microns  *)
R2=0.008;                    (*  final radius of cylinder in microns  *)
Rinc=0.0005;            (*  increment of radius  *)
r1 = 0.2;                  (*  initial radius in microns for calculation  *)
r2 = 0.2;                  (*  final radius in microns for calculation  *)
rinc = 0.0;              (*  increment of radius for calculation  *)
nout = 1;                   (*  refractive index of air around wire  *)
ng=1.1;                    (*  initial guess for SP refractive index  *)

Off[FindRoot::lstol];
λ1 = 0.532;                  (*  initial wavelength  *)
λ2 = 0.532;                  (*  final wavelength  *)
nin = -6;                   (* refractive index of gold wire  *)  
kz1= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];
nin = -7;                   (* refractive index of silver wire  *)  
kz2= cylinder[λ1, λ2, λinc, r1, r2, rinc,R1,R2,Rinc, nin, nout,ng];

ex1a=ListLinePlot[{kz1[[All,4]],kz2[[All,4]]}, PlotLegend -> {"Au", "Ag"},PlotLabel->"532 nm", LegendPosition -> {-0.7, 0.3}, LegendSize -> {0.25,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Propagation Distance (μm)"}, Evaluate[option]]
ex1aBW=ListLinePlot[{kz1[[All,4]],kz2[[All,4]]}, PlotLegend -> {"Au", "Ag"},PlotLabel->"532 nm", LegendPosition -> {-0.4, 0.3}, LegendSize -> {0.25,0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Propagation Distance (μm)"}, Evaluate[optB]];

Export[FileNameJoin[{base,"Fig8.ex1a.JPG"}],ex1aBW,ImageResolution->600];
On[FindRoot::lstol];

Off[FindRoot::lstol];
λ1 = 0.820;                  (*  initial wavelength  *)
λ2 = 0.820;                  (*  final wavelength  *)
nin = -6;                   (* refractive index of gold wire  *)  
kz1 = cylinder[λ1, λ2, λinc, r1, r2, rinc, R1, R2, Rinc, nin, nout, ng];
nin = -7;                   (* refractive index of silver wire  *)  
kz2 = cylinder[λ1, λ2, λinc, r1, r2, rinc, R1, R2, Rinc, nin, nout, ng];
ex1b = ListLinePlot[{kz1[[All, 4]], kz2[[All, 4]]}, PlotLegend -> {"Au", "Ag"}, PlotLabel -> "820 nm",
  LegendPosition -> {-0.7, 0.3}, LegendSize -> {0.25, 0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Propagation Distance (μm)"},  Evaluate[option]]
ex1bBW = ListLinePlot[{kz1[[All, 4]], kz2[[All, 4]]}, PlotLegend -> {"Au", "Ag"}, PlotLabel -> "820 nm",
  LegendPosition -> {-0.4, 0.3}, LegendSize -> {0.25, 0.2}, LegendShadow -> {.01, -.01}, FrameLabel -> {"R (μm)", "Propagation Distance (μm)"},  Evaluate[optB]];
On[FindRoot::lstol];

Export[FileNameJoin[{base,"Fig8.ex1b.JPG"}],ex1bBW,ImageResolution->600];
Clear[λ1, λ2, λinc, r1, r2, rinc, R1, R2, Rinc, nin, nout, ng, kz1, kz2];

λ1 = 0.633;                 (*  initial wavelength  *)
λ2 = 0.633;                  (*  final wavelength  *)
λinc = 0;                     (*  increment of wavelength  *)
R1 =0.5;                       (*  inner radius of cylinder in microns  *)
R2 =0.55;                       (*  outer radius of cylinder in microns  *)
r1 = 0.001;                  (*  initial radius in microns for calculation  *)
r2 = 1;                  (*  final radius in microns for calculation  *)
rinc = 0.001;              (*  increment of radius for calculation  *)
n1 = 1.5;               (* refractive index of region inside tube  *)  
n2 = -2;                  (* refractive index of silver tube  *)  
n3 = 1;                   (*  refractive index of air outside tube  *)
ng=1.6;                   (*  initial guess for SP refractive index  *)
fields1= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,+1]; 
ng=1.1;                   (*  initial guess for SP refractive index  *)
fields2= tube[λ1,λ2,λinc,r1,r2,rinc,R1,R2,n1,n2,n3,ng,+1];

ex2a=ListLinePlot[{fields1[[All,1]],fields1[[All,2]],fields1[[All,3]]}, PlotLegend -> {"Re((Subscript[E, z]))","Im((Subscript[E, ρ]))","|E|"}, LegendPosition -> {0.05, 0.15}, LegendSize -> {0.5,0.35}, LegendShadow -> {.01, -.01},FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]];
ex2aBW=ListLinePlot[{fields1[[All,1]],fields1[[All,2]],fields1[[All,3]]}, PlotLegend -> {"Re((Subscript[E, z]))","Im((Subscript[E, ρ]))","|E|"}, LegendPosition -> {0.1, 0.2}, LegendSize -> {0.5,0.35}, LegendShadow -> {.01, -.01},FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]];

ex2b=ListLinePlot[{fields2[[All,1]],fields2[[All,2]],fields2[[All,3]]}, FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[option]];
ex2bBW=ListLinePlot[{fields2[[All,1]],fields2[[All,2]],fields2[[All,3]]}, FrameLabel -> {"R (μm)", "Field Amplitudes"}, GridLines -> Automatic, Frame -> True, PlotRange->All, Evaluate[optB]];
ex2=GraphicsRow[{ex2a,ex2b},ImageSize->gr2,Spacings->{-30,30},Frame->True]
ex2BW=GraphicsRow[{ex2aBW,ex2bBW},ImageSize->gr2,Spacings->{-30,30},Frame->True];

Export[FileNameJoin[{base,"Fig8.ex2.JPG"}],ex2BW,ImageResolution->600];
Clear[λ1,λ2,λinc,R1,R2,r1,r2,rinc,n1,n2,n3,ng,fields1,fields2];

Figures

8.2  Propagating surface plasmons on metallic wires of circular cross section

8.3  Propagating surface plasmons on metallic wires of noncircular cross section

8.4  Propagating surface plasmons on hollow cylindrical waveguides

8.5  Propagating surface plasmons on hollow cylindrical shells

8.6  Excitation of surface plasmons on nanowires with plane waves

8.6.1  General solution

8.6.2  Nonpropagating surface plasmons

8.6.3  Scattering and absorption coefficients

8.7  Nanowires with noncircular cross sections

Exercises

1. How does the propagation distance for a silver nanowire compare to that for a gold nanowire?  Plot the propagation distance as a function of cylinder radius from 10 to 100 nm at a wavelength of 532 nm and 820 nm.

ex1a

ex1b

2. An example was given (Fig. 8.29) of two n = 0 modes for a cylindrical tube at the same frequency, one of which exhibited the same surface charge polarity on the inner and outer surfaces and one of which exhibited opposite surface charge polarity.  Consider the example of a silver tube with an inner diameter of 1 μm and an outer diameter of 1.1 μm at a free space wavelength of 633 nm.  Plot the fields for the first two modes as a function of radial position (hint: for one mode try an initial SP effective index of 1.6 and for the other mode try an initial SP index of 1.1).  Are the surface charges of the same sign or opposite sign on the inner and outer surfaces of the tube for each mode?

ex2

References

[1] C. A. Pfeiffer, E. N. Economou and K. L. Ngai.  Surface polaritons in a circularly cylindrical interface: surface plasmons.  Phys. Rev. B 10 (1974) 3038.
[2] J. Takahara, S. Yamagishi, H. Taki, A. Morimoto and T. Kobayashi.  Guiding of a one-dimensional optical beam with nanometer diameter.  Opt. Lett. 22 (1997) 475.
[3] H. Khosravi, D. R. Tilley and R. Loudon.  Surface polaritons in cylindrical optical fibers.   J. Opt. Soc. Am A 8 (1991) 112.
[4] U. Schröter and A. Dereux.  Surface plasmon polaritions on metal cylinders with dielectric core.  Phys. Rev. B 64 (2001) 125420.
[5] C. Miziumski.  Utilization of a cylindrical geometry to promote radiative interaction with slow surface excitations.  Phys. Lett. A 40 (1972) 187.
[6] P. Berini.  Plasmon-polariton modes guided by a metal film of finite width.  Opt. Lett. 24 (1999) 1011.
[7] P. Berini.  Plasmon-polariton waves guided by thin lossy metal films of finite width:  bound modes of symmetric structures.  Phys. Rev. B 61 (2000) 10484.
[8] I. Breukelaar, R. Charbonneau and P. Berini.  Long-range surface plasmon-polariton mode cutoff and radiation in embedded strip waveguides.  J. Appl. Phys. 100 (2006) 043104.
[9] R. Zia, M. D. Selker and M. L. Brongersma.  Leaky and bound modes of surface plasmon waveguides.  Phys. Rev. B 71 (2005) 165431.
[10] R. M. Dickson and L. A. Lyon.  Unidirectional Plasmon Propagation in Metallic Nanowires.  J. Phys. Chem. B 104 (2000) 6095.
[11] J-C. Weeber, A. Dereux, C. Girard, J. R. Krenn and J-P. Goudonnet.  Plasmon polaritons of metallic nanowires for controlling submicron propagation of light.  Phys. Rev. B 60 (1999) 9061.
[12] J. R. Krenn, B. Lamprecht, H. Ditlbacher, G. Schider, M. Salerno, A. Leitner and F. R. Aussenegg.  Non-diffraction-limited light transport by gold nanowires.  Europhys. Lett. 60 (2002) 663.
[13] B. Lamprecht, J. R. Krenn, H. Ditlbacher, M. Salerno, N. Felidj, A. Leitner, F. R. Aussenegg and J. C. Weeber.  Surface plasmon propagation in microscale metal stripes.  Appl. Phys. Lett. 79 (2001) 51.
[14] H. Ditlbacher, A. Hohenau, D. Wagner, U. Kriebig, M. Rogers, F. Hofer, F. R. Aussenegg and J. R. Krenn.  Silver nanowires as surface plasmon resonators.  Phys. Rev. Lett. 95 (2005) 257403.
[15] R. Zia, J. A. Schuller and M. L. Brongersma.  Near-field characterization of guided polariton propagation and cutoff in surface plasmon waveguides.  Phys. Rev. B 74 (2006) 165415.
[16] T. Nikolajsen, K. Leosson, I. Salakhutdinov and S. I. Bozhevolnyi.  Polymer-based surface-plasmon-polariton stripe waveguides at telecommunication wavelengths.  Appl. Phys. Lett. 82 (2003) 668.
[17] T. Nikolajsen, K. Leosson and S. Bozhevolni.  Surface plasmon polariton based modulators and switches operating at telecom wavelengths.  Appl. Phys. Lett. 85 (2004) 5833.
[18] N. C. J. van der Valk and P. C. M. Planken.  Effect of dielectric coating on terahertz surface plasmon polaritons on metal wires.  Appl. Phys. Lett. 87 (2005) 071106.
[19] B. Steinberger, A. Hohenau, H. Ditlbacher, A. L. Stepanov, A. Drezet, F. R. Aussenegg, A. Leitner and J. R. Krenn.  Dielectric stripes on gold as surface plasmon waveguides.  J. Appl. Phys. 88 (2006) 094104.  
[20] J. Q. Lu and A. A. Maradudin.  Channel plasmons.  Phys. Rev. B 42 (1990) 11159.
[21] I. V. Novikov and A. A. Maradudin.  Channel polaritons.  Phys. Rev. B 66 (2002) 035403.
[22] D. F. P. Pile and D. K. Gramotnev.  Plasmonic subwavelength waveguides:  next to zero losses at sharp bends.  Opt. Lett. 30 (2005) 1186.
[23] S. Bozhevolyni, V. S. Volkov, E. Devaux, J-Y. Laluet and T. W. Ebbesen.  Channel plasmon subwavelength waveguide components including interferometers and ring resonators.  Nature 440 (2006) 508.
[24] E. C. Jordan and K. G. Balmain.  Electromagnetic Waves and Radiating Systems (Englewood Cliffs, Prentice-Hall, 1968) Chap. 8.
[25] L.Novotny and C. Hafner.  Light propagation in a cylindrical waveguide with a complex, metallic, dielectric function.  Phys. Rev. E 50 (1994) 4094.
[26] H. Shin, P. B. Catrysse and S. Fan.  Effect of the plasmonic dispersion relation on the transmission properties of subwavelength cylindrical holes.  Phys. Rev. B 72 (2005) 85436.
[27] M. Abramowitz and I. A. Stegun.  Handbook of Mathematical Functions (New York, Dover, 1972).
[28] C. F. Bohren and D. R. Huffman.  Absorption and Scattering of Light by Small Particles (New York, Wiley, 1998) Chap. 8.
[29] R. D. Birkhoff, J. C. Ashley, H. H. Hubbell, Jr., and L. C. Emerson.  Light scattering from micron-size fibers.   J. Opt. Soc. Am. 67 (1977) 564.
[30] J. P. Kottmann, O. J. F. Martin, D. R. Smith and S. Schultz.  Dramatic localized electromagnetic enhancement in plasmon resonant nanowires. Chem. Phys. Lett. 341 (2001) 1.
[31] J. P. Kottmann, O. J. F. Martin, D. R. Smith and S. Schultz.  Plasmon resonances of silver nanowires with a nonregular cross section.  Phys. Rev. B 64 (2001) 235402.
[32] J. P. Kottmann and O. J. F. Martin.  Influence of the cross section and the permittivity on the plasmon-resonance spectrum of silver nanowires.  Appl. Phys. B 73 (2001) 299.
[33] Q. Xu, J. Bao, F. Capasso and G. M. Whitesides.  Surface plasmon resonances of free-standing gold nanowires fabricated by nanoskiving.  Angew. Chem. 118 (2006) 3713.
[34] V. Giannini and J. A. Sánchez-Gil.  Calculations of light scattering from isolated and interacting metallic nanowires of arbitrary cross section by means of Green’s theorem surface integral equations in parametric form.  J. Opt. Soc. Am. A 24 (2007) 2822.