% Scattering of surface plasmon polaritons by metal nanoparticles placed on
% a metal surface. This program is based on the theory described in the
% article "Phys. Rev. B, 71, 134304, Bozhevolnyi and Evlyukhin".
clear all


%% MISCELLANEOUS INITIALIZATIONS ------------------------------------------

ii = sqrt(-1);
nm = 10^-9;

e0 = 8.8541878176*10^-12;
em = -16.5;
er = 1.00;
ep = em;
a  = sqrt(er/-em);

scatrad = 50*nm;	% radius of scatterers, used in calculation of polarizabilities
effscatrad = 0*nm;	% "radius" of scatterers (points), used when calculating intensity maps

wavelength = 700*nm; % 700 nm, em = -16.5;
k0=2*pi/(wavelength);
kspp = k0*(er*em/(er+em))^(1/2);       % propagation constant

gaussianbeam = 1;           % 1 for gaussian beam, 0 for plane wave
beamwaist = 500*nm;
beamincomingy = 00*nm;
	
obsheight = scatrad; % observation height

StructureType = '';
LatticeConstant = 0;
RotAngle = 0;
StructureX = 0;
StructureY = 0;
NumScatterersX = 0;
NumScatterersY = 0;
MissingRows = 0;
LatticeConstant = 445*nm;

DyadXX = [[1 0 0]; [0 0 0]; [0 0 0]]; DyadXY = [[0 1 0]; [0 0 0]; [0 0 0]]; DyadXZ = [[0 0 1]; [0 0 0]; [0 0 0]];
DyadYX = [[0 0 0]; [1 0 0]; [0 0 0]]; DyadYY = [[0 0 0]; [0 1 0]; [0 0 0]]; DyadYZ = [[0 0 0]; [0 0 1]; [0 0 0]];
DyadZX = [[0 0 0]; [0 0 0]; [1 0 0]]; DyadZY = [[0 0 0]; [0 0 0]; [0 1 0]]; DyadZZ = [[0 0 0]; [0 0 0]; [0 0 1]];
DyadUnit = DyadXX + DyadYY + DyadZZ;

XSizePx = 201; % x and y sizes in pixels
YSizePx = 201;
XSizeNm = 30000*nm; % physical x and y size
YSizeNm = 10000*nm;
XOffsetNm = 0*nm;
YOffsetNm = 0;

WriteDataToFile = 0;
FilePath = 'Images\ParaChain\';
FileName = 'ParaChain';
FileExtension = 'bmp';


%% FILE NAMING SETUP ------------------------------------------------------

% This part is just for getting the files named correctly and not
% overwriting old images.
Counter = 1;
fn = sprintf('%s%s%u.%s',FilePath,FileName,Counter,FileExtension);
while(fopen(fn) ~= -1)
	Counter = Counter + 1;
	fn = sprintf('%s%s%u.%s',FilePath,FileName,Counter,FileExtension);
end
fclose('all');


%% POSITIONING OF SCATTERERS ----------------------------------------------

ScatStruct.x=0;
ScatStruct.y=0;
ScatStruct.z=0;
ScatStruct.rad=0;	% radius
ScatStruct.perm=0;	% permittivity

% Single Particle
N=1;
Scats=repmat(ScatStruct,2,1);
Scats(1).x=0*nm;
Scats(1).y=0*nm;
Scats(1).z=scatrad;
Scats(1).rad=scatrad;
Scats(1).perm=ep;

% Nanoparticle Parabolic Chain
parabheight = 20000*nm;
parabwidth  = 5000*nm;
paraba = -4*parabheight / parabwidth^2;
parabystep = 100*nm;
parabc = 0;
parabb = 10000*nm;

yy = -parabwidth/2;
steps = parabwidth/parabystep;
N = floor(steps)+1;
Scats=repmat(ScatStruct,N,1);
for i=1:N
	Scats(i).x = paraba*yy^2+parabb*yy+parabc;
	Scats(i).y = yy;
	Scats(i).z = scatrad;
	Scats(i).rad = scatrad;
	Scats(i).perm = ep;
	yy = yy + parabystep;
end



%% CALCULATIONS -----------------------------------------------------------

Mat= [];

