double integral and plot of physics function.

Hi.

I'm trying to calculate a function, which is supposed to give me the conductivity for certain functions, defined as a double integral with limits on kx and ky (in my code, I put them as 'x' and 'y'). The problem is how do I perform the double integral of sigman (the function to integrate)? I have tried various ways, but I have no idea how to do it.

Is for my Thesis work.

Thanks.

My code from EDITOR window:

x = linspace(-1.6*pi,1.6*pi,100);
y = linspace(-1.6*pi,1.6*pi,100);
z = linspace(-1.6*pi,1.6*pi,100);
format long
es = 1; 
%% Previous functions for the final integral
deltaE = 0.3; %meV/nm
Ebar = 0.1; %meV/nm
omega = 0;
hbar = 4.135; %meV
e=1 % in eV
dx = sin(y); 
dy = sin(x);
dz = 2 - cos(x) - cos (y) - es; 
d  = sqrt(6 + es.*es + 2.*(es.*cos(x) + es.*cos(y) + cos(kx).*cos(y)) - 4.*(es + cos(x) + cos(y)) );
theta = acos(dz./d);
dvarphikx = (cos(x).*sin(y))./(sin(x).*sin(x) + sin(y).*sin(y));
dvarphiky = (sin(x).*cot(y).*csc(y))./(sin(x).*csc(y) + 1);
dthetakx = sin(x).*((2 - es).*cos(x) - cos(y).^(2) + cos(x).^(2) -2)./ (sqrt((2-cos(x).^(2) - cos(y).*cos(x))./(es.*es + 2.*cos(y).*(es + cos(x) - 2) + 2.*(es-2).*cos(x) - 4.*es + 6)).*(es.*es + 2.*cos(y).*(es + cos(x) -2) + 2.*(es-2).*cos(x) - 4.*es + 6).^(1.5));
dthetaky = sin(y).*((2 - es).*cos(y) - cos(x).^(2) + cos(y).^(2) -2)./ (sqrt((2-cos(x).^(2) - cos(y).*cos(x))./(es.*es + 2.*cos(x).*(es + cos(y) - 2) + 2.*(es-2).*cos(y) - 4.*es + 6)).*(es.*es + 2.*cos(x).*(es + cos(y) -2) + 2.*(es-2).*cos(y) - 4.*es + 6).^(1.5)); 
gc = - e*Ebar.*(sin(theta).*dvarphikx +1i.*dthetakx)/2;
gcc = - e*Ebar.*(sin(theta).*dvarphikx -1i.*dthetakx)/2;
gq = - e*deltaE.*(sin(theta).*dvarphikx +1i.*dthetakx)/2;
gqc= - e*deltaE.*(sin(theta).*dvarphikx -1i.*dthetakx)/2;
gq2 = e^(2)*deltaE^(2).*(dthetakx.^(2) + dvarphikx.^(2).*sin(theta).^(2))/4 ;
gc2 = e^(2)*Ebar^(2).*(dthetakx.^(2) + dvarphikx.^(2).*sin(theta).^(2))/4 ;
Ecplus = sqrt((d - hbar*omega/2).^(2) + gc2 );
Delta = 2*Ecplus - hbar*omega ;
betac = atan(dvarphikx./(sin(theta).*dvarphikx))
alphac = acos((2.*d - hbar*omega)./(sqrt((2.*d -hbar*omega).^(2) + 4.*gc2)));
cos2alphac2 = (sqrt((2.*d - hbar*omega).^(2) + 4.*gc2) + 2.*d - hbar*omega)./(2.*sqrt((2.*d - hbar*omega).*(2.*d - hbar*omega) + 4.*gc2)) ;
sin2alphac2 = (sqrt((2.*d - hbar*omega).^(2) + 4.*gc2) - 2.*d + hbar*omega)./(2.*sqrt((2.*d - hbar*omega).*(2.*d - hbar*omega) + 4.*gc2)) ;
eta = -gq.*cos2alphac2.*exp(-1i.*betac) + gqc.*sin2alphac2.*exp(1i.*betac);
dcosn = abs(eta).*abs(eta)./(Delta.*Delta + 4.*abs(eta).^(2)).^(1.5);
Berrycurv = (cos(x) + cos(y) + cos(x).*cos(y).*(es -2))./((6 + es.*es + 2.*(es.*cos(x) + es.*cos(y) + cos(x).*cos(y)) - 4.*(es + cos(x) + cos(y))).^(1.5));
%Function to integrate from -1.5π, to 0.1 in kx and ky (INTERVAL OF FBZ)
sigman =  (2*Ebar.*dcosn.*Berrycurv )/(4*pi^(2));
%The integral looks like this (in units of e^2/hbar)...

Answer 1

[Rohit Pappu]

A plausible approach would be -

• Add Line no. 4 till the last line in a function

function sigman = BerryCurve(x,y)
    %% All the necessary code
    %% Remove all occurences of kx and ky as it will cause an error
end

• Integrate that function using integral2()

sol = integral2(BerryCurve,xmin,xmax,ymin,ymax);