func_nf.m

%near field absorption cross-section conversion for fixed distance from the surface vs. wavenumber (starting + 400)
%disp ('etot func_nf');
clear global
nf_mode = 1;

disp('...');disp('0 = default    1 = gold');disp('2 = platinum   3 = silver');disp('4 = silicon    5 = silicon fit');disp('6 = silicon nitride');disp('7 = coated tip');
disp('...');material1 = input ('Enter the tip material code:  ');
if material1 == 7,
	%disp('...');disp('0 = default    1 = gold');disp('2 = platinum   3 = silver');disp('4 = silicon    5 = silicon fit');disp('6 = silicon nitride   7 = iron carbonyl');
	disp('...');material_0 = input ('Enter tip core material code:  ');
	material_1 = input ('Enter tip coating material code (except for iron carbonyl):  ');
	ad = input ('Enter tip core radius (10 to 100), nm:  ');
	a = ad * 1.000E-7 + input ('Enter the tip coating thickness (25), nm:  ') * 1.000E-7;
   	%a = ad * 1.000E-7 + 25 * 1.000E-7;
   Rparticle1 = a;
else
   %type in the end of tip radius of curvature (cm)
	Rparticle1 = input ('Enter the end-of-tip radius of curvature (10 to 120), nm:  ') * 1.000E-7;
end

disp('...');disp('0 = default    1 = gold');disp('2 = platinum   3 = silver');disp('4 = silicon    5 = silicon fit');disp('6 = silicon nitride');disp('7 = reference material');
disp('...');material2 = input ('Enter sample material code:  ');
if material2 == 7,
	material2 = input ('Enter the sample material code (0 for default, 7 for iron carbonyl):  ');
   
   if material2 == 7,
		ic_shapka
	else
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   REFparticle = inf;
   %vFermi = input ('Enter the sample vFermi (e.g., 1.4E8 for noble metals):  ');
   vFermi = 1.4E8;
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
   resonance_type = 2;
   %omega_r = 6E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   omega_r = 6E10 * pi * 1815;
   %REFgamma = 6E10 * pi * input ('Enter the reference gamma, cm-1:  ');
   REFgamma = 6E10 * pi * 15;
   Mif = 1E30 * input ('Enter absorption strength (1 to 1000):  ');
   %sigma_eps =  1E15 * input ('Enter sigma, typically much less than 1:  ');
   sigma_eps = 0;
   trial_number = input('Enter the model number:  ');
	end
end

%type in the particle radius (cm)
%Rparticle2 = input ('Enter sample particle radius (1 to inf), nm:  ');
	Rparticle2 = inf;

%type in the tip-surface distance (cm)
d = 1.000E-7 * input ('Enter tip-sample distance, nm:   ');
	%d = 0;

%nu = input ('Starting wavenumber:  ');
%nu = (100 * round (omega_r/(6.000E12 * pi)) - 100);

poglot = input ('Choose 0 for epsilon plot, 1 for NF absorption, 2 for NF scattering, 3 for normal reflectance:  ');
if poglot == 0,
   kakmany = input ('Enter 0 for combination plot, 1 for ref, 2 for kappa*nu:  ');
end

%figure %plot function vs. nu
hold on

nomer = 1;
plot_function

done = 1;%input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
if poglot == 0,
   kakmany = input ('Enter 0 for combination plot, 1 for ref, 2 for kappa*nu:  ');
end
   if material2 == 7,
   	ic_shapka   
   else
   	trial_number = input('Enter the model number:  ');
   	Mif = input ('Enter absorption strength (1 to 1000):  ');
   end
   nomer = 2;
   plot_function
   clear nu

done = 1;%input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   if material2 == 7,
   	ic_shapka   
   else
   	trial_number = input('Enter the model number:  ');
   	Mif = input ('Enter absorption strength (1 to 1000):  ');
   end
   nomer = 3;
	plot_function
   clear nu

done = 1;%input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   if material2 == 7,
   	ic_shapka   
   else
   	trial_number = input('Enter the model number:  ');
   	Mif = input ('Enter absorption strength (1 to 1000):  ');
   end
   nomer = 4;
	plot_function
end
end
end
hold off






plot_function.m

%function vs. wavenumber
stpn = 0;
nu = 1700;%00;
%nu = (100 * round (omega_r/(6.000E12 * pi)) - 100);
%hold on
%   epsilonrej = zeros (1, 2000);
%   epsilonimj = zeros (1, 2000);
   %epsilonimj2 = zeros (1, 2000);
   %epsilonimj3 = zeros (1, 2000);
   %epsilonimj4 = zeros (1, 2000);
%   nuj = zeros (1, 2000);
%   refj = zeros (1, 2000);
%   kappaj = zeros (1, 2000);
%   responsej = zeros (1, 2000);
if nomer == 1,
	for j=1:351;%2000;
eps_ic
      funkcia
%responsej(j) = (abs ((sqrt(epsilonre+i*epsilonim) - 1.33)/(sqrt(epsilonre+i*epsilonim) + 1.33))) ^ 2;
%responsej(j) = 0.1+70*(((real(sqrt(epsilonre+i*epsilonim)) - 2)/(real(sqrt(epsilonre+i*epsilonim)) + 2))) ^ 2;
%self-homodyned%responsej(j) = abs(((real(sqrt(epsilonre+i*epsilonim)) - 1)/(real(sqrt(epsilonre+i*epsilonim)) + 1)));
%      epsilonrej(j)=epsilonre;
%      epsilonimj(j)=epsilonim;
      %epsilonimj2(j)=epsilonim2;
      %epsilonimj3(j)=epsilonim3;
      %epsilonimj4(j)=epsilonim4;
%      kappa=epsilonim;
%      kappa = imag(sqrt(epsilonre+i*epsilonim))%;%sqrt(sqrt((epsilonre/2)^2+(epsilonim/2)^2)-epsilonre/2);
      %ref=epsilonre;
%      ref = real(sqrt(epsilonre+i*epsilonim));%epsilonim/(2*kappa);%
%      refj(j) = -LOG10(1-(((ref-1)/(ref+1))^2+bgnd_ref));
%      kappa = 4*pi*kappa*nu;
%      kappaj(j) = LOG10(kappa*nu*0.008/7.86);
%   vv=ref^2+kappa^2+1;
responsej(j) = response;%(vv-2*ref)/(vv+2*ref);

      nuj(j) = nu;
%   	clear funkcia
   	nu = nu + 1;%0.1;
   end
   if poglot ~= 0,
hold on
       plot (nuj, responsej, 'k.', 'markersize', 4);
%      figure
%      plot (nuj, refj, 'k.', 'markersize', 5);
%      hold on
%      plot (nuj, kappaj, 'k.', 'markersize', 3);
   	if poglot == 1,
			xlabel('wavenumber, [cm-1]'); ylabel('absorption response, [a.u.]');
   	end
   	if poglot == 2,
			xlabel('wavenumber, [cm-1]'); ylabel('scattering response, [a.u.]');
   	end
%   	if poglot == 3,
%			xlabel('wavenumber, [cm-1]'); ylabel('Reflection coefficient');
%   	end
   else
     	if kakmany == 0,
         plot (nuj, epsilonrej, 'k.-', 'markersize', 5);
         %axis ([1700, 2300, -2, 2]);  
         hold on  
      	plot (nuj, epsilonimj, 'k.--', 'markersize', 3);
      	%plot (nuj, epsilonimj2, 'k*', 'markersize', 3);
      	%plot (nuj, epsilonimj3, 'k*', 'markersize', 4);
      	%plot (nuj, epsilonimj4, 'k*', 'markersize', 5);
         %hold on
         xlabel('wavenumber, [cm-1]'); ylabel('epsilon, Re - and Im --');
     	end
      if kakmany == 1,
         plot (nuj, refj, 'k.', 'markersize', 5);%epsilonre
         hold on
			xlabel('wavenumber, [cm-1]'); ylabel('reflectance');
   	end
     	if kakmany == 2,
         plot (nuj, kappaj, 'k--', 'markersize', 5);%epsilonim
         hold on  
			xlabel('wavenumber, [cm-1]'); ylabel('absorbance');
      end
   end
	axis square;
end

if nomer == 2,
	for j=1:2000;
   	funkcia
      epsilonrej(j)=epsilonre;
      epsilonimj(j)=epsilonim;
      kappa = imag(sqrt(epsilonre+i*epsilonim));%sqrt(sqrt((epsilonre/2)^2+(epsilonim/2)^2)-epsilonre/2);
      ref = real(sqrt(epsilonre+i*epsilonim));%epsilonim/(2*kappa);%
      refj(j) = -log10(1-(((ref-1)/(ref+1))^2+bgnd_ref));
      kappaj(j) = kappa*nu*0.008/7.86;
      nuj(j) = nu;
   	clear funkcia
		nu = nu + 0.3;
	end
   if poglot ~= 0,
      plot (nu, response, '.', 'markersize', 6);
   else
      if kakmany == 0,
      	plot (nuj, epsilonrej, '.-', 'markersize', 6);
      	plot (nuj, epsilonimj, '.--', 'markersize', 4);
      end
      if kakmany == 1,
      	plot (nuj, refj, '.', 'markersize', 6);
		end
      if kakmany == 2,
      	plot (nuj, kappaj, '--', 'markersize', 6);
      end
   end
end

if nomer == 3,
	for j=1:2000;
		funkcia
      epsilonrej(j)=epsilonre;
      epsilonimj(j)=epsilonim;
      kappa = imag(sqrt(epsilonre+i*epsilonim));%sqrt(sqrt((epsilonre/2)^2+(epsilonim/2)^2)-epsilonre/2);
      ref = real(sqrt(epsilonre+i*epsilonim));%epsilonim/(2*kappa);%
      refj(j) = -log10(1-(((ref-1)/(ref+1))^2+bgnd_ref));
      kappaj(j) = kappa*nu*0.008/7.86;
      nuj(j) = nu;
   	clear funkcia
		nu = nu + 0.3;
	end
   if poglot ~= 0,
      plot (nu, response, 'm.', 'markersize', 9);
   else
      if kakmany == 0,
      	plot (nuj, epsilonrej, 'm.-', 'markersize', 7);
      	plot (nuj, epsilonimj, 'm.--', 'markersize', 5);
      end
      if kakmany == 1,
         plot (nuj, refj, 'm.', 'markersize', 9);
      end
      if kakmany == 2,
      	plot (nuj, kappaj, 'm--', 'markersize', 9);
      end
   end
end

if nomer == 4,
	for j=1:2000;
		funkcia
      epsilonrej(j)=epsilonre;
      epsilonimj(j)=epsilonim;
      kappa = imag(sqrt(epsilonre+i*epsilonim));%sqrt(sqrt((epsilonre/2)^2+(epsilonim/2)^2)-epsilonre/2);
      ref = real(sqrt(epsilonre+i*epsilonim));%epsilonim/(2*kappa);%
      refj(j) = -log10(1-(((ref-1)/(ref+1))^2+bgnd_ref));
      kappaj(j) = kappa*nu*0.008/7.86;
      nuj(j) = nu;
   	clear response funkcia
		nu = nu + 0.3;
	end
   if poglot ~= 0,
      plot (nu, response, 'c*', 'markersize', 7);
   else
      if kakmany == 0,
      	plot (nuj, epsilonrej, 'c.-', 'markersize', 8);
      	plot (nuj, epsilonimj, 'c.--', 'markersize', 6);
      end
      if kakmany == 1,
      	plot (nuj, refj, 'c*', 'markersize', 7);
   	end
      if kakmany == 2,
      	plot (nuj, kappaj, 'c--', 'markersize', 7);
      end
   end
end






nf.m

%near field signal plots
clear global
if nf_mode ~= 0,
	nu = input ('Enter wavenumber (800-1950), cm-1:  ');
	k = 2 * pi * nu;
	omega = k * 3.000E10;
end

