I want to plot SPR curve for Reflectance Vs Wavelength in place of angle of incidence.can anyone guide me ?
18 views (last 30 days)
Show older comments
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
1 Comment
Rehaan
on 10 Nov 2024 at 8:07
@Ananya Banerjee,Hi ,are you able to resolve the issue regarding wavelength plot,If yes then can your share the things with me.
Answers (1)
Abhinaya Kennedy
on 16 May 2024
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:
- Loop over wavelengths ('lambda') instead of angles ('theta').
- Update calculations that depend on wavelength (wavevector k, z-component for each layer 'kz1' to 'kz5').
- 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.
2 Comments
Rehaan
on 10 Nov 2024 at 8:06
@Abhinaya Kennedy Thanks for helping out the in research SPR sensor simulation
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!