Structural Characterisation and Optical Properties of Nanoporous Gold Films

Structural Characterisation and Optical Properties of Nanoporous Gold Films

Anton Segerkvist December 3, 2014

Abstract

Nanoporous metal films have many applications in a great variety of scientific fields, and especially nanoporous gold films have many applica- tions in green nanotechnology. Hence, structural and optical properties of such materials were investigated. Local density functions and local per- colation functions were calculated by using scaning electron micrographs and the optical properties of the films were calculated using the Hilfer equation. The results are presented in the report as graphs and show how the materials optical properties depend on the structure of the gold films.

Introduction

Nanoporous metal films have great potential in many scientific fields. Especially gold nanoporous films have many applications in green nanotechnology [5, 6], in devices such as gas sensors, solar cells, light-emiting diodes and in detection of cancer and many more areas [5].

In this report, investigations of thin nanoporous gold films with properties that goes from island structure to homogeneous films were done. These films have been deposited using DC magnetron sputtering. The images used in the report have been taken with scanning electron microscopy. A LEO 1550 FEG with inlense detection was used, and the images were processed with MATLAB image analysis, including an improved version of the image precessing tool [6].

The optical properties of these nanopurous materials depend partially on their size and form, but also on variations in the porosity. These variations can be analysed with the help of local porosity theory in terms of local density functions and local connectivity functions.

E↵ective medium theory is a way of describing the optical and electrical properties of nanoparticles mixed with another material, and by using the Hilfer equation (17) one can include local denstiy functions and local connectivity functions into the relation [1, 2].

1

Theory

Topological Structure

The nanoporous gold films that are investigated are all placed upon glass films substrates of given thickness. The medium sorrounding the glass and the gold is air and the small holes that might appear in the gold films are hence also filled with air.

Optical Properties

The response in a material, when an electric field is applied, can be described by the dielectric permittivity, equation (1), where ✏

and ✏

are defined by equation (2) and (3) respectively.

✏ = ✏

+ i✏

(1)

✏

= n

k

(2)

✏

= 2nk (3)

The optical quantities, n (refractive index) and k (extinction coefficient), are related to the optical properties which is described by transmitance, reflectance and absorbtion. From experimental data of transmitance and reflectance, the optical quantities n and k are derivable, using Fresnel formalism. Equally, from n and k the permittivity ✏ is obtainable.

In this work, however, the permittivity ✏

has also been calculated using the Hilfer equation (17), where this is then compared to the experimental values of

✏

.

In order to obtain the theoretical optical properties, the calculated permit- tivities need to be converted to reflection, transmission and absorbtion. Due to the time-consuming computations only the transmission is computed. The results are obtained from Maxwell’s equations, but only the results are stated[4].

First of all the refractive index is needed and is given by equation (4) and (5) for the real and imaginary part respectively.

n =

s p ✏

+ ✏

+ ✏

2 (4)

k =

s p ✏

+ ✏

✏

2 (5)

From this the transmission of the entire air, gold and glass compund can

be calculated. The equation below requires a bit more notation and we start

by defining them. Most quantities below are already described, but now index

notation is introduced to distinguish between the di↵erent media. An index 0

refers to the air, an index 1 refers to the gold film and an index 2 refers to the glass on which the gold is placed. The d

refers to the thickness of the gold film and d

refers to the thickness of the glass on which the gold resides. The other quantities are derived quantites, and their physical relevance is not essensial for the computations.

r

= (n

n

) + (k

k

)

(n

+ n

)

+ (k

+ k

)

(6) r

= 2(n

k

n

k

)

(n

+ n

)

+ (k

+ k

)

(7) R

= (r

)

+ (r

)

(8)

j

= 4⇡ n

d

(9)

↵

= 4⇡ k

(10)

x

= exp( ↵

d

) (11)

h = ((1 + r

)

+ (r

)

)((1 + r

)

+ (r

)

)((1 + r

)

+ (r

)

)x

x

(12)

u = 1 + R

R

x

R

R

x

x

R

R

x

+ 2(r

r

(1 R

x

) r

r

(1 + R

x

))x

cos(

) + 2(r

r

(1 R

x

) + r

r

(1 R

x

))x

sin(

)

(13) Now the transmission is given from equation (14) [9].

T = h

u (14)

Local Porosity Theory

Local porosity theory is used to calculate the permittivity of the nanoporous gold films. It’s assumed that that the area fraction of the film can be used to approximate the volume fraction of the films. In order to calculate the properties of the gold films, some structural data is needed. Hence, two quantities called the local density function and the local percolation function are used. The local density function describes fluctuations of the porosity the films, while the local percolation function describes how well the di↵erent dense sections of the films are connected. These two quantities can be combined in the Hilfer Equation (17) to calculate the e↵ective permittivity [1, 2, 3].

3

Converting Electroscopy Images to Black-White Images

MATLAB is the software that is used for all the calculations, but before any calculations can be done on the electron microscopy images, they need to be converted to black and white images. The converted images in this report are given. There is some image processing behind it and the code in the appendix o↵eres limited support for converting these. The way it’s done in the appendix is that the images are processed by a linear smoothing filter and then thresholded with an appropriate value. This, however, means that di↵erent sections of the films will not be equally intensified due to the initial conditions of the images.

Also, the images given, are slightly under thresholded, which will show up in the result section. A better method for processing these images are needed.

However, with some fiddling with the numbers in the code one can find certain sections of the films that have the appropriate properties.

Calculating Local Density Functions and Local Perculation Functions

The local density is strightforward to calculate. We divide the image into small squares of equal side length. The number of white pixels, which corresponds to the gold particles, in each of the squares are counted and a distribution of the number of squares with a number of white pixel are created. Then we divide the distribution by the number of small boxes. This distribution is the local density function. It’s possible from this distribution to calculate the total area fraction of the film. This is done by calculating the expectation fraction, and since the system is discrete the sum is taken instead, as is seen in equation (15) where f denotes the mentioned quantity.

f = Z

0

f µ(f )df ! f = X

f

f µ(f ) (15)

The local density function is calculated in a range of di↵erent box sizes in order to find the box size which gives the maximal information entropy. The box size with maximal information entropy is expected to give the best description of the optical properties.

After that it’s possible to calculate the local percolation function and it’s cal- culated using a depth-first search algorithm. In each of the squares it’s checked whether the opposing side of each box is connected the other side of the box trough the conducting medium. Then statistical data is created from this. A 1 means that every side is connected to it’s other side trough the conducting medium and a 0 means that no side is connected to it’s opposite side.

Information Entropy

The information entropy is required in order to choose square size. The infor-

mation entropy is defined through the Gibbs-Shannon entropy and is defined

through that of equation (16) and is stright forward to calculate [1, 2].

H(L) = 1 log(L

+ 1)

X

k

µ(f

; L)log(µ(f

; L)) (16)

Hilfer Equation

The Hilfer equation is used to calculate the e↵ective permittivity of the film.

However, equation is complex in both it’s domain and range. This means a more advanced method for solving it is needed. The permittivities are frequency de- pendent and thus a whole set of equations needs to be solved. In equation (17) f denotes the area fraction of the gold and is an integration variable, varying from zero to one. The di↵erent permittivities are the e↵ective permittivity,

✏

(!), the conducting material permittivity, ✏

(f ; !), and the non-conducting material permittivity, ✏

(f ; !). The di↵erent perculation functions are the conducting material perculation function,

(f ), and the non-conducting mate- rial perculation function,

(f ). Finally the local density function is denoted by µ(f ) [1, 2, 3].

Z

[

(f ) ✏

(f ; !) ✏

(!)

✏

(f ; !) + 2✏

(!) +

(f ) ✏

(f ; !) ✏

(!)

✏

(f ; !) + 2✏

(!) ]µ(f )df = 0 (17) The filling factor dependent permittivities are given by equations (18) and (19) where the permittivities on the right hand side are the usually tabulated frequency dependent permittivity values [1, 8].

✏

(f ; !) = ✏

(!) · (1 1 f

1 1 ^{✏air (!)}

✏gold(!)

1

3

f ) (18)

✏

(f ; !) = ✏

(!) · (1 f

1 1 ^✏gold(!)

✏air (!)

1

3

(1 f ) ) (19)

The values for ✏

and ✏

need to take into account the increased grain boundary scattering in thin films, which leads to bulk conductivities consider- ably lower than values tabulated for polycrystalline gold [6]. For these values, experimental data from a homogenous gold film, using DC magnetron sputtering were produced. The actuall data for these can be found in the appendix.

Muller’s Method

This method is capable of generating complex solutions, and is based on find- ing zeroth roots to a second, or higher, order polynomial, and is the method that is used to solve the Hilfer equation. However, since a second order polyno- mial equation is the simplest to solve, this type of polynomial is choosen. The method is based on choosing three initial guess values of the function which we want to solve, and then interpolating between these. Then the zeroth roots are

5

calculated of the interpolated curve and the one root that makes the change in the recurence value the smallest is choosen. The resulting recurence equation is that of equation (22) where the denominator is choosen such that it’s absolute value is as large as possible [7]. The functions, g, with brackets after them are their divided di↵erence defined by equation (20) and (21).

