Bala's Website
Research and day-to-day Ramblings
Electric field at localized plasmon resonance using MEEP
23 months ago
This article is about simulating localized plasmon resonances in metal nanospheres using MEEP package. Generally, I am interested in solving three problems in LSPR systems:
- Calculate the extinction, scattering, absorption spectra of metal nanoparticle
- Calculating the electric field enhancement spatially as function of wavelength
This involves taking electric field distributions with a particle in time domain and taking FFT of them. Also to be noted is that the electric fields near the particle should be normalized with electric fields with no nanoparticle. This has to be done by an external program like octave/matlab. One has to do this for three components of electric field in two planes. This was the tough part to pull off. - Calculate electric field at plasmon resonance as function of time
This was very easy. Get the electric fields at certain wavelength and then multiply it with a time varying sin function.
The procedure for doing this is very similar to the method I mentioned here.
Below are some results for a silver sphere of 25 nm radius. I am happy with these results. They seem to match the analytical results. Some of the animations are large in size and may take time to load on your browser. If you are interested in code, please contact me.
Update: This method also works well for gold nano-particles of 50 nm radius. See updated figure.
Building my new desktop system
06 months ago
by Bala Krishna Juluri
in Ramblings
As I had access to great workstations at my previous labs, I never needed a personal desktop system to do heavy duty operations such as running simulations, photo/video editing etc. But now that I moved out of graduate school and have some money to spend, I am finding myself more in need for a good desktop system. My plan with a new system is to make more tutorials on plasmonics and do more photo editing (now that I own a Nikon D5100). I currently have three old laptops, but all of them are just good enough for browsing but not any of things I want to do. I have a great system at my current workspace, but I do not want to use it for my personal adventures. 
I could have bought a pre-installed system, but there is no fun in it . In some cases, buying a pre-configured system may be a little bit cheaper, but at an expense of no flexibility and customization. I have never assembled a system before, so I thought it will be a good learning process. I decided to build a personal system couple of months ago and I wanted to get decent stuff at a good price and time was no factor as I am in no hurry.
Justin (my colleague) referred me to a great website called pc partpicker. Here one can make their own builds based on user ratings, prices. One can set alerts for any price drops. Alternatively, we would also check for any deals that would pop-up on slickdeals website.
Currently I have the following parts:
- Power supply: PSU SEASONIC X650 GOLD SS-650KM R – $96 (Current price is $140, got a 40$ discount)
- Motherboard: ASRock Z77 Extreme4 ATX Intel Motherboard- $115 (Current price is 124$+taxes)
- Memory: (8Gb) GEIL GEV38GB1600C9DC R – Free with motherboard (Current price is$40+taxes)
- Case: Antec Nine Hundred Mid Tower ATX Computer Case- ~$50 (Current price $100+taxes, got a discount on fry’s)
- Processor and Storage: Intel I7-3770K + 240 GB intel SSD 520 Series : $355 (Current price is $586+taxes, but got it from a friend who works for Intel at a reduced price )
- Operating System: Ubuntu 12.10: Free
- I paid a total price of ~$630.
For time being, I have decieded not to overclock the processor.
I have learned that the key is to wait for good deals and explore all options (find friends if anyone works for companies such as intel or amd). I will assemble all the parts soon and post a picture of my build. You can find my final configuration at http://pcpartpicker.com/b/yUT
Parallelization in Octave using parcellfun/pararrayfun
01 year ago
My computer has many processors and I would like to run some octave scripts so that all the processors are being used.
One can use octave function called “pararrayfun” for this purpose. This function is part of “general” package on octave-forge.
On my ubuntu 11.10, I used “sudo apt-get install octave-general” to install this package and ran the following script
1; # this is kept to Prevent Octave from thinking that this is a function file:
close all;
clear all;
function y=test(a,b)
y=sin(a)+cos(b)
endfunction
num_process=8
a_test_inputs=[0:3.14/20:3.14];
b_test_inputs=[0:3.14/20:3.14]*2;
tic ();
tt_par= pararrayfun(num_process,@test,(a_test_inputs),(b_test_inputs));
parallel_elapsed_time = toc ()
tic ();
tt_ser= test((a_test_inputs),(b_test_inputs));
serial_elapsed_time = toc ()
plot(a_test_inputs,tt_par,'-o',a_test_inputs,tt_ser,'-s');
legend('Parallel result','Serial result');
disp(sprintf('Elapsed time during serial computation is %e',serial_elapsed_time))
disp(sprintf('Elapsed time during parallel computation is %e',parallel_elapsed_time))
close all;
clear all;
function y=test(a,b)
y=sin(a)+cos(b)
endfunction
num_process=8
a_test_inputs=[0:3.14/20:3.14];
b_test_inputs=[0:3.14/20:3.14]*2;
tic ();
tt_par= pararrayfun(num_process,@test,(a_test_inputs),(b_test_inputs));
parallel_elapsed_time = toc ()
tic ();
tt_ser= test((a_test_inputs),(b_test_inputs));
serial_elapsed_time = toc ()
plot(a_test_inputs,tt_par,'-o',a_test_inputs,tt_ser,'-s');
legend('Parallel result','Serial result');
disp(sprintf('Elapsed time during serial computation is %e',serial_elapsed_time))
disp(sprintf('Elapsed time during parallel computation is %e',parallel_elapsed_time))
Spaces in strings in matlab/octave
01 year ago
by Bala Krishna Juluri
in Matlab
To get spaces in the strings to work in matlab or octave, use
t1={‘test test’}
Result is
t1 =
{
[1,1] = test test
}
t2=strcat({‘test test ‘},{‘blah blah’})
Result is
t2 =
{
[1,1] = test test blah blah
}
you can use this string in your figures by
plot([1:4])
title(t2{})
van der Pauw correction factor
01 year ago
The van der Pauw Method is a technique commonly used to measure the Resistivity and the Hall Coefficient of a sample. A correction factor goes into calculating the resistivity as described in van der Pauw paper. A iterative method is generally used to calculate the correction factor and this correction factor is plotted in Figure 5 of van der Pauw paper
I reproduced the same figure below using fsolve function in octave.
This figure was produced by the octave code shown below. The raw data is here.
%This octave/matlab code calculates the correction factor,f as a function of Rmnop/Rnopm. This correcton factor will be used in calculation of sheet resistance/resistivity of thin films.
% Fore more details see the paper L. J. van der Pauw, A method of measuring the resistivity and hall coefficients of Lamellae of arbitrary shape, Philips Technical Review, Vol 26, 220.
clear all
close all
x = logspace(0,3,100);
for i=1:length(x)
Lp=@(f)cosh(((x(i)-1)./(x(i)+1))*(log(2)/f))-0.5*exp((log(2)/f));
F(i) = fsolve(Lp,0.5);
end
semilogx(x,F,'-o')
xlabel('R_{RM,OP}/R_{NO,PM}')
ylabel('Correction Factor (f)')
grid('minor');
title('van der Pauw Correction factors for resistivity measurements');
print('Van_der_Paw_correction_factor.png','-dpng','-r200','-S600,600');
% Fore more details see the paper L. J. van der Pauw, A method of measuring the resistivity and hall coefficients of Lamellae of arbitrary shape, Philips Technical Review, Vol 26, 220.
clear all
close all
x = logspace(0,3,100);
for i=1:length(x)
Lp=@(f)cosh(((x(i)-1)./(x(i)+1))*(log(2)/f))-0.5*exp((log(2)/f));
F(i) = fsolve(Lp,0.5);
end
semilogx(x,F,'-o')
xlabel('R_{RM,OP}/R_{NO,PM}')
ylabel('Correction Factor (f)')
grid('minor');
title('van der Pauw Correction factors for resistivity measurements');
print('Van_der_Paw_correction_factor.png','-dpng','-r200','-S600,600');
Surface plasmon dispersion relation for thin metal films
21 year ago
A thin metal film in dielectric (also known as dielectric-metal-dielectric configuration) can support surface plasmons that are different in nature to the ones observed in thick metal-dielectric interfaces. Unlike, a single mode that is observed in thick metal film, thin metal films exhibit two types of modes for the same wavevector due to excitation and interaction of surface plasmons on both sides of the film. One mode (L+) is at higher energy and other (L-) is at a lower energy. The high energy has anti-symmetric field distribution whereas the low energy one has symmetric field distribution. The dispersion relations of these modes can be obtained by applying appropriate boundary conditions and solving Maxwell equations (page 25 of Raether book on Surface plasmons).
The dispersion relations (energy of different modes as function of complex wavevector, 
) are given by the solutions of :
![\omega^+: L+= \epsilon_1k_{z2}+\epsilon_2k_{z1} \tanh\left [ \frac{-ik_{z1}d}{2}\right ]=0](http://juluribk.com/wp-content/ql-cache/quicklatex.com-6644705a10a239aa0cfd17a04af5a24a_l3.png "Rendered by QuickLaTeX.com")
![\omega^-: L-= \epsilon_1k_{z2}+\epsilon_2k_{z1} \coth\left [ \frac{-ik_{z1}d}{2}\right ]=0](http://juluribk.com/wp-content/ql-cache/quicklatex.com-539adc384d3a2fbfed20118c86c952a3_l3.png "Rendered by QuickLaTeX.com")
,
Where 
and 
represent the dielectric functions of dielectric and metal respectively, 
is the thickness of metal film, 
is frequency of high energy mode, 
is frequency of low energy mode, ![k_{z1,z2}=\sqrt{\epsilon_{1,2}\left[ \frac{\omega}{c}\right ]^2-(k_x^{'}(\omega)+ik_x^{''}(\omega))^2}](http://juluribk.com/wp-content/ql-cache/quicklatex.com-d08477f9efc46edbaa551422c010fe27_l3.png "Rendered by QuickLaTeX.com")
.
What does solving these equations mean? if one considers one mode, say 
mode, for a certain , there exists a particular 
and 
that will make the above left hand side of L+ equation equal to zero. How do we get those special 
? Unfortunately, there is no exact solution to this, to solve them, one has to use numerical techniques such as Nelder-Mead minimization algorithm that does unconstrained minimization. The algorithm works cleverly by searching for 
at a particular 
, such that and 
will be minimum (in this case as close to zero as possible).
I have used “fmins” function in octave to solve these dispersion relations. Fmins function in octave can be obtained by installing octave-optmin package. On Ubuntu use “sudo apt-get install octave-optim” in your terminal. Documentation of fmins function is given here.
The numerical solution of dispersion relation for a sample configuration is shown below and qualitatively matches with a figure 5 in this paper. I have not plotted the dispersion relation with imaginary part of wave-vector to simplify the case.
Numerical solution to thin films surface plasmon dispersion obtained by unconstrained minimization algorithm
Below is the octave code written for octave3.2 and has not been checked in matlab.
%% Octave script to solve the dispersion relations of surface
%% plasmons propagating on a thin film embedded in dielectric media.
%% Written by Bala Krishna Juluri (www.juluribk.com), Dec-6-2011.
%% Requires octave-optim package and octave3.2 or greater.
%% A note on units. The code is written in normalized units.
%% All length units are normalized to lambda_p (frequency corresponding to plasma frequency).
%%By normalizing to lambda_p, the dispersion curves become independent of material.
clear all;
close all;
c=1;% velocity of light in vacuum in normalized units.
wp=1;% plasma frequency of metal in normalized units, this will also make kp=1, as kp=wp/c;
d=0.3;% thickness of metal film in normalized units.
% One normalized length unit is equal to 212.99 nm for silver and 377 nm for gold. This is because,
% 1) for silver, plasma frequency is 8.85e15 1/sec , so the plasma wavelength = 212.99 nm
% 2) for copper, plasma frequency is 5e15 1/sec, so the plasma wavelength = 377 nm
% so if d (thickness of the film) is 0.3 in normalized units, it means that in real units,
%% it is equal to 0.3*212.99nm of silver or 0.3*377nm for gold.
epsd=1;# dielectric of surrounding.
tol=1e-10;# Relative size of simplex for fmins function
%creates a non-unifom density of points of w
a=linspace(0.01,3,100);b=1-exp(-a); w=b/max(b)*1.25; % creates a non-uniform density of frequencies
%%from 0 to ~1.3 at which real and complex wavevectors will be calculated.
%%finds the solutions using fmins
for j=1:length(w)
epsm=1-(wp/w(j))^2; # The code uses Drude model for material and assumes no losses. plasma_frequency (omega_p) =1,
#Define first mode
Lp=@(kx)abs(epsm*(sqrt(epsd*(w(j)/c)^2-(kx(1)+i*kx(2))^2))+epsd*(sqrt(epsm*(w(j)/c)^2-(kx(1)+i*kx(2))^2))*tanh(-0.5*i*(sqrt(epsm*(w(j)/c)^2-(kx(1)+i*kx(2))^2))*d));
#solve for first mode using Nelder-Mead technique
[kx_Lp(j,:)] = fmins(Lp,[1.5; 0.25],[0,tol]); # 1.5 and 0.25 are the initial quesses for kx(1) and kx(2)
#Define other mode
Lm=@(kx)abs(epsm*(sqrt(epsd*(w(j)/c)^2-(kx(1)+i*kx(2))^2))+epsd*(sqrt(epsm*(w(j)/c)^2-(kx(1)+i*kx(2))^2))*coth(-0.5*i*(sqrt(epsm*(w(j)/c)^2-(kx(1)+i*kx(2))^2))*d));
#solve for other mode using Nelder-Mead technique
[kx_Lm(j,:)] = fmins(Lm,[1.5; 0.25],[0,tol]);# 1.5 and 0.25 are the initial quesses for kx(1) and kx(2)
end
%Plotting and printing section
putenv("GNUTERM",'wxt')
plot(abs(kx_Lp(:,1)),w,'ro',abs(kx_Lm(:,1)),w,'rs');
xlabel('Real(kx) /k_p')
ylabel('w/w_p')
legend('L+','L-','location','SouthEast')
title('Analytical solution to metal thin-film surface plasmon dispersion relation')
xlim([0,3]);
print('dmd_analytical_web.png','-dpng','-r200','-S600,600');
%% plasmons propagating on a thin film embedded in dielectric media.
%% Written by Bala Krishna Juluri (www.juluribk.com), Dec-6-2011.
%% Requires octave-optim package and octave3.2 or greater.
%% A note on units. The code is written in normalized units.
%% All length units are normalized to lambda_p (frequency corresponding to plasma frequency).
%%By normalizing to lambda_p, the dispersion curves become independent of material.
clear all;
close all;
c=1;% velocity of light in vacuum in normalized units.
wp=1;% plasma frequency of metal in normalized units, this will also make kp=1, as kp=wp/c;
d=0.3;% thickness of metal film in normalized units.
% One normalized length unit is equal to 212.99 nm for silver and 377 nm for gold. This is because,
% 1) for silver, plasma frequency is 8.85e15 1/sec , so the plasma wavelength = 212.99 nm
% 2) for copper, plasma frequency is 5e15 1/sec, so the plasma wavelength = 377 nm
% so if d (thickness of the film) is 0.3 in normalized units, it means that in real units,
%% it is equal to 0.3*212.99nm of silver or 0.3*377nm for gold.
epsd=1;# dielectric of surrounding.
tol=1e-10;# Relative size of simplex for fmins function
%creates a non-unifom density of points of w
a=linspace(0.01,3,100);b=1-exp(-a); w=b/max(b)*1.25; % creates a non-uniform density of frequencies
%%from 0 to ~1.3 at which real and complex wavevectors will be calculated.
%%finds the solutions using fmins
for j=1:length(w)
epsm=1-(wp/w(j))^2; # The code uses Drude model for material and assumes no losses. plasma_frequency (omega_p) =1,
#Define first mode
Lp=@(kx)abs(epsm*(sqrt(epsd*(w(j)/c)^2-(kx(1)+i*kx(2))^2))+epsd*(sqrt(epsm*(w(j)/c)^2-(kx(1)+i*kx(2))^2))*tanh(-0.5*i*(sqrt(epsm*(w(j)/c)^2-(kx(1)+i*kx(2))^2))*d));
#solve for first mode using Nelder-Mead technique
[kx_Lp(j,:)] = fmins(Lp,[1.5; 0.25],[0,tol]); # 1.5 and 0.25 are the initial quesses for kx(1) and kx(2)
#Define other mode
Lm=@(kx)abs(epsm*(sqrt(epsd*(w(j)/c)^2-(kx(1)+i*kx(2))^2))+epsd*(sqrt(epsm*(w(j)/c)^2-(kx(1)+i*kx(2))^2))*coth(-0.5*i*(sqrt(epsm*(w(j)/c)^2-(kx(1)+i*kx(2))^2))*d));
#solve for other mode using Nelder-Mead technique
[kx_Lm(j,:)] = fmins(Lm,[1.5; 0.25],[0,tol]);# 1.5 and 0.25 are the initial quesses for kx(1) and kx(2)
end
%Plotting and printing section
putenv("GNUTERM",'wxt')
plot(abs(kx_Lp(:,1)),w,'ro',abs(kx_Lm(:,1)),w,'rs');
xlabel('Real(kx) /k_p')
ylabel('w/w_p')
legend('L+','L-','location','SouthEast')
title('Analytical solution to metal thin-film surface plasmon dispersion relation')
xlim([0,3]);
print('dmd_analytical_web.png','-dpng','-r200','-S600,600');
All entries of array except certain indices in octave/matlab
01 year ago
In Octave or Matlab, some times one needs to eliminate certain elements in an array.
For example, if
a=[10,20,30,40,50,60];
and suppose I want to create a matrix “b” such that it has all the elements of “a” except 20 and 40. This can be achieved by:
b=a(1:end~=2&1:end~=4);
The result is:
b =
10 30 50 60
10 30 50 60
Scattering / extinction / absorption cross-sections of silver nanowires (infinite cylinders) using meep
91 year ago
Particles scatter and absorb electromagnetic radiation. One often needs to compare the amount of scattering/absorption/extinction for particles of different shapes, composition, sizes and incident light properties (polarization, frequency and angle). In this regard, the concept of cross-sections comes into picture. There are three types of cross-sections, 1) scattering 2) absorption and 3) extinction. All of them have units of area, $m^2$, and provide a measure to quantify scattering/absorption process. Here using MEEP I calculate the crossections of silver nanowires and compare them with numerical solution (code from Bohren and Hauffman book).
Comparison of meep results with analytical results for silver nanowires
To achieve this, I wrote a meep code that performs a 2D simulation (x-y) with the cylinder axis along z axis with sources and monitors places as shown below. The source is a line source which travels along y direction and has polarization with electric field along x axis (along radius). I also use mirror-symmetries anti-symetry along X direction (which reduces the simulation time by half). PML layers are used on all sides. Calculation of cross-sections involves creating multiple 1-d flux monitors and running multiple simulations as shown below.
Here is my project file.. You would need ubuntu like system with meep and octave installed. Shell script will do all the work.
Gnuplot rgb color schemes in matlab as colormaps
22 years ago
I like the color schemes that are used as palettes in gnuplot’s pm3d plots. I wanted similar color schemes that can be used as colormaps in matlab.
After reading this article, I found that one can easily incorporate the traditional rgbpallette schemes in matlab
To start, there are several schemes in gnuplot that can be used in pm3d plots and they are:
7,5,15 ... traditional pm3d (black-blue-red-yellow) 3,11,6 ... green-red-violet 23,28,3 ... ocean (green-blue-white); try also all other permutations 21,22,23 ... hot (black-red-yellow-white) 30,31,32 ... color printable on gray (black-blue-violet-yellow-white) 33,13,10 ... rainbow (blue-green-yellow-red) 34,35,36 ... AFM hot (black-red-yellow-white)
The numbers represent different mapping functions for red, green and blue and the functions are
0: 0 1: 0.5 2: 1 3: x 4: x^2 5: x^3 6: x^4 7: sqrt(x) 8: sqrt(sqrt(x)) 9: sin(90x) 10: cos(90x) 11: |x-0.5| 12: (2x-1)^2 13: sin(180x) 14: |cos(180x)| 15: sin(360x) 16: cos(360x) 17: |sin(360x)| 18: |cos(360x)| 19: |sin(720x)| 20: |cos(720x)| 21: 3x 22: 3x-1 23: 3x-2 24: |3x-1| 25: |3x-2| 26: (3x-1)/2 27: (3x-2)/2 28: |(3x-1)/2| 29: |(3x-2)/2| 30: x/0.32-0.78125 31: 2*x-0.84 32: 4x;1;-2x+1.84;x/0.08-11.5 33: |2*x - 0.5| 34: 2*x 35: 2*x - 0.5 36: 2*x - 1 * negative numbers mean inverted=negative colour component * thus the ranges in `set pm3d rgbformulae' are -36..36
For example, if we choose, 7,5,15 , ie, traditional pm3d (black-blue-red-yellow), the three numbers correspond to functions used for red, green and blue respectively. So in this case, red(x)=sqrt(x), green(x)=x^3 and blue(x)=sin(2*pi*x), where x has range from 0 to 1. Note that red(x)/green(x)/blue(x) should not be negative numbers, so if they are negative make then zeros.
With that knowledge, here is a matlab code that creates 7,5,15 and 3,11,6 colormaps and uses it in the display of an image.
clear all
x=linspace(0,1,128);
r = sqrt(x);
g = x.^3.0;
b=sin(x*2*pi);
b(b<0)=0;
pm3d7515_bbry=[r;g;b]'
figure
load mandrill;
imagesc(X);
colormap(pm3d7515_bbry);
colorbar
r = x;
g = abs(x-0.5);
b=x.^4;
pm3d3116_grv=[r;g;b]'
figure
load mandrill;
imagesc(X);
colormap(pm3d3116_grv);
colorbar
x=linspace(0,1,128);
r = sqrt(x);
g = x.^3.0;
b=sin(x*2*pi);
b(b<0)=0;
pm3d7515_bbry=[r;g;b]'
figure
load mandrill;
imagesc(X);
colormap(pm3d7515_bbry);
colorbar
r = x;
g = abs(x-0.5);
b=x.^4;
pm3d3116_grv=[r;g;b]'
figure
load mandrill;
imagesc(X);
colormap(pm3d3116_grv);
colorbar
This results are shown below. Other color schemes can be similarly incorporated into matlab codes.
Arbitrary 2d shapes in MEEP
22 years ago
In MEEP (1.1.1), dielectric structures are often created by constructive geometry (adding and subtracting primitive shapes). The primitive shapes that are allowed are blocks, cylinders, ellipsoids and cones. To create a complex shape, one has to decompose the geometry into these primitive shapes. Over the weekend, I was wondering if it was possible to somehow create any complex shape in 2d without figuring out the exact positions and operations with the available primitive shapes. Here I report how I solve this problem.
The first thing I figured out was to create a 2d triangle with known vertices using a certain primitive shape. One can cut a cone to create a triangle, but that would limit me to isosceles triangle. My very smart colleague (Mengqian Lu) suggested me to cut a block/brick object with non-orthogonal axes such that you get a triangle. In meep, block object requires, 1) the lengths of the block edges along each of its three axes, 2) the directions of the axes of the block and 3)center. A triangle with known (x1,y1,0), (x2,y2,0) and (x3,y3,0) can be represented with a block (with possible non-orthogonality) element by modifications as shown in the figure. I pick a random point “E” which is not in the plane of my triangle and assume it as origin of the block. Then the vectors corresponding to axes of block are given by EA, EB, EC vectors, sizes are given by |EA|,|EB|,|EC| and the center of the block is OE+0.5*(EA+EB+EC).
Triangle in meep using block primitive object
To test this, I wrote a octave code which produces a meep code that that can be executed to create a triangle with vertices of (-10,0,0), (10,0,0) and (0,10,0). The result of compiling the meep code and exporting the dielectric distribution is shown below. The length of the simulation domain is 20 units in both x and y directions.
Triangle shape implemented in meep using block primitive object
With this information of incorporating a single triangle in meep, I can now obtain complex 2d shapes in meep. I start by writing a file that contains information regarding points that make up my structure and do a constrained delaunay triangulation. One can use matlab (versions >2009) to do such constrained delaunay triangulation. However, many of us (poor graduate students 
) do not have access to matlab. I figured out that constrained delaunay triangulation can be done by using Triangle package, which is free. After I compile this code, I see two executables 1) Triangle and 2) Showme in my source folder. Triangle does the triangulation and showme does the visualization before and after triangulation.
First, I write a “.poly” file, to mention my structure. ”.poly” file contains information of the points that make up the structure, one has to also mention segments (the border of the polygon) and any holes (if present). More information can be found here. Once a poly file is written, it can visualized using showme executable . If everything is ok, then do a constrained delaunay triangulation using Triangle executable (beware of concavities). I again visualize the structure after triangulation using showme executable to check whether triangulation has happened correctly or not. Once triangulated, I use the method above to convert each triangle to a corresponding meep block element (with suitable centers, sizes and axes) and stitch all of them together (once again use a octave code to do that automatically).
Using the poly file for shape “A” obtained from here I create a dielectric structure in MEEP with each triangle representing a random dielectric constant. Below is the final result. This method can also be extended to a group of arbitrary shapes in one simulation each with different dielectric functions.
Screenshot of showme executable showing the triangulation of 'A' shape