v3 = input('Change size parameters? 1=Yes 0=No:  ');
if v3 == 1,
	if material1 == 7,
     	ad = input ('Enter tip core radius (10 to 100), nm:  ');
     	a = ad * 1.000E-7 + input ('Enter tip coating thickness (25), nm:  ') * 1.000E-7;
      Rparticle1 = a;
   else
      Rparticle1 = input ('Enter the end-of-tip radius of curvature (10 to 120), nm: ') * 1.000E-7;
   end
else
	if v3 == 1,
		if material1 == 7,
     		ad = input ('Enter tip core radius (10 to 100), nm:  ');
     		a = ad * 1.000E-7 + input ('Enter tip coating thickness (25), nm:  ') * 1.000E-7;
      	Rparticle1 = a;
   	else
      	Rparticle1 = input ('Enter the end-of-tip radius of curvature (10 to 120), nm: ') * 1.000E-7;
   	end
   end
end

if v3 == 1,
   Rparticle2 = input ('Enter sample particle radius (1 to inf), nm:  ');
   dmin = 1.000E-7 * input ('Enter tip-sample distance, nm:   ');
   d = dmin;
	ao = 1.000E-7 * input ('Enter amplitude of oscillation, nm:   ');
end

v1 = input('Change tip material? 1=Yes 0=No:  ');
if v1 == 1,
	disp('...');disp('0 = default   1 = gold');disp('2 = platinum   3 = silver');disp('4 = silicon   5 = silicon fit');disp('6 = silicon nitride');disp('7 = coated tip');
	disp('...');material1 = input ('Enter the tip material code:  ');
   
   if material1 == 0,
	%material1 = input ('Enter the sample material code (8 for default, 0 for iron carbonyl):  ');
   	%double-check
   	if material1 == 0,
         ic_shapka
		else
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   REFparticle = inf;
   %vFermi = input ('Enter the sample vFermi (e.g., 1.4E8 for noble metals):  ');
   vFermi = 1.4E8;
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
   resonance_type = 2;
   %omega_r = 6E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   omega_r = 6E10 * pi * 1815;
   %REFgamma = 6E10 * pi * input ('Enter the reference gamma, cm-1:  ');
   REFgamma = 6E10 * pi * 15;
   Mif = 1E30 * input ('Enter absorption strength (1 to 1000):  ');
   %sigma_eps =  1E15 * input ('Enter sigma, typically much less than 1:  ');
   sigma_eps = 0;
   trial_number = input('Enter the model number:  ');
		end
	end

	if material1 == 7,
		disp('...');disp('0 = default   1 = gold');disp('2 = platinum   3 = silver');disp('4 = silicon   5 = silicon fit');disp('6 = silicon nitride   7 = iron carbonyl');
		disp('...');material_0 = input ('Enter tip core material code:  ');
		material_1 = input ('Enter tip coating material code (except for iron carbonyl):  ');
      
      
		if material_0 == 7,
	%material1 = input ('Enter the sample material code (8 for default, 0 for iron carbonyl):  ');
			if material1 == 0,
				ic_shapka
			else
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   REFparticle = inf;
   %vFermi = input ('Enter the sample vFermi (e.g., 1.4E8 for noble metals):  ');
   vFermi = 1.4E8;
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
   resonance_type = 2;
   %omega_r = 6E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   omega_r = 6E10 * pi * 1815;
   %REFgamma = 6E10 * pi * input ('Enter the reference gamma, cm-1:  ');
   REFgamma = 6E10 * pi * 15;
   Mif = 1E30 * input ('Enter absorption strength (1 to 1000):  ');
   %sigma_eps =  1E15 * input ('Enter sigma, typically much less than 1:  ');
   sigma_eps = 0;
   trial_number = input('Enter the model number:  ');
			end
		end
      
      
      
   end
end

v2 = input('Change sample material? 1=Yes 0=No:  ');
if v2 == 1,
	disp('...');disp('0 = default   1 = gold');disp('2 = platinum   3 = silver');disp('4 = silicon   5 = silicon fit');disp('6 = silicon nitride');disp('7 = reference material');
	disp('...');material2 = input ('Enter sample material code:  ');
	if material2 == 7,
      
      
      
		material2 = input ('Enter the sample material code (0 for default, 7 for iron carbonyl):  ');
		if material2 == 7,
         eps_inf= input ('Enter epsilon inf for iron carbonyl:  ');
			ic_shapka
		else
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   REFparticle = inf;
   %vFermi = input ('Enter the sample vFermi (e.g., 1.4E8 for noble metals):  ');
   vFermi = 1.4E8;
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
   resonance_type = 2;
   %omega_r = 6E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   omega_r = 6E10 * pi * 1815;
   %REFgamma = 6E10 * pi * input ('Enter the reference gamma, cm-1:  ');
   REFgamma = 6E10 * pi * 15;
   Mif = 1E30 * input ('Enter absorption strength (1 to 1000):  ');
   %sigma_eps =  1E15 * input ('Enter sigma, typically much less than 1:  ');
   sigma_eps = 0;
   trial_number = input('Enter the model number:  ');
		end

      
      
      %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
		%omega_r = 6.000E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   	%REFgamma = input ('Enter the reference gamma, cm-1:  ');
   	%REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
      %vFermi = input ('Enter the sample vFermi (1.4000E8 for noble metals):  ');
   end
end






ic_shapka.m

   eps_inf=1.8;%input ('Enter epsilon inf (recommended 1.8):  ');
%   bgnd_abs=input ('Enter background absorption (recommended 0.15 to 2):  ');%0.55;%
%   bgnd1_labs=0;%input ('Enter first local absorption strength (recommended 0.04 to 10):  ');%0.04;%10
%   bgnd2_labs=0;%.051;%input ('Enter second local absorption strength (recommended 0 to 1.5):  ');%1.5
   bgnd_ref=0.5;%input ('Enter background reflection (recommended 0.5, %):  ')/100;
   %eps_inf=2;
   sigma_eps=0;

   %omega...=omega/60*pi to convert into wavenumbers
	omegac=10.25/30;
   scale_coeff = 45 - eps_inf; %32; for eps_inf = 1 
   eps_ratio=4;
	gamma1=7;
	omega1=sqrt (omegac^2+omegac*gamma1*eps_ratio);
   
   Mif1=scale_coeff * ((omega1^2-omegac^2)^2+omegac^2*gamma1^2) / (omega1^2-omegac^2);
   
	%e_charge =1.6021773349E-19*2997924580; %unit charge, statC (cgs or Gaussian system)
	%m_au = 1.66958E-24;
	%m_C = 12.011*m_au;
	%m_O = 15.9994*m_au;
	%m_eff = m_C*m_O /(m_C+m_O);
	%m_Fe = 55.847*m_au;
	%n = 1E21/0.0673478925113312; %mol. concentration, 1/cm^3 @density ~9g/cm^3 + one electron per molecule
	%Mif = omegap = 4*pi*(e_charge)^2*(3*n)/(m_eff);

   Mif2 = 0.029;%input ('Enter the oscillator strength (recommended 0.02):  ');
   %iron carbonyl absorption peak: 1815 to 1830, cm-1
   omega2=1826;%input ('Enter iron carbonyl absorption peak: 1817 to 1826 (cm-1):  ');

    omega_r = 6.000E10 * pi * omega2;
   gamma2=17.5;%input ('Enter 1826 peak width: 15 to 60 (cm-1):  ');%35-45(7.5-15@100K)

%   Mif2 = Mif2 * gamma2/15;
   
if scale_coeff > 8 * (Mif2 - 1.5), %Mif2 < 5.5,
   eps_ratio=eps_ratio/scale_coeff;
   scale_coeff = scale_coeff - 8*Mif2;
   eps_ratio=eps_ratio*scale_coeff;
	omega1 = sqrt (omegac^2+omegac*gamma1*eps_ratio);
	Mif1 = scale_coeff * ((omega1^2-omegac^2)^2+omegac^2*gamma1^2) / (omega1^2-omegac^2);      
else
	Mif1 = omegac*gamma1*8;%scale_coeff/eps_ratio;      
   scale_coeff = 0;
   eps_ratio=0;
end

%effective IR resonances, cm-1
	omega0=599;%input('600:  ');%598;%650;
	gamma0=15;%input('600 width (15):  ');%13;%25;

	omega3=2012;%input('2012:  ');%2012;%(2011);%2020 @ 100K
	gamma3=11;%input('10:  ');%10;%input ('Enter 2018 peak width: 9 to 60 (cm-1):  ');%30(9-21@100K)

	omega4=2086;%2088 @ 100K
    gamma4=12.5;%150;%input ('Enter 2085 peak width: 6 to 40 (cm-1):  ');%6-20
   
	omega5=2018;%input('2018:  ');%2018;%input('adjust frequency (2450):  ');%2350;%2000;%2500;
    gamma5=14;            %input('10:  ');%10;%input('band width (400):  ');%400;%10000;%1100;

	omega6=1E4;
    gamma6=1E3;

	omega7=454;%1730;%input ('Enter first absorption peak : 1650 to 1750 (cm-1):  ');%1650;%1700;
    gamma7=16;%100;%input ('Enter first absorption peak width: 200 to 500 (cm-1):  ');%200;%500;
   
	omega8=675;%2030;%input ('Enter second absorption peak : 2052 to 2088 (cm-1):  ');%2054;
    gamma8=38;%30;%input ('Enter second absorption peak width: 200 to 500 (cm-1):  ');%200;%500;
   
if eps_inf==1,
   ratio2=0;
   ratio3=0;
   ratio4=0;
   ratio5=0;
   ratio0=0;%input('ratio0 (600):  ');
   ratio7=0;%input('ratio7 (454):  ');
   ratio8=0;%input('ratio8 (675):  ');
else
   if eps_inf==1.2,
   	ratio2=0.0037;
   	ratio3=0.0028;
   	ratio4=0.0007;
   	ratio5=0.01;
   	ratio0=0.011;
   	ratio7=0.002;
   	ratio8=0.003;
   else
   	if eps_inf==1.5,
   		ratio2=0.013;
   		ratio3=0.0053;
   		ratio4=0.005;
   		ratio5=0.025;
   ratio0=0.041;%input('ratio0 (600):  ');
   ratio7=0.0075;%input('ratio7 (454):  ');
   ratio8=0.011;%input('ratio8 (675):  ');
    else
   		if eps_inf==1.8,
   			ratio2=0.017;
   			ratio3=0.006;
   			ratio4=0.0021;
   			ratio5=0.05;
   			ratio0=0.042;
   			ratio7=0.009;
   			ratio8=0.01;
      	else
   		if eps_inf==1.9,
   			ratio2=0.0185;
   			ratio3=0.006;
   			ratio4=0.0022;
   			ratio5=0.0606;
   			ratio0=0.0722;
   			ratio7=0.014;
   			ratio8=0.0189;
      	else
   			if eps_inf==2,
   				ratio2=0.02;
   				ratio3=0.006;
   				ratio4=0.0024;
   				ratio5=0.07;
   ratio0=0.11;%input('ratio0 (600):  ');
   ratio7=0.02;%input('ratio7 (454):  ');
   ratio8=0.03;%input('ratio8 (675):  ');
      		else
    			if eps_inf==2.1,
   				ratio2=0.0216;
   				ratio3=0.00603;
   				ratio4=0.00265;
   				ratio5=0.0754;
   ratio0=0.121;%input('ratio0 (600):  ');
   ratio7=0.022;%input('ratio7 (454):  ');
   ratio8=0.033;%input('ratio8 (675):  ');
                else
                    if eps_inf==2.2,
   				ratio2=0.0233;
   				ratio3=0.0061;
   				ratio4=0.0029;
   				ratio5=0.08;
   ratio0=0.132;%input('ratio0 (600):  ');
   ratio7=0.024;%input('ratio7 (454):  ');
   ratio8=0.036;%input('ratio8 (675):  ');
                    else
                        if eps_inf==2.4,
   				ratio2=0.0269;
   				ratio3=0.00633;
   				ratio4=0.0033;
   				ratio5=0.09;
   ratio0=0.154;%input('ratio0 (600):  ');
   ratio7=0.028;%input('ratio7 (454):  ');
   ratio8=0.042;%input('ratio8 (675):  ');
                        else
                            if eps_inf==3,
   					ratio2=0.039;
   					ratio3=0.008;
   					ratio4=0.0042;
                    ratio5=0.12;
   ratio0=0.22;%input('ratio0 (600):  ');
   ratio7=0.04;%input('ratio7 (454):  ');
   ratio8=0.06;%input('ratio8 (675):  ');
                            else
   ratio2=input('ratio2 (1826):  ');
   ratio3=input('ratio3 (2012):  ');
   ratio4=input('ratio4 (2086):  ');
   ratio5=input('ratio5 (2018):  ');
   ratio0=input('ratio0 (600):  ');
   ratio7=input('ratio7 (454):  ');
   ratio8=input('ratio8 (675):  ');
                            end
                        end
                    end
                end
            end
        end
        end
    end
   end
