how to fix this code in calling function
3 次查看(过去 30 天)
显示 更早的评论
%% Call the transfer matrix R, T, A calculator function and build spectra
%%%may be need to % Convert angle to radians
[n1]=xlsread('SiO2_n_k.xlsx');
[n2]=xlsread('TiO2_n_k.xlsx');
L= 20;
i=1:L;
A = n1(i,2); %%%%% SiO2 refractive coefficient
B = n2(i,2); %%%%% TiO2 refractive coefficient
l=n2(i,1); %%%% wavelength
% Preallocate memory
R = zeros(1,L);
T = zeros(1,L);
A = zeros(1,L);
for m = 1:L
% Reflection and transmission coefficients, r, t, not used, so
% replace output with ~. Can add back in if needed.
% [r(m),t(m),R(m),T(m),A(m)]=jreftran_rt(wl(m),d,n(m,:),t0,polarization);
% full form of jreftran_rt
[r,t,R,T,A]=reftran_rt(l,[NaN,100,500,2000],[1,B,A,1.5],t0,polarization)
end
figure
plot(l,R);
figure
plot(l,A)
%%%%%%%%%%%%%%%% the function
function [r,t,R,T,A]=reftran_rt(l,d,n,t0,polarization)
%l = free space wavelength, nm
%d = layer thickness vector, nm
%n = layer complex refractive index vector
%t0= angle of incidence
%polarization should be 0 for TE (s-polarized), otherwise TM (p-polarized)
Z0=376.730313; %impedance of free space, Ohms
%the line below had mistakenly been a Z0/n instead of a n/Z0 in version 1!
Y=n./Z0; %admittance in terms of impedance of free space and refractive index, assuming non-magnetic media
g=1i*2*pi*n/l; %propagation constant in terms of free space wavelength and refractive index
%all of the calculations rely on cosine of the complex angle, but we can
%only find the sine of the complex angle from snells law. So we use the
%fact that cos(asin(x))=sqrt(1-x^2)
%t=asin(n(1)./n*sin(t0)), complex theta for each layer
ct=sqrt(1-(n(1)./n*sin(t0)).^2); %cosine theta
if polarization==0
eta=Y.*ct; %tilted admittance, TE case
else
eta=Y./ct; %tilted admittance, TM case
end
delta=1i*g.*d.*ct;
M=zeros(2,2,length(d));
for j=1:length(d)
M(1,1,j)=cos(delta(j));
M(1,2,j)=1i./eta(j).*sin(delta(j));
M(2,1,j)=1i*eta(j).*sin(delta(j));
M(2,2,j)=cos(delta(j));
end
M_t=[1,0;0,1]; %M total
for j=2:(length(d)-1)
M_t=M_t*M(:,:,j);
end
r=(eta(1)*(M_t(1,1)+M_t(1,2)*eta(end))-(M_t(2,1)+M_t(2,2)*eta(end)))/(eta(1)*(M_t(1,1)+M_t(1,2)*eta(end))+(M_t(2,1)+M_t(2,2)*eta(end)));
t=2*eta(1)/(eta(1)*(M_t(1,1)+M_t(1,2)*eta(end))+(M_t(2,1)+M_t(2,2)*eta(end)));
R=abs(r)^2;
T=real(eta(end)/eta(1))*abs(t)^2;
A=(4*eta(1)*real((M_t(1,1)+M_t(1,2)*eta(end))*conj(M_t(2,1)+M_t(2,2)*eta(end))-eta(end)))/abs(eta(1)*(M_t(1,1)+M_t(1,2)*eta(end))+(M_t(2,1)+M_t(2,2)*eta(end)))^2;
end
2 个评论
Walter Roberson
2020-11-12
you did not post your data and you did not indicate which line the error is on.
采纳的回答
Walter Roberson
2020-11-12
i=1:L;
That is a vector.
A = n1(i,2); %%%%% SiO2 refractive coefficient
B = n2(i,2); %%%%% TiO2 refractive coefficient
You are indexing a 2D array with a vector of i values. The results in A and B is going to be a column vector of length L.
A = zeros(1,L);
You throw away everything you wrote into A, and replace it with a row vector of length L. not a column vector.
[r,t,R,T,A]=reftran_rt(l,[NaN,100,500,2000],[1,B,A,1.5],t0,polarization)
In the sub-expression [1,B,A,1.5], the 1 is a scalar. In order for [] (horzcat) to work, everything after that 1 in the [] has to have one row. However, your B is a column vector of length L. Your A is a row vector of length L, so the A,1.5 part is okay, but the column B cannot go there.
If somehow B did fit, such as if you had transposed B into a row vector earlier, then the term would be a vector of length 1 + L + L + 1 = 2*L+2 . Are you sure that is appropriate?
I was going to say that you do the same calculation for every iteration of the for m loop, but then I noticed that you are updating A and it is the updated A that is going back into the next iteration, so potentially that is reasonable.
... But I can't help but think that what you really want is
[r,t,R,T,A(m)]=reftran_rt(l,[NaN,100,500,2000],[1,B(m),A(m),1.5],t0,polarization);
0 个评论
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!