How to do double numerical integration and go from finding transmission coefficient in terms of energy to wave vector?

I got the following code in which I want to find my transmission coefficient using the wave vector instead of energy and also I want to compute the current density using the double integration over k_parallel and k_perpendicular values. I have tried doing the integration but I am not sure about it. I need some suggestions.

clc;
clear;
close all;


x1 = -10.847e-10; %Zone B Starting point
x7 = 0e-10; %Zone B Ending Point
h = (x7-x1)/300;
hc = 1.054e-34; %hbar
q = 1.6e-19; % electron charge
N = 303; %Square Matrix dimensions

etaMax = -500;
etaMin = 500;

W = 4.74*q; %work function for Cu (111) flat
%W = 4.46*q; %work function for Cu cylinder
%W = 3.67*q; %work function for Li/Cu cylinder

Ef = 7.05*q; %Fermi Energy, calculate by hand
%Ef=0;
m = 9.11e-31; %mass of electron
pi = 3.14159265359;

n = 10; %No. of iterations
E = linspace(0.00000001*q,W+Ef,n); % Electron Energy (Upper and Lower Limit for Integral)

%F0=0;
F0 = linspace(1e9,10e9,n); %Applied Field, generates a row vector of evenly spaced points between 1e9 and 10e9

    
%% Zone B Potential (Curve fit equation) Sum of sine curve fitting method used here of order 4


data = readmatrix('Cu111.txt');

Error using readmatrix
Unable to find or open 'Cu111.txt'. Check the path and filename or file permissions.

t = data(:,1);
col1 = data(:,2);
col2 = data(:,3);

samples = linspace(t(1),t(291),301);

V1 = spline(t,col1,samples);

V = V1 .* q;
%V=7.2e-19 * ones(1, 301);

data1 = readmatrix('kdistr');
k_parallel = data1(:,1);
k_perp = data1(:,2);
dis = data1(:,3);
k_pa_up = linspace(min(k_parallel), max(k_parallel),n);
k_per_up = linspace(0, max(k_perp),n);
dis_up = linspace(0, max(dis),n);


%        V1 = V./q; 

kB = zeros(301,10);
eta = zeros(1,10);
kA = zeros(1,10);
A = zeros(N,N);
TC = zeros(1,n);
TT = zeros(1,n);
RR = zeros(1,n);
RC = zeros(1,n);
J_abc_100 = zeros(1,n);

b_f_1 = @(kA) exp(1i*kA*x1);
b_f_3 = @(kA) -1i*kA*h*exp(1i*kA*x1);

for dk = 1:n
    Ez1(dk) =sqrt(E(dk)^2 - (hc^2*k_pa_up(dk)^2/(2*m)^2));
end

%Ez1 = linspace(0.00000001,(Ef+W)/q,n); %Range of energy for integral
 kz1=sqrt(2.*m.*Ez1)./hc;
 Ef1 = Ef/q; %Fermi Level Tungsten
 c = 2.086412e12; %constant qmkT/4pihbar^3
 dd=(q^2*6.582119569e-16) / (2 * pi^2 * m);
 temp = 300; %temp in kelvin
 kbolt = 8.617333262145e-5; %bolt constant eV/K
 kTinv = 1/(temp*kbolt); % 1/kT


  %% Wave vectors (kB) using Zone B potential at 301 discrete points     
for f = 1:n
for i = 1:length(E)
    eta(i) = -((E(i)-W-Ef)/q/F0(f))*(2*q*m*F0(f)/hc^2)^(0.3333);
    ke = (2*q*F0(f)*m/hc^2)^(1/3);
    kA(i) = ((2*m*(E(i)-V(1)))^0.5)/hc;
    for j = 1:301
        kB(j,i) = sqrt(2*m*(E(i)-V(j)))/hc;
    end
end

% Building the matrix of equations
           for j = 1:n
               %%% A %%%
               A(1,1) = 1;
               A(1,N-1) = -exp(-1i*kA(j)*x1);
               A(1,N) = 0;

               A(2,N-2) = 1;
               A(2,N) = -(airy(0,eta(j))-1i*airy(2,eta(j))); %airy function

               A(3,1) = 1;
               A(3,2) = -1;
               A(3,N-1) = -1i*h*kA(j)*exp(-1i*kA(j)*x1);

               A(4,N-3) = 1;
               A(4,N-2) = -1;
               A(4,N) = -h*((2*q*m*F0(f)/hc^2)^(1/3))*(airy(1,eta(j))-1i*airy(3,eta(j))); %airy function

               for ii = 5:N
                   A(ii,ii-4) = 1+(h^2*kB(ii-4,j)^2)/12;
                   A(ii,ii-3) = -2*(1-(5/12)*h^2*kB(ii-3,j)^2);
                   A(ii,ii-2) = 1+(h^2*kB(ii-2,j)^2)/12;
               end

               %%% A Matrix %%%

               %%% b Matrix %%%
               b(1,:) = b_f_1(kA(j));
               b(2,:) = 0;
               b(3,:) = b_f_3(kA(j));
               for i = 4:N
                   b(i,:) = 0;
               end
               %%% b %%%
               x11 = A\b; % variables = bA^-1
               R = x11(302);
               T = x11(303);
               x = abs(A\b);
               TC (j) = ((x(303)^2)/pi)*ke/kA(j); %transmission coefficient
               TT (j) = x(303);
               RR (j) = x(302);
               RC (j) = RR(j)^2; %reflection coefficient

           end
          
