% Elastic (in plane) scattering of surface plasmons from particles on the
% surface. The program is based on the theory described in "Elastic
% scattering of surface plasmon polaritons: Modelling and experiment" by
% Bozhevolnyi and Coello, Physical Review B, 58, 10899
%
% The incident field is always incident from the negative x direction. The
% imaginary unit is called z throughout this program.

clear Mat ES E0S N alpha e1 e2 k0 beta FieldIntensity;

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

nm = 10^-9;

z = sqrt(-1);
alpha = 3;      % Polarizability (as given in Bozevolnyi and Coello)
e1=-16+z;
e2=1.00;
wavelength = 400*nm;
k0=2*pi/(wavelength);
beta = k0*(e1*e2/(e1+e2))^(1/2);       % propagation constant

XSizePx = 1200; % x and y sizes in pixels
YSizePx = 600;
XSizeNm = 12000*nm; % physical x and y size
YSizeNm = 6000*nm;
XOffsetNm = 0;
YOffsetNm = 0;

MaxIntensity = 0;
MeanIntensity = 0;

ScatRadius = 2*nm;

% Beam parameters
GaussianBeam = 0;           % 1 for gaussian beam, 0 for plane wave
BeamWaist = 1000*nm;
BeamIncomingY = 00*nm;


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

% This part is just for getting the files named correctly and not
% overwriting old images.
FilePath = 'Images\';
FileName = 'Chain';
FileExtension = 'bmp';
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
%----------------------------------------------------------------------

% N=1;
% ScatPosX(1)=8000*nm;
% ScatPosY(1)=0;

% N=300;
% k=1;
% for i=1:30
% 	for j=1:10
% 		ScatPosX(k)=(j-1)*300*nm-4000*nm;
% 		ScatPosY(k)=(i-1)*300*nm-4500*nm;		
% 		if(i > 15)
% 		%	ScatPosY(k)=ScatPosY(k)+300*3*nm;
% 		end
% 		k=k+1;
% 	end
% end

% % Mirror
% k=1;
% N=250;
% angle = pi/180*60;
% scatdist = 220*nm;
% rowdistance = 350*nm;
% for i=1:5
% 	for j=1:50
% 		ScatPosX(k)=cos(angle)*scatdist*(j-1) + rowdistance*(i-1) - 1000*nm;
% 		ScatPosY(k)=sin(angle)*scatdist*(j-1) - 8000*nm;
% 		k=k+1;
% 	end
% end

% Parabola - Nanoparticle Chain
width=1400*nm;
c = 5000*nm;
a = -4*abs(c)/width^2;
posx = -c;
N=55;
for i=1:N
	ScatPosY(i) = -width/2 + (i-1)*(width)/(N-1);
	ScatPosX(i) = a*ScatPosY(i)^2 + c + posx;
end

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

% First we setup the matrix describing the relation between the incident
% field [E0] and the field at the position of the scatterers [E(ri)]
for i=1:N
    for j=1:N
        if(i == j)
            Mat(i,j) = 1;
        else
            p = beta * ( (ScatPosX(i) - ScatPosX(j))^2 + (ScatPosY(i) - ScatPosY(j))^2)^(1/2);
            Mat(i,j) = -alpha*(1/4)*besselh(0,1,p);
        end
    end
end

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

% Initiate the vector of incident fields (at the positions of the
% scatterers) which we have to multiply by the (inverted) matrix Mat
for i=1:N
    if(GaussianBeam)
        E0S(i,1) = exp(z*beta*ScatPosX(i)-(ScatPosY(i)-BeamIncomingY)^2/BeamWaist^2); % E0S: E0 Scatterer
    else
        E0S(i,1) = exp(z*beta*ScatPosX(i)); % E0S: E0 Scatterer
    end
end