end
   Mif2 = (ratio2+Mif2)*omega2^2;%((omega2^2-omegac^2)^2+omegac^2*gamma2^2) / (omega2^2-omegac^2);
   Mif0 = (ratio0+0.075)*omega0^2;%Mif2 * ((134+120+798+763)/(682+899))/2;
%   	dolya=input('ratio (1):  ');%0.6-0.8*gamma4/gamma3;

Mif3 = 1.2*(ratio3+0.0195)*omega3^2;%*Mif2 * ((1322+1527+1747+2210)/(682+899)+3)/2;%input('main ratio (1):  ')*...%0.36-0.5

   Mif4 = (ratio4+0.003)*omega4^2;%dolya*Mif3;

Mif5 = 2.4*0.175*(ratio5+0.0157)*omega5^2;%bgnd_abs*
   
   Mif6 = (eps_inf-1)*omega6^2;
   Mif7 = (ratio7+0.017)*omega7^2;%bgnd1_labs*omega7^2;
   Mif8 = (ratio8+0.015)*omega8^2;%bgnd2_labs*omega8^2;






eps_ic.m

%epsilonre2=Mif2*(omega2^2-nu^2)/((omega2^2-nu^2)^2+nu^2*gamma2^2);
%epsilonim2=Mif2*nu*gamma2/((omega2^2-nu^2)^2+nu^2*gamma2^2);


epsilonre0=Mif0*(omega0^2-nu^2)/((omega0^2-nu^2)^2+nu^2*gamma0^2);
epsilonre1=Mif1*(omega1^2-nu^2)/((omega1^2-nu^2)^2+nu^2*gamma1^2);
epsilonre2=Mif2*(omega2^2-nu^2)/((omega2^2-nu^2)^2+nu^2*gamma2^2);
epsilonre3=Mif3*(omega3^2-nu^2)/((omega3^2-nu^2)^2+nu^2*gamma3^2);
epsilonre4=Mif4*(omega4^2-nu^2)/((omega4^2-nu^2)^2+nu^2*gamma4^2);
epsilonre5=Mif5*(omega5^2-nu^2)/((omega5^2-nu^2)^2+nu^2*gamma5^2);%complementary
epsilonre6=Mif6*(omega6^2-nu^2)/((omega6^2-nu^2)^2+nu^2*gamma6^2);
epsilonre7=Mif7*(omega7^2-nu^2)/((omega7^2-nu^2)^2+nu^2*gamma7^2);%complementary
epsilonre8=Mif8*(omega8^2-nu^2)/((omega8^2-nu^2)^2+nu^2*gamma8^2);%complementary
epsilonre=eps_inf+epsilonre0+epsilonre1+epsilonre2+epsilonre3+epsilonre4+epsilonre5+epsilonre6+epsilonre7+epsilonre8;
%epsL = eps_inf + Mif0/omega0^2 + Mif1/omega1^2 + Mif2/omega2^2 + Mif3/omega3^2;

epsilonim0=Mif0*nu*gamma0/((omega0^2-nu^2)^2+nu^2*gamma0^2);
epsilonim1=Mif1*nu*gamma1/((omega1^2-nu^2)^2+nu^2*gamma1^2);
epsilonim2=Mif2*nu*gamma2/((omega2^2-nu^2)^2+nu^2*gamma2^2);
epsilonim3=Mif3*nu*gamma3/((omega3^2-nu^2)^2+nu^2*gamma3^2);
epsilonim4=Mif4*nu*gamma4/((omega4^2-nu^2)^2+nu^2*gamma4^2);
epsilonim5=Mif5*nu*gamma5/((omega5^2-nu^2)^2+nu^2*gamma5^2);
epsilonim6=Mif6*nu*gamma6/((omega6^2-nu^2)^2+nu^2*gamma6^2);%complementary
epsilonim7=Mif7*nu*gamma7/((omega7^2-nu^2)^2+nu^2*gamma7^2);
epsilonim8=Mif8*nu*gamma8/((omega8^2-nu^2)^2+nu^2*gamma8^2);
%epsilonim = 4*pi*sigma_eps + epsilonim0 + epsilonim1 + epsilonim2 + epsilonim3;
epsilonim=epsilonim0+epsilonim1+epsilonim2+epsilonim3+epsilonim4+epsilonim5+epsilonim6+epsilonim7+epsilonim8;






funkcia.m

omega = 6E10 * pi * nu;
if material2 == 7,
   %get iron carbonyl epsilon
   eps_ic
else
   %get model epsilon
   %abs2eps
   
   sepsilon
end

if poglot ~= 0,
	es = epsilonre + i * epsilonim;
	beta = (es - 1) / (es + 1);
	%alpha eff for tip over the flat sample surface
	if material1 == 7,
		et	%coated tip
	else
   	tipepsilon	%etip=effective tip epsilon (omega)
		alpha = 4 * pi * Rparticle1 ^ 3 * (etip - 1) / (etip + 2);
   end
	alphae = alpha * (1 + beta) / (1 - alpha * beta / (16 * pi * (Rparticle1 + d) ^ 3));

	k = 2 * pi * nu;
   if poglot == 1,
      response = 4 * pi * k * imag(alphae);
   end
   if poglot == 2,
      response = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
   end
end






abswn.m

%absorption for fixed distance from the surface vs. wavenumber (800-1950 cm-1)
clear global
nf_mode = 0;

material1 = input ('Enter the tip material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for coated tip):  ');
if material1 == 7,
	material_0 = input ('Enter tip core material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for coated tip):  ');
	material_1 = input ('Enter tip coating material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride):  ');
	ad = input ('Enter tip core radius (10 to 100), nm:  ');
	a = ad * 1.000E-7 + input ('Enter the tip coating thickness (25), nm:  ') * 1.000E-7;
else
   %type in the end of tip radius of curvature (cm)
	Rparticle1 = input ('Enter the end-of-tip radius of curvature (10 to 120), nm: ') * 1.000E-7;
end

material2 = input ('Enter sample material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for reference):  ');
if material2 == 7,
   resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
	omega_r = 6.000E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   REFgamma = input ('Enter the reference gamma, cm-1:  ');
   REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   vFermi = input ('Enter the sample vFermi (1.4000E8 for noble metals):  ');
end

%type in the particle radius (cm)
Rparticle2 = input ('Enter particle radius (1 to inf), nm:  ');

%type in the minimum tip-surface distance (cm)
dmin = 1.000E-7 * input ('Enter tip-sample distance, nm:   ');

%type in the amplitude of cantilever oscillation (cm)
ao = 1.000E-7 * input ('Enter amplitude of oscillation, nm:   ');

stpn=0;

figure		%plot Cabs vs. nu
hold on

nomer = 1;
plot_abswn

done = input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   nf
	nomer = 2;
	plot_abswn
end

done = input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   nf
	nomer = 3;
	plot_abswn
end
hold off






csca.m

%scattering for single wavelength and distance from the surface

%parameters: k, alphae (omega)
%alphae=alpha eff for single wavelength

%type in alfa (for homogeneous tip) or et (for coated tip)
alfa

%scattering cross-section
C_sca = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;






eps_low.m

ic_shapka

x1=8.5/30;	y1=48.3; z1=0.86;
x2=9/30;		y2=53.4; z2=1.55;
x3=9.5/30;	y3=42.7; z3=0.45;
x4=10/30;	y4=41.1; z4=1.95;
x5=10.5/30;	y5=49.8; z5=2.19;
x6=11/30;	y6=42.6; z6=2.71;
x7=11.5/30;	y7=42.8; z7=1.72;
x8=12/30;	y8=44.3; z8=1.02;

figure
hold on

subplot (2, 2, 1)
plot (x1, y1, 'k*', x2, y2, 'k*', x3, y3, 'k*', x4, y4, 'k*', x5, y5, 'k*', x6, y6, 'k*', x7, y7, 'k*', x8, y8, 'k*')
hold on

y1=52.7;
y2=55.4;
y3=43.5;
y4=42.7;
y5=47.1;
y6=45.0;
y7=47.5;
y8=43.7;

plot (x1, y1, 'k*', x2, y2, 'k*', x3, y3, 'k*', x4, y4, 'k*', x5, y5, 'k*', x6, y6, 'k*', x7, y7, 'k*', x8, y8, 'k*')
hold on

E = zeros (1, 150);
E1 = zeros (1, 150);
nu_m = zeros (1, 150);
for j=1:1500,	%nu=.2:.002:.5,
   nu=.2+.002*j;
   nu_m(1,j)=nu;
   eps_ic
   E(1,j)=epsilonre;
	E1(1,j)=epsilonre1;
%plot (nu, epsilonre1, 'k', 'markersize', 5);
%plot (nu, epsilonre, 'k', 'markersize', 5);
end
plot (nu_m, E1, 'k--', 'markersize', 5);
plot (nu_m, E, 'k-.', 'markersize', 5);
axis ([.2 .5 0 100])
title ('epsilon Re')
clear E

subplot (2, 2, 2)
plot (x1, z1, 'k*', x2, z2, 'k*', x3, z3, 'k*', x4, z4, 'k*', x5, z5, 'k*', x6, z6, 'k*', x7, z7, 'k*', x8, z8, 'k*')
hold on

z1=0.93;
z2=1.73;
z3=0.17;
z4=1.51;
z5=1.98;
z6=3.06;
z7=1.54;
z8=1.16;

plot (x1, z1, 'k*', x2, z2, 'k*', x3, z3, 'k*', x4, z4, 'k*', x5, z5, 'k*', x6, z6, 'k*', x7, z7, 'k*', x8, z8, 'k*')
hold on

for j=1:1500,	%nu=.2:.002:.5,
   nu=.2+.002*j;
   nu_m(1,j)=nu;
   eps_ic
%plot (nu, epsilonim1, 'k', 'markersize', 5)
	E1(1,j)=epsilonim1;
end
plot (nu_m, E1, 'k-.', 'markersize', 5)
axis ([.2 .5 0 100])
title ('epsilon Im')
clear E E1 nu_m

x1=8.5/30;	y1=48.3; z1=0.86;
x2=9/30;		y2=53.4; z2=1.55;
x3=9.5/30;	y3=42.7; z3=0.45;
x4=10/30;	y4=41.1; z4=1.95;
x5=10.5/30;	y5=49.8; z5=2.19;
x6=11/30;	y6=42.6; z6=2.71;
x7=11.5/30;	y7=42.8; z7=1.72;
x8=12/30;	y8=44.3; z8=1.02;

subplot (2, 2, 3)
semilogx (x1, y1, 'k*', x2, y2, 'k*', x3, y3, 'k*', x4, y4, 'k*', x5, y5, 'k*', x6, y6, 'k*', x7, y7, 'k*', x8, y8, 'k*')
hold on

y1=52.7;
y2=55.4;
y3=43.5;
y4=42.7;
y5=47.1;
y6=45.0;
y7=47.5;
y8=43.7;