g[x

, x

] = g(x

) g(x

) x

x

(20)

g[x

, x

, x

] = g[x

, x

] g[x

, x

] x

x

(21)

x

= x

2g(x

)

w ± p

w

4g(x

)g[x

, x

, x

] (22) w = g[x

, x

] + g[x

, x

] g[x

, x

] (23) The Hilfer equation (17) is approximated with a linear approximation, where the distance between the di↵erent points, can be neglected, since that won’t af- fect where the roots are found. In this case the di↵erent e↵ective permittivities are the di↵erent points. This yields the equation (24) where the di↵erent quan- tities are already defined.

g(✏

(!)) = X

i

[

(f

) ✏

(f

; !) ✏

(!)

✏

(f

; !) + 2✏

(!) +

(f

) ✏

(f

; !) ✏

(!)

✏

(f

; !) + 2✏

(!) ]µ(f

) (24)

The method has a convergence order of about 1.9, but the method does

not garuantee convergence. Neither is there a good way to choose the initial

guessing values, but these are taken to be around the origin with a given search

radius and an imaginary part larger than zero.

Results

Scanning Electron Micrographs

The images investigated are shown in figure 1.

(a) Gold film 1: d

= 1.4nm (b) Gold film 2: d

= 3.0nm

(c) Gold film 3: d

= 3.8nm (d) Gold film 4: d

= 4.6nm

Figure 1: Scanning electron micrographs. About 110 pixels corresponds to 200nm as given by the equipment specifications.

The films are orded in increasing dense order. It’s seen from these images that the gold particles are forming connected threads as the density increases and in the last few images they become more homogenous. The white sections denote the gold and dark sections denote the air.

Local Density Function, Local Perculation Function and Information Entropy Function

From the images in figure 1 the simulation results in the local density functions in figure 2. From these graphs it’s possible to estimate the average density of the entire film and compare to the experimental values. These values can be found in table 1.

7

(a) Gold film 1

(b) Gold film 2

(c) Gold film 3

(d) Gold film 4

Figure 2: Local density analysis of the di↵erent films at maximal information

entropy. The horozontal axis represents the filling factor and the vertical rep-

resents the density fraction.

It’s clear that the peak of the curve keeps shifting to the right in the images, which is expected as the density keeps increasing. Table 1 shows the total density of the entire films as determined from experimental procedures, denoted by f

, agains the total film density as calculated from the images using equation (15), denoted by f . It’s seen from there that the images have been converted with a slightly to low threshold value.

f

f

Film 1 0.363 0.423 Film 2 0.443 0.480 Film 3 0.563 0.604 Film 4 0.728 0.770

Table 1: Experimental [6] vs calculated total film density

Density functions have been calculated at the maximum information entropy.

The information entropy as a function of the box side length can be seen in figure 3.

9

(a) Gold film 1

(b) Gold film 2

(c) Gold film 3

(d) Gold film 4

Figure 3: Information entropy analysis of the di↵erent films. The horizontal axis

denotes the box width and the vertical axis denotes the information entropy.

It’s seen from figure 3 that the curve maintains it’s shape. The peak, how- ever, does not maintain it’s position. The box size resulting in maximal infor- mation entropy for the di↵erent films are documented in table 2.

Optimal Box Size [px] Optimal Box Size [nm]

Film 1 12 21.82

Film 2 8 14.55

Film 3 9 16.36

Film 4 11 20.00

Table 2: Box length resulting in maximal information entropy

The percolation functions for the di↵erent films are also calculated at the box width resulting in maximal information entropy. These are found in figure 4 and show how well the these films are connected.

11

(a) Gold film 1

(b) Gold film 2

(c) Gold film 3

(d) Gold film 4

Figure 4: Perculation function for the di↵erent films at maximal information

entropy. The horizontal axis denotes the filling factor and the vertical axis

denotes the percolation.

As the density of the film increases, a more gradual increase in the perco- lation is observed. That means, as the film density increases, the perculation threshold increases slightly.

E↵ective Permittivity and Transmismittance

From this point on the calculations get more heavy. That means that the data dealt with here are not necessarily complete. These simulations can take days to complete and do not necessarily need to finish what so ever, meaning if the initial conditions are not set properly, convergence may never happen and the simulation continuous for ever.

The simulations here are done on film two and three from figure 1. E↵ec- tive permittivity is showed first. This information comes dirrectly from solving equation (17).

13

(a) Gold film 2

(b) Gold film 3

Figure 5: Real part of permittivity for two films where filled lines are experi-

mental [6] data and dotted lines are calculated data.

(a) Gold film 2

(b) Gold film 3

Figure 6: Imaginary part of permittivity for two films where filled lines are experimental data and dotted lines are calculated data.

It’s seen from figure 5 and figure 6 that we have quite good agreement at short wavelengths. At longer wavelengths, however, a di↵erence becomes more clear. It’s worth noting that the code might generate false solutions due to the convergence condition. However, this is highly unlikely but might be a source of errors in some simulations.

From these curves we can use equation (14) to calculate the transmittance curves and those are shown in figure 7.

15

(a) Gold film 2

(b) Gold film 3

Figure 7: Transmittance curves for two gold films where dotted points are cal- culated values and filled lines are experimental values [6]

Some consistency in the experimental curve and the simulated values are observed. If the films are more carefully thresholded, one would expect smaller transmittance values and thus further consistency. Although the experimen- tal and theoretical values di↵er, the method shows promise, but need further development.

Improvements

The method used to solve the Hilfer equation is very slow and takes days to

give reasonable amount of points to compare to data. Hence a better method

is needed to give more results in a shorter time interval. There are lots of

room for parallel computations in the simulation and this could be done using

threading or alternatively distributed systems. This could significanly increase

performance in the computations. Another method for aquiring more data in the structural analysis, is to shift the boxes that the electromicroscopy images is divided into. This could potentially give less statistical errors.

Conclusions

The structural characterisation in the project is showing good agreement to what is expected. Solving the Hilfer equation, however, turns out to give some strange results. From the looks of it, there is good correspondence at high frequencies but less at lower frequencies. This could be due to one of many problems, one which has not been established. There could be something in the theory that is not taken in to consideration that is present in the experimental data. There could be something wrong in the code. There could be multiple solutions to the equation and for some reason a branch becomes the prefered solution. Further analysis of the problem is definitely needed.

References

[1] C. Andraud* and J. Lafait. Entropic model for the optical properties of heterogeneous media: Validation on granular gold films near percolation.

PHYSICAL REVIEW B, 57(20), 1998.

[2] C. Andraud*, J. Lafaita, A. Beghdadi, and J. Peiro. Experimental validation on granular metallic films of an optical model based on the entropic analysis of their morphology. PHYSICA A, pages 133–137, 1997.

[3] C. Andrauda, A. Beghdadi, E. Haslundc, R. Hilferc’d’*, J. Lafaita, and B. Virgin. Local entropy characterization of correlated random microstruc- tures. PHYSICA A, pages 307–318, 1997.

[4] Philip Hofmann. Solid state physics: an introduction. Wiley-VCH, page 175, 2008.

[5] Pia Lans˚ aker. Gold-based nanoparticles and thin films: Applications to green nanotechnology. 2012.

[6] Pia Lans˚ aker, E. Tuncer, I. Valyukh, H. Arwin, C. G. Granqvist, and G. A.

Niklasson. Spectral density analysis of thin gold films: Thickness and struc- ture dependence of the optical properties. pages 443–447, 2013.

[7] David E. Muller. A method for solving algebraic equations using an auto- matic computer. JSTOR, pages 208–215, 1956.

[8] T. Rage R. Hilfer and B. Virgin. Local perculation probabilities for a natural sandstone. arXiv, 1996.

[9] R. Swanpoel. Transmission and reflection of an absorbing thin film on an absorbing substrate. S.-Afr. Tydskr., pages 148–156, 1989.

17

Appendix

Dielectric permittivity as a function of wavelength for a homogenous gold film with thickness d

= 7.9

lamda(nm) Eps1 Eps2

250 -2,04904393 4,02671161

251,01 -2,03822497 4,0832348

252,03 -2,02637399 4,13924459

253,06 -2,01350611 4,19469669

254,1 -1,999638 4,2495481

255,14 -1,98478776 4,30375706

256,2 -1,96897605 4,35727926

257,26 -1,95222426 4,41007423

258,33 -1,9345564 4,46209907

259,41 -1,91599797 4,51331231

260,5 -1,89657725 4,56367043

261,6 -1,87632323 4,61313368

262,71 -1,85526722 4,66166148

263,83 -1,83344306 4,70921267

264,96 -1,81088561 4,75574869

266,09 -1,7876333 4,80122857

267,24 -1,76372409 4,84561693

268,4 -1,73920058 4,88887428

269,57 -1,71410597 4,93096483

270,74 -1,68848558 4,97185363

271,93 -1,66238745 5,01150589

273,13 -1,63586183 5,04988797

274,34 -1,60896066 5,08696824

275,56 -1,58173761 5,12271708

276,79 -1,55425028 5,15710416

278,03 -1,52655696 5,19010284

279,28 -1,49871882 5,22168749

