Terahertz Opto-Electronics Analyzation

This project is to analyze the relationship between the chemical potential 𝜇 and the doping concentration N and the relationship between an auxiliary value ∏02 and the doping concentration N. Here is the report:

The MATLAB code is given below:

clear all;
%% Constant Parameters setting
%Wave-Vector
q=1e8;
%Frequency
w=1e8;
%Mass
m_0=9.10953e-31;
m=0.067*m_0;
%Electron Temperature
T_e=1;
%Boltzmann Constant
k_B=1.3806599e-23;
%Plank Constant
h=1.05459e-34;
%Electron Charge
e=1.6e-19;
%Others
alpha=12;
miu_0=2400;
%% Variable Parameters setting
%B
B=[3.5 1.5 0.5 0.05];
%nsy if for calculate miu-N
%l_B
l_B=sqrt(h./(e*B));
nsy=[10 30 50 300];
%nsy if for calculate Pai-N
nspai=[10 10 30 80];
%Initialize miu
miu=2e-22:0.1e-22:3.2e-21;
%Initialize the result
y=zeros(4,length(miu));
pais=zeros(4,length(miu));
for t=1:4
    %wc
    wc=e*B(t)/m;
    %gama
    gama=sqrt(2*e*h^2*wc*alpha/(pi*m*miu_0));
    %% Calculate N with miu
    %coefficient
    s=1/(pi^2*l_B(t)^2*gama)*sqrt(pi/2);
    % For every miu
    for j=1:length(miu)
        pp=0;
        %Initialize n
        n=0:nsy(t);
        xs=zeros(1,length(n));
        % Travel n from 0 to nsy
        for k=1:length(n)
            % Calculate Landau energy level
            xs(k)=h*(n(k)+0.5)*wc;
            % Set the bound of x
            minx=min([miu(j) xs(k)]);
            maxx=max([miu(j) xs(k)]);
            x=minx-3e-21:2.5e-24:maxx+3e-21;
            % Initialize the result
            zz=zeros(1,length(x));
            % Travel x from lower bound to upper bound
            for i=1:length(x)
                zz(i)=(1/(1+exp((x(i)-miu(j))/(k_B*T_e))))*exp(-((x(i)-xs(k))^2)/(2*(gama)^2));
                pp=pp+zz(i)*2.5e-24;
            end
        end
        % Calculate N
        y(t,j)=s*pp;
    end
    % Clear
    clear j zz x minx maxx;
    %% Calculate Pai with miu
    temp=l_B(t)^2*q^2/2;
    for j=1:length(miu)
        %Initialize n and n'
        n1=0:nspai(t);
        n2=0:nspai(t);
        pai=zeros(length(n1),length(n2));
        L=zeros(length(n1),length(n2));
        C=zeros(length(n1),length(n2));
        % Travel n
        for k=1:length(n1)
            % Travel n'
            for p=max([1 k-15]):min([length(n2) k+15])
                %% Calculate Pai_2
                % Calculate Landau energy level
                xs1=h*(n1(k)+0.5)*wc;
                xs2=h*(n2(p)+0.5)*wc;
                % Set the bound of x
                mins=[miu(j) xs1 xs2];
                minx=min(mins)-h*w;
                maxx=max(mins);
                x=minx-3e-21:1e-24:maxx+3e-21;
                % Initialize the result
                zz=zeros(1,length(x));
                % Travel x from lower bound to upper bound
                for i=1:length(x)
                    zz(i)=(1/(1+exp((x(i)-miu(j))/(k_B*T_e)))-1/(1+exp((x(i)-miu(j)+h*w)/(k_B*T_e))))*exp(-((x(i)-xs1)^2+(x(i)-xs2+h*w)^2)/(2*(gama)^2));
                    pai(k,p)=pai(k,p)+zz(i)*1e-24;
                end
                pai(k,p)=-pai(k,p)/(gama^2);
              %% Calculate L and C
                nmax=max([n1(k) n2(p)]);
                nmin=min([n1(k) n2(p)]);
                % Calculate L
                for s=0:nmin
                    L(k,p)=L(k,p)+(-1)^s*(factorial(nmax)*temp^s)/(factorial(nmax-nmin+s)*factorial(nmin-s)*factorial(s));
                end
                % Calculate C
                C(k,p)=(factorial(nmin)/factorial(nmax))*temp^(nmax-nmin)*exp(-temp)*L(k,p)^2;
                pais(t,j)=pais(t,j)+C(k,p)*pai(k,p);
            end
        end
    end
    pais(t,:)=pais(t,:)/(2*pi*l_B(t)^2);
end
%% Plot Pi-N
figure;
semilogy(y(1,:),-pais(1,:),'b','LineWidth',2);
gtext('B=3.5T', 'color', 'b')
xlabel('N(m^-^2)'),ylabel('\Pi_0_2');
xlim([1e15 5e15])
ylim([1e28 1e34])
title('\Pi_0_2-N');
hold on;
semilogy(y(2,:),-pais(2,:),'r','LineWidth',2);
gtext('B=1.5T', 'color', 'r')
hold on;
semilogy(y(3,:),-pais(3,:),'g','LineWidth',2);
gtext('B=0.5T', 'color', 'g')
hold on;
% Since B=0.05T is strange, the plot is now showed
%semilogy(y(4,:),-pais(4,:),'y','LineWidth',2);
%gtext('B=0.05T', 'color', 'y')
%% Plot miu-N
figure;
plot(y(1,:),miu/1.6e-22,'b','LineWidth',2);
gtext('B=3.5T', 'color', 'b')
xlabel('N(m^-^2)'),ylabel('\mu(meV)');
xlim([1e15 5e15])
title('\mu-N');
hold on
plot(y(2,:),miu/1.6e-22,'r','LineWidth',2);
gtext('B=1.5T', 'color', 'r')
hold on
plot(y(3,:),miu/1.6e-22,'g','LineWidth',2);
gtext('B=0.5T', 'color', 'g')
hold on
plot(y(4,:),miu/1.6e-22,'y','LineWidth',2);
gtext('B=0.05T', 'color', 'y')
hold on