%============================================================
% Tony Hyun Kim
% 6.152 MOSCAP analysis
%============================================================

filename = '500nosint_5to-5_sd500m_lightoff_n2off.txt';

close all;

dimension = str2double(filename(1:3));
A = (dimension*1e-4)^2; % Convert to centimeters

% Format: V, C, G
data = load(filename);
data = sortrows(data,1); % Go from left to right
N = length(data(:,1));
data(:,2) = data(:,2)/A; % Now we have capacitance per area

plot(data(:,1),data(:,2),'linewidth',1); hold on;
plot(data(:,1),data(:,2),'o','markersize',2);
xlabel('Bias voltage (V)');
ylabel('Capacitance per area (F/cm^2)');
set(0,'defaulttextinterpreter','none');
title(filename);
set(0,'defaulttextinterpreter','tex');
grid on;

Cox = max(data(:,2));
Cinv = min(data(1:35,2));

%------------------------------------------------------------
% Computations
%------------------------------------------------------------
q = 1.6e-19;
Vth = 0.026; % Thermal voltage, Vth = kT/q
es = 1.05e-12;
eox = 3.45e-13;
Nc = 2.86e19;
ni = 1e10;

% A. Oxide thickness
t = eox/Cox;

% B. Bulk dopant concentration
Cs = Cinv/(Cox-Cinv)*Cox; % Semiconductor capacitance

f = @(Nd) abs(log(Nd/ni)) - q*es*Nd/(4*Cs^2*Vth);
Nd = fzero(f,1e15);

% C. Fixed interface charge

% Flatband voltage
Cfb = (1/Cox + sqrt(Vth/(q*es*Nd)))^-1;
for i = 1:N
    if (data(i,2)>Cfb) 
        break;
    end
end
Vfb = 1/2*(data(i,1)+data(i-1,1));

% Built-in potential
phi_MS = Vth*log(Nd/Nc);

Qf = Cox*(phi_MS - Vfb);

%------------------------------------------------------------
% Print results
%------------------------------------------------------------
fprintf(1,'\n========================================\n');
fprintf(1,filename);
fprintf(1,'\n----------------------------------------');
fprintf(1,'\nC_{ox}*  = %.2e F/cm^2\n',Cox);
fprintf(1,'C_{inv}* = %.2e F/cm^2\n',Cinv);
fprintf(1,'V_{fb}  = %.2f V\n\n',Vfb);
fprintf(1,'t_{ox}  = %.1f Angstroms\n',t*1e8);
fprintf(1,'N_d     = %.2e #/cm^3\n',Nd);
fprintf(1,'Q_f     = %.2e C/cm^2\n',Qf);
fprintf(1,'Q_f/e   = %.2e #/cm^2\n',abs(Qf/q));
fprintf(1,'========================================\n\n');

% Draw on the CV diagram
Xbounds = [min(data(:,1)) max(data(:,1))];
Ybounds = [min(data(:,2)) max(data(:,2))];
plot(Xbounds,[Cox Cox],'r');
plot(Xbounds,[Cinv Cinv],'r');
plot([Vfb Vfb],Ybounds,'r');