280,54 -1,4707989 5,25183491

281,82 -1,44286371 5,28052277

283,11 -1,4149811 5,30773212

284,4 -1,38722189 5,33344587

285,71 -1,35965829 5,3576505

287,04 -1,33236601 5,38033436

288,37 -1,30542262 5,40148936

289,72 -1,2789081 5,4211108

291,08 -1,25290472 5,43919763

292,45 -1,22749652 5,45575306

293,84 -1,20277016 5,47078429

295,24 -1,17881424 5,48430322

296,65 -1,1557183 5,49632725

298,08 -1,13357396 5,50687883 299,52 -1,11247447 5,51598612 300,97 -1,09251363 5,52368366 302,44 -1,0737855 5,53001273 303,92 -1,0563843 5,53502156 305,42 -1,0404039 5,53876575 306,93 -1,02593762 5,54130856 308,46 -1,01307613 5,5427216 310 -1,00190831 5,54308471 311,56 -0,992519365 5,54248646 313,13 -0,984990307 5,54102421 314,72 -0,979397105 5,5388043 316,33 -0,975808975 5,5359418 317,95 -0,97428807 5,53256083 319,59 -0,974887581 5,52879376 321,24 -0,977650896 5,52478127 322,92 -0,982610055 5,52067164 324,61 -0,989784531 5,51661994 326,32 -0,999180197 5,51278709 328,04 -1,01078742 5,50933901 329,79 -1,0245804 5,50644498 331,55 -1,04051604 5,50427625 333,33 -1,05853228 5,50300434 335,14 -1,07854831 5,50279896 336,96 -1,10046269 5,503826 338,8 -1,12415364 5,50624506 340,66 -1,14947916 5,51020711 342,54 -1,17627518 5,51585153 344,44 -1,2043589 5,5233038 346,37 -1,23352756 5,53267247 348,31 -1,26355868 5,54404581 350,28 -1,29421494 5,55749018 352,27 -1,32524193 5,5730457 354,29 -1,35637387 5,59072502 356,32 -1,38733439 5,61051002 358,38 -1,41783968 5,63234943 360,47 -1,44760376 5,65615824 362,57 -1,47633979 5,68181432 364,71 -1,50376598 5,70915866 366,86 -1,52960928 5,73799412 369,05 -1,55360976 5,76808499 371,26 -1,57552693 5,79915943 373,49 -1,59514164 5,83090687 375,76 -1,61226387 5,86298368 378,05 -1,62673544 5,8950128 380,37 -1,63843509 5,92658737

19

382,72 -1,64728287 5,95727492 385,09 -1,65324365 5,98662084 387,5 -1,65633056 6,01415313 389,94 -1,65660787 6,03938785 392,41 -1,65419312 6,06183451 394,9 -1,64925859 6,08100258 397,44 -1,64203234 6,09640636 400 -1,63279766 6,10757314 402,6 -1,62189285 6,11404805 405,23 -1,60970938 6,11540119 407,89 -1,59668926 6,11123373 410,6 -1,58332179 6,10118335 413,33 -1,57014038 6,08493048 416,11 -1,55771681 6,06220268 418,92 -1,54665629 6,03277874 421,77 -1,53759182 5,99649269 424,66 -1,53117801 5,95323774 427,59 -1,52808394 5,90296562 430,56 -1,52898677 5,84569244 433,57 -1,53456448 5,78149318 436,62 -1,54548896 5,71050658 439,72 -1,5624194 5,63292978 442,86 -1,58599555 5,54901817 446,04 -1,61683184 5,45908111 449,28 -1,65551058 5,36348096 452,55 -1,70257898 5,26262297 455,88 -1,75854251 5,15695586 459,26 -1,82386112 5,04696509 462,69 -1,89894811 4,93316375 466,17 -1,98416746 4,81608778 469,7 -2,07982985 4,696293 473,28 -2,18619554 4,57434379 476,92 -2,30347542 4,45080675 480,62 -2,43182393 4,32625282 484,38 -2,57135172 4,20123982 488,19 -2,72212194 4,07631353 492,06 -2,88415494 3,952001

496 -3,0574319 3,82880597

500 -3,24189658 3,70720633

504,07 -3,43746411 3,58764766

508,2 -3,64402365 3,47054198

512,4 -3,86144068 3,35626694

516,67 -4,08956955 3,24516017

521,01 -4,328251 3,13752255

525,42 -4,57731949 3,03361636

529,91 -4,83661083 2,93366445

534,48 -5,1059612 2,83785323 539,13 -5,38522005 2,74633068 543,86 -5,67425162 2,65920947 548,67 -5,9729266 2,5765723 553,57 -6,28114945 2,4984673 558,56 -6,59883737 2,42491736 563,64 -6,92594431 2,35591746 568,81 -7,26244149 2,29144248 574,07 -7,6083444 2,23144495 579,44 -7,96369225 2,17586243 584,91 -8,32856179 2,12461793 590,48 -8,70306504 2,07762345 596,15 -9,08735647 2,03478202 601,94 -9,48161693 1,99599213 607,84 -9,88607352 1,96114802 613,86 -10,3009931 1,93014284 620 -10,7266755 1,90287121 626,26 -11,1634632 1,87923039 632,65 -11,6117343 1,85912234 639,18 -12,0719091 1,84245495 645,83 -12,5444534 1,82914334 652,63 -13,0298609 1,81911153 659,57 -13,5286738 1,81229276 666,67 -14,0414756 1,80863061 673,91 -14,5688919 1,80807971 681,32 -15,1115899 1,81060628 688,89 -15,6702754 1,81618855 696,63 -16,2457151 1,82481735 704,55 -16,8387121 1,83649629 712,64 -17,4501215 1,85124214 720,93 -18,0808478 1,86908503 729,41 -18,7318613 1,89006919 738,1 -19,4041777 1,91425256 746,99 -20,0988864 1,94170794 756,1 -20,817134 1,97252271 765,43 -21,5601383 2,00679957

775 -22,329199 2,04465753

784,81 -23,1256859 2,08623162 794,87 -23,9510659 2,13167481 805,19 -24,8068957 2,18115821 815,79 -25,6948218 2,23487166 826,67 -26,6166229 2,29302702 837,84 -27,5741616 2,35585614 849,32 -28,5694657 2,423617 861,11 -29,6046766 2,49659162 873,24 -30,6820815 2,57508942

21

885,71 -31,8041396 2,65945065 898,55 -32,973491 2,75004883 911,76 -34,1929677 2,84729356 925,37 -35,4655858 2,951632 939,39 -36,7946216 3,06355744 953,85 -38,1835892 3,1836104 968,75 -39,6362646 3,31238366 984,13 -41,1567402 3,45053065 1000 -42,7494347 3,59877039 1016,4 -44,4191138 3,75789385 1033,3 -46,1709426 3,92877399 1050,8 -48,0105424 4,11237712 1069 -49,9439895 4,30976944 1087,7 -51,9779095 4,52213401 1107,1 -54,1195039 4,75078208 1127,3 -56,3766449 4,99717306 1148,1 -58,7579117 5,26292966 1169,8 -61,2727093 5,54986402 1192,3 -63,9313287 5,85999956 1215,7 -66,7450603 6,19560119 1240 -69,726353 6,5592144 1265,3 -72,8888801 6,95369654

MATLAB Code

clear ; close all ;

% % I n p u t data thresh = 100;

% % Load the i m a g e

% i m a g e = i m r e a d ( ’ pic1 . tif ’) ;

% i m a g e = i m r e a d ( ’ pic2 . tif ’) ;

% i m a g e = i m r e a d ( ’ pic3 . tif ’) ;

% i m a g e = i m r e a d ( ’ pic4 . tif ’) ;

% i m a g e = i m r e a d ( ’ pic5 . tif ’) ;

% i m a g e = i m r e a d ( ’ pic6 . tif ’) ;

% i m a g e = i m r e a d ( ’ I Y 1 _ 0 2 _ a n a l y s i s . tif ’) ;

% i m a g e = i m r e a d ( ’ I Y 2 _ 0 4 _ a n a l y s i s . tif ’) ;

% i m a g e = i m r e a d ( ’ I Y 3 _ 1 5 _ a n a l y s i s . tif ’) ; image = imread ( ’ I Y 4 _ 0 8 _ a n a l y s i s . tif ’) ;

% i m a g e = i m r e a d ( ’ I Y 5 _ 0 6 _ a n a l y s i s . tif ’) ; [ height , width ] = size ( image ) ;

% i m s h o w ( i m a g e ) ;

% % Crop i m a g e

croped = imcrop ( image , [0 , 0 , width , height ]) ; [ cheight , cwidth ] = size ( croped ) ;

% i m s h o w ( c r o p e d ) ;

% % F i l t e r i m a g e

htf = f s p e c i a l ( ’ a v e r a g e ’) ; c r o p e d 2 = i m f i l t e r ( croped , htf ) ;

% { hold on ; s u b p l o t (1 ,2 ,1) ; imshow ( croped ) ; s u b p l o t (1 ,2 ,2) ; imshow ( c r o p e d 2 ) ; hold off ;

% }