% First we setup the matrix describing the relation between the incident
% field [E0] and the field at the position of the scatterers [E(ri)]
adenom = 2*(1-a^2)*(1-a^4);
for i=1:N
	for j=1:N
		if(i == j)
			Mat{i,j} = [1,0,0; 0,1,0; 0,0,1];
		else
			xi = (er - em)*(Scats(j).perm - er)/((er + em)*(Scats(j).perm + 2*er));
			beta = (Scats(j).rad / (2*Scats(j).z))^3;
			a0 = 4*e0*er*Scats(j).rad^3*pi*(Scats(j).perm-er)/(ep+2*er);
			alpha = a0 * (1/(1+xi*beta)*(DyadXX + DyadYY) + 1/(1+2*xi*beta)*DyadZZ);
			
			rx = Scats(i).x - Scats(j).x; % rx: relative x
			ry = Scats(i).y - Scats(j).y;
			rz = Scats(i).z - Scats(j).z;
			dxy = (rx^2+ry^2)^(1/2);
			cosphi = rx/dxy;
			sinphi = ry/dxy;
			distance = (rx^2+ry^2+rz^2)^(1/2);
			
			%h01 = besselh(0,1,distance*kspp);
			param = distance*kspp;
			h01 = sqrt(2/(pi*param))*exp(ii*(param-pi/4));
			
			Gspp = ii*a*kspp*h01*exp(-a*kspp*(Scats(i).z+Scats(j).z))/adenom*[(a^2*cosphi^2 - ii*a*cosphi), (a^2*cosphi*sinphi), 0;
				(a^2*sinphi*cosphi - ii*a*sinphi), a^2*sinphi^2, 0;
				(ii*a*cosphi), ii*a*sinphi, 1];

			Mat{i,j} = k0^2/e0 * Gspp*alpha;
		end
	end
end

invMat=[];
for i=1:N
	invMat=vertcat(invMat,horzcat([],Mat{i,1:N}));
end

if(N > 0)
	invMat = inv(invMat);
end

% Initiate the vector of incident fields (at the positions of the
% scatterers) which we have to multiply by the (inverted) matrix invMat
E0S = [];
for i=1:N
% 	f = 1 + ii*2*kspp*Scats(i).x/(kspp*beamwaist)^2;
% 	kz = sqrt(k0^2*er - kspp^2);
% 	fdenom = (kspp*beamwaist)^2*f;
% 	E0S = vertcat(E0S, exp(ii*kspp*Scats(i).x + ii*kz*Scats(i).z -Scats(i).y^2/(beamwaist^2*f))/sqrt(f)*[(-kz/kspp*(1+1/fdenom-2*(kspp*Scats(i).y)^2/fdenom^2)); (-ii*2*Scats(i).y*kz/fdenom); 1]);
	E0S = vertcat(E0S, exp(ii*kspp*Scats(i).x - a*kspp*Scats(i).z)*[-i*a; 0; 1]);
end

if(N > 0)
    ES = invMat*E0S; % ES: E Scatterer, self consistent field at scatterer positions
end

% Now that we know the fields at the positions of the scatterers 
% we can calculate the field at any position
FieldIntensity = zeros(YSizePx,XSizePx);
for y=1:YSizePx
    yy = YSizeNm / YSizePx * (y-1) - YSizeNm / 2 + YOffsetNm;
    for x=1:XSizePx
        xx = XSizeNm / XSizePx * (x-1) - XSizeNm / 2 + XOffsetNm;
		zz = obsheight;
        
		E = exp(ii*kspp*xx-a*kspp*zz)*[(-ii*a); 0; 1];
		E = 0;
		
