I want to figure out the transmittance intensity of metasurfaces using equation in MATLAB code

Question

I am working on a paper related to metasurfaces. I want to figure out the transmittance intensity using the equation they provide in MATLAB code. But I can't figure out the exact plot. Can anyone help me regarding this issue? Here I upload the equation and desired figure for your reference:

equation transmittance intensity

Here is the code that I tried-

nd=3.45; % refractive index of si
% nd=1.66;
c0=3e8; % light velocity

lambda_el=1340e-9; % electric dipole spectral
lambda_mag=1340e-9; % magnetic dipole spectral

gamma_mag=30e12; % mag dipole reso %6.8
gamma_el=30e12; %ele dipole resonance
lc_a=666e-9; %lattice parameter

w_e=2*pi*c0/lambda_el;
w_m=2*pi*c0/lambda_mag;

lambda = linspace(900*1e-9,1700*1e-9,100);
r = linspace(120*1e-9,320*1e-9,100); % metasurfaces radius

h=220*10^-9;
%----------------------------------
T=zeros(length(lambda),length(r));
% t=zeros(100,1);

%---------------------------
for j=1:length(lambda)
    for i=1:length(r)

a(i) = 1.38* 2 * r(i);
% A(i)=2*pi*r(i)*(h+r(i));
% A(i)=2*pi*r(i)*(h);
A(i)=2*pi*r(i).^2; % area
w(j)=2*pi*c0./lambda(j);
kd(j)=nd*w(j)./c0; % wavevector


alfa0_e(i)=4*A(i)*c0*gamma_el/nd; % electric_polarizability amplitudes
alfa0_m(i)=4*A(i)*c0*gamma_mag/nd; % magnetic_polarizability amplitudes

alfa_e(i,j) = alfa0_e(i)./((w_e.^2)-(w(j).^2)-2i*gamma_el*(w(j))); % electric polarizability

alfa_m(i,j) = alfa0_m(i)./((w_m.^2)-(w(j).^2)-2i*gamma_mag*(w(j))); % magnetic polarizability

  if (w(j)==w_e || w(j) == w_m)
      T(i,j) = 0
  end

T(i,j)= 1 + (alfa_m(i,j) + alfa_e(i,j)) * (1i*kd(j)./(2*A(i))); % transmission intensity

T1(i,j) = (abs(T(i,j))).^2;

t(j) = 1 + 2i*gamma_el*w(j)./((w_e.^2)-(w(j).^2)-2i*gamma_el*(w(j)))+ 2i*gamma_mag*w(j)./((w_m.^2)-(w(j).^2)-2i*gamma_mag*(w(j)));
   
    end
end


imagesc(r, lambda, T1)
colormap jet
colorbar

% zLimits = [0, 1]; % Set your desired limits
% caxis(zLimits);
% zlim([0 1])
% imagesc(r,lambda,t)
% imagesc(r,lambda,abs(alfa_m))
% plot(lambda, t)

My plot look like a red solid backgroud- (https://i.stack.imgur.com/ALdov.png)