% % R e m o v e u n w a n t e d data croped = c r o p e d 2 ; clear image ; clear c r o p e d 2 ;

% % C r e a t e b l a c k and w h i t e i m a g e imgbw = im2bw ( croped , thresh /255) ; imshow ( imgbw ) ;

% % S e l e c t p o r t i o n of i m a g e imshow ( imgbw ) ;

rect = g e t r e c t ;

imgbw = imcrop ( imgbw , rect ) ;

% % C r e a t e the r e q u i r e d grid

% [ grid ] = c r e a t e _ g r i d (1 , 1 , cwidth , cheight , b o x W i d t h ) ;

% % Plot the grid { hold on ;

% {

for i =1: length ( grid )

plot ([ grid (i ,1) grid (i ,3) ] , [ grid (i ,2) grid (i ,2) ] , ’b - ’) ; plot ([ grid (i ,1) grid (i ,1) ] , [ grid (i ,2) grid (i ,4) ] , ’b - ’) ; plot ([ grid (i ,1) grid (i ,3) ] , [ grid (i ,4) grid (i ,4) ] , ’b - ’) ; plot ([ grid (i ,3) grid (i ,3) ] , [ grid (i ,2) grid (i ,4) ] , ’b - ’) ; end

hold off ;

% }

% % C a l c u l a t e the l o c a l d e n s i t y p r o b a b i l i t y f u n c t i o n sy = [];

sx = [];

for i =2:50 b o x W i d t h = i

[~ , s ] = c a l c u l a t e _ l o c a l _ d a t a _ 1 ( imgbw , b o x W i d t h ) ; sx (1+ end ) = i ;

sy (1+ end ) = s ; end

% % C a l c u l a t e the data for o p t i m a l i n f o r m a t i o n e n t r o p y [~ , index ] = max ( sy ) ;

b o x W i d t h = sx ( index ) ;

f = l i n s p a c e (0 , 1 , b o x W i d t h ^2+1) ;

[ mu , la , s ] = c a l c u l a t e _ l o c a l _ d a t a ( imgbw , b o x W i d t h ) ;

% % C a l c u l a t e f i l l i n g f a c t o r for w h o l e film f_fem = sum ( f .* mu )

% % C a l c u l a t e e f f i c i e n t p e r m i t t i v i t y . cerror = 0.01;

s r a d i u s = 30;

i t e r a t i o n s = 50;

[ wavelength , ~] = p e r m i t t i v i t y _ a u () ; done = zeros (1 , length ( w a v e l e n g t h ) ) ; c o n v e r g e = zeros (1 , length ( w a v e l e n g t h ) ) ;

% %

while sum ( c o n v e r g e ) ~= length ( w a v e l e n g t h )

23

list = find (~ c o n v e r g e ) ; l e n g t h _ w a v e = length ( list ) ; values = zeros (4 , l e n g t h _ w a v e ) ;

values (2:4 ,:) = ( rand (3 , l e n g t h _ w a v e ) - 0.5) * 2 * s r a d iu s ...

+ rand (3 , l e n g t h _ w a v e ) * s r a d i u s * 1 i ; sum (~ c o n v e r g e )

% Use m u l l e r s m e t h o d and try to find c o n v e r g e n c e for i =1: i t e r a t i o n s

[ y0 ] = f z e r o _ f u n c ( values (4 ,:) , f , la , mu , list ) ; [ y1 ] = f z e r o _ f u n c ( values (3 ,:) , f , la , mu , list ) ; [ y2 ] = f z e r o _ f u n c ( values (2 ,:) , f , la , mu , list ) ; [ values (1 ,:) ] = m u l l e r _ r e c u r r e n c e ( y0 , y1 , y2 , ...

values (4 ,:) , ...

values (3 ,:) , ...

values (2 ,:) ) ; values (4 ,:) = values (3 ,:) ;

values (3 ,:) = values (2 ,:) ; values (2 ,:) = values (1 ,:) ; end

% C h e c k for c o n v e r g e n c e for i =1: l e n g t h _ w a v e

if imag ( values (1 , i ) ) >= 0 && ...

abs ( values (2 , i ) - values (3 , i ) ) < cerror c o n v e r g e ( list ( i ) ) = 1;

done ( list ( i ) ) = values (1 , i ) ; end

end end

% % Plot P e r m i t t i v i t y and O p t i c s [~ , temp ] = p e r m i t t i v i t y _ a i r () ; [~ , temp2 ] = r e f r a c t i o n _ g l a s s () ; wp = w a v e l e n g t h ( find ( c o n v e r g e ) ) ; e10 = real ( temp ( find ( c o n v e r g e ) ) ) ; e20 = imag ( temp ( find ( c o n v e r g e ) ) ) ; e1a = real ( done ( find ( c o n v e r g e ) ) ) ; e2a = imag ( done ( find ( c o n v e r g e ) ) ) ;

d (1) = 6.7 e -9; % gold t h i c k n e s s d (2) = 1e -3; % gl a s s t h i c k n e s s n = zeros (3 , length ( e10 ) ) ; k = zeros (3 , length ( e10 ) ) ;

n (1 ,:) = sqrt (( sqrt ( e10 .^2 + e20 .^2) + e10 ) / 2) ; k (1 ,:) = sqrt (( sqrt ( e10 .^2 + e20 .^2) - e10 ) / 2) ; n (2 ,:) = sqrt (( sqrt ( e1a .^2 + e2a .^2) + e1a ) / 2) ; k (2 ,:) = sqrt (( sqrt ( e1a .^2 + e2a .^2) - e1a ) / 2) ; n (3 ,:) = real ( temp2 ( find ( c o n v e r g e ) ) ) ;

k (3 ,:) = - imag ( temp2 ( find ( c o n v e r g e ) ) ) ; Tp = [];

for ki = 1: length ( wp ) r = zeros (3 ,3) ; for i =1:3

for j =1:3

r (i , j ) = (( n (i , ki ) - 1 i * k (i , ki ) ) - ( n (j , ki ) - 1 i * k (j , ki ) ) ) ...

./ (( n (i , ki ) - 1 i * k (i , ki ) ) + ( n (j , ki ) - 1 i * k (j , ki ) ) ) ; end

end

R = zeros (3 ,3) ;

for i =1:3 for j =1:3

R (i , j ) = ( real ( r (i , j ) ) ) .^2 + ( imag ( r (i , j ) ) ) .^2;

end end

delta (1:2) = 4 * pi * n (2:3 , ki ) ’ .* d (1:2) / wp ( ki ) ; alpha (1:2) = 4 * pi * k (2:3 , ki ) ’ / wp ( ki ) ;

x (1:2) = exp ( - alpha (1:2) .* d (1:2) ) ; u = 1 ...

+ R (1 ,2) .* R (2 ,3) .* x (1) .* x (1) ...

- R (1 ,2) .* R (3 ,1) .* x (1) .^2 .* x (2) .^2 ...

- R (2 ,3) .* R (3 ,1) .* x (2) .* x (2) ...

+ 2 * ( real ( r (1 ,2) ) .* real ( r (2 ,3) ) .* (1 - R (3 ,1) .* x (2) .^2) ...

- ( imag ( r (1 ,2) ) .* imag ( r (2 ,3) ) .* (1 + R (3 ,1) .* x (2) .^2) ) ) ...

.* x (1) .* cos ( delta (1) ) ...

+ 2 * ( real ( r (1 ,2) ) .* imag ( r (2 ,3) ) .* (1 + R (3 ,1) .* x (2) .^2) ...

+ ( imag ( r (1 ,2) ) .* real ( r (2 ,3) ) .* (1 - R (3 ,1) .* x (2) .^2) ) ) ...

.* x (1) .* sin ( delta (1) ) ;

h = ((1 + real ( r (1 ,2) ) ) .^2 + imag ( r (1 ,2) ) .^2) ...

.* ((1 + real ( r (2 ,3) ) ) .^2 + imag ( r (2 ,3) ) .^2) ...

.* ((1 + real ( r (3 ,1) ) ) .^2 + imag ( r (3 ,1) ) .^2) ...

.* x (1) .* x (2) ; Tp ( end +1) = h / u ; end

hold on ;

% T r a n s m i t a n c e 1 = load ( ’ T r a n s m i t a n c e s / AU1T . SP ’) ;

% T r a n s m i t a n c e 2 = load ( ’ T r a n s m i t a n c e s / AU2T . SP ’) ;

% T r a n s m i t a n c e 3 = load ( ’ T r a n s m i t a n c e s / AU3T . SP ’) ; T r a n s m i t a n c e 4 = load ( ’ T r a n s m i t a n c e s / AU4T . SP ’) ;

% T r a n s m i t a n c e 5 = load ( ’ T r a n s m i t a n c e s / AU5T . SP ’) ;

% T r a n s m i t a n c e 6 = load ( ’ T r a n s m i t a n c e s / AU6T . SP ’) ; plot ( wp *1 e9 , Tp *100 , ’* b ’) ;

ylabel ( ’ T r a n s m i s s i o n [%] ’) ; xlabel ( ’ Wave Length [ nm ] ’) ;

% plot ( T r a n s m i t a n c e 1 (: ,1) , T r a n s m i t a n c e 1 (: ,2) ,’r ’) ;

% plot ( T r a n s m i t a n c e 2 (: ,1) , T r a n s m i t a n c e 2 (: ,2) ,’g ’) ;

% plot ( T r a n s m i t a n c e 3 (: ,1) , T r a n s m i t a n c e 3 (: ,2) ,’b ’) ; plot ( T r a n s m i t a n c e 4 (: ,1) , T r a n s m i t a n c e 4 (: ,2) , ’b ’) ;

% plot ( T r a n s m i t a n c e 5 (: ,1) , T r a n s m i t a n c e 5 (: ,2) ,’c ’) ;

% plot ( T r a n s m i t a n c e 6 (: ,1) , T r a n s m i t a n c e 6 (: ,2) ,’m ’) ;

% % hold on ;

p e r m i t t i v i t y 3 = load ( ’ p e r m i t t i v i t y 3 . nk ’) ;

plot ( p e r m i t t i v i t y 3 (: ,1) , p e r m i t t i v i t y 3 (: , 3) , ’ - ’) ; plot ( wp *1 e9 , e2a , ’* ’) ;

ylabel ( ’ C om p l e x Part P e r m i t t i v i t y [] ’) ; xlabel ( ’ Wave Length [ nm ] ’) ;

hold off ;

% % Plot the s t u f f !

% hold on ;

% plot ( f (1: end ) , log ( mu (1: end ) ) , ’* ’) ;

% y l a b e l ( ’ l o c a l d e n s i t y p r o b a b i l i t y [] ’) ;

% x l a b e l ( ’ v o l u m e f r a c t i o n [] ’) ;

% hold off ;

% { hold on ; s u b p l o t (2 ,2 ,1) ;

25

imshow ( imgbw ) ;

plot ( f (1: end ) , mu (1: end ) , ’* ’) ;

ylabel ( ’ local d e n s i t y p r o b a b i l i t y [] ’) ; xlabel ( ’ volume f r a c t i o n [] ’) ;

% }

plot (f , la , ’* ’) ;

ylabel ( ’ local p e r c o l a t i o n [] ’) ; xlabel ( ’ volume f r a c t i o n [] ’) ;

% {

plot ( sx , sy , ’* ’) ;

ylabel ( ’ i n f o r m a t i o n e n t r o p y [] ’) ; xlabel ( ’ number of pixels [] ’) ; hold off ;

% }

f u n c t i o n [ ...

local_den sity , ...

l o c a l _ p e r c u l a t i o n , ...

e n t r o p y ] = c a l c u l a t e _ l o c a l _ d a t a ( image , b o x W i d t h )

% The f u n c t i o n c a l c u l a t e s the l o c a l i m a g e data of a g i v e n s p e c t r o s c o p y

% i m a g e and r e t u r n s the v a l u e .

% P A R A M E T E R S :

% i m a g e - A s p e c t r o s c o p y i m a g e c o n v e r t e d into b l a c k w h i t e c o l o r .

% b o x W i d t h - The size of e v e r y box in the grid .

% R E T U R N :

% l o c a l _ d e n s i t y - The c a l c u l a t e d l o c a l d e n s i t y p r o b a b i l i t y f u n c t i o n .

% l o c a l _ p e r c u l a t i o n - The c a l c u l a t e d l o c a l p e r c u l a t i o n p r o b a b i l i t y

% f u n c t i o n .

% e n t r o p y - The c a l c u l a t e d i n f o r m a t i o n e n t r o p y . l o c a l _ d e n s i t y = zeros (1 , b o x W i d t h ^2 + 1) ;

l o c a l _ p e r c u l a t i o n = zeros (1 , b o x W i d t h ^2 + 1) ; e n t r o p y = 0;

% % P a r t i t i o n i m a g e into a grid [h , w ] = size ( image ) ;

grid = c r e a t e _ g r i d (1 , 1 , w , h , b o x W i d t h ) ;

% % C a l c u l a t e the l o c a l d e n s i t y p r o b a b i l i t y and p e r c u l a t i o n f u n c t i o n for i =1: length ( grid )

wCount = 0;

exists = 0;

croped = imcrop ( image , [ grid (i ,1) , grid (i ,2) , boxWidth -1 , boxWidth -1]) ; [ th , tw ] = size ( croped ) ;

for j =1: th for k =1: tw

if croped (j , k ) == 1 wCount = wCount + 1;

end end

[ tmp ] = f i n d _ p a t h ( croped , [1 j ] , [ tw j ]) ; exists = exists + tmp ;

[ tmp ] = f i n d _ p a t h ( croped , [ j 1] , [ j th ]) ; exists = exists + tmp ;

end

l o c a l _ p e r c u l a t i o n ( wCount +1) = l o c a l _ p e r c u l a t i o n ( wCount +1) + ...

exists / ( th + tw ) ;

l o c a l _ d e n s i t y ( wCount +1) = l o c a l _ d e n s i t y ( wCount +1) + 1;

end

for i =1: length ( l o c a l _ p e r c u l a t i o n )

l o c a l _ p e r c u l a t i o n ( i ) = l o c a l _ p e r c u l a t i o n ( i ) / l o c a l _ d e n s i t y ( i ) ; end

l o c a l _ d e n s i t y = l o c a l _ d e n s i t y / length ( grid ) ;

% % C a l c u l a t e the i n f o r m a t i o n e n t r o p y for i =1: length ( l o c a l _ d e n s i t y )

if l o c a l _ d e n s i t y ( i ) ~= 0

e n t r o p y = e n t r o py - l o c a l _ d e n s i t y ( i ) .* log ( l o c a l _ d e n s i t y ( i ) ) ; end

end

e n t r o p y = e n t r op y / log ( b o x W i d t h ^2 + 1) ; end

f u n c t i o n [ ...

local_dens ity , ...

e n t r o p y ] = c a l c u l a t e _ l o c a l _ d a t a _ 1 ( image , b o x W i d t h )

% The f u n c t i o n c a l c u l a t e s the l o c a l i m a g e data of a g i v e n s p e c t r o s c o p y

% i m a g e and r e t u r n s the v a l u e .

% P A R A M E T E R S :

% i m a g e - A s p e c t r o s c o p y i m a g e c o n v e r t e d into b l a c k w h i t e c o l o r .

% b o x W i d t h - The size of e v e r y box in the grid .

% R E T U R N :

% l o c a l _ d e n s i t y - The c a l c u l a t e d l o c a l d e n s i t y p r o b a b i l i t y f u n c t i o n .

% e n t r o p y - The c a l c u l a t e d i n f o r m a t i o n e n t r o p y . l o c a l _ d e n s i t y = zeros (1 , b o x W i d t h ^2 + 1) ;

e n t r o p y = 0;

f = l i n s p a c e (0 , 1 , b o x W i d t h ^2+1) ;

% % P a r t i t i o n i m a g e into a grid [h , w ] = size ( image ) ;

grid = c r e a t e _ g r i d (1 , 1 , w , h , b o x W i d t h ) ;

% % C a l c u l a t e the l o c a l d e n s i t y p r o b a b i l i t y for i =1: length ( grid )

wCount = 0;

croped = imcrop ( image , [ grid (i ,1) , grid (i ,2) , boxWidth -1 , boxWidth -1]) ; [ th , tw ] = size ( croped ) ;

for j =1: th for k =1: tw

if croped (j , k ) == 1 wCount = wCount + 1;

end end end

l o c a l _ d e n s i t y ( wCount +1) = l o c a l _ d e n s i t y ( wCount +1) + 1;

end

l o c a l _ d e n s i t y = l o c a l _ d e n s i t y / length ( grid ) ;

% % C a l c u l a t e the i n f o r m a t i o n e n t r o p y for i =1: length ( l o c a l _ d e n s i t y )

if l o c a l _ d e n s i t y ( i ) ~= 0

e n t r o p y = e n t r o py - l o c a l _ d e n s i t y ( i ) .* log ( l o c a l _ d e n s i t y ( i ) ) ; end

end

e n t r o p y = e n t r op y / log ( b o x W i d t h ^2 + 1) ; end

f u n c t i o n [ value ] = d i v i d e d _ d i f f e r e n c e 2 ( y0 , y1 , x0 , x1 )

% R e t u r n s the d i v i d e d d i f f e r e n c e of the g i v e n p o i n t s .

% P A R A M E T E R S

27

% y0 - The f i r s t e v a l u a t e d f u n c t i o n p o i n t .

% y1 - The s e c o n d e v a l u a t e d f u n c t i o n p o i n t s .

% x0 - F i r s t e v a l u a t i o n p o i n t s .

% x1 - S e c o n d e v a l u a t i o n p o i n t s .

% R E T U R N S

% v a l u e - The d i v i d e d d i f f e r e n c e . value = ( y1 - y0 ) ./ ( x1 - x0 ) ; end

f u n c t i o n [ value ] = d i v i d e d _ d i f f e r e n c e 3 ( y0 , y1 , y2 , x0 , x1 , x2 )

% R e t u r n s the d i v i d e d d i f f e r e n c e of the g i v e n p o i n t s .

% P A R A M E T E R S

% y0 - The f i r s t e v a l u a t e d f u n c t i o n p o i n t s .

% y1 - The s e c o n d e v a l u a t e d f u n c t i o n p o i n t s .

% y2 - The t h i r d e v a l u a t e d f u n c t i o n p o i n t s .

% x0 - F i r s t e v a l u a t i o n p o i n t s .

% x1 - S e c o n d e v a l u a t i o n p o i n t s .

% x2 - T h i r d e v a l u a t i o n p o i n t s .

% R E T U R N S

% v a l u e - The d i v i d e d d i f f e r e n c e .

value = ( y2 - y1 ) ./ (( x2 - x1 ) .* ( x2 - x0 ) ) ...

- ( y1 - y0 ) ./ (( x1 - x0 ) .* ( x2 - x0 ) ) ; end

f u n c t i o n [ grid ] = c r e a t e _ g r i d ( xmin , ymin , width , height , b o x W i d t h )

% C r e a t e a grid of the g i v e n i n p u t data and r e t u r n s the r e c t a n g l e

% c o o r d i n a t e s .

% P A R A M E T E R S :

% xmin - The m i n i m u m x v a l u e that the grid c o n t a i n s .

% ymin - The m i n i m u m y v a l u e that the grid c o n t a i n s .

% w i d t h - The grid w i d t h in p i x e l s .

% h e i g h t - The grid h e i g h t in p i x e l s .

% b o x W i d t h - The w i d t h of the s q u a r e s in the grid .

% R E T U R N :

% grid - The s p e c i f i e d grid . grid = [0 0 0 0];

index = 1;

xi = xmin ; yi = ymin ;

while yi + b o x W i d t h < height while xi + boxWidth -1 < width

grid ( index , 1:4) = [ xi , yi , xi + boxWidth -1 , yi + boxWidth -1];

index = index + 1;

xi = xi + b o x W i d t h ; end

xi = xmin ;

yi = yi + b o x W i d t h ; end

end

f u n c t i o n [ exists ] = f i n d _ p a t h ( graph , start , stop )

% L o o k s for a path in the g i v e n g r a p h b e t w e e n two g i v e n p o i n t s .

% P A R A M E T E R S :

% g r a p h - A m a t r i x with n o d e s r e p r e s e n t e d by 1 ’ s and e v e r y node is

% c o n n e c t d to e v e r y n e i g h b o u r i n g node .

% s t a r t - The i n d e x of the s t a r t node [ x y ].

% stop - The i n d e x of the stop node [ x y ].

% R E T U R N :

% e x i s t s - True if t h e r e is a path b e t w e e n the two n o d e s .

% %

data = zeros ( size ( graph ) ) ; exists = 0;

stack = start ; [h , w ] = size ( graph ) ;

data ( start (1 , 2) , start (1 , 1) ) = 1;

if data ( start (1 , 2) , start (1 , 1) ) == 0 return ;

end

% % S t a r t depth - f i r s t s e a r c h with c o n n e c t i v i t y 4 while i s e m p t y ( stack ) == 0 ...

&& ( stack ( end , 1) ~= stop (1) || stack ( end , 2) ~= stop (2) ) x = stack ( end , 1) ;

y = stack ( end , 2) ;

if x < w && graph (y , x +1) == 1 && data (y , x +1) < 1 stack (1+ end , 1:2) = [ x +1 , y ];

data (y , x +1) = 1;

elseif y < h && graph ( y +1 , x ) == 1 && data ( y +1 , x ) < 1 stack (1+ end , 1:2) = [x , y +1];

data ( y +1 , x ) = 1;

elseif x > 1 && graph (y , x -1) == 1 && data (y , x -1) < 1 stack (1+ end , 1:2) = [x -1 , y ];

data (y , x -1) = 1;

elseif y > 1 && graph (y -1 , x ) == 1 && data (y -1 , x ) < 1 stack (1+ end , 1:2) = [x , y -1];

data (y -1 , x ) = 1;

else

stack ( end , :) = [];

end end

if i s e m p t y ( stack ) == 0 exists = 1;

end end

f u n c t i o n [ ret ] = f z e r o _ f u n c ( value , ff , connectivity , density , list )

% The s o l u t i o n to the e q u a t i o n f z e r o _ f u n c ( v a l u e ) = 0 is the r e q u i r e d

% e f f e c t i v e p e r m i t t i v i t y of the s a m p l e .

% P A R A M E T E R S

% v a l u e - The v a l u e of the e f f e c t i v e p e r m i t t i v i t y .

% ff - The f i l l i n g f a c t o r c o r r e s p o n d i n g to the data p o i n t s .

% c o n n e c t i v i t y - The c o n n e c t i v i t y f u n c t i o n of the s a m p l e .

% d e n s i t y - The d e n s i t y d i s t r o b u t i o n of the s a m p l e .

% list - List of v a l u s to use

% R E T U R N S

% ret - The e v a l u a t e d f u n c t i o n at all the s p e c i f i e d w a v e l e n g t h s . ret = zeros (1 , length ( value ) ) ;

for i =1: length ( ff )

[~ , ec , enc ] = p e r m i t t i v i t y _ a l l ( ff ( i ) ) ; ret = ret + ...

( c o n n e c t i v i t y ( i ) * ( ec ( list ) - value ) ./ ( ec ( list ) + 2 * value ) + ...

(1 - c o n n e c t i v i t y ( i ) ) * ( enc ( list ) - value ) ./ ( enc ( list ) + 2 * value ) ) * ...

d e n s i t y ( i ) ; end

end

f u n c t i o n [ value ] = m u l l e r _ r e c u r r e n c e ( y0 , y1 , y2 , x0 , x1 , x2 )

% C a l c u l a t e s the next step in the m u l l e r r e c u r r e n c e r e l a t i o n .

% P A R A M E T E R S

29

% y2 - v a l u e s at p o i n t two .

% y1 - v a l u e s at p o i n t one .

% y0 - v a l u e s at p o i n t zero .

% x2 - p o i n t two .

% x1 - p o i n t one .

% x0 - p o i n t zero .

% R E T U R N S

% v a l u e - The next i t e r a t i o n of the m u l l e r r e c u r r e n c e . value = zeros (1 , length ( y0 ) ) ;

w = d i v i d e d _ d i f f e r e n c e 2 ( y2 , y1 , x2 , x1 ) ...

+ d i v i d e d _ d i f f e r e n c e 2 ( y2 , y0 , x2 , x0 ) ...

- d i v i d e d _ d i f f e r e n c e 2 ( y1 , y0 , x1 , x0 ) ;

d1 = sqrt ( w .* w - 4 * y2 .* d i v i d e d _ d i f f e r e n c e 3 ( y2 , y1 , y0 , x2 , x1 , x0 ) ) ; d2 = w - d1 ;

d1 = w + d1 ; for i =1: length ( y0 )

if abs ( d1 ( i ) ) > abs ( d2 ( i ) )

value ( i ) = x2 ( i ) - 2 * y2 ( i ) ./ d1 ( i ) ; else

value ( i ) = x2 ( i ) - 2 * y2 ( i ) ./ d2 ( i ) ; end

end end

f u n c t i o n [ wavelength , p e r m i t t i v i t y ] = p e r m i t t i v i t y _ a i r ()

% R e t u r n s the p e r m i t t i v i t y of air .

% R E T U R N

% w a v e l e n g t h - The w a v e l e n g t h s at w h i c h the p e r m i t i v i t y is g i v e n .

% p e r m i t t i v i t y - The c o m p l e x p e r m i t t i v i t y of air at the g i v e n r a n g e of

% w a v e l e n g t h s .

[ wavelength , ~] = p e r m i t t i v i t y _ a u () ; p e r m i t t i v i t y = ones (1 , length ( w a v e l e n g t h ) ) ; end

f u n c t i o n [ wavelength , p e r m i t t iv i t y _ c , p e r m i t t i v i t y _ n c ] = p e r m i t t i v i t y _ a l l ( f )

% This f u n c t i o n c a l c u l a t e s the p e r m i t t i v i t y of the c o n d u c t i n g and the

% n o n c o n d u c t i n g c e l l s .

% P A R A M E T E R S

% f - The f i l l i n g f a c t o r .

% R E T U R N S

% w a v e l e n g t h - A v e c t o r with all the w a v e l e n g t h s at w h i c h the

% p e r m i t t i v i t i e s are s p e c i f i e d .

% p e r m i t t i v i t y _ c - The p e r m i t t i v i t y of the c o n d u c t i n g c e l l s .

% p e r m i t t i v i t y _ n c - The p e r m i t t i v i t y of the n o n c o n d u c t i n g c e l l s . [ wavelength , eg ] = p e r m i t t i v i t y _ a u () ;

[~ , ea ] = p e r m i t t i v i t y _ a i r () ;

% p e r m i t t i v i t y _ c = eg .* (1 - (1 - f ) ./ (1 ./ (1 - ea ./ eg ) - 1 / 3 * f ) ) ;

% p e r m i t t i v i t y _ n c = ea .* (1 - f ./ (1 ./ (1 - eg ./ ea ) - 1 / 3 * (1 - f ) ) ) ; p e r m i t t i v i t y _ c = (2/3 * eg .* f .* ( eg - ea ) + eg .* ea ) ./ ...

( eg - 1/3 .* f .* ( eg - ea ) ) ;

p e r m i t t i v i t y _ n c = ( ea .^ 2 - 2/3 * ea .* f .* ( ea - eg ) - 1/3 * ( ea - eg ) .*

ea ) ./ ...

( ea - 1/3 * (1 - f ) .* ( ea - eg ) ) ; end

f u n c t i o n [ wavelength , p e r m i t t i v i t y ] = p e r m i t t i v i t y _ a u ()

% R e t u r n s the p e r m i t t i v i t y of gold .

% R E T U R N

% w a v e l e n g t h - The w a v e l e n g t h s at w h i c h the p e r m i t i v i t y is g i v e n .

% p e r m i t t i v i t y - The c o m p l e x p e r m i t t i v i t y of gold at the g i v e n r a n g e of

% w a v e l e n g t h s .

w a v e l e n g t h = 1e -9 * [ 2 5 0 ; 2 5 2 . 0 3 ; 2 5 4 . 1 ; ...

2 5 6 . 2 ; 2 5 8 . 3 3 ; 2 6 0 . 5 ; ...

2 6 2 . 7 1 ; 2 6 4 . 9 6 ; 2 6 7 . 2 4 ; ...

2 6 9 . 5 7 ; 2 7 1 . 9 3 ; 2 7 4 . 3 4 ; ...

2 7 6 . 7 9 ; 2 7 9 . 2 8 ; 2 8 1 . 8 2 ; ...

2 8 4 . 4 ; 2 8 7 . 0 4 ; 2 8 9 . 7 2 ; ...

2 9 2 . 4 5 ; 2 9 5 . 2 4 ; 2 9 8 . 0 8 ; ...

3 0 0 . 9 7 ; 3 0 3 . 9 2 ; 3 0 6 . 9 3 ; ...

3 1 0 ; 3 1 3 . 1 3 ; 3 1 6 . 3 3 ; ...

3 1 9 . 5 9 ; 3 2 2 . 9 2 ; 3 2 6 . 3 2 ; ...

3 2 9 . 7 9 ; 3 3 3 . 3 3 ; 3 3 6 . 9 6 ; ...

3 4 0 . 6 6 ; 3 4 4 . 4 4 ; 3 4 8 . 3 1 ; ...

3 5 2 . 2 7 ; 3 5 6 . 3 2 ; 3 6 0 . 4 7 ; ...

3 6 4 . 7 1 ; 3 6 9 . 0 5 ; 3 7 3 . 4 9 ; ...

3 7 8 . 0 5 ; 3 8 2 . 7 2 ; 3 8 7 . 5 ; ...

3 9 2 . 4 1 ; 3 9 7 . 4 4 ; 4 0 2 . 6 ; ...

4 0 7 . 8 9 ; 4 1 3 . 3 3 ; 4 1 8 . 9 2 ; ...

4 2 4 . 6 6 ; 4 3 0 . 5 6 ; 4 3 6 . 6 2 ; ...

4 4 2 . 8 6 ; 4 4 9 . 2 8 ; 4 5 5 . 8 8 ; ...

4 6 2 . 6 9 ; 4 6 9 . 7 ; 4 7 6 . 9 2 ; ...

4 8 4 . 3 8 ; 4 9 2 . 0 6 ; 5 0 0 ; 5 0 8 . 2 ; ...

5 1 6 . 6 7 ; 5 2 5 . 4 2 ; 5 3 4 . 4 8 ; ...

5 4 3 . 8 6 ; 5 5 3 . 5 7 ; 5 6 3 . 6 4 ; ...

5 7 4 . 0 7 ; 5 8 4 . 9 1 ; 5 9 6 . 1 5 ; ...

6 0 7 . 8 4 ; 6 2 0 ; 6 3 2 . 6 5 ; ...

6 4 5 . 8 3 ; 6 5 9 . 5 7 ; 6 7 3 . 9 1 ; ...

6 8 8 . 8 9 ; 7 0 4 . 5 5 ; 7 2 0 . 9 3 ; ...

7 3 8 . 1 ; 7 5 6 . 1 ; 7 7 5 ; ...

7 9 4 . 8 7 ; 8 1 5 . 7 9 ; 8 3 7 . 8 4 ; ...

8 6 1 . 1 1 ; 8 8 5 . 7 1 ; 9 1 1 . 7 6 ; ...

9 3 9 . 3 9 ; 9 6 8 . 7 5 ; 1 0 0 0 ; ...

1 0 3 3 . 3 ; 1 0 6 9 ; 1 1 0 7 . 1 ; ...

1 1 4 8 . 1 ; 1 1 9 2 . 3 ; 1 2 4 0 ; ...

1 2 6 5 . 3 ; 1 3 1 9 . 1 ; 1 4 7 6 . 2 ; ...

1 5 5 0 ; 1 6 3 1 . 6 ; 1 7 2 2 . 2 ; ...

1 8 2 3 . 5 ; 1 9 3 7 . 5 ; ] ’ ; e1 = [ - 2 . 0 4 9 0 4 3 9 3 ; - 2 . 0 2 6 3 7 3 9 9 ; ...

- 1 . 9 9 9 6 3 8 ; - 1 . 9 6 8 9 7 6 0 5 ; ...

- 1 . 9 3 4 5 5 6 4 ; - 1 . 8 9 6 5 7 7 2 5 ; ...

- 1 . 8 5 5 2 6 7 2 2 ; - 1 . 8 1 0 8 8 5 6 1 ; ...

- 1 . 7 6 3 7 2 4 0 9 ; - 1 . 7 1 4 1 0 5 9 7 ; ...

- 1 . 6 6 2 3 8 7 4 5 ; - 1 . 6 0 8 9 6 0 6 6 ; ...

- 1 . 5 5 4 2 5 0 2 8 ; - 1 . 4 9 8 7 1 8 8 2 ; ...

- 1 . 4 4 2 8 6 3 7 1 ; - 1 . 3 8 7 2 2 1 8 9 ; ...

- 1 . 3 3 2 3 6 6 0 1 ; - 1 . 2 7 8 9 0 8 1 ; ...

- 1 . 2 2 7 4 9 6 5 2 ; - 1 . 1 7 8 8 1 4 2 4 ; ...

- 1 . 1 3 3 5 7 3 9 6 ; - 1 . 0 9 2 5 1 3 6 3 ; ...

- 1 . 0 5 6 3 8 4 3 ; - 1 . 0 2 5 9 3 7 6 2 ; ...

- 1 . 0 0 1 9 0 8 3 1 ; - 0 . 9 8 4 9 9 0 3 0 7 ; ...

-0.975808 975; ...

- 0 . 9 7 4 8 8 7 5 8 1 ; - 0 . 9 8 2 6 1 0 0 5 5 ; ...

- 0 . 9 8 9 7 8 4 5 3 1 ; - 1 . 0 1 0 7 8 7 4 2 ; ...

- 1 . 0 4 0 5 1 6 0 4 ; - 1 . 0 7 8 5 4 8 3 1 ; ...

- 1 . 1 2 4 1 5 3 6 4 ; - 1 . 1 7 6 2 7 5 1 8 ; ...

- 1 . 2 3 3 5 2 7 5 6 ; - 1 . 2 9 4 2 1 4 9 4 ; ...

- 1 . 3 5 6 3 7 3 8 7 ; - 1 . 4 1 7 8 3 9 6 8 ; ...

- 1 . 4 7 6 3 3 9 7 9 ; - 1 . 5 2 9 6 0 9 2 8 ; ...

- 1 . 5 7 5 5 2 6 9 3 ; - 1 . 6 1 2 2 6 3 8 7 ; ...

- 1 . 6 3 8 4 3 5 0 9 ; - 1 . 6 5 3 2 4 3 6 5 ; ...

- 1 . 6 5 6 6 0 7 8 7 ; - 1 . 6 4 9 2 5 8 5 9 ; ...

- 1 . 6 3 2 7 9 7 6 6 ; - 1 . 6 0 9 7 0 9 3 8 ; ...

- 1 . 5 8 3 3 2 1 7 9 ; - 1 . 5 5 7 7 1 6 8 1 ; ...

31

- 1 . 5 3 7 5 9 1 8 2 ; - 1 . 5 2 8 0 8 3 9 4 ; ...

- 1 . 5 3 4 5 6 4 4 8 ; - 1 . 5 6 2 4 1 9 4 ; ...

- 1 . 6 1 6 8 3 1 8 4 ; - 1 . 7 0 2 5 7 8 9 8 ; ...

- 1 . 8 2 3 8 6 1 1 2 ; - 1 . 9 8 4 1 6 7 4 6 ; ...

- 2 . 1 8 6 1 9 5 5 4 ; - 2 . 4 3 1 8 2 3 9 3 ; ...

- 2 . 7 2 2 1 2 1 9 4 ; - 3 . 0 5 7 4 3 1 9 ; ...

- 3 . 4 3 7 4 6 4 1 1 ; - 3 . 8 6 1 4 4 0 6 8 ; ...

- 4 . 3 2 8 2 5 1 ; - 4 . 8 3 6 6 1 0 8 3 ; ...

- 5 . 3 8 5 2 2 0 0 5 ; - 5 . 9 7 2 9 2 6 6 ; ...

- 6 . 5 9 8 8 3 7 3 7 ; - 7 . 2 6 2 4 4 1 4 9 ; ...

- 7 . 9 6 3 6 9 2 2 5 ; - 8 . 7 0 3 0 6 5 0 4 ; ...

- 9 . 4 8 1 6 1 6 9 3 ; - 1 0 . 3 0 0 9 9 3 1 ; ...

- 1 1 . 1 6 3 4 6 3 2 ; - 1 2 . 0 7 1 9 0 9 1 ; ...

- 1 3 . 0 2 9 8 6 0 9 ; - 1 4 . 0 4 1 4 7 5 6 ; ...

- 1 5 . 1 1 1 5 8 9 9 ; - 1 6 . 2 4 5 7 1 5 1 ; ...

- 1 7 . 4 5 0 1 2 1 5 ; - 1 8 . 7 3 1 8 6 1 3 ; ...

- 2 0 . 0 9 8 8 8 6 4 ; - 2 1 . 5 6 0 1 3 8 3 ; ...

- 2 3 . 1 2 5 6 8 5 9 ; - 2 4 . 8 0 6 8 9 5 7 ; ...

- 2 6 . 6 1 6 6 2 2 9 ; - 2 8 . 5 6 9 4 6 5 7 ; ...

- 3 0 . 6 8 2 0 8 1 5 ; - 3 2 . 9 7 3 4 9 1 ; ...

- 3 5 . 4 6 5 5 8 5 8 ; - 3 8 . 1 8 3 5 8 9 2 ; ...

- 4 1 . 1 5 6 7 4 0 2 ; - 4 4 . 4 1 9 1 1 3 8 ; ...

- 4 8 . 0 1 0 5 4 2 4 ; - 5 1 . 9 7 7 9 0 9 5 ; ...

- 5 6 . 3 7 6 6 4 4 9 ; - 6 1 . 2 7 2 7 0 9 3 ; ...

- 6 6 . 7 4 5 0 6 0 3 ; - 7 2 . 8 8 8 8 8 0 1 ; ...

- 7 9 . 8 1 9 6 5 5 6 ; - 1 0 1 . 6 0 0 5 0 7 ; ...

- 1 1 2 . 6 3 9 7 7 1 ; - 1 2 5 . 4 3 4 1 7 6 ; ...

- 1 4 0 . 3 7 6 2 9 7 ; - 1 5 7 . 9 7 4 1 0 2 ; ...

-178.894047;] ’;

e2 = [ 4 . 0 2 6 7 1 1 6 1 ; 4 . 1 3 9 2 4 4 5 9 ; ...

4 . 2 4 9 5 4 8 1 ; 4 . 3 5 7 2 7 9 2 6 ; ...

4 . 4 6 2 0 9 9 0 7 ; 4 . 5 6 3 6 7 0 4 3 ; ...

4 . 6 6 1 6 6 1 4 8 ; 4 . 7 5 5 7 4 8 6 9 ; ...

4 . 8 4 5 6 1 6 9 3 ; 4 . 9 3 0 9 6 4 8 3 ; ...

5 . 0 1 1 5 0 5 8 9 ; 5 . 0 8 6 9 6 8 2 4 ; ...

5 . 1 5 7 1 0 4 1 6 ; 5 . 2 2 1 6 8 7 4 9 ; ...

5 . 2 8 0 5 2 2 7 7 ; 5 . 3 3 3 4 4 5 8 7 ; ...

5 . 3 8 0 3 3 4 3 6 ; 5 . 4 2 1 1 1 0 8 ; ...

5 . 4 5 5 7 5 3 0 6 ; 5 . 4 8 4 3 0 3 2 2 ; ...

5 . 5 0 6 8 7 8 8 3 ; 5 . 5 2 3 6 8 3 6 6 ; ...

5 . 5 3 5 0 2 1 5 6 ; 5 . 5 4 1 3 0 8 5 6 ; ...

5 . 5 4 3 0 8 4 7 1 ; 5 . 5 4 1 0 2 4 2 1 ; ...

5 . 5 3 5 9 4 1 8 ; 5 . 5 2 8 7 9 3 7 6 ; ...

5 . 5 2 0 6 7 1 6 4 ; 5 . 5 1 2 7 8 7 0 9 ; ...

5 . 5 0 6 4 4 4 9 8 ; 5 . 5 0 3 0 0 4 3 4 ; ...

5 . 5 0 3 8 2 6 ; 5 . 5 1 0 2 0 7 1 1 ; ...

5 . 5 2 3 3 0 3 8 ; 5 . 5 4 4 0 4 5 8 1 ; ...

5 . 5 7 3 0 4 5 7 ; 5 . 6 1 0 5 1 0 0 2 ; ...

5 . 6 5 6 1 5 8 2 4 ; 5 . 7 0 9 1 5 8 6 6 ; ...

5 . 7 6 8 0 8 4 9 9 ; 5 . 8 3 0 9 0 6 8 7 ; ...

5 . 8 9 5 0 1 2 8 ; 5 . 9 5 7 2 7 4 9 2 ; ...

6 . 0 1 4 1 5 3 1 3 ; 6 . 0 6 1 8 3 4 5 1 ; ...

6 . 0 9 6 4 0 6 3 6 ; 6 . 1 1 4 0 4 8 0 5 ; ...

6 . 1 1 1 2 3 3 7 3 ; 6 . 0 8 4 9 3 0 4 8 ; ...

6 . 0 3 2 7 7 8 7 4 ; 5 . 9 5 3 2 3 7 7 4 ; ...

5 . 8 4 5 6 9 2 4 4 ; 5 . 7 1 0 5 0 6 5 8 ; ...

5 . 5 4 9 0 1 8 1 7 ; 5 . 3 6 3 4 8 0 9 6 ; ...

5 . 1 5 6 9 5 5 8 6 ; 4 . 9 3 3 1 6 3 7 5 ; ...

4 . 6 9 6 2 9 3 ; 4 . 4 5 0 8 0 6 7 5 ; ...

4 . 2 0 1 2 3 9 8 2 ; 3 . 9 5 2 0 0 1 ; ...

3 . 7 0 7 2 0 6 3 3 ; 3 . 4 7 0 5 4 1 9 8 ; ...

3 . 2 4 5 1 6 0 1 7 ; 3 . 0 3 3 6 1 6 3 6 ; ...

2 . 8 3 7 8 5 3 2 3 ; 2 . 6 5 9 2 0 9 4 7 ; ...

2 . 4 9 8 4 6 7 3 ; 2 . 3 5 5 9 1 7 4 6 ; ...

2 . 2 3 1 4 4 4 9 5 ; 2 . 1 2 4 6 1 7 9 3 ; ...

2 . 0 3 4 7 8 2 0 2 ; 1 . 9 6 1 1 4 8 0 2 ; ...

1 . 9 0 2 8 7 1 2 1 ; 1 . 8 5 9 1 2 2 3 4 ; ...

1 . 8 2 9 1 4 3 3 4 ; 1 . 8 1 2 2 9 2 7 6 ; ...

1 . 8 0 8 0 7 9 7 1 ; 1 . 8 1 6 1 8 8 5 5 ; ...

1 . 8 3 6 4 9 6 2 9 ; 1 . 8 6 9 0 8 5 0 3 ; ...

1 . 9 1 4 2 5 2 5 6 ; 1 . 9 7 2 5 2 2 7 1 ; ...

2 . 0 4 4 6 5 7 5 3 ; 2 . 1 3 1 6 7 4 8 1 ; ...

2 . 2 3 4 8 7 1 6 6 ; 2 . 3 5 5 8 5 6 1 4 ; ...

2 . 4 9 6 5 9 1 6 2 ; 2 . 6 5 9 4 5 0 6 5 ; ...

2 . 8 4 7 2 9 3 5 6 ; 3 . 0 6 3 5 5 7 4 4 ; ...

3 . 3 1 2 3 8 3 6 6 ; 3 . 5 9 8 7 7 0 3 9 ; ...

3 . 9 2 8 7 7 3 9 9 ; 4 . 3 0 9 7 6 9 4 4 ; ...

4 . 7 5 0 7 8 2 0 8 ; 5 . 2 6 2 9 2 9 6 6 ; ...

5 . 8 5 9 9 9 9 5 6 ; 6 . 5 5 9 2 1 4 4 ; ...

7 . 3 8 2 2 6 9 5 1 ; 8 . 3 5 6 7 2 9 3 7 ; ...

1 1 . 7 1 4 1 1 1 3 ; 1 3 . 5 7 5 6 0 8 1 ; ...

1 5 . 8 5 9 4 9 0 7 ; 1 8 . 6 9 1 3 4 4 3 ; ...

2 2 . 2 4 4 3 4 0 5 ; 2 6 . 7 6 1 5 9 3 2 ] ’ ; p e r m i t t i v i t y = e1 + 1 i * e2 ; end

33