%% Integral using the associated transmission coefficent
for p = 1:n
    sum = 0;
    
    for f = 1:n
        % Calculate the integrand for the current density
        integrand =dd*(k_pa_up(p) * k_per_up(f) * dis_up(f) * TC(f));
        
        % Add the integrand to the sum
        sum = sum + integrand;
    end
    
    % Store the result in the J_abc_100 array
    J_abc_100(p) = sum;
end

%J_abc_100(f) = c*trapz(Ez1,TC.*log(1+exp(-kTinv.*(Ez1-Ef1)))); % Current Density

end


F0_conv = F0/1e9;
        J_abc_conv = J_abc_100/1e12;
        size = 16;
        semilogy(F0_conv,abs(J_abc_conv),'o-','Color',[0 0.447 0.741],'linewidth',2.8,'markersize',15)
        xlabel('E-Field (GV.m)', 'FontSize',size,'FontWeight', 'bold')
        ylabel('J (A/m^2)', 'FontSize', size,'FontWeight', 'bold')
        legend('Current Density', 'FontSize', size, 'FontWeight', 'bold', 'Location', 'northwest')
        grid on; grid minor
        figure;
        %plot(Ez1,TC,'bo-','linewidth',3,'markersize',15)
        %ylabel('Transmission Coefficient')
        %xlabel('Energy (eV)')
        %grid on
        E2 = E/1.6e-19;
        plot(E2,abs(TC),'r','linewidth',3) %converted energy to eV
        hold off
        set(gca,'FontSize',size,'fontweight','bold')
        xlabel('Energy (eV)','FontSize',size,'fontweight','bold')
        ylabel('Transmission Coefficient','FontSize',size,'fontweight','bold');
        legend('Transmission Coefficient', 'FontSize', size, 'FontWeight', 'bold','Location', 'northwest')
        grid on; grid minor;
        trans = [E2(:), TC(:)];
        curde = [F0_conv(:), J_abc_conv(:)];
        dlmwrite('ztrans.txt',trans, 'delimiter', '\t');
        dlmwrite('zcurde.txt',curde, 'delimiter', '\t');
        
        figure;

                E2 = E/1.6e-19;
        plot(E2,abs(RC),'r','linewidth',3) %converted energy to eV
        hold off
        set(gca,'FontSize',size,'fontweight','bold')
        xlabel('Energy (eV)','FontSize',size,'fontweight','bold')
        ylabel('Reflection Coefficient','FontSize',size,'fontweight','bold');
        legend('Reflection Coefficient', 'FontSize', size, 'FontWeight', 'bold','Location', 'northwest')
        grid on; grid minor;
        trans = [E2(:), RC(:)];
        curde = [F0_conv(:), J_abc_conv(:)];
        dlmwrite('ztrans.txt',trans, 'delimiter', '\t');
        dlmwrite('zcurde.txt',curde, 'delimiter', '\t');
       
        
 

Answer 1

[Walter Roberson]

If you have a 2D numeric array, then usually you would integrate over the whole array by using two trapz() calls. Make sure you pass in the coordinates or unit spacing.

[Yasir Iqbal]

I tried trapz but I am not sure about it, right now I have two arrays in the form of k_parallel and k_perpendicular, what I tried was trapz(k_perpendicular, k_parallel, function, limits of both) but I did not get a good results.

[Walter Roberson]

rowmargins = 0:.2:1;
colmargins = 0:.25:2;
nrow = length(rowmargins);
ncol = length(colmargins);
data_to_integrate = rand(nrow, ncol);

Should we trapz over rows or columns first?

trapz(colmargins, trapz(rowmargins, data_to_integrate))

ans = 0.9257

trapz(rowmargins, trapz(colmargins, data_to_integrate))

Error using trapz
Point spacing must be a scalar specifying uniform spacing or a vector of x-coordinates for each data point.

So the trapz() against rowmargins should be down first, and the result should be passed to trapz with column margins.