how to fix this code in calling function
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
zina shadidi
le 12 Nov 2020
Réponse apportée : Walter Roberson
le 12 Nov 2020
%% 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 commentaires
Walter Roberson
le 12 Nov 2020
you did not post your data and you did not indicate which line the error is on.
Réponse acceptée
Walter Roberson
le 12 Nov 2020
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 commentaires
Plus de réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!