In this code, I’m using the Abeles matrix method and Bruggeman’s effective medium approximation to calculate reflection and absorption. The problem is that I’m obtaining negative absorption values at high incident angles above the critical angle. To address this, I considered only the real part of the angle in the transmitted medium, which partially resolved the issue. However, at normal incidence, I’m now getting transmission values exceeding 100%, which is physically impossible. I suspect this might be due to numerical precision issues, as MATLAB may struggle to compute exponentials with very small values. Any ideas on how to fix this? Thank you
clear all
lambda = 500*10^-3:10*10^-3:2500*10^-3; % in microns
A1= 1.4182;
B1= 0.021304;
n_pc=sqrt(1+(A1.*lambda.^2)./(lambda.^2-B1));
wp= 15.92; % eV
f0 = 0.096;
Gamma0 = 0.048; % eV
f1 = 0.100;
Gamma1 = 4.511; % eV
omega1 = 0.174; % eV
f2 = 0.135;
Gamma2 = 1.334; % eV
omega2 = 0.582; % eV
f3 = 0.106;
Gamma3 = 2.178; % eV
omega3 = 1.597; % eV
f4 = 0.729;
Gamma4 = 6.292; % eV
omega4 = 6.089; % eV
OmegaP = sqrt(f0) * wp; % eV
eV = 4.13566733e-1 * 2.99792458 ./ lambda;
epsilon = 1 – OmegaP^2 ./ (eV .* (eV + 1i * Gamma0));
epsilon = epsilon + f1 * wp^2 ./ ((omega1^2 – eV.^2) – 1i * eV * Gamma1);
epsilon = epsilon + f2 * wp^2 ./ ((omega2^2 – eV.^2) – 1i * eV * Gamma2);
epsilon = epsilon + f3 * wp^2 ./ ((omega3^2 – eV.^2) – 1i * eV * Gamma3);
epsilon = epsilon + f4 * wp^2 ./ ((omega4^2 – eV.^2) – 1i * eV * Gamma4);
n = real(sqrt(epsilon));
k = imag(sqrt(epsilon));
nMetal = n+1i.*k;
%NWATER
NWATER = 0.0738.*lambda.^6 – 0.6168.*lambda.^5 + 2.0263.*lambda.^4 – 3.3315.*lambda.^3 + 2.8708.*lambda.^2 – 1.2367.*lambda + 1.5411;
nAir=1+0.05792105./(238.0185-(lambda.^(-2)))+0.00167917./(57.362-(lambda.^(-2)));
% permetivity
epsMetal = nMetal.^2;
epsWATER=NWATER.^2;
%parmenter for each axes
%1 specific volume
ni=n_pc;
nt=NWATER;
%angels
thetaI =0; % incidence angle
thetaI = thetaI/180*pi;
beta=20;
beta = beta/180*pi;
fi=0;
fi =fi/180*pi;
%%%%%
a=30*10.^-9;
b=5*10.^-9;
d=2*a*cos(beta);
f=0.3;
%555
Rpp = [];
Rsp = [];
Rps = [];
Rss = [];
Tpp = [];
Tsp = [];
Tps = [];
Tss = [];
for i=1:size(lambda,2)
%snel law §
thetaT=asin(ni(i)*sin(thetaI)/nt(i));
k0=2*pi/(lambda(i)*10^-6);
VX=ni(i)*sin(thetaI);
% epsilon for compost materals SHAPE Factor
e= sqrt(1-(b/a)^2);
la=((1-e^.2)/e^2)*((log((1+e)/(1-e)))/2*e^2 -1);
lb=0.5*(1-la);
% la=((1-e^.2)/e^2)*((log((1+e)/(1-e)))/2*e^2 -1);
% lb=0.5*(1-la);
epsi=epsMetal(i);
epsh=epsWATER(i);
eps_x = epsh .* (epsh + (la.* (1 – f) + f) .* (epsi – epsh)) ./ (epsh + lb .* (1 – f) .* (epsi – epsh));
eps_y = epsh .* (epsh + (la.* (1 – f) + f) .* (epsi – epsh)) ./ (epsh + lb .* (1 – f) .* (epsi – epsh));
eps_z = epsh .* (epsh + (lb.* (1 – f) + f) .* (epsi – epsh)) ./ (epsh + la .* (1 – f) .* (epsi – epsh));
eps1=eps_x;
eps2=eps_y;
eps3=eps_z;
% eps_xx=eps_2+(eps_1*cos(beta)*cos(beta) +eps_3*sin(beta)*sin(beta) -eps_2)*cos(fi)*cos(fi);
% eps_xy=0.5*(eps_1*cos(beta)*cos(beta)+eps_3*sin(beta)*sin(beta) -eps_2)*sin(2*fi);
% eps_xz=0.5*(eps_3 -eps_1)*sin(2*beta)*cos(fi);
% eps_yx=0.5.*(eps_1.*cos(beta)*cos(beta) +eps_3.*sin(beta).*sin(beta) -eps_2).*sin(2.*fi);
% eps_yy=eps_2 +(eps_1*cos(beta)*cos(beta) +eps_3.*sin(beta).*sin(beta) -eps_2)*sin(fi)*sin(fi);
% eps_yz=0.5*(eps_3-eps_1).*sin(2.*beta).*sin(fi);
% eps_zx=0.5.*(eps_3 – eps_1).*sin(2.*beta).*cos(fi);
% eps_zy=0.5.*(eps_3 – eps_1).*sin(2.*beta).*sin(fi);
% eps_zz=eps_1+(eps_3 – eps_1).*cos(beta)*cos(beta);
% delt11= -(V_x.*eps_zx)./eps_zz;
% delt12=1 – (V_x.*V_x)./eps_zz;
% delt13=-V_x.*eps_zy./eps_zz;
% delt14=0;
% delt21=eps_xx – (eps_xz.*eps_zx)./eps_zz;
% delt22=-V_x.*eps_xz./eps_zz;
% delt23=eps_xy – (eps_xz.*eps_zy)./eps_zz;
% delt24=0;
% delt31=0;
% delt32=0;
% delt33=0;
% delt34=1;
% delt41=eps_yx- (eps_yz*eps_zx) /eps_zz;
% delt42=-V_x*eps_yz /eps_zz;
% delt43=eps_yy -(V_x.^2-eps_yz*eps_zy)./eps_zz;
% delt44=0;
% delt=[delt11 delt12 delt13 delt14 ;delt21 delt22 delt23 delt24 ;delt31 delt32 delt33 delt34;delt41 delt42 delt43 delt44];
P=exp(1)^(1i.*k0.*d.*delt);
a1 = ni(i).*(nt(i).*P(1,2)-cos(thetaT).*P(2,2))+cos(thetaI).*(nt(i).*P(1,1)-cos(thetaT).*P(2,1));
a2 = ni(i).*(nt(i).*P(1,2)-cos(thetaT)*P(2,2))-cos(thetaI).*(nt(i).*P(1,1)-cos(thetaT).*P(2,1));
a3 = (nt(i).*P(1,3)-cos(thetaT).*P(2,3))+ni(i).*cos(thetaI).*(nt(i).*P(1,4)-cos(thetaT).*P(2,4));
a4 = (nt(i).*P(1,3)-cos(thetaT).*P(2,3))-ni(i).*cos(thetaI).*(nt(i).*P(1,4)-cos(thetaT).*P(2,4));
a5 = ni(i).*(nt(i).*cos(thetaT).*P(3,2 )-P(4,2)) +cos(thetaI).*(nt(i).*cos(thetaT).*P(3,1)-P(4,1));
a6 = ni(i).*(nt(i).*cos(thetaT).*P(3,2) -P(4,2))-cos(thetaI).*(nt(i).*cos(thetaT).*P(3,1)-P(4,1));
a7 = (nt(i).*cos(thetaT).*P(3,3)-P(4,3))+ni(i).*cos(thetaI).*(nt(i).*P(3,4).*cos(thetaT)-P(4,4));
a8 = (nt(i).*cos(thetaT).*P(3,3)-P(4,3))-ni(i).*cos(thetaI).*(nt(i).*P(3,4).*cos(thetaT)-P(4,4));
b1 = (ni(i).*P(2,2)+cos(thetaI).*P(2,1))./nt(i);
b2 = (ni(i).*P(2,2)-cos(thetaI).*P(2,1))./nt(i);
b3 = (P(2,3)-ni(i).*cos(thetaI).*P(2,4))./nt(i);
b4 = (P(2,3)+ni(i).*cos(thetaI).*P(2,4))./nt(i);
b5 = ni(i).*P(3,2)+cos(thetaI).*P(3,1);
b6 = ni(i).*P(3,2)-cos(thetaI).*P(3,1);
b7 = P(3,3)-ni(i).*cos(thetaI).*P(3,4);
b8 = P(3,3)+ni(i).*cos(thetaI).*P(3,4);
rpp=(a1.*a8-a4.*a5)./(a4.*a6-a2.*a8);
rps=(a3.*a8-a4.*a7)./(a4.*a6-a2.*a8);
rsp=(a2.*a5-a1.*a6)./(a4.*a6-a2.*a8);
rss=(a2.*a7-a6.*a3)./(a4.*a6-a2.*a8);
tpp=b1+b2.*rpp+b3.*rsp;
tps=b4+b2.*rps+b3.*rss;
tsp=b5+b6.*rpp+b7.*rsp;
tss=b8+b6.*rps+b7.*rss;
Rpp = [Rpp (abs(rpp))^2];
Rsp = [Rsp (abs(rsp))^2];
Rps = [Rps (abs(rps))^2];
Rss = [Rss (abs(rss))^2];
Tpp = [Tpp (abs(tpp .*tpp.* (nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI))))];
Tsp = [Tsp (abs(tsp .*tsp.* nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI)))];
Tps = [Tps (abs(tps .*tps.* nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI)))];
Tss = [Tss (abs(tss .*tss.* nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI)))];
T=(Tpp+Tsp+Tps+Tss)/2;
R=(Rpp+Rsp+Rps+Rss)/2;
A=1-R-T;
end
figure(1)
hold on
plot(lambda*10^3,A*100,’LineWidth’, 3)
hold on
xlabel(‘Wavelength (nm)’,’fontsize’,16)
hold on
ylabel(‘Absorbtion’,’fontsize’, 16)In this code, I’m using the Abeles matrix method and Bruggeman’s effective medium approximation to calculate reflection and absorption. The problem is that I’m obtaining negative absorption values at high incident angles above the critical angle. To address this, I considered only the real part of the angle in the transmitted medium, which partially resolved the issue. However, at normal incidence, I’m now getting transmission values exceeding 100%, which is physically impossible. I suspect this might be due to numerical precision issues, as MATLAB may struggle to compute exponentials with very small values. Any ideas on how to fix this? Thank you
clear all
lambda = 500*10^-3:10*10^-3:2500*10^-3; % in microns
A1= 1.4182;
B1= 0.021304;
n_pc=sqrt(1+(A1.*lambda.^2)./(lambda.^2-B1));
wp= 15.92; % eV
f0 = 0.096;
Gamma0 = 0.048; % eV
f1 = 0.100;
Gamma1 = 4.511; % eV
omega1 = 0.174; % eV
f2 = 0.135;
Gamma2 = 1.334; % eV
omega2 = 0.582; % eV
f3 = 0.106;
Gamma3 = 2.178; % eV
omega3 = 1.597; % eV
f4 = 0.729;
Gamma4 = 6.292; % eV
omega4 = 6.089; % eV
OmegaP = sqrt(f0) * wp; % eV
eV = 4.13566733e-1 * 2.99792458 ./ lambda;
epsilon = 1 – OmegaP^2 ./ (eV .* (eV + 1i * Gamma0));
epsilon = epsilon + f1 * wp^2 ./ ((omega1^2 – eV.^2) – 1i * eV * Gamma1);
epsilon = epsilon + f2 * wp^2 ./ ((omega2^2 – eV.^2) – 1i * eV * Gamma2);
epsilon = epsilon + f3 * wp^2 ./ ((omega3^2 – eV.^2) – 1i * eV * Gamma3);
epsilon = epsilon + f4 * wp^2 ./ ((omega4^2 – eV.^2) – 1i * eV * Gamma4);
n = real(sqrt(epsilon));
k = imag(sqrt(epsilon));
nMetal = n+1i.*k;
%NWATER
NWATER = 0.0738.*lambda.^6 – 0.6168.*lambda.^5 + 2.0263.*lambda.^4 – 3.3315.*lambda.^3 + 2.8708.*lambda.^2 – 1.2367.*lambda + 1.5411;
nAir=1+0.05792105./(238.0185-(lambda.^(-2)))+0.00167917./(57.362-(lambda.^(-2)));
% permetivity
epsMetal = nMetal.^2;
epsWATER=NWATER.^2;
%parmenter for each axes
%1 specific volume
ni=n_pc;
nt=NWATER;
%angels
thetaI =0; % incidence angle
thetaI = thetaI/180*pi;
beta=20;
beta = beta/180*pi;
fi=0;
fi =fi/180*pi;
%%%%%
a=30*10.^-9;
b=5*10.^-9;
d=2*a*cos(beta);
f=0.3;
%555
Rpp = [];
Rsp = [];
Rps = [];
Rss = [];
Tpp = [];
Tsp = [];
Tps = [];
Tss = [];
for i=1:size(lambda,2)
%snel law §
thetaT=asin(ni(i)*sin(thetaI)/nt(i));
k0=2*pi/(lambda(i)*10^-6);
VX=ni(i)*sin(thetaI);
% epsilon for compost materals SHAPE Factor
e= sqrt(1-(b/a)^2);
la=((1-e^.2)/e^2)*((log((1+e)/(1-e)))/2*e^2 -1);
lb=0.5*(1-la);
% la=((1-e^.2)/e^2)*((log((1+e)/(1-e)))/2*e^2 -1);
% lb=0.5*(1-la);
epsi=epsMetal(i);
epsh=epsWATER(i);
eps_x = epsh .* (epsh + (la.* (1 – f) + f) .* (epsi – epsh)) ./ (epsh + lb .* (1 – f) .* (epsi – epsh));
eps_y = epsh .* (epsh + (la.* (1 – f) + f) .* (epsi – epsh)) ./ (epsh + lb .* (1 – f) .* (epsi – epsh));
eps_z = epsh .* (epsh + (lb.* (1 – f) + f) .* (epsi – epsh)) ./ (epsh + la .* (1 – f) .* (epsi – epsh));
eps1=eps_x;
eps2=eps_y;
eps3=eps_z;
% eps_xx=eps_2+(eps_1*cos(beta)*cos(beta) +eps_3*sin(beta)*sin(beta) -eps_2)*cos(fi)*cos(fi);
% eps_xy=0.5*(eps_1*cos(beta)*cos(beta)+eps_3*sin(beta)*sin(beta) -eps_2)*sin(2*fi);
% eps_xz=0.5*(eps_3 -eps_1)*sin(2*beta)*cos(fi);
% eps_yx=0.5.*(eps_1.*cos(beta)*cos(beta) +eps_3.*sin(beta).*sin(beta) -eps_2).*sin(2.*fi);
% eps_yy=eps_2 +(eps_1*cos(beta)*cos(beta) +eps_3.*sin(beta).*sin(beta) -eps_2)*sin(fi)*sin(fi);
% eps_yz=0.5*(eps_3-eps_1).*sin(2.*beta).*sin(fi);
% eps_zx=0.5.*(eps_3 – eps_1).*sin(2.*beta).*cos(fi);
% eps_zy=0.5.*(eps_3 – eps_1).*sin(2.*beta).*sin(fi);
% eps_zz=eps_1+(eps_3 – eps_1).*cos(beta)*cos(beta);
% delt11= -(V_x.*eps_zx)./eps_zz;
% delt12=1 – (V_x.*V_x)./eps_zz;
% delt13=-V_x.*eps_zy./eps_zz;
% delt14=0;
% delt21=eps_xx – (eps_xz.*eps_zx)./eps_zz;
% delt22=-V_x.*eps_xz./eps_zz;
% delt23=eps_xy – (eps_xz.*eps_zy)./eps_zz;
% delt24=0;
% delt31=0;
% delt32=0;
% delt33=0;
% delt34=1;
% delt41=eps_yx- (eps_yz*eps_zx) /eps_zz;
% delt42=-V_x*eps_yz /eps_zz;
% delt43=eps_yy -(V_x.^2-eps_yz*eps_zy)./eps_zz;
% delt44=0;
% delt=[delt11 delt12 delt13 delt14 ;delt21 delt22 delt23 delt24 ;delt31 delt32 delt33 delt34;delt41 delt42 delt43 delt44];
P=exp(1)^(1i.*k0.*d.*delt);
a1 = ni(i).*(nt(i).*P(1,2)-cos(thetaT).*P(2,2))+cos(thetaI).*(nt(i).*P(1,1)-cos(thetaT).*P(2,1));
a2 = ni(i).*(nt(i).*P(1,2)-cos(thetaT)*P(2,2))-cos(thetaI).*(nt(i).*P(1,1)-cos(thetaT).*P(2,1));
a3 = (nt(i).*P(1,3)-cos(thetaT).*P(2,3))+ni(i).*cos(thetaI).*(nt(i).*P(1,4)-cos(thetaT).*P(2,4));
a4 = (nt(i).*P(1,3)-cos(thetaT).*P(2,3))-ni(i).*cos(thetaI).*(nt(i).*P(1,4)-cos(thetaT).*P(2,4));
a5 = ni(i).*(nt(i).*cos(thetaT).*P(3,2 )-P(4,2)) +cos(thetaI).*(nt(i).*cos(thetaT).*P(3,1)-P(4,1));
a6 = ni(i).*(nt(i).*cos(thetaT).*P(3,2) -P(4,2))-cos(thetaI).*(nt(i).*cos(thetaT).*P(3,1)-P(4,1));
a7 = (nt(i).*cos(thetaT).*P(3,3)-P(4,3))+ni(i).*cos(thetaI).*(nt(i).*P(3,4).*cos(thetaT)-P(4,4));
a8 = (nt(i).*cos(thetaT).*P(3,3)-P(4,3))-ni(i).*cos(thetaI).*(nt(i).*P(3,4).*cos(thetaT)-P(4,4));
b1 = (ni(i).*P(2,2)+cos(thetaI).*P(2,1))./nt(i);
b2 = (ni(i).*P(2,2)-cos(thetaI).*P(2,1))./nt(i);
b3 = (P(2,3)-ni(i).*cos(thetaI).*P(2,4))./nt(i);
b4 = (P(2,3)+ni(i).*cos(thetaI).*P(2,4))./nt(i);
b5 = ni(i).*P(3,2)+cos(thetaI).*P(3,1);
b6 = ni(i).*P(3,2)-cos(thetaI).*P(3,1);
b7 = P(3,3)-ni(i).*cos(thetaI).*P(3,4);
b8 = P(3,3)+ni(i).*cos(thetaI).*P(3,4);
rpp=(a1.*a8-a4.*a5)./(a4.*a6-a2.*a8);
rps=(a3.*a8-a4.*a7)./(a4.*a6-a2.*a8);
rsp=(a2.*a5-a1.*a6)./(a4.*a6-a2.*a8);
rss=(a2.*a7-a6.*a3)./(a4.*a6-a2.*a8);
tpp=b1+b2.*rpp+b3.*rsp;
tps=b4+b2.*rps+b3.*rss;
tsp=b5+b6.*rpp+b7.*rsp;
tss=b8+b6.*rps+b7.*rss;
Rpp = [Rpp (abs(rpp))^2];
Rsp = [Rsp (abs(rsp))^2];
Rps = [Rps (abs(rps))^2];
Rss = [Rss (abs(rss))^2];
Tpp = [Tpp (abs(tpp .*tpp.* (nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI))))];
Tsp = [Tsp (abs(tsp .*tsp.* nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI)))];
Tps = [Tps (abs(tps .*tps.* nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI)))];
Tss = [Tss (abs(tss .*tss.* nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI)))];
T=(Tpp+Tsp+Tps+Tss)/2;
R=(Rpp+Rsp+Rps+Rss)/2;
A=1-R-T;
end
figure(1)
hold on
plot(lambda*10^3,A*100,’LineWidth’, 3)
hold on
xlabel(‘Wavelength (nm)’,’fontsize’,16)
hold on
ylabel(‘Absorbtion’,’fontsize’, 16) In this code, I’m using the Abeles matrix method and Bruggeman’s effective medium approximation to calculate reflection and absorption. The problem is that I’m obtaining negative absorption values at high incident angles above the critical angle. To address this, I considered only the real part of the angle in the transmitted medium, which partially resolved the issue. However, at normal incidence, I’m now getting transmission values exceeding 100%, which is physically impossible. I suspect this might be due to numerical precision issues, as MATLAB may struggle to compute exponentials with very small values. Any ideas on how to fix this? Thank you
clear all
lambda = 500*10^-3:10*10^-3:2500*10^-3; % in microns
A1= 1.4182;
B1= 0.021304;
n_pc=sqrt(1+(A1.*lambda.^2)./(lambda.^2-B1));
wp= 15.92; % eV
f0 = 0.096;
Gamma0 = 0.048; % eV
f1 = 0.100;
Gamma1 = 4.511; % eV
omega1 = 0.174; % eV
f2 = 0.135;
Gamma2 = 1.334; % eV
omega2 = 0.582; % eV
f3 = 0.106;
Gamma3 = 2.178; % eV
omega3 = 1.597; % eV
f4 = 0.729;
Gamma4 = 6.292; % eV
omega4 = 6.089; % eV
OmegaP = sqrt(f0) * wp; % eV
eV = 4.13566733e-1 * 2.99792458 ./ lambda;
epsilon = 1 – OmegaP^2 ./ (eV .* (eV + 1i * Gamma0));
epsilon = epsilon + f1 * wp^2 ./ ((omega1^2 – eV.^2) – 1i * eV * Gamma1);
epsilon = epsilon + f2 * wp^2 ./ ((omega2^2 – eV.^2) – 1i * eV * Gamma2);
epsilon = epsilon + f3 * wp^2 ./ ((omega3^2 – eV.^2) – 1i * eV * Gamma3);
epsilon = epsilon + f4 * wp^2 ./ ((omega4^2 – eV.^2) – 1i * eV * Gamma4);
n = real(sqrt(epsilon));
k = imag(sqrt(epsilon));
nMetal = n+1i.*k;
%NWATER
NWATER = 0.0738.*lambda.^6 – 0.6168.*lambda.^5 + 2.0263.*lambda.^4 – 3.3315.*lambda.^3 + 2.8708.*lambda.^2 – 1.2367.*lambda + 1.5411;
nAir=1+0.05792105./(238.0185-(lambda.^(-2)))+0.00167917./(57.362-(lambda.^(-2)));
% permetivity
epsMetal = nMetal.^2;
epsWATER=NWATER.^2;
%parmenter for each axes
%1 specific volume
ni=n_pc;
nt=NWATER;
%angels
thetaI =0; % incidence angle
thetaI = thetaI/180*pi;
beta=20;
beta = beta/180*pi;
fi=0;
fi =fi/180*pi;
%%%%%
a=30*10.^-9;
b=5*10.^-9;
d=2*a*cos(beta);
f=0.3;
%555
Rpp = [];
Rsp = [];
Rps = [];
Rss = [];
Tpp = [];
Tsp = [];
Tps = [];
Tss = [];
for i=1:size(lambda,2)
%snel law §
thetaT=asin(ni(i)*sin(thetaI)/nt(i));
k0=2*pi/(lambda(i)*10^-6);
VX=ni(i)*sin(thetaI);
% epsilon for compost materals SHAPE Factor
e= sqrt(1-(b/a)^2);
la=((1-e^.2)/e^2)*((log((1+e)/(1-e)))/2*e^2 -1);
lb=0.5*(1-la);
% la=((1-e^.2)/e^2)*((log((1+e)/(1-e)))/2*e^2 -1);
% lb=0.5*(1-la);
epsi=epsMetal(i);
epsh=epsWATER(i);
eps_x = epsh .* (epsh + (la.* (1 – f) + f) .* (epsi – epsh)) ./ (epsh + lb .* (1 – f) .* (epsi – epsh));
eps_y = epsh .* (epsh + (la.* (1 – f) + f) .* (epsi – epsh)) ./ (epsh + lb .* (1 – f) .* (epsi – epsh));
eps_z = epsh .* (epsh + (lb.* (1 – f) + f) .* (epsi – epsh)) ./ (epsh + la .* (1 – f) .* (epsi – epsh));
eps1=eps_x;
eps2=eps_y;
eps3=eps_z;
% eps_xx=eps_2+(eps_1*cos(beta)*cos(beta) +eps_3*sin(beta)*sin(beta) -eps_2)*cos(fi)*cos(fi);
% eps_xy=0.5*(eps_1*cos(beta)*cos(beta)+eps_3*sin(beta)*sin(beta) -eps_2)*sin(2*fi);
% eps_xz=0.5*(eps_3 -eps_1)*sin(2*beta)*cos(fi);
% eps_yx=0.5.*(eps_1.*cos(beta)*cos(beta) +eps_3.*sin(beta).*sin(beta) -eps_2).*sin(2.*fi);
% eps_yy=eps_2 +(eps_1*cos(beta)*cos(beta) +eps_3.*sin(beta).*sin(beta) -eps_2)*sin(fi)*sin(fi);
% eps_yz=0.5*(eps_3-eps_1).*sin(2.*beta).*sin(fi);
% eps_zx=0.5.*(eps_3 – eps_1).*sin(2.*beta).*cos(fi);
% eps_zy=0.5.*(eps_3 – eps_1).*sin(2.*beta).*sin(fi);
% eps_zz=eps_1+(eps_3 – eps_1).*cos(beta)*cos(beta);
% delt11= -(V_x.*eps_zx)./eps_zz;
% delt12=1 – (V_x.*V_x)./eps_zz;
% delt13=-V_x.*eps_zy./eps_zz;
% delt14=0;
% delt21=eps_xx – (eps_xz.*eps_zx)./eps_zz;
% delt22=-V_x.*eps_xz./eps_zz;
% delt23=eps_xy – (eps_xz.*eps_zy)./eps_zz;
% delt24=0;
% delt31=0;
% delt32=0;
% delt33=0;
% delt34=1;
% delt41=eps_yx- (eps_yz*eps_zx) /eps_zz;
% delt42=-V_x*eps_yz /eps_zz;
% delt43=eps_yy -(V_x.^2-eps_yz*eps_zy)./eps_zz;
% delt44=0;
% delt=[delt11 delt12 delt13 delt14 ;delt21 delt22 delt23 delt24 ;delt31 delt32 delt33 delt34;delt41 delt42 delt43 delt44];
P=exp(1)^(1i.*k0.*d.*delt);
a1 = ni(i).*(nt(i).*P(1,2)-cos(thetaT).*P(2,2))+cos(thetaI).*(nt(i).*P(1,1)-cos(thetaT).*P(2,1));
a2 = ni(i).*(nt(i).*P(1,2)-cos(thetaT)*P(2,2))-cos(thetaI).*(nt(i).*P(1,1)-cos(thetaT).*P(2,1));
a3 = (nt(i).*P(1,3)-cos(thetaT).*P(2,3))+ni(i).*cos(thetaI).*(nt(i).*P(1,4)-cos(thetaT).*P(2,4));
a4 = (nt(i).*P(1,3)-cos(thetaT).*P(2,3))-ni(i).*cos(thetaI).*(nt(i).*P(1,4)-cos(thetaT).*P(2,4));
a5 = ni(i).*(nt(i).*cos(thetaT).*P(3,2 )-P(4,2)) +cos(thetaI).*(nt(i).*cos(thetaT).*P(3,1)-P(4,1));
a6 = ni(i).*(nt(i).*cos(thetaT).*P(3,2) -P(4,2))-cos(thetaI).*(nt(i).*cos(thetaT).*P(3,1)-P(4,1));
a7 = (nt(i).*cos(thetaT).*P(3,3)-P(4,3))+ni(i).*cos(thetaI).*(nt(i).*P(3,4).*cos(thetaT)-P(4,4));
a8 = (nt(i).*cos(thetaT).*P(3,3)-P(4,3))-ni(i).*cos(thetaI).*(nt(i).*P(3,4).*cos(thetaT)-P(4,4));
b1 = (ni(i).*P(2,2)+cos(thetaI).*P(2,1))./nt(i);
b2 = (ni(i).*P(2,2)-cos(thetaI).*P(2,1))./nt(i);
b3 = (P(2,3)-ni(i).*cos(thetaI).*P(2,4))./nt(i);
b4 = (P(2,3)+ni(i).*cos(thetaI).*P(2,4))./nt(i);
b5 = ni(i).*P(3,2)+cos(thetaI).*P(3,1);
b6 = ni(i).*P(3,2)-cos(thetaI).*P(3,1);
b7 = P(3,3)-ni(i).*cos(thetaI).*P(3,4);
b8 = P(3,3)+ni(i).*cos(thetaI).*P(3,4);
rpp=(a1.*a8-a4.*a5)./(a4.*a6-a2.*a8);
rps=(a3.*a8-a4.*a7)./(a4.*a6-a2.*a8);
rsp=(a2.*a5-a1.*a6)./(a4.*a6-a2.*a8);
rss=(a2.*a7-a6.*a3)./(a4.*a6-a2.*a8);
tpp=b1+b2.*rpp+b3.*rsp;
tps=b4+b2.*rps+b3.*rss;
tsp=b5+b6.*rpp+b7.*rsp;
tss=b8+b6.*rps+b7.*rss;
Rpp = [Rpp (abs(rpp))^2];
Rsp = [Rsp (abs(rsp))^2];
Rps = [Rps (abs(rps))^2];
Rss = [Rss (abs(rss))^2];
Tpp = [Tpp (abs(tpp .*tpp.* (nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI))))];
Tsp = [Tsp (abs(tsp .*tsp.* nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI)))];
Tps = [Tps (abs(tps .*tps.* nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI)))];
Tss = [Tss (abs(tss .*tss.* nt(i) .* cos(real(thetaT)) ./ni(i) .* cos(thetaI)))];
T=(Tpp+Tsp+Tps+Tss)/2;
R=(Rpp+Rsp+Rps+Rss)/2;
A=1-R-T;
end
figure(1)
hold on
plot(lambda*10^3,A*100,’LineWidth’, 3)
hold on
xlabel(‘Wavelength (nm)’,’fontsize’,16)
hold on
ylabel(‘Absorbtion’,’fontsize’, 16) 4×4 matrix, abeles matrix, negative absorption val, matlab MATLAB Answers — New Questions

​