I want to plot SPR curve for Reflectance Vs Wavelength in place of angle of incidence.can anyone guide me ?

function SPR_matrix_profile_comparison_v2

%variables

n_Ma=1.493; %målt ved 589nm @ 20 C

n_buffer=1.33156; %HEPES buffer

d_Ma=1E-9;

%vinkelting

theta_start=50.01; %GRADER! !

theta_slut=60.00; %GRADER! !

thetasteps=300;

theta_start=theta_start/360*2*pi ;

theta_slut=theta_slut/360*2*pi ;

deltatheta=(theta_slut-theta_start)/(thetasteps-1);

for t =1:1:thetasteps;

theta=theta_start+deltatheta*(t-1);

vinkel(t)=theta/(2* pi)*360 ;

end

figure(1) ;

hold on

for l =2:8

d_Au=l*10E-9;

[R_list]=spr_calc(n_Ma,n_buffer,d_Ma,d_Au,vinkel,theta_start,theta_slut,deltatheta,thetasteps,theta);

% [y_data]=data_import(vinkel);

if l==2

plot( vinkel , R_list ,' k ')

end

if l==3

plot ( vinkel , R_list ,' g ')

end

if l==4

plot ( vinkel , R_list ,'m')

end

if l==5

plot ( vinkel , R_list ,'c')

end

if l==6

plot ( vinkel , R_list ,' b ')

end

if l==7

plot ( vinkel , R_list , ' y ' )

end

if l==8

plot ( vinkel , R_list , ' r ' )

end

end

%axis ( [ theta_start *180/ pi theta_slut *180/ pi+1e-3 0 1 ] )

legend ( '20nm' , ' 30nm' , ' 40nm' , ' 50nm' , ' 60nm' , ' 70nm' , ' 80nm' ,'Location','Southheast');

