surface conductivity of graphene afunction of frequency in THz

function [ sigma ] = conduction( f , u ) % f is frequency in THZ % u is chemical potential in eV % omega, T and gamma are in the unit of eV % The following equations are for converting unit to eV (to use in the eq.(13) of the reference) % delta is choosed to be zero in the formula omega=1.055e-34*2*pi*1e12*f/1.6e-19; T=4.14e-23/1.6e-19; gamma=0.43e-3; % dw is the parameter for discretization %m is used to change the integal bound and also the bound of the parameter omega dw=1e-4; m=1e6; w=-m*dw:dw:m*dw; n1=1./(exp((w-u)./T)+1); %for fremi distribution (nf(w)) in eq.(13) n2=1./(exp((-w-u)./T)+1); %for nf(-w)in eq(13) d=n1-n2; X=(w./(dw)).*(diff([d d(length(d))])); % the first integral parameter Y=((n2-n1)./((omega+2i*gamma)^2-(4.*w.^2))); % the second integral parameter sigma1=(-2i*(omega+2i*gamma)).*((1/((omega+2i*gamma)^2)).*trapz(dw.*(0:1:m),X((m+1):1:(2*m+1)))); % first term of the equation(13) (trapz function is used for intrgral) sigma2=(-2i*(omega+2i*gamma)).*(-trapz(dw.*(0:1:m),Y((m+1):1:(2*m+1)))); % second term of the equation(13) sigma= 3.86e-2*(sigma1+sigma2); end