% 		f = 1 + ii*2*kspp*xx/(kspp*beamwaist)^2;
% 		kz = sqrt(k0^2*er - kspp^2);
% 		fdenom = (kspp*beamwaist)^2*f;
% 		E = exp(ii*kspp*xx + ii*kz*zz -yy^2/(beamwaist^2*f))/sqrt(f)*[(-kz/kspp*(1+1/fdenom-2*(kspp*yy)^2/fdenom^2)); (-ii*2*yy*kz/fdenom); 1];
 		
		
        for l=1:N
            distance = ((Scats(l).x-xx)^2 + (Scats(l).y-yy)^2 + (Scats(l).z-zz)^2)^(1/2);
            
			xi = (er - em)*(Scats(l).perm - er)/((er + em)*(Scats(l).perm + 2*er));
			beta = (Scats(l).rad / (2*Scats(l).z))^3;
			a0 = 4*e0*er*Scats(l).rad^3*pi*(Scats(l).perm-er)/(ep+2*er);
			alpha = a0 * (1/(1+xi*beta)*(DyadXX + DyadYY) + 1/(1+2*xi*beta)*DyadZZ);

			rx = Scats(l).x - xx;
			ry = Scats(l).y - yy;
			dxy = (rx^2 + ry^2)^(1/2);
			cosphi = rx/dxy;
			sinphi = ry/dxy;
			
			%h01 = besselh(0,1,distance*kspp)
			param = distance*kspp;
			h01 = sqrt(2/(pi*param))*exp(ii*(param-pi/4));
			Gspp = ii*a*kspp*h01*exp(-a*kspp*(Scats(l).z + zz))/adenom*[(a^2*cosphi^2 - ii*a*cosphi), (a^2*cosphi*sinphi), 0;
				 (a^2*sinphi*cosphi - ii*a*sinphi), a^2*sinphi^2, 0;
				 (ii*a*cosphi), ii*a*sinphi, 1];

			EScurr = [ES((l-1)*3+1); ES((l-1)*3+2); ES((l-1)*3+3)];
            E = E + k0^2/e0*Gspp*alpha*EScurr;
		end        
        FieldIntensity(y,x)=abs(E(3))^2;
	end
	y/YSizePx
end


%% DATA PROCESSING --------------------------------------------------------

MaxIntensity=0;
for y=1:YSizePx
	for x=1:XSizePx		
		if(FieldIntensity(y,x)>MaxIntensity)
			MaxIntensity=FieldIntensity(y,x);
		end
	end
end
%MaxIntensity=2.0;

%% MARK SCATTERER POSITIONS -----------------------------------------------

for i=1:N
	x=int16(Scats(i).x*XSizePx/XSizeNm+XSizePx/2);
	y=int16(Scats(i).y*YSizePx/YSizeNm+YSizePx/2);
	%FieldIntensity(y-1,x-1)=255; 	FieldIntensity(y-1,x)=255;	FieldIntensity(y-1,x+1)=255;
	%FieldIntensity(y,x-1)=255;		FieldIntensity(y,x)=255;	FieldIntensity(y,x+1)=255;
	%FieldIntensity(y+1,x-1)=255; 	FieldIntensity(y+1,x)=255;	FieldIntensity(y+1,x+1)=255;
	FieldIntensity(y,x)=MaxIntensity;
end

%% OUTPUT TO FILE ---------------------------------------------------------

colormap(hot(256));
fn = sprintf('%s%s%u.%s',FilePath,FileName,Counter,FileExtension);
imagesc(FieldIntensity,[0,MaxIntensity]);
colorbar;

if(WriteDataToFile==1)	
	imwrite(FieldIntensity/MaxIntensity*255,colormap, fn, FileExtension);
	
	fn = sprintf('%s%s%uFull.%s',FilePath,FileName,Counter,'eps');
	print('-depsc',fn);

	fn = sprintf('%s%s%u.%s',FilePath,FileName,Counter,'txt');
	file = fopen(fn,'w+');

	s='StructureType: %s\r\n';

	numparams=22;
	for i=1:numparams
		s = sprintf('%s%s',s,'%s: %e\r\n');
	end

	text = sprintf(s,StructureType,'gaussianbeam',gaussianbeam,'beamwaist',beamwaist/nm, ...
		'beamincomingy',beamincomingy/nm,'XSizeNm',XSizeNm/nm,'YSizeNm',YSizeNm/nm,'XSizePx',XSizePx, ...
		'YSizePx',YSizePx,'XOffsetNm',XOffsetNm,'YOffsetNm',YOffsetNm,'real(em)',real(em),'imag(em)',imag(em),'er',er, ...
		'wavelength',wavelength/nm, 'LatticeConstant', LatticeConstant, 'RotAngle', RotAngle, ...
		'NumScatterersX', NumScatterersX, 'NumScatterersY', NumScatterersY, 'StructureX', StructureX, ...
		'StructureY', StructureY, 'effscatrad',effscatrad, 'MissingRows', MissingRows);
	fwrite(file,text);
	fclose(file);
end