%l egend ( 'n_3=1.3300 ' , 'n_3=1.3340 ' , 'n_3=1.3380 ' , 'n_3=1.3420 ' , 'n_3=1.3460 ' , 'n_3=1.3500 ' , ' Locatxlabel ( ' Angle o f incidence ' )

ylabel ( ' Reflectance(R) ' );

xlabel('Angle of incidence');

load SPR_data.txt ;

[ sx , sy ]= size (SPR_data) ;

for i =1:300

y_data ( i )=SPR_data( i +150 ,3)/7E4 ;

%x_data ( i )=i ;

end

% options = fitoptions ( ' NonlinearLeastSquares ' ,method ) ;

%

%

% [ fresult , gof , fout ]= fit (vinkel, y_data , spr_calc , options ) ;

%figure(1) ;

%hold on

%plot( vinkel , R_list )

%hold on ;

%figure ( 1 ) ;

%plot(vinkel , R_list )

%figure (2) ;

%plot(vinkel , y_data ) ;

%figure(2) ;

%plot(vinkel , R_list-y_data ) ;

end

function [R_list] = spr_calc(n_Ma,n_buffer,d_Ma,d_Au,vinkel,theta_start,theta_slut,deltatheta,thetasteps,theta)

lambda_0=780E-9;

n_glass =1.77;

n_Au=0.17-4.93i;

n_Cr=3.11982+3.44408i ; % http : / / refractive index . info/? group=METALS&material=Chromium d.11n_Au=0.17-4.93i ; %Stenberg etal . 1991

n_Ma=n_Ma;

n_buffer=n_buffer ;

d_Cr=1E-9;

d_Au=d_Au;

d_Ma=d_Ma;

% theta_start=51; %GRADER! !

% theta_slut=55; %GRADER! !

% thetasteps=250;

n1=n_glass ;

n2=n_Cr ;

n3=n_Au;

n4=n_Ma;

n5=n_buffer;

eps1=n1^2;

eps2=n2^2;

eps3=n3^2;

eps4=n4^2;

eps5=n5^2;

k=2*pi/lambda_0 ;

% theta_start=theta_start /360*2* pi ;

% theta_slut=the ta_slut /360*2* pi ;

% deltatheta=( theta_slut-theta_start ) / (thetasteps-1);

for t =1:1: thetasteps

theta=theta_start+deltatheta *( t-1);

kx=k* n_glass * sin ( theta ) ;

kz1=sqrt (-(kx^2)+k^2* eps1 ) ;

if imag ( kz1 )>0 ,

kz1=-kz1 ;

end

kz2=sqrt (-(kx^2)+k^2* eps2 ) ;

if imag ( kz2 )>0 ,

kz2=-kz2 ;

end

kz3=sqrt (-(kx^2)+k^2* eps3 ) ;

if imag ( kz3 )>0 ,

kz3=-kz3 ;

end

kz4=sqrt (-(kx^2)+k^2* eps4 ) ;

if imag ( kz4 )>0 ,

kz4=-kz4 ;

end

kz5=sqrt (-(kx^2)+k^2* eps5 ) ;

if imag ( kz5 )>0 ,

kz5=-kz5 ;

end

beta_Cr=kz2*d_Cr ; %2*pi /lambda_0 *(n_Cr*d_Cr* cos ( theta ) ) ;

beta_Au=kz3*d_Au; %2*pi /lambda_0 *(n_Au*d_Au* cos ( theta ) ) ;

beta_Ma=kz4*d_Ma; %2*pi /lambda_0 *(n_Ma*d_Ma* cos ( theta ) ) ;

L_Cr=zeros(2) ;

L_Cr(1 ,1)=exp(-1i *beta_Cr ) ;

L_Cr(2 ,2)=exp (1i *beta_Cr ) ;

L_Au=zeros(2) ;

L_Au(1 ,1)=exp(-1i *beta_Au ) ;

L_Au(2 ,2)=exp (1i *beta_Au ) ;

L_Ma=zeros(2) ;

L_Ma(1 ,1)=exp(-1i *beta_Ma ) ;

L_Ma(2 ,2)=exp (1i *beta_Ma ) ;

b1=(n1/n2 )^2*( kz2 /kz1 ) ;

b2=(n2/n3 )^2*( kz3 /kz2 ) ;

b3=(n3/n4 )^2*( kz4 /kz3 ) ;

b4=(n4/n5 )^2*( kz5 /kz4 ) ;

rho_glass_Cr=(1-b1)/(1+b1 ) ;

rho_Cr_Au=(1-b2)/(1+b2 ) ;

rho_Au_Ma=(1-b3)/(1+b3 ) ;

rho_Ma_buffer=(1-b4)/(1+b4 ) ;

tau_glass_Cr=(2*(n1/n2 ))/(1+b1 ) ;

tau_Cr_Au=(2*(n2/n3 ))/(1+b2 ) ;

tau_Au_Ma=(2*(n3/n4 ))/(1+b3 ) ;

tau_Ma_buffer=(2*(n4/n5 ))/(1+b4 ) ;

H_glass_Cr=ones (2)/ tau_glass_Cr ;

H_glass_Cr(1,2)=rho_glass_Cr / tau_glass_Cr ; %rho i s r e f l e c t i o n c o e f f i c i e n t s

H_glass_Cr(2,1)=rho_glass_Cr / tau_glass_Cr ; %rho i s r e f l e c t i o n c o e f f i c i e n t s

H_Cr_Au=ones(2)/tau_Cr_Au ;

H_Cr_Au(1,2)=rho_Cr_Au/tau_Cr_Au ; %rho i s r e f l e c t i o n c o e f f i c i e n t s

H_Cr_Au(2,1)=rho_Cr_Au/tau_Cr_Au ; %rho i s r e f l e c t i o n c o e f f i c i e n t s

H_Au_Ma=ones (2)/tau_Au_Ma ;

H_Au_Ma(1,2)=rho_Au_Ma/tau_Au_Ma ; %rho i s r e f l e c t i o n c o e f f i c i e n t s

H_Au_Ma(2,1)=rho_Au_Ma/tau_Au_Ma ; %rho i s r e f l e c t i o n c o e f f i c i e n t s

H_Ma_buffer=ones (2)/ tau_Ma_buffer ;

H_Ma_buffer(1,2)=rho_Ma_buffer/tau_Ma_buffer ; %rho i s r e f l e c t i o n c o e f f i c i e n t s

H_Ma_buffer(2,1)=rho_Ma_buffer/tau_Ma_buffer ; %rho i s r e f l e c t i o n c o e f f i c i e n t s

S_glass_buffer=H_glass_Cr*L_Cr*H_Cr_Au*L_Au*H_Au_Ma*L_Ma*H_Ma_buffer ;

R_list(1,t )=( abs(S_glass_buffer ( 1 , 2 ) / S_glass_buffer ( 2 , 2 ) ) )^2 ;

% vinkel (t)=theta /(2* pi ) * 3 6 0 ;

end

% figure(1) ;

% hold on

% plot(vinkel , R_list )

% [ sy , sx ]= s i z e ( R_l i s t ) ;

%

% for i =1:(sx-5)

% theta=the ta_start+deltatheta *(i -1);

% diff_R_list(i)=(-R_list(i)+R_list(i+5))/(5) ;

% diff_vinkel(i)=theta/(2* pi)* 360 ;

% end

%

% for i =1:(sx-6)

% theta=theta_start+deltatheta *( i-1);

% if sign ( diff_R_list ( i ))==sign ( diff_R_list (i +1))

% signchange (i)=0;

% else

% sign_note=theta /(2* pi )*360

% end

% end

%

% figure(2)

% plot ( diff_vinkel , diff_R_list )

end

function [ y_data ] = data_import(vinkel)

load SPR_data.txt ;

[sx,sy]= size(SPR_data) ;

for i =1:250

y_data( i )=SPR_data( i +150 ,2)/7E4 ;

%x_data( i )=i ;

end

x_data=vinkel ;

% hold on ;

%

% figure(2) ;

% plot( x_data , y_data )

end

Answer 1

[Abhinaya Kennedy]

Hi Ananya,

Your code simulates SPR for different gold thicknesses but plots reflectance vs. incident angle. To generate a reflectance vs. wavelength plot, modify the 'spr_calc' function:

1. Loop over wavelengths ('lambda') instead of angles ('theta').
2. Update calculations that depend on wavelength (wavevector k, z-component for each layer 'kz1' to 'kz5').
3. Reflectance calculation remains similar using the updated 'kz' values.

See the modified 'spr_calc' function for details.

function [R_list] = spr_calc(n_Ma, n_buffer, d_Ma, d_Au, lambda_range, theta_start, theta_slut, deltatheta, thetasteps)

% Define material properties (permittivities can be replaced with functions for wavelength dependence)
n_glass = 1.77;
n_Cr = 3.11982 + 3.44408i;
n_Au = 0.17 - 4.93i;
n_Ma = n_Ma;
n_buffer = n_buffer;

% Loop over wavelengths
R_list = zeros(length(lambda_range), thetasteps);  % Pre-allocate reflectance matrix

for l = 1:length(lambda_range)
  lambda = lambda_range(l);
  
  % Calculate wavevector
  k = 2*pi/lambda;

  % Calculate z-component of wavevector in each layer
  for t = 1:thetasteps
    theta = theta_start + deltatheta*(t-1);
    kx = k*n_glass*sin(theta);
    kz1 = sqrt(-(kx^2) + k^2*n_glass^2);
    if imag(kz1) > 0, kz1 = -kz1; end
    kz2 = sqrt(-(kx^2) + k^2*n_Cr^2);
    if imag(kz2) > 0, kz2 = -kz2; end
    kz3 = sqrt(-(kx^2) + k^2*n_Au^2);
    if imag(kz3) > 0, kz3 = -kz3; end
    kz4 = sqrt(-(kx^2) + k^2*n_Ma^2);
    if imag(kz4) > 0, kz4 = -kz4; end
    kz5 = sqrt(-(kx^2) + k^2*n_buffer^2);
    if imag(kz5) > 0, kz5 = -kz5; end
  
    % ... (rest of the calculations using kz1 to kz5 for reflectance)

    R_list(l, t) = (calculated reflectance value);  % Update reflectance based on kz values
  end
end
end

This function now returns a reflectance matrix where rows represent reflectance for a gold layer thickness and columns represent reflectance for a specific wavelength. You can then plot the desired SPR curve.