semilogx (x1, y1, 'k*', x2, y2, 'k*', x3, y3, 'k*', x4, y4, 'k*', x5, y5, 'k*', x6, y6, 'k*', x7, y7, 'k*', x8, y8, 'k*')
hold on

for schet=1:10000,
	nu=0.2*10^(schet/3000);
   eps_ic
	semilogx (nu, epsilonre, 'k')
end
for schet=10001:60000,
	nu=200*10^(1/3+(schet-10000)/50000);
   eps_ic
	semilogx (nu, epsilonre, 'k')
end
for schet=60001:62500,
	nu=2000*10^(1/3+(schet-60000)/5000);
   eps_ic
	semilogx (nu, epsilonre, 'k')
end
%axis ([.1 10000 30 50])
title ('overall Re')

%subplot (2, 2, 4)
%semilogx (x1, z1, 'k*', x2, z2, 'k*', x3, z3, 'k*', x4, z4, 'k*', x5, z5, 'k*', x6, z6, 'k*', x7, z7, 'k*', x8, z8, 'k*')
%hold on

z1=0.93;
z2=1.73;
z3=0.17;
z4=1.51;
z5=1.98;
z6=3.06;
z7=1.54;
z8=1.16;

%semilogx (x1, z1, 'k*', x2, z2, 'k*', x3, z3, 'k*', x4, z4, 'k*', x5, z5, 'k*', x6, z6, 'k*', x7, z7, 'k*', x8, z8, 'k*')

subplot (2, 2, 4)
%for schet=1:10000,
%nu=0.2*10^(schet/3000);
%epsilonim0=Mif0*nu*gamma0/((omega0^2-nu^2)^2+nu^2*gamma0^2);
%epsilonim1=Mif1*nu*gamma1/((omega1^2-nu^2)^2+nu^2*gamma1^2);
%epsilonim2=Mif2*nu*gamma2/((omega2^2-nu^2)^2+nu^2*gamma2^2);
%epsilonim3=Mif3*nu*gamma3/((omega3^2-nu^2)^2+nu^2*gamma3^2);
%epsilonim = 4*pi*sigma_eps + epsilonim0 + epsilonim1 + epsilonim2 + epsilonim3;
%semilogx (nu, epsilonim, 'k')
%hold on
%axis ([.1 10000 0 50])
%end
for schet=1:10000,
	nu=200*10^(schet/20000);
   eps_ic
	semilogx (nu, epsilonim, 'k')
hold on
end
for schet=10001:90000,
	nu=200*10^(0.5+(schet-10000)/100000);
   eps_ic
	semilogx (nu, epsilonim, 'k')
end
%axis ([.1 10000 0 50])
title ('overall Im in IR')

hold off






epslow.m

sigma_eps=0;

%omega=omega/60*pi
omegac=10.25/30;
%omegac=0.02;
scale_coeff=32;
eps_ratio=4;
%eps_ratio=70;
gamma1=7; %+ 1.39E8 / ad;
omega1 = sqrt (omegac^2+omegac*gamma1*eps_ratio);

Mif1 = scale_coeff * ((omega1^2-omegac^2)^2+omegac^2*gamma1^2) / (omega1^2-omegac^2);

%e_charge =1.6021773349E-19*2997924580; %unit charge, statC (cgs or Gaussian system)
%m_au = 1.66958E-24;
%m_C = 12.011*m_au;
%m_O = 15.9994*m_au;
%m_eff = m_C*m_O /(m_C+m_O);
%m_Fe = 55.847*m_au;
%n = 1E21/0.0673478925113312; %density ~9g/cm^3
%Mif = omegap = 4*pi*(e_charge)^2*n/(m_eff);
Mif2 = input ('Enter the IR resonance factor:  ');

	if Mif2 < 4.5,
      %scale_coeff = 32 + 8*(1.5-Mif2);
      %eps_ratio=scale_coeff*7/30;
      eps_ratio=eps_ratio/scale_coeff;
      scale_coeff = scale_coeff + 8*(1.5-Mif2);
      eps_ratio=eps_ratio*scale_coeff;
		omega1 = sqrt (omegac^2+omegac*gamma1*eps_ratio)
		Mif1 = scale_coeff * ((omega1^2-omegac^2)^2+omegac^2*gamma1^2) / (omega1^2-omegac^2);
   else
      scale_coeff = 0;
      eps_ratio=0;
		%Mif1 = 2*omegac*gamma1*30/7;
		Mif1 = 2*omegac*gamma1*scale_coeff/eps_ratio;
   end
   
%effective IR resonances, cm-1
omega0=650;
gamma0=25; %+ 1.39E8 / ad;

omega2=1815;
gamma2=15; %+ 1.39E8 / ad;

omega3=2016;
gamma3=9; %+ 1.39E8 / ad;

%e_charge =1.6021773349E-19*2997924580; %unit charge, statC (cgs or Gaussian system)
%m_au = 1.66958E-24;
%m_C = 12.011*m_au;
%m_O = 15.9994*m_au;
%m_eff = m_C*m_O/(m_C+m_O);
%m_Fe = 55.847*m_au;
%n = 1E21/0.0673478925113312; %density ~9g/cm^3
%Mif = omegap = 4*pi*(e_charge)^2*n/(m_eff);
Mif2 = Mif2 * ((omega2^2-omegac^2)^2+omegac^2*gamma2^2) / (omega2^2-omegac^2);
Mif0 = Mif2 * ((134+120+798+763)/(682+899))/2;
Mif3 = Mif2 * ((1322+1527+1747+2210)/(682+899)+3)/2;

%epsilonre2=Mif2*(omega2^2-omega^2)/((omega2^2-omega^2)^2+omega^2*gamma2^2);
%epsilonim2=Mif2*omega*gamma2/((omega2^2-omega^2)^2+omega^2*gamma2^2);

x1=8.5/30;	y1=48.3; z1=0.86;
x2=9/30;		y2=53.4; z2=1.55;
x3=9.5/30;	y3=42.7; z3=0.45;
x4=10/30;	y4=41.1; z4=1.95;
x5=10.5/30;	y5=49.8; z5=2.19;
x6=11/30;	y6=42.6; z6=2.71;
x7=11.5/30;	y7=42.8; z7=1.72;
x8=12/30;	y8=44.3; z8=1.02;

figure
hold on

subplot (2, 2, 1)
plot (x1, y1, 'k*', x2, y2, 'k*', x3, y3, 'k*', x4, y4, 'k*', x5, y5, 'k*', x6, y6, 'k*', x7, y7, 'k*', x8, y8, 'k*')
hold on

y1=52.7;
y2=55.4;
y3=43.5;
y4=42.7;
y5=47.1;
y6=45.0;
y7=47.5;
y8=43.7;

plot (x1, y1, 'k*', x2, y2, 'k*', x3, y3, 'k*', x4, y4, 'k*', x5, y5, 'k*', x6, y6, 'k*', x7, y7, 'k*', x8, y8, 'k*')
hold on

E = zeros (1, 150);
E1 = zeros (1, 150);
omega_m = zeros (1, 150);
for j=1:1500,	%omega=.2:.002:.5,
   omega=.2+.002*j;
   omega_m(:,j)=omega;
epsilonre0=Mif0*(omega0^2-omega^2)/((omega0^2-omega^2)^2+omega^2*gamma0^2);
epsilonre1= 1 + Mif1*(omega1^2-omega^2)/((omega1^2-omega^2)^2+omega^2*gamma1^2);
epsilonre2=Mif2*(omega2^2-omega^2)/((omega2^2-omega^2)^2+omega^2*gamma2^2);
epsilonre3=Mif3*(omega3^2-omega^2)/((omega3^2-omega^2)^2+omega^2*gamma3^2);
epsilonre = epsilonre0 + epsilonre1 + epsilonre2 + epsilonre3;
	E(:,j)=epsilonre;
	E1(:,j)=epsilonre1;
%plot (omega, epsilonre1, 'k', 'markersize', 5);
%plot (omega, epsilonre, 'k', 'markersize', 5);
end
plot (omega_m, E1, 'k--', 'markersize', 5);
plot (omega_m, E, 'k-.', 'markersize', 5);
axis ([.2 .5 0 100])
title ('epsilon Re')
clear E

subplot (2, 2, 2)
plot (x1, z1, 'k*', x2, z2, 'k*', x3, z3, 'k*', x4, z4, 'k*', x5, z5, 'k*', x6, z6, 'k*', x7, z7, 'k*', x8, z8, 'k*')
hold on

z1=0.93;
z2=1.73;
z3=0.17;
z4=1.51;
z5=1.98;
z6=3.06;
z7=1.54;
z8=1.16;

plot (x1, z1, 'k*', x2, z2, 'k*', x3, z3, 'k*', x4, z4, 'k*', x5, z5, 'k*', x6, z6, 'k*', x7, z7, 'k*', x8, z8, 'k*')
hold on

for j=1:1500,	%omega=.2:.002:.5,
   omega=.2+.002*j;
   omega_m(:,j)=omega;
epsilonim1= 4*pi*sigma_eps + Mif1*omega*gamma1/((omega1^2-omega^2)^2+omega^2*gamma1^2);
%plot (omega, epsilonim1, 'k', 'markersize', 5)
	E1(:,j)=epsilonim1;
end
plot (omega_m, E1, 'k-.', 'markersize', 5)
axis ([.2 .5 0 100])
title ('epsilon Im')
clear E E1 omega_m

x1=8.5/30;	y1=48.3; z1=0.86;
x2=9/30;		y2=53.4; z2=1.55;
x3=9.5/30;	y3=42.7; z3=0.45;
x4=10/30;	y4=41.1; z4=1.95;
x5=10.5/30;	y5=49.8; z5=2.19;
x6=11/30;	y6=42.6; z6=2.71;
x7=11.5/30;	y7=42.8; z7=1.72;
x8=12/30;	y8=44.3; z8=1.02;

subplot (2, 2, 3)
semilogx (x1, y1, 'k*', x2, y2, 'k*', x3, y3, 'k*', x4, y4, 'k*', x5, y5, 'k*', x6, y6, 'k*', x7, y7, 'k*', x8, y8, 'k*')
hold on

y1=52.7;
y2=55.4;
y3=43.5;
y4=42.7;
y5=47.1;
y6=45.0;
y7=47.5;
y8=43.7;

semilogx (x1, y1, 'k*', x2, y2, 'k*', x3, y3, 'k*', x4, y4, 'k*', x5, y5, 'k*', x6, y6, 'k*', x7, y7, 'k*', x8, y8, 'k*')
hold on

for schet=1:10000,
omega=0.2*10^(schet/3000);
epsilonre0=Mif0*(omega0^2-omega^2)/((omega0^2-omega^2)^2+omega^2*gamma0^2);
epsilonre1= 1 + Mif1*(omega1^2-omega^2)/((omega1^2-omega^2)^2+omega^2*gamma1^2);
epsilonre2=Mif2*(omega2^2-omega^2)/((omega2^2-omega^2)^2+omega^2*gamma2^2);
epsilonre3=Mif3*(omega3^2-omega^2)/((omega3^2-omega^2)^2+omega^2*gamma3^2);
epsilonre = epsilonre0 + epsilonre1 + epsilonre2 + epsilonre3;
semilogx (omega, epsilonre, 'k')
end
for schet=10001:60000,
omega=200*10^(1/3+(schet-10000)/50000);
epsilonre0=Mif0*(omega0^2-omega^2)/((omega0^2-omega^2)^2+omega^2*gamma0^2);
epsilonre1= 1 + Mif1*(omega1^2-omega^2)/((omega1^2-omega^2)^2+omega^2*gamma1^2);
epsilonre2=Mif2*(omega2^2-omega^2)/((omega2^2-omega^2)^2+omega^2*gamma2^2);
epsilonre3=Mif3*(omega3^2-omega^2)/((omega3^2-omega^2)^2+omega^2*gamma3^2);
epsilonre = epsilonre0 + epsilonre1 + epsilonre2 + epsilonre3;
semilogx (omega, epsilonre, 'k')
end
for schet=60001:62500,
omega=2000*10^(1/3+(schet-60000)/5000);
epsilonre0=Mif0*(omega0^2-omega^2)/((omega0^2-omega^2)^2+omega^2*gamma0^2);
epsilonre1= 1 + Mif1*(omega1^2-omega^2)/((omega1^2-omega^2)^2+omega^2*gamma1^2);
epsilonre2=Mif2*(omega2^2-omega^2)/((omega2^2-omega^2)^2+omega^2*gamma2^2);
epsilonre3=Mif3*(omega3^2-omega^2)/((omega3^2-omega^2)^2+omega^2*gamma3^2);
epsilonre = epsilonre0 + epsilonre1 + epsilonre2 + epsilonre3;
semilogx (omega, epsilonre, 'k')
end
%axis ([.1 10000 30 50])
title ('overall Re')

