filebase="Graphene";
filesuffix=["_CQplot-dos-NSP.txt" "_CQplot-dos-SP.txt" "_dos-cv.eps"];
area=5.246e-16;
for l=1:2
    if l == 1
        load DOSCARgraphenesheetuc22aGGAPAWNSP % load the whole data to Workspace for VASP
        dosE = double(DOSCARgraphenesheetuc22aGGAPAWNSP(:,1))--1.23947321; % take 1st column as E and shift to E-E_F
        dosD = double(DOSCARgraphenesheetuc22aGGAPAWNSP(:,2)); % take 2nd column as DOS(E). in 1/eV --> atom/eV. Total number of electrons: \SUM dosDi*(dosEi-dosEj)
        sumD=0.0;
        shiftedtoFermi = 0.0;
        [value,index] = min(abs(dosE - shiftedtoFermi));
        for i = 2:length(dosD)
            sumD=sumD+dosD(i)*(dosE(i)-dosE(i-1));
            if i == index
                x=sprintf('%s: Total number of electrons up to Fermi level: %f',filebase,sumD);
                disp(x);
            end
        end
        x=sprintf('%s: Total number of electrons: %f',filebase,sumD);
        disp(x);
    end
    if l == 2
        load DOSCARgraphenesheetuc22aGGAPAWSP % load the whole data to Workspace for VASP
        dosE = double(DOSCARgraphenesheetuc22aGGAPAWSP(:,1))--1.24006904; % take 1st column as E and shift to E-E_F
        dosDup = double(DOSCARgraphenesheetuc22aGGAPAWSP(:,2)); % take 2nd column as up spin DOS(E)
        dosDdown = double(DOSCARgraphenesheetuc22aGGAPAWSP(:,3)); % take 3rd column as down spin DOS(E)
        dosD=(dosDup+dosDdown); % Total DOS(E). in 1/eV --> atom/eV. Total number of electrons: \SUM dosDi*(dosEi-dosEj)
        sumDup=0.0;
        shiftedtoFermi = 0.0;
        [value,index] = min(abs(dosE - shiftedtoFermi));
        for i = 2:length(dosDup)
            sumDup=sumDup+dosDup(i)*(dosE(i)-dosE(i-1));
            if i == index
                x=sprintf('%s: Total number of "spin up" electrons up to Fermi level: %f',filebase,sumDup);
                disp(x);
            end
        end
        x=sprintf('%s: Total number of "spin up" electrons: %f',filebase,sumDup);
        disp(x);
        %
        sumDdown=0.0;
        for i = 2:length(dosDup)
            sumDdown=sumDdown+dosDdown(i)*(dosE(i)-dosE(i-1));
            if i == index
                x=sprintf('%s: Total number of "spin down" electrons up to Fermi level: %f',filebase,sumDdown);
                disp(x);
            end
        end
        x=sprintf('%s: Total number of "spin down" electrons: %f',filebase,sumDdown);
        disp(x);
    end
    titlesize=18;
    k=1;
    for j=-1.5:0.01:1.5
        CQ(k)=0.0;
        v(k)=j; % v as x axis for the plot
        for i = 2:length(dosE)
            FT = sech((dosE(i)-j)/(2*8.617333262145e-5*300))^2; 
            CQ(k) = CQ(k)+dosD(i)*FT*(dosE(i)-dosE(i-1)); % multiply FT with DOS  and \Delta E : 1/eV*UNITLESS*eV
        end
        e_sqr = (1.602e-19)^2; % Coulombs^2
        CQ(k) = e_sqr*(1/(4*(1.380649e-23)*300))*CQ(k); % CQ multiply by e square.  C^2/J=Farad --> unit is Farad!!
        CQ(k) = (CQ(k)*1e6)/(2*area); % \mu F/cm^2--> unit is microFarad/cm2!!
        x=sprintf('CQ(%d)=%e v(%d)=%f',k,CQ(k),k,v(k));
%         disp(x);
        k=k+1;
    end
    %
    file=strcat(filebase,filesuffix(l));
    fileID = fopen(file,'w');
    k=1;
    for j=-1.5:0.01:1.5
        v(k)=j;
        fprintf(fileID,' %5.2f %12.8f\n',v(k),CQ(k));
        k=k+1;
    end
    fclose(fileID);
    %
    if l == 1
        plot (v,CQ)
        grid on
        hold on;
    end
    if l == 2
        plot (v,CQ)
        xlabel('\phi (Volt)','fontsize',18,'fontweight','b');      % x label
        ylabel('C_Q (\muF/cm^2)','fontsize',18,'fontweight','b');      % y label
        title(filebase,'fontsize',titlesize,'fontweight','b');
        legend('NSP','SP');
        set(0,'defaultfigureposition',[300 300 600 450 ]); 
        filename=strcat(filebase,filesuffix(3)); print(filename,'-r250','-depsc2'); 
        % clear all;
        grid on
        hold off;
    end
end