if(N > 0)
    ES = Mat*E0S; % ES: E Scatterer
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;
        
        if(GaussianBeam)
            E = exp(z*beta*xx-((yy-BeamIncomingY)/BeamWaist)^2);
        else
            E = exp(z*beta*xx);
        end

        for l=1:N
            distance = ( (xx - ScatPosX(l))^2 + (yy - ScatPosY(l))^2)^(1/2);
            
            % This is to avoid calculating the fields at the position of
            % the scatterer itself. Doing so leads to extremely high
            % intensities at the positions of the scatterer. ScatRadius is
            % an effective radius and not a real radius because the
            % scatterers are considered to be point particles in this
            % model.
            if(distance < ScatRadius)
                continue;
            end
                
            p = beta * distance;
            %E = E + alpha*ES(l)*(z/4)*besselh(0,1,p);
            E = E + alpha*ES(l)*z*sqrt(1/(8*pi*p))*exp(z*(p-pi/4));
        end        
        FieldIntensity(y,x)=abs(E)^2;
	end
	y/YSizePx
end


%----------------------------------------------------------------------
% DATA PROCESSING AND OUTPUT TO FILE
%----------------------------------------------------------------------


% To avoid extremely high and non-physical field intensities around the
% scatterers we average over an area around each scatterer.
AreaDimensionNm = 100*nm;
PxPerNmX = XSizePx / XSizeNm;
PxPerNmY = YSizePx / YSizeNm;
AreaDimensionPxX = AreaDimensionNm*PxPerNmX;
AreaDimensionPxY = AreaDimensionNm*PxPerNmY;
for i=1:N
	PosPxX1 = uint16(ScatPosX(i)*PxPerNmX - AreaDimensionPxX/2 + PxPerNmX*XSizeNm/2);
	PosPxY1 = uint16(ScatPosY(i)*PxPerNmY - AreaDimensionPxY/2 + PxPerNmY*YSizeNm/2);
	PosPxX2 = uint16(PosPxX1 + AreaDimensionPxX);
	PosPxY2 = uint16(PosPxY1 + AreaDimensionPxY);
	
	temp=FieldIntensity;
	for y=PosPxY1:PosPxY2
		for x=PosPxX1:PosPxX2
			temp(y,x) = (FieldIntensity(y-1,x-1) + FieldIntensity(y-1,x) + FieldIntensity(y-1,x+1) + FieldIntensity(y,x-1) + FieldIntensity(y,x) + FieldIntensity(y,x+1) + FieldIntensity(y+1,x-1) + FieldIntensity(y+1,x) + FieldIntensity(y+1,x+1))/9;
		end
	end
	%FieldIntensity=temp;
end

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

for y=1:YSizePx
    for x=1:XSizePx
        FieldIntensity(y,x) = 255 / (MaxIntensity) * FieldIntensity(y,x);
    end
end

colormap(hot(256));
fn = sprintf('%s%s%u.%s',FilePath,FileName,Counter,FileExtension);
imwrite(FieldIntensity,colormap, fn, 'bmp');
image(imread(fn, 'bmp'));

fn = sprintf('%s%s%u.%s',FilePath,FileName,Counter,'txt');
file = fopen(fn,'w+');
text = sprintf('%s: %e\r\n%s: %e\r\n%s: %e\r\n%s: %e\r\n%s: %e\r\n%s: %e\r\n%s: %e\r\n%s: %e\r\n%s: %e\r\n%s: %e\r\n%s: %e\r\n%s: %e\r\n','GaussianBeam',GaussianBeam,'BeamWaist',BeamWaist/nm,'BeamIncomingY',BeamIncomingY/nm,'XSizeNm',XSizeNm/nm,'YSizeNm',YSizeNm/nm,'XSizePx',XSizePx,'YSizePx',YSizePx,'real(e1)',real(e1),'imag(e1)',imag(e1),'e2',e2,'alpha',alpha,'wavelength',wavelength/nm);
fwrite(file,text);
fclose(file);








