Installing Meep 1.2 on ubuntu


Pre-compiled Meep binaries for meep1.1 exist for Ubuntu distribution. This makes it very easy to install meep on ubuntu using “apt-get install” command or from the ubuntu software center. However recently, Meep developers have release meep1.2 which has more functions compared to meep1.1. I have recently installed meep1.2 from source on ubuntu 12.04 using the instructions shown at

I have root access to my computer, so I installed all the libraries/bin files in their default location (i.e, libraries go in /usr/local/lib, programs in /usr/local/bin, etc)

These are the steps I followed:

1) To avoid any complications, I uninstalled meep1.1 and meep-openmpi, which were installed earlier on my ubuntu system from the software center.

2) Setting the environmental variables:
export LDFLAGS=”-L/usr/local/lib”
export CPPFLAGS=”-I/usr/local/include”
export LD_LIBRARY_PATH=”/usr/local/lib:$LD_LIBRARY_PATH”

3) I installed G77, gfortran, f77, g++ and make packages on my system. I installed g++ and make using “sudo apt-get install build-essential”. Some of the packages might have been pre-installed.

4)Installing blas and lapack

As said at, Openblas at should install both blas and lapack libraries. I downloaded the tar file from the website and extracted it. Then I cd’ed into that folder and issued a “make” command. I had some problems with “make” command as it was using f77 for compiling and was throwing an error. After some search, I found that “make FC=gfortran NO_AFFINITY=1 USE_OPENMP=1” solves this problem. This resulted in “libopenblas.a” in /usr/local/lib folder.

5) Installing Harminv
I downloaded the tar file from Extracted the tar file and cd’ed into it. I installed the harminv using ./configure, make and sudo make commands as said at ” title=””> This resulted in “libharminv.a”. Harminv required blas and lapack installed.

6) Installing Guile
This was simple. I used the command “sudo apt-get install guile-2.0” and “sudo apt-get install guile-2.0-dev”

7) Installing libctl
Downloaded the tar file, extracted it, cd’ed into it. I followed ./configure, make and sudo make commands as told on at and after the process I saw libctl.a in the /usr/local/lib folder.

8) Installing hdf5
I downloaded the hdf5 tar file from hdf5-1.8.8 source page. I used 1.8.8 version (no special reason for this version). Extracted the tar file. Cd’ed into the folder and use “./configure, make, make install” . It took a while for this to complete. At the end, I saw all the libraries, executables and include files were placed in “hdf5” sub-folder. I saw “bin”,”include”,”share”,”lib” inside this sub-folder. I figured they should be placed at default locations so I issued following commands inside the hdf5-1.8.6 folder
sudo cp hdf5/bin/* /usr/local/bin/
sudo cp hdf5/lib/* /usr/local/lib/
sudo cp hdf5/include/* /usr/local/include/
sudo cp hdf5/share/* /usr/local/share/
By the end of the process, I saw libhdf5.a in the /usr/local/lib. I did not enable MPI parallel I/O support in hdf5.

9) Installing meep
I downloaded the tar file from meep downloaded page and then extracted it. Then I cd’ed into that folder and issued “./configure” command. I had no errors as guile, libctl, harminv, hdf5 libraries were all available. See my configure_log_file. Then “make” and then “sudo make install” commands were issued. I did not enabled parallel version for this post. This resulted in meep executable in “/usr/local/bin” folder. I renamed it to meep12 (sudo mv /usr/local/bin/meep meep12) so that I know that it is running meep1.2 and not to confuse with meep1.1.  Now I opened a new terminal and issued “/usr/local/bin/meep12” to run meep1.2. you should see:
“GNU Guile 2.0.5-deb+1-1
Copyright (C) 1995-2012 Free Software Foundation, Inc.
Guile comes with ABSOLUTELY NO WARRANTY; for details type `,show w’.
This program is free software, and you are welcome to redistribute it
under certain conditions; type `,show c’ for details.
Enter `,help’ for help.

In future, I will post some meep examples that use the added functionalities of meep1.2. Stay tuned!

Maximum electric field intensity near silver nanoparticle as function of wavelength

Electric field at localized plasmon resonance using MEEP


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:

      1. Calculate the extinction, scattering, absorption spectra of metal nanoparticle

The procedure for doing this is very similar to the method I mentioned here.

      1. 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.
      2. Calculate electric field at plasmon resonance as function of timeThis was very easy. Get the electric fields at certain wavelength and then multiply it with a time varying sin function.

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.



Update1: You can find my project file HERE. Enjoy!
(please remember that the project file is designed to run on ubuntu like system with meep and octave installed. The shell script does everything needed to get the results, make sure to “chmod” the shell file, so that you can run the file. In addition, also make sure to increase the resl (resolution) from 10 to 50 for decent results. This means it will take longer time to run the simulation. The results you see in this post were generated by resl=100 and it took 2 days to run on my computer. If you find any bugs, please let me know.)


Update2: This method also works well for gold nano-particles of 50 nm radius. See the last figure in above album.




Building my new desktop system


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

Parallelization in Octave using parcellfun/pararrayfun


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)



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 ()

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


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

Van der Pauw correction factor

van der Pauw correction factor


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)
F(i) = fsolve(Lp,0.5);
ylabel('Correction Factor (f)')
title('van der Pauw Correction factors for resistivity measurements');
Analytically solution to metal thin film surface plasmon dispersion relation

Surface plasmon dispersion relation for thin metal films


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, k_x^{'}(\omega)+ik_x^{''}(\omega)) are given by the solutions of :

\omega^+: L+= \epsilon_1k_{z2}+\epsilon_2k_{z1} \tanh\left [  \frac{-ik_{z1}d}{2}\right ]=0

\omega^-: L-= \epsilon_1k_{z2}+\epsilon_2k_{z1} \coth\left [  \frac{-ik_{z1}d}{2}\right ]=0,

Where \epsilon_1 and \epsilon_2 represent the dielectric functions of dielectric and metal respectively, d is the thickness of metal film, \omega^+ is frequency of high energy mode, \omega^- 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}.

What does solving these equations mean? if one considers one mode, say L+ mode, for a certain \omega^+, there exists a particular k_x^{'} and k_x^{''} that will make the above left hand side of L+ equation equal to zero. How do we get those special k? 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 k_{z1,z2} at a particular \omega, such that L+ and L- 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 (, 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
#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
#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)

%Plotting and printing section
xlabel('Real(kx) /k_p')
title('Analytical solution to metal thin-film surface plasmon dispersion relation')

All entries of array except certain indices in octave/matlab


In Octave or Matlab, some times one needs to eliminate certain elements in an array.

For example, if


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:


The result is:

b =

   10   30   50   60
Comparision of meep results with analytical results for silver nanowires

Scattering / extinction / absorption cross-sections of silver nanowires (infinite cylinders) using meep


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.

I have also shared the project at github.

colormap of pm3d 3,11,6 (green-red-violet)

Gnuplot rgb color schemes in matlab as colormaps


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
r = sqrt(x);
g = x.^3.0;
load mandrill;

r = x;
g = abs(x-0.5);
load mandrill;

This results are shown below. Other color schemes can be similarly incorporated into matlab codes.

Go to Top