%subplot (2, 2, 4)
%semilogx (x1, z1, 'k*', x2, z2, 'k*', x3, z3, 'k*', x4, z4, 'k*', x5, z5, 'k*', x6, z6, 'k*', x7, z7, 'k*', x8, z8, 'k*')
%hold on

z1=0.93;
z2=1.73;
z3=0.17;
z4=1.51;
z5=1.98;
z6=3.06;
z7=1.54;
z8=1.16;

%semilogx (x1, z1, 'k*', x2, z2, 'k*', x3, z3, 'k*', x4, z4, 'k*', x5, z5, 'k*', x6, z6, 'k*', x7, z7, 'k*', x8, z8, 'k*')

subplot (2, 2, 4)
%for schet=1:10000,
%omega=0.2*10^(schet/3000);
%epsilonim0=Mif0*omega*gamma0/((omega0^2-omega^2)^2+omega^2*gamma0^2);
%epsilonim1=Mif1*omega*gamma1/((omega1^2-omega^2)^2+omega^2*gamma1^2);
%epsilonim2=Mif2*omega*gamma2/((omega2^2-omega^2)^2+omega^2*gamma2^2);
%epsilonim3=Mif3*omega*gamma3/((omega3^2-omega^2)^2+omega^2*gamma3^2);
%epsilonim = 4*pi*sigma_eps + epsilonim0 + epsilonim1 + epsilonim2 + epsilonim3;
%semilogx (omega, epsilonim, 'k')
%hold on
%axis ([.1 10000 0 50])
%end
for schet=1:10000,
omega=200*10^(schet/20000);
epsilonim0=Mif0*omega*gamma0/((omega0^2-omega^2)^2+omega^2*gamma0^2);
epsilonim1= 4*pi*sigma_eps + Mif1*omega*gamma1/((omega1^2-omega^2)^2+omega^2*gamma1^2);
epsilonim2=Mif2*omega*gamma2/((omega2^2-omega^2)^2+omega^2*gamma2^2);
epsilonim3=Mif3*omega*gamma3/((omega3^2-omega^2)^2+omega^2*gamma3^2);
epsilonim = epsilonim0 + epsilonim1 + epsilonim2 + epsilonim3;
semilogx (omega, epsilonim, 'k')
hold on
end
for schet=10001:90000,
omega=200*10^(0.5+(schet-10000)/100000);
epsilonim0=Mif0*omega*gamma0/((omega0^2-omega^2)^2+omega^2*gamma0^2);
epsilonim1= 4*pi*sigma_eps + Mif1*omega*gamma1/((omega1^2-omega^2)^2+omega^2*gamma1^2);
epsilonim2=Mif2*omega*gamma2/((omega2^2-omega^2)^2+omega^2*gamma2^2);
epsilonim3=Mif3*omega*gamma3/((omega3^2-omega^2)^2+omega^2*gamma3^2);
epsilonim = epsilonim0 + epsilonim1 + epsilonim2 + epsilonim3;
semilogx (omega, epsilonim, 'k')
end
%axis ([.1 10000 0 50])
title ('overall Im in IR')

hold off






et.m

%alpha eff for single wavelength, for coated tip

%parameters:
%a=end of tip radius; ad=tip core radius; d=tip to surface distance
%ed=tip core epsilon; em=tip coating epsilon; f=relative core volume

f = ad ^ 3 / a ^ 3;

Rparticle = a;
%type in the tip coating material under this line, e.g. gold, platinum, silver, silicon
if material_1 == 0,
default
end

if material_1 == 1,
gold
end

if material_1 == 2,
platinum
end

if material_1 == 3,
silver
end

if material_1 == 4,
silicon
end

if material_1 == 5,
%silicon fit
siepsilon
end

if material_1 == 6,
%silicon nitride
sinit
end

em = epsilonre + i * epsilonim;
clear Rparticle

Rparticle = ad;
%type in the tip core material under this line, e.g. gold, platinum, silver, silicon, sin
if material_0 == 0,
default
end

if material_0 == 1,
gold
end

if material_0 == 2,
platinum
end

if material_0 == 3,
silver
end

if material_0 == 4,
silicon
end

if material_0 == 5,
%silicon fit
siepsilon
end

if material_0 == 6,
%silicon nitride
sinit
end

if material_0 == 7,
%iron carbonyl
eps_ic
end

ed = epsilonre + i * epsilonim;

alpha = 4 * pi * a ^ 3 * ((em - 1) * (ed + 2 * em) + f * (ed - em) * (1 + 2 * em)) / ((em + 2) * (ed + 2 * em) + 2 * f * (ed - em) * (em - 1));

clear epsilonre epsilonim ed em Rparticle

%es=sample surface epsilon
%if material2 ~= 7,
	sepsilon
	beta = (es - 1) / (es + 1);

	%alpha eff for single wavelength, for coated tip
	alphae = alpha * (1 + beta) / (1 - alpha * beta / (16 * pi * (a + d) ^ 3));
   clear Rparticle
%end






gold.m

%gold IR parameters

if omega < 3E10 * 2 * pi * 750,
	omega = 3E10 * 2 * pi * 750;
	gold
end

if omega > 3E10 * 2 * pi * 2000,
	omega = 3E10 * 2 * pi * 2000;
	gold
end

N = 5.9600E+22;
%omegap2 = 4 * pi * N * (4.80321276758E-10) ^ 2 / (9.10938188E-28);
omegap2 = 1.3000E+32;
epsilonb = 0.25;
epsiloninf = epsilonb + 1;
%vFermi = 1.4000E8;
OMEGA2 = 5 * 1.4000E8 ^ 2 / (3 * Rparticle ^ 2);
%gammab = 9.000E13;
gamma = 9.000E13 + 1.4000E8 / Rparticle;

epsilonre = epsiloninf + omegap2 * (OMEGA2 - omega ^ 2) / ((OMEGA2 - omega ^ 2) ^ 2 + omega ^ 2 * gamma ^ 2);
epsilonim = omegap2 * omega * gamma / ((OMEGA2 - omega ^ 2) ^ 2 + omega ^ 2 * gamma ^ 2);






platinum.m

%platinum IR parameters
%pre-requisites: Rparticle, omega

if omega < 3E10 * 2 * pi * 750,
	omega = 3.000E10 * 2 * pi * 750;
	platinum
end

if omega > 3E10 * 2 * pi * 2000,
	omega = 3.000E10 * 2 * pi * 2000;
	platinum
end

N = 5.9600E+22;
%omegap2 = 4 * pi * N * (4.80321276758E-10) ^ 2 / (9.10938188E-28);
omegap2 = 6.600E+31;
%epsilonb = 0.25;
epsiloninf = 0.25 + 1;
%vFermi = 1.400E8;
OMEGA2 = 5 * 1.400E8 ^ 2 / (3 * Rparticle ^ 2);
%gammab = 1.0500E14;
gamma = 1.0500E14 + 1.400E8 / Rparticle;

epsilonre = epsiloninf + omegap2 * (OMEGA2 - omega ^ 2) / ((OMEGA2 - omega ^ 2) ^ 2 + omega ^ 2 * gamma ^ 2);
epsilonim = omegap2 * omega * gamma / ((OMEGA2 - omega ^ 2) ^ 2 + omega ^ 2 * gamma ^ 2);






plot_abs.m

%absorption cross-section vs. wavenumber
stpn=0;

if nomer == 1,
for j=1:1150;
	nu = 800 + stpn;
	k = 2 * pi * nu;
	omega = k * 3.000E10;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
   plot(nu, C_abs, 'k.', 'markersize', 4);
   clear C_abs alfa alphae
	stpn=stpn+1;
end
	xlabel('wavenumber, [cm-1]'); ylabel('Cabs, [a.u.]');
	axis square;
end

if nomer == 2,
for j=1:1150;
	nu = 800 + stpn;
	k = 2 * pi * nu;
	omega = k * 3.000E10;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
   plot(nu, C_abs, '.', 'markersize', 5);
   clear C_abs alfa alphae
	stpn=stpn+1;
end
end

if nomer == 3,
for j=1:1150;
	nu = 800 + stpn;
	k = 2 * pi * nu;
	omega = k * 3.000E10;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
   plot(nu, C_abs, 'm.', 'markersize', 9);
   clear C_abs alfa alphae
	stpn=stpn+1;
end
end






plot_absamp.m

%absorption vs. amplitude
stpn=0;

if nomer == 1,
for j=1:1000;
	ao = stpn * 1.000E-7; %amplitude of cantilever oscillation (cm)
	d = dmin;
   alfa
	Cabs_d = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + 2 * ao;
	alfa
	Cabs_a = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + ao;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
	clear alfa d
	Cabs_1 = 0.5 * (Cabs_d - C_abs);
	Cabs_2 = 0.25 * (Cabs_d + Cabs_a - 2 * C_abs);
	ao = ao / Rparticle1;
   plot(ao, Cabs_2, 'k.', 'markersize', 4);
   clear alfa alphae C_abs Cabs_2 Cabs_d Cabs_ao ao
	stpn=stpn+1;
end
   xlabel('amplitude/tip radius'); ylabel('Cabs**, [a.u.]');
   axis square;
end

if nomer == 2,
for j=1:1000;
	ao = stpn * 1.000E-7; %amplitude of cantilever oscillation (cm)
	d = dmin;
   alfa
	Cabs_d = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + 2 * ao;
	alfa
	Cabs_a = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + ao;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
	clear alfa d
	Cabs_1 = 0.5 * (Cabs_d - C_abs);
	Cabs_2 = 0.25 * (Cabs_d + Cabs_a - 2 * C_abs);
	ao = ao / Rparticle1;
   plot(ao, Cabs_2, '.', 'markersize', 5);
   clear alfa alphae C_abs Cabs_2 Cabs_d Cabs_ao ao
	stpn=stpn+1;
end
end

if nomer == 3,
for j=1:1000;
	ao = stpn * 1.000E-7; %amplitude of cantilever oscillation (cm)
	d = dmin;
   alfa
	Cabs_d = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + 2 * ao;
	alfa
	Cabs_a = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + ao;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
	clear alfa d
	Cabs_1 = 0.5 * (Cabs_d - C_abs);
	Cabs_2 = 0.25 * (Cabs_d + Cabs_a - 2 * C_abs);
	ao = ao / Rparticle1;
   plot(ao, Cabs_2, 'm.', 'markersize', 9);
   clear alfa alphae C_abs Cabs_2 Cabs_d Cabs_ao ao
	stpn=stpn+1;
end
end






plot_absd.m

%absorption vs. distance
stpn=0;

if nomer == 1,
for j=1:500;
	d = stpn * 1.000E-7;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
   plot(stpn, C_abs, 'k.', 'markersize', 4);
   clear C_abs alfa alphae
	stpn=stpn+1;
end
   xlabel('d, [nm]'); ylabel('Cabs, [a.u.]');
   axis square;
end

if nomer == 2,
for j=1:500;
	d = stpn * 1.000E-7;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
   plot(stpn, C_abs, '.', 'markersize', 5);
   clear C_abs alfa alphae
	stpn=stpn+1;
end
end

if nomer == 3,
for j=1:500;
	d = stpn * 1.000E-7;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
   plot(stpn, C_abs, 'm.', 'markersize', 9);
   clear C_abs alfa alphae
	stpn=stpn+1;
end
end






plot_abswn.m

%absorption vs. wavenumber
stpn=0;

if nomer == 1,
for j=1:1150;
	nu = 800 + stpn;
   k = 2 * pi * nu;
	omega = k * 3.000E10;
	d = dmin;
  	alfa
	Cabs_d = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + 2 * ao;
	alfa
	Cabs_ao = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + ao;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
	clear alfa d
	Cabs_1 = 0.5 * (Cabs_d - C_abs);
	Cabs_2 = 0.25 * (Cabs_d + Cabs_ao - 2 * C_abs);
   plot(nu, Cabs_2, 'k.', 'markersize', 4);
   clear alphae C_abs Cabs_2 Cabs_d Cabs_ao nu k
	stpn=stpn+1;
end
   xlabel('wavenumber, [cm-1]'); ylabel('Cabs**, [a.u.]');
   axis square;
end

if nomer == 2,
for j=1:1150;
	nu = 800 + stpn;
   k = 2 * pi * nu;
	omega = k * 3.000E10;
	d = dmin;
  	alfa
	Cabs_d = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + 2 * ao;
	alfa
	Cabs_ao = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + ao;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
	clear alfa d
	Cabs_1 = 0.5 * (Cabs_d - C_abs);
	Cabs_2 = 0.25 * (Cabs_d + Cabs_ao - 2 * C_abs);
   plot(nu, Cabs_2, '.', 'markersize', 5);
   clear alphae C_abs Cabs_2 Cabs_d Cabs_ao nu k
	stpn=stpn+1;
end
end

if nomer == 3,
for j=1:1150;
	nu = 800 + stpn;
   k = 2 * pi * nu;
	omega = k * 3.000E10;
	d = dmin;
  	alfa
	Cabs_d = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + 2 * ao;
	alfa
	Cabs_ao = 4 * pi * k * imag(alphae);
	clear alfa d
	d = dmin + ao;
	alfa
	C_abs = 4 * pi * k * imag(alphae);
	clear alfa d
	Cabs_1 = 0.5 * (Cabs_d - C_abs);
	Cabs_2 = 0.25 * (Cabs_d + Cabs_ao - 2 * C_abs);
   plot(nu, Cabs_2, 'm.', 'markersize', 9);
   clear alphae C_abs Cabs_2 Cabs_d Cabs_ao nu k
	stpn=stpn+1;
end
end






plot_alfapf..m

%scattering vs. wavenumber
nu=1800;

Rparticle=0.4*1E-7;
Rparticle1=35E-7;


dff=100;
dff=dff*1E-7;

material1=2;
material2=1; Rparticle2=0.4*1E-7;
		for j=1:101;
   		k = 2 * pi * nu;
			omega = k * 3.000E10;
            d = 0;
   		alpha_new
			Csca_d = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
			d = dff;
   		alpha_new
			Csca_ff = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
			Csca_e = 0.5 * (abs(Csca_ff) - abs(Csca_d));
			Csca_p = 0.5 * (abs(Csca_ff)^2 - abs(Csca_d)^2);
   		xj (j) = nu;
   		yj (j) = abs(Csca_e);
        yyj(j) = abs(Csca_p);
			nu = nu + 1;
        end
%			figure
			hold on
%            plot(xj,yj);

material2=7; Rparticle2=inf; ic_shapka; eps_ic;
nu=1800;

		for j=1:101;
   		k = 2 * pi * nu;
			omega = k * 3.000E10;
            d = 0;
   		alpha_new
			Csca_d = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
			d = dff;
   		alpha_new
			Csca_ff = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
			Csca_e = 0.5 * (abs(Csca_ff) - abs(Csca_d));
			Csca_p = 0.5 * (abs(Csca_ff)^2 - abs(Csca_d)^2);
   		xj (j) = nu;
   		zj (j) = abs(Csca_e);
        zzj(j) = abs(Csca_p);
			nu = nu + 1;
        end
%			plot(xj,zj);
			hold on
            
            sj=zj./yj;
%            figure; hold on; plot(xj,sj);
            
            yj=yyj; zj=zzj; ssj=zj./yj;
%            figure; hold on; plot(xj,yj); plot(xj,zj);
%            figure; hold on; plot(xj,sj);






plot_polar.m

%polarizability angle vs. wavenumber
func_nf
nu = 1700;
hold on
   ugol = zeros (1, 2000);
   nuj = zeros (1, 2000);
	for j=1:2000;
omega = 3E10 * 2 * pi * nu;
      alfa_polar
      ugol1(j)=angle(alphae);
      nuj(j) = nu;
   	nu = nu + 0.3;
    end
func_nf
nu = 1700;
hold on
   ugol = zeros (1, 2000);
   nuj = zeros (1, 2000);
	for j=1:2000;
omega = 3E10 * 2 * pi * nu;
      alfa_polar
      ugol2(j)=angle(alphae);
      nuj(j) = nu;
   	nu = nu + 0.3;
    end
plot (nuj, (ugol2-ugol1), 'k.', 'markersize', 5);
hold on
xlabel('wavenumber, [cm-1]'); ylabel('phase');
axis square;






plot_ref.m

%ic reflectance (albedo)
%nu = (100 * round (omega_r/(6.000E12 * pi)) - 100);

%eps_inf = input ('Enter epsilon inf for iron carbonyl (1 to 10):  ');
ic_shapka
   dirre = zeros (1, 150);
   nuj = zeros (1, 150);
for j=1:150;
   nuj (j) = j + 1750;
   nu = j + 1750;
   eps_ic
   %nre = real(sqrt(0.3*(epsilonre+i*epsilonim)+0.7*1.5));
   nre = real(sqrt((epsilonre+i*epsilonim)));
      %dirre(j) = -log10(1-(2-(((nre-1)/(nre+1))^2))*(((nre-1)/(nre+1))^2)-bgnd_ref);
      dirre(j) = -log10(1-(2-(((nre-1.4)/(nre+1.4))^2))*(((nre-1.4)/(nre+1.4))^2)-bgnd_ref);
end
plot (nuj, dirre, 'c.-', 'markersize', 5);
hold on
xlabel('wavenumber, [cm-1]'); ylabel('reflectance');
axis square;
%hold off






plot_sca.m

%absorption cross-section vs. wavenumber
stpn=0;

if nomer == 1,
for j=1:1150;
	nu = 800 + stpn;
	k = 2 * pi * nu;
	omega = k * 3.000E10;
	alfa
   C_sca = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
   plot(nu, C_sca, 'k.', 'markersize', 4);
   clear C_sca alfa alphae
	stpn=stpn+1;
end
	xlabel('wavenumber, [cm-1]'); ylabel('Csca, [a.u.]');
	axis square;
end

if nomer == 2,
for j=1:1150;
	nu = 800 + stpn;
	k = 2 * pi * nu;
	omega = k * 3.000E10;
	alfa
   C_sca = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
   plot(nu, C_sca, '.', 'markersize', 5);
   clear C_sca alfa alphae
	stpn=stpn+1;
end
end

if nomer == 3,
for j=1:1150;
	nu = 800 + stpn;
	k = 2 * pi * nu;
	omega = k * 3.000E10;
	alfa
   C_sca = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
   plot(nu, C_sca, 'm.', 'markersize', 9);
   clear C_sca alfa alphae
	stpn=stpn+1;
end
end






plot_scad.m

%scattering cross-section vs. distance
stpn=0;
dj = zeros (500);
C_sca = zeros(500);
if nomer == 1,
for j=1:500;
	d = stpn * 1.000E-7;
   alfa
   C_sca (j) = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
   dj (j) = stpn;
	stpn=stpn+1;
end
	plot (dj, C_sca, 'k-.', 'markersize', 4);
	set (gca, 'yscale', 'log');
	hold on
   xlabel ('d, [nm]'); ylabel('Csca, [a.u.]');
   axis square;
   clear alfa C_sca dj alphae
end

if nomer == 2,
for j=1:500;
	d = stpn * 1.000E-7;
   alfa
   C_sca (j) = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
   dj (j) = stpn;
	stpn=stpn+1;
end
   plot (dj, C_sca, '-', 'markersize', 5);
   clear alfa C_sca dj alphae
end

if nomer == 3,
for j=1:500;
	d = stpn * 1.000E-7;
   alfa
   C_sca (j) = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
   dj (j) = stpn;
	stpn=stpn+1;
end
   plot (dj, C_sca, 'm.', 'markersize', 9);
   clear alfa C_sca dj alphae
end






plot_scatamp

%scattering vs. amplitude
stpn=0;

if nomer == 1,
for j=1:1000;
	ao = stpn * 1.000E-7; %amplitude of cantilever oscillation (cm)
	d = dmin;
   alfa
	Csca_d = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	d = dmin + 2 * ao;
	alfa
	Csca_a = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	d = dmin + ao;
	alfa
	C_sca = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	Csca_1 = 0.5 * (Csca_d - C_sca);
	Csca_2 = 0.25 * (Csca_d + Csca_a - 2 * C_sca);
	ao = ao / Rparticle1;
   plot(ao, Csca_2, 'k.', 'markersize', 4);
   clear alfa alphae C_sca Csca_2 Csca_d Csca_ao ao
	stpn=stpn+1;
end
   xlabel('amplitude/tip radius'); ylabel('Csca**, [a.u.]');
   axis square;
end

if nomer == 2,
for j=1:1000;
	ao = stpn * 1.000E-7; %amplitude of cantilever oscillation (cm)
	d = dmin;
   alfa
	Csca_d = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	d = dmin + 2 * ao;
	alfa
	Csca_a = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	d = dmin + ao;
	alfa
	C_sca = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	Csca_1 = 0.5 * (Csca_d - C_sca);
	Csca_2 = 0.25 * (Csca_d + Csca_a - 2 * C_sca);
	ao = ao / Rparticle1;
   plot(ao, Csca_2, '.', 'markersize', 5);
   clear alfa alphae C_sca Csca_2 Csca_d Csca_ao ao
	stpn=stpn+1;
end
end

if nomer == 3,
for j=1:1000;
	ao = stpn * 1.000E-7; %amplitude of cantilever oscillation (cm)
	d = dmin;
   alfa
	Csca_d = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	d = dmin + 2 * ao;
	alfa
	Csca_a = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	d = dmin + ao;
	alfa
	C_sca = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	Csca_1 = 0.5 * (Csca_d - C_sca);
	Csca_2 = 0.25 * (Csca_d + Csca_a - 2 * C_sca);
	ao = ao / Rparticle1;
   plot(ao, Csca_2, 'm.', 'markersize', 9);
   clear alfa alphae C_sca Csca_2 Csca_d Csca_ao ao
	stpn=stpn+1;
end
end






plot_scatwn.m

%scattering vs. wavenumber
stpn=0;

if nomer == 1,
	chto = input ('Enter 0=load; 1=load, plot; 2=plot; 3=plot, save; 4=save; 5=re-calculate:  ');   
   if chto < 2,
      load fname1r_03_35
   else
   	xj = zeros(1, 2000);
	yj = ones(1, 2000);
		for j=1:2000;
			nu = 1700 + 0.1*stpn;
   		k = 2 * pi * nu;
			omega = k * 3.000E10;
			d = dmin;
   		alfa
			Csca_d = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
			clear alfa d
			d = dmin + 2 * ao;
			alfa
			Csca_ao = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
			clear alfa d
			d = dmin + ao;
			alfa
			C_sca = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
			clear alfa d
			Csca_1 = 0.5 * (sqrt(Csca_ao) - sqrt(Csca_d));
   		Csca_2 = 0.25 * (Csca_d + Csca_ao - 2 * C_sca);
   		xj (j) = nu;
   		yj (j) = abs(Csca_1);
			stpn=stpn+1;
        end
      if chto > 2,
         if chto < 5,
            save fname1r_03_35 xj yj
         end
      end
   end
   if chto > 0, 
      if chto < 4,
          %plot(xj,yj,'--rs','LineWidth',2,...
          %      'MarkerEdgeColor','k',...
          %      'MarkerFaceColor','g',...
          %      'MarkerSize',10)
			plot(xj, yj, 'k-', 'markersize', 5);
			hold on
      	clear alphae C_sca Csca_2 Csca_d Csca_ao nu k
      end
   end
	xlabel('wavenumber, [cm-1]'); ylabel('Csca**, [a.u.]');
   axis square;
end

if nomer == 2,
	ynewj = zeros(1, 2000);
	ynew = zeros(1, 2000);
for j=1:2000;
	nu = 1700 + 0.1*stpn;
   k = 2 * pi * nu;
	omega = k * 3.000E10;
	d = dmin;
   alfa
	Csca_d = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	d = dmin + 2 * ao;
	alfa
	Csca_ao = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	d = dmin + ao;
	alfa
	C_sca = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	Csca_1 = 0.5 * (sqrt(Csca_d) - sqrt(Csca_ao));
	Csca_2 = 0.25 * (Csca_d + Csca_ao - 2 * C_sca);
   xj (j) = nu;
   ynewj (j) = abs(Csca_1);
	stpn=stpn+1;
end
if chto > 3,
   ynew = ynewj./yj;
   plot(xj, ynew, '.-', 'markersize', 7);
   hold on
else
   if chto == 0,
   	ynew = ynewj./yj;
            save fname1r18_03_35 xj ynew
   	plot(xj, ynew, '.-', 'markersize', 7);
   	hold on
	else
   	plot(xj, ynewj, '.-', 'markersize', 7);
	end   
end
   clear alphae C_sca Csca_2 Csca_d Csca_ao nu k
end

if nomer == 3,
   xyj=0;
for j=1:200;
plot(1815, xyj); plot(1835, xyj); plot(1857,xyj); plot(1881,xyj);   
%	nu = 1700 + 0.1*stpn;
%	k = 2 * pi * nu;
%	omega = k * 3.000E10;
%	d = dmin;
%	alfa
%	Csca_d = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
%	clear alfa d
%	d = dmin + 2 * ao;
%	alfa
%	Csca_ao = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
%	clear alfa d
%	d = dmin + ao;
%	alfa
%	C_sca = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
%	clear alfa d
%	Csca_1 = 0.5 * (Csca_d - C_sca);
%	Csca_2 = 0.25 * (Csca_d + Csca_ao - 2 * C_sca);
%	xj (j) = nu;
%	yj (j) = abs(Csca_2);
%	stpn=stpn+1;
xyj=xyj+0.01;
end
%	plot(xj, yj, 'm.', 'markersize', 9);
%	clear alphae C_sca Csca_2 Csca_d Csca_ao nu k xj yj ynewj ynew chto
end

if nomer == 4,
for j=1:200;
	nu = 1700 + 0.1*stpn;
	k = 2 * pi * nu;
	omega = k * 3.000E10;
	d = dmin;
	alfa
	Csca_d = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	d = dmin + 2 * ao;
	alfa
	Csca_ao = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	d = dmin + ao;
	alfa
	C_sca = 8 * pi * k ^ 4 * abs(alphae) ^ 2 / 3;
	clear alfa d
	Csca_1 = 0.5 * (Csca_d - C_sca);
	Csca_2 = 0.25 * (Csca_d + Csca_ao - 2 * C_sca);
	xj (j) = nu;
	yj (j) = abs(Csca_2);
	stpn=stpn+1;
end
	plot(xj, yj, 'm.', 'markersize', 9);
	clear alphae C_sca Csca_2 Csca_d Csca_ao nu k xj yj ynewj ynew chto
end






plot_vs.m

%ic direct vs kramers-kronig
%nu = (100 * round (omega_r/(6.000E12 * pi)) - 100);

%eps_inf= input ('Enter epsilon inf for iron carbonyl (1 to 10):  ');
ic_shapka
   dirre = zeros (1, 600);
   krkr = zeros (1, 600);
   nuj = zeros (1, 600);
for j=1:600;
   nuj (j) = j + 1700;
   nuc = j + 1700;
   nu = nuc;
   %abs5eps
   abs3eps
   dirre (j) = eps_re;
   krkr (j) = epsilonre;
end
plot (nuj, krkr, '.', 'markersize', 2);
hold on
plot (nuj, dirre, 'k.', 'markersize', 2);
%save fic1_05 dirre krkr
xlabel('wavenumber, [cm-1]'); ylabel('epsilon Re');
axis square;
%hold off






sca.m

%scattering cross-section for fixed distance from the surface vs. wavenumber
clear global
nf_mode = 4;

material1 = input ('Enter the tip material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for coated tip):  ');
if material1 == 7,
	material_0 = input ('Enter tip core material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for coated tip):  ');
	material_1 = input ('Enter tip coating material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride):  ');
	ad = input ('Enter tip core radius (10 to 100), nm:  ');
	a = ad * 1.000E-7 + input ('Enter the tip coating thickness (25), nm:  ') * 1.000E-7;
else
   %type in the end of tip radius of curvature (cm)
	Rparticle1 = input ('Enter the end-of-tip radius of curvature (10 to 120), nm: ') * 1.000E-7;
end

material2 = input ('Enter sample material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for reference):  ');
if material2 == 7,
   resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
	omega_r = 6.000E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   REFgamma = input ('Enter the reference gamma, cm-1:  ');
   REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   vFermi = input ('Enter the sample vFermi (1.4000E8 for noble metals):  ');
end

%type in the particle radius (cm)
Rparticle2 = input ('Enter sample particle radius (1 to inf), nm:  ');

%type in the tip-surface distance (cm)
d = 1.000E-7 * input ('Enter tip-sample distance, nm:   ');

figure %plot Csca vs. nu
hold on

nomer = 1;
plot_sca

done = input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   nf
	nomer = 2;
	plot_sca
end

done = input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   nf
	nomer = 3;
	plot_sca
end
hold off






scad.m

%scattering for single wavenumber vs. distance
clear global
nf_mode = 6;

material1 = input ('Enter the tip material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for coated tip):  ');
if material1 == 0,
   %material1 = input ('Enter the sample material code (8 for default, 0 for iron carbonyl):  ');
   %double-check
	if material1 == 0,
		ic_shapka
	else
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   REFparticle = inf;
   %vFermi = input ('Enter the sample vFermi (e.g., 1.4E8 for noble metals):  ');
   vFermi = 1.4E8;
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
   resonance_type = 2;
   %omega_r = 6E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   omega_r = 6E10 * pi * 1815;
   %REFgamma = 6E10 * pi * input ('Enter the reference gamma, cm-1:  ');
   REFgamma = 6E10 * pi * 15;
   Mif = 1E30 * input ('Enter absorption strength (1 to 1000):  ');
   %sigma_eps =  1E15 * input ('Enter sigma, typically much less than 1:  ');
   sigma_eps = 0;
   trial_number = input('Enter the model number:  ');
	end
end

if material1 == 7,
	ad = input ('Enter tip core radius (10 to 100), nm:  ');
	a = ad * 1.000E-7 + input ('Enter the tip coating thickness (25), nm:  ') * 1.000E-7;
	material_0 = input ('Enter tip core material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for iron carbonyl):  ');
   
   
   
   if material_0 == 7,
		ic_shapka
	else
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
%   REFparticle = inf;
   %vFermi = input ('Enter the sample vFermi (e.g., 1.4E8 for noble metals):  ');
%   vFermi = 1.4E8;
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
%   resonance_type = 2;
   %omega_r = 6E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
%   omega_r = 6E10 * pi * 1815;
   %REFgamma = 6E10 * pi * input ('Enter the reference gamma, cm-1:  ');
%   REFgamma = 6E10 * pi * 15;
%   Mif = 1E30 * input ('Enter absorption strength (1 to 1000):  ');
   %sigma_eps =  1E15 * input ('Enter sigma, typically much less than 1:  ');
%   sigma_eps = 0;
%   trial_number = input('Enter the model number:  ');
	end
   
   
   
   material_1 = input ('Enter tip coating material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride):  ');
else
   %type in the end of tip radius of curvature (cm)
	Rparticle1 = input ('Enter the end-of-tip radius of curvature (10 to 120), nm: ') * 1.000E-7;
end

%type in the particle radius (cm)
Rparticle2 = input ('Enter sample particle radius (1 to inf), nm:  ');

material2 = input ('Enter sample material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for reference):  ');
if material2 == 0,
	ic_shapka
end

if material2 == 7,
   
   
   
	material2 = input ('Enter the sample material code (0 for default, 7 for iron carbonyl):  ');
	if material2 == 7,
		ic_shapka
	else
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   REFparticle = inf;
   %vFermi = input ('Enter the sample vFermi (e.g., 1.4E8 for noble metals):  ');
   vFermi = 1.4E8;
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
   resonance_type = 2;
   %omega_r = 6E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   omega_r = 6E10 * pi * 1815;
   %REFgamma = 6E10 * pi * input ('Enter the reference gamma, cm-1:  ');
   REFgamma = 6E10 * pi * 15;
   Mif = 1E30 * input ('Enter absorption strength (1 to 1000):  ');
   %sigma_eps =  1E15 * input ('Enter sigma, typically much less than 1:  ');
   sigma_eps = 0;
   trial_number = input('Enter the model number:  ');
	end
  
   
   
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
	%omega_r = 6.000E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   %REFgamma = input ('Enter the reference gamma, cm-1:  ');
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   %vFermi = input ('Enter the sample vFermi (1.4000E8 for noble metals):  ');
end

%type in the wavenumber below
nu = input ('Enter wavenumber (800-1950), cm-1:  ');

k = 2 * pi * nu;
omega = k * 3.000E10;

stpn=0;

figure %plot Csca vs. d
hold on

nomer = 1;
plot_scad

done = input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   nf
	nomer = 2;
	plot_scad
end

done = input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   nf
	nomer = 3;
	plot_scad
end
hold off






scatamp.m

%scattering for single wavenumber vs. amplitude of cantilever oscillation
clear global
nf_mode = 5;

material1 = input ('Enter the tip material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for coated tip):  ');
if material1 == 7,
	material_0 = input ('Enter tip core material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for coated tip):  ');
	material_1 = input ('Enter tip coating material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride):  ');
	ad = input ('Enter tip core radius (10 to 100), nm:  ');
	a = ad * 1.000E-7 + input ('Enter the tip coating thickness (25), nm:  ') * 1.000E-7;
else
   %type in the end of tip radius of curvature (cm)
	Rparticle1 = input ('Enter the end-of-tip radius of curvature (10 to 120), nm: ') * 1.000E-7;
end

material2 = input ('Enter sample material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for reference):  ');
if material2 == 7,
   resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
	omega_r = 6.000E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   REFgamma = input ('Enter the reference gamma, cm-1:  ');
   REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   vFermi = input ('Enter the sample vFermi (1.4000E8 for noble metals):  ');
end

%type in the particle radius (cm)
Rparticle2 = input ('Enter particle radius (1 to inf), nm:  ');

%type in the minimum tip-surface distance (cm)
dmin = 1.000E-7 * input ('Enter minimum tip-sample distance, nm:   ');

%type in the wavenumber below
nu = input ('Enter wavenumber (800-1950), cm-1:  ');

k = 2 * pi * nu;
omega = k * 3.000E10;

stpn=0;

figure %plot Csca vs. amplitude
hold on

nomer = 1;
plot_scatamp

done = input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   nf
	nomer = 2;
	plot_scatamp
end

done = input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   nf
	nomer = 3;
	plot_scatamp
end
hold off






scawn.m

%scattering for fixed distance from the surface vs. wavenumber (800-1950 cm-1)
clear global
nf_mode = 0;

material1 = input ('Enter the tip material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for coated tip):  ');
if material1 == 0,
   %material1 = input ('Enter the sample material code (8 for default, 0 for iron carbonyl):  ');
   %double-check
	if material1 == 0,
		ic_shapka
	else
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   REFparticle = inf;
   %vFermi = input ('Enter the sample vFermi (e.g., 1.4E8 for noble metals):  ');
   vFermi = 1.4E8;
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
   resonance_type = 2;
   %omega_r = 6E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   omega_r = 6E10 * pi * 1815;
   %REFgamma = 6E10 * pi * input ('Enter the reference gamma, cm-1:  ');
   REFgamma = 6E10 * pi * 15;
   Mif = 1E30 * input ('Enter absorption strength (1 to 1000):  ');
   %sigma_eps =  1E15 * input ('Enter sigma, typically much less than 1:  ');
   sigma_eps = 0;
   trial_number = input('Enter the model number:  ');
	end
end

if material1 == 7,
	ad = input ('Enter tip core radius (10 to 100), nm:  ');
	a = ad * 1.000E-7 + input ('Enter the tip coating thickness (25), nm:  ') * 1.000E-7;
	material_0 = input ('Enter tip core material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for iron carbonyl):  ');
   
   
   
   if material_0 == 7,
		ic_shapka
	else
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
%   REFparticle = inf;
   %vFermi = input ('Enter the sample vFermi (e.g., 1.4E8 for noble metals):  ');
%   vFermi = 1.4E8;
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
%   resonance_type = 2;
   %omega_r = 6E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
%   omega_r = 6E10 * pi * 1815;
   %REFgamma = 6E10 * pi * input ('Enter the reference gamma, cm-1:  ');
%   REFgamma = 6E10 * pi * 15;
%   Mif = 1E30 * input ('Enter absorption strength (1 to 1000):  ');
   %sigma_eps =  1E15 * input ('Enter sigma, typically much less than 1:  ');
%   sigma_eps = 0;
%   trial_number = input('Enter the model number:  ');
	end
   
   
   
   material_1 = input ('Enter tip coating material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride):  ');
else
   %type in the end of tip radius of curvature (cm)
	Rparticle1 = input ('Enter the end-of-tip radius of curvature (10 to 120), nm: ') * 1.000E-7;
end

%type in the particle radius (cm)
Rparticle2 = input ('Enter sample particle radius (1 to inf), nm:  ') * 1.000E-7;

material2 = input ('Enter sample material code (0 for default, 1 for gold, 2 for platinum, 3 for silver, 4 for silicon, 5 for silicon fit, 6 for silicon nitride, 7 for reference):  ');
if material2 == 0,
   material2 = 7;
end
if material2 == 7,
   
   
   
	material2 = input ('Enter the sample material code (0 for default, 7 for iron carbonyl):  ');
	if material2 == 7,
		ic_shapka
	else
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   REFparticle = inf;
   %vFermi = input ('Enter the sample vFermi (e.g., 1.4E8 for noble metals):  ');
   vFermi = 1.4E8;
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
   resonance_type = 2;
   %omega_r = 6E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   omega_r = 6E10 * pi * 1815;
   %REFgamma = 6E10 * pi * input ('Enter the reference gamma, cm-1:  ');
   REFgamma = 6E10 * pi * 15;
   Mif = 1E30 * input ('Enter absorption strength (1 to 1000):  ');
   %sigma_eps =  1E15 * input ('Enter sigma, typically much less than 1:  ');
   sigma_eps = 0;
   trial_number = input('Enter the model number:  ');
	end
  
   
   
   %resonance_type = input ('Enter the type of resonance (1 for real, 2 for imaginary):  ');
	%omega_r = 6.000E10 * pi * input ('Enter the reference resonance line, cm-1:  ');
   %REFgamma = input ('Enter the reference gamma, cm-1:  ');
   %REFparticle = input ('Enter the reference particle radius, cm (inf for very large):  ');
   %vFermi = input ('Enter the sample vFermi (1.4000E8 for noble metals):  ');
end

%type in the minimum tip-surface distance (cm)
dmin = 1.000E-7 * input ('Enter tip-sample distance, nm:   ');

%type in the amplitude of cantilever oscillation (cm)
ao = 1.000E-7 * input ('Enter amplitude of oscillation, nm:   ');

stpn=0;

%figure		%plot Csca vs. nu
hold on

nomer = 1;
plot_scatwn

done = input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   nf
	nomer = 2;
	plot_scatwn
end

done = input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   nf
	nomer = 3;
	plot_scatwn
end

done = input('Are you done? Enter 1 for Yes, 0 for No:  ');
if done == 0,
   nf
	nomer = 4;
	plot_scatwn
end
hold off






sepsilon.m

%sample epsilon

%type in the particle radius (cm)
Rparticle = Rparticle2;
%type in the sample material under this line, e.g. gold, platinum, silver, silicon
if material2 == 0,
	default
end

if material2 == 1,
	gold
end

if material2 == 2,
	platinum
end

if material2 == 3,
	silver
end

if material2 == 4,
	silicon
end

if material2 == 5,
	siepsilon	%silicon fit
end

if material2 == 6,
	sinit		%silicon nitride
end

if material2 == 7,
   eps_ic	%iron carbonyl
   %abs1eps		%reference data
end

es = epsilonre + i * epsilonim;






siepsilon.m

%silicon epsilon for 750 to 2000 wavenumbers

epsilonb = 10.7;
epsilonre = epsilonb + 1;

if omega < 3E10 * 2 * pi * 750,
	omega = 3.000E10 * 2 * pi * 750;
	siepsilon
end

if omega > 3E10 * 2 * pi * 2000,
	omega = 3.000E10 * 2 * pi * 2000;
	siepsilon
end

if omega == 3E10 * 2 * pi * 750,
	epsilonim = 2 * 3.4209 * 0.000259;
end

if omega > 3E10 * 2 * pi * 750, 
	if omega < 3E10 * 2 * pi * 1000, 
		t1 = 3.000E10 * 2 * pi * 750; e1 = 2 * 3.4209 * 0.000259;		
		t2 = 3.000E10 * 2 * pi * 1000; e2 = 2 * 3.4215 * 0.0000676;
		epsilonim = e1 + (e2 - e1) * (omega - t1) / (t2 - t1);
	end
end

if omega == 3E10 * 2 * pi * 1000,
	epsilonim = 2 * 3.4215 * 0.0000676;
end

if omega > 3E10 * 2 * pi * 1000, 
	if omega < 3E10 * 2 * pi * 1110, 
		t1 = 3.000E10 * 2 * pi * 1000; e1 = 2 * 3.4215 * 0.0000676;		
		t2 = 3.000E10 * 2 * pi * 1110; e2 = 2 * 3.4219 * 0.0000846;
		epsilonim = e1 + (e2 - e1) * (omega - t1) / (t2 - t1);
	end
end

if omega == 3E10 * 2 * pi * 1110,
	epsilonim = 2 * 3.4219 * 0.0000846;
end

if omega > 3E10 * 2 * pi * 1110, 
	if omega < 3E10 * 2 * pi * 1428, 
		t1 = 3.000E10 * 2 * pi * 1110; e1 = 2 * 3.4219 * 0.0000846;		
		t2 = 3.000E10 * 2 * pi * 1428; e2 = 2 * 3.4231 * 2.4800E-5;
		epsilonim = e1 + (e2 - e1) * (omega - t1) / (t2 - t1);
	end
end

if omega == 3E10 * 2 * pi * 1428,
	epsilonim = 2 * 3.4231 * 2.4800E-5;
end

if omega > 3E10 * 2 * pi * 1428, 
	if omega < 3E10 * 2 * pi * 1667, 
		t1 = 3.000E10 * 2 * pi * 1428; e1 = 2 * 3.4231 * 2.4800E-5;		
		t2 = 3.000E10 * 2 * pi * 1667; e2 = 2 * 3.4242 * 5.1400E-7;
		epsilonim = e1 + (e2 - e1) * (omega - t1) / (t2 - t1);
	end
end

if omega == 3E10 * 2 * pi * 1667,
	epsilonim = 2 * 3.4242 * 5.1400E-7;
end

if omega > 3E10 * 2 * pi * 1667, 
	if omega < 3E10 * 2 * pi * 2000, 
		t1 = 3.000E10 * 2 * pi * 1667; e1 = 2 * 3.4242 * 5.1400E-7;		
		t2 = 3.000E10 * 2 * pi * 2000; e2 = 2 * 3.4261 * 0.0000002;
		epsilonim = e1 + (e2 - e1) * (omega - t1) / (t2 - t1);
	end
end

if omega == 3E10 * 2 * pi * 2000,
	epsilonim = 2 * 3.4261 * 0.0000002;
end






silicon.m

%silicon IR parameters
%pre-requisites: Rparticle, omega

if omega < 3E10 * 2 * pi * 750,
	omega = 3E10 * 2 * pi * 750;
	silicon
end

if omega > 3E10 * 2 * pi * 2000,
	omega = 3E10 * 2 * pi * 2000;
	silicon
end

N = 5.85E+22;
%omegap2 = 4 * pi * N * (4.80321276758E-10) ^ 2 / (9.10938188E-28);
omegap2 = 1.3E+25;
epsilonb = 10.7;
epsiloninf = epsilonb + 1;
%vFermi = 1.39E8;
OMEGA2 = 5 * 1.39E8 ^ 2 / (3 * Rparticle ^ 2);
%gammab = 1E14;
gamma = 1E14 + 1.39E8 / Rparticle;

epsilonre = epsiloninf + omegap2 * (OMEGA2 - omega ^ 2) / ((OMEGA2 - omega ^ 2) ^ 2 + omega ^ 2 * gamma ^ 2);
epsilonim = omegap2 * omega * gamma / ((OMEGA2 - omega ^ 2) ^ 2 + omega ^ 2 * gamma ^ 2);






silver.m

%silver IR parameters
%pre-requisites: Rparticle, omega

if omega < 3E10 * 2 * pi * 750,
	omega = 3.000E10 * 2 * pi * 750;
	silver
end

if omega > 3E10 * 2 * pi * 2000,
	omega = 3.000E10 * 2 * pi * 2000;
	silver
end

N = 5.8500E+22;
%omegap2 = 4 * pi * N * (4.80321276758E-10) ^ 2 / (9.10938188E-28);
omegap2 = 1.3000E+32;
epsilonb = 3.66;
epsiloninf = epsilonb + 1;
%vFermi = 1.3900E8;
OMEGA2 = 5 * 1.3900E8 ^ 2 / (3 * Rparticle ^ 2);
%gammab = 1.000E14;
gamma = 1.000E14 + 1.3900E8 / Rparticle;

epsilonre = epsiloninf + omegap2 * (OMEGA2 - omega ^ 2) / ((OMEGA2 - omega ^ 2) ^ 2 + omega ^ 2 * gamma ^ 2);
epsilonim = omegap2 * omega * gamma / ((OMEGA2 - omega ^ 2) ^ 2 + omega ^ 2 * gamma ^ 2);

clear N omegap2 epsilonb epsiloninf OMEGA2 gamma






sin.m

%silicon nitride IR parameters

epsilonb = 3;
epsilonre = epsilonb + 1;
epsilonim = 0;






sinit.m

%silicon nitride IR parameters

epsilonb = 3;
epsilonre = epsilonb + 1;
epsilonim = 0;






tipepsilon.m

%tip epsilon

%type in the particle radius (cm)
Rparticle = Rparticle1;
%type in the tip material under this line, e.g. gold, platinum, silver, silicon
if material1 == 0,
	default
end

if material1 == 1,
	gold
end

if material1 == 2,
	platinum
end

if material1 == 3,
	silver
end

if material1 == 4,
	silicon
end

if material1 == 5,
	siepsilon	%silicon fit
end

if material1 == 6,
	sinit		%silicon nitride
end

etip = epsilonre + i * epsilonim;