Skip to content

Commit 651b57c

Browse files
CerilogCerilog
authored
Add files via upload
0 parents  commit 651b57c

4 files changed

+439
-0
lines changed

‎DrudeFitting.m

+137
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,137 @@
1+
%%
2+
% Provide Drude model Fitting with n,k data
3+
% By comparing complex conductivity with classical EM theory
4+
% Date : 2020/03/03
5+
% Check epi_inf at line 12
6+
% Initial guess at line 47
7+
% Revised : 2020/04/10
8+
%% Data I/O and Basic Treatment
9+
D = fscanf(fopen('ThinFilm-RIX-RTA-1e14.txt','r'),'%f %f %f',[3,inf]);
10+
global epi0
11+
epi0 = 8.85e-12;
12+
epi_inf = 10.89;
13+
fTHz = D(1,:);
14+
n = D(2,:);
15+
k = D(3,:);
16+
Range = and(fTHz>=0.2,fTHz<=1.2);
17+
fTHz = fTHz(Range);
18+
n = n(Range);
19+
k = k(Range);
20+
f = fTHz*1e+12;
21+
w = 2*pi*f;
22+
23+
figure(1);
24+
title('Complex Refractive Index');
25+
yyaxis left
26+
plot(fTHz,n,'Linewidth',0.9);
27+
xlabel('Frequency(THz)');
28+
ylabel('Real Refractive Index');
29+
yyaxis right
30+
plot(fTHz,k,'Linewidth',0.9);
31+
ylabel('Imag Refractive Index');
32+
%% Complex Conductivity from Classical EM Wave Theory
33+
ReCond_EM = 2*n.*k.*w*epi0;
34+
ImCond_EM = epi0*w.*(epi_inf-n.^2+k.^2);
35+
ComplexCond_EM = complex(ReCond_EM,ImCond_EM);
36+
37+
figure(2);
38+
title('Complex Conductivity');
39+
yyaxis left
40+
plot(fTHz,ReCond_EM,'Linewidth',0.9);
41+
xlabel('Frequency(THz)');
42+
ylabel('Real Conductivity');
43+
yyaxis right
44+
plot(fTHz,ImCond_EM,'Linewidth',0.9);
45+
ylabel('Imag Conductivity');
46+
%% Drude Model Fitting
47+
iniguess = [3e+14,10e-15];
48+
% Initial Guess for plasma freq
49+
g = @(x) sum(abs(CondDrude(epi0,w,x(1),x(2)) - ComplexCond_EM));
50+
[result,residue] = fminsearch(g,iniguess);
51+
wp = result(1);
52+
tau = result(2);
53+
fprintf('Plasma Frequency = %g (rad*THz)\n',wp/1e+12);
54+
fprintf('Collision Time = %g (fs)\n',tau*1e+15);
55+
%% Comparison between Experimental Result and Drude Model Fitting
56+
ReCondFit = real(CondDrude(epi0,w,wp,tau));
57+
ImCondFit = imag(CondDrude(epi0,w,wp,tau));
58+
CondFit = complex(ReCondFit,ImCondFit);
59+
60+
figure(3);
61+
sgtitle('Comparsion of Experiment and Drude Fitting');
62+
subplot(2,1,1);
63+
plot(fTHz,ReCond_EM,'Linewidth',0.9);
64+
xlabel('Frequency(THz)');
65+
ylabel('Real Conductivity');
66+
hold on
67+
plot(fTHz,ReCondFit,'o');
68+
hold off
69+
legend('Experimental','Drude Fitting');
70+
subplot(2,1,2);
71+
plot(fTHz,ImCond_EM,'Linewidth',0.9);
72+
xlabel('Frequency(THz)');
73+
ylabel('Imag Conductivity');
74+
hold on
75+
plot(fTHz,ImCondFit,'o');
76+
hold off
77+
legend('Experimental','Drude Fitting');
78+
%% Error Contour
79+
% wp_err = linspace(0.5*wp,1.5*wp,1000);
80+
% tau_err = linspace(0.5*tau,1.5*tau,1000);
81+
% err = zeros(1000);
82+
%
83+
% for a = 1:1000
84+
% for b = 1:1000
85+
% err(a,b) = abs(g([wp_err(a),tau_err(b)]));
86+
% end
87+
% end
88+
%
89+
% figure(4);
90+
% contourf(wp_err,tau_err,err,50);
91+
% title('Residue Contour of Drude Fitting');
92+
% xlabel('Variation of Plasma Frequency');
93+
% ylabel('Variation of Collision Time');
94+
% colormap (jet(51))
95+
% colorbar
96+
%% Derivation of Carrier Concentration, Mobility and DC Conductivity
97+
q = 1.6e-19;
98+
m0 = 9.1e-31;
99+
effm = 0.063*m0;
100+
Ne = epi0*wp^2*effm/q^2;
101+
mobility = q*tau/effm;
102+
DCCond = epi0*wp^2*tau;
103+
fprintf('Carrier Concentration = %g (cm^-3)\n',Ne/1e+6);
104+
fprintf('Mobility = %g (cm^2/V/s)\n',mobility*1e+4);
105+
fprintf('DC Conductivity = %g (S/cm)\n',DCCond/100);
106+
%% Generation of Comparison Plots of Complex RIX
107+
epi = epi_inf + 1i*CondFit./w/epi0;
108+
n_pre = real(sqrt(epi));
109+
k_pre = imag(sqrt(epi));
110+
111+
figure(4);
112+
sgtitle('Comparsion of Refractive Index');
113+
subplot(2,1,1);
114+
plot(fTHz,n,'o');
115+
xlabel('Frequency(THz)');
116+
ylabel('n');
117+
title('Real Refractive Index');
118+
hold on
119+
plot(fTHz,n_pre,'Linewidth',0.9);
120+
legend('Experimental','Predicted');
121+
hold off
122+
subplot(2,1,2);
123+
plot(fTHz,k,'o');
124+
xlabel('Frequency(THz)');
125+
ylabel('k');
126+
title('Imag Refractive Index');
127+
hold on
128+
plot(fTHz,k_pre,'Linewidth',0.9);
129+
legend('Experimental','Predicted');
130+
hold off
131+
fclose('all');
132+
%% Function Body (Do Not Modify)
133+
function CompCond = CondDrude(epi0,w,wp,tau)
134+
numer = epi0*wp^2*tau;
135+
denom = 1-1i*w*tau;
136+
CompCond = numer./denom;
137+
end

‎TDTS_SingleLayer_Analytic.m

+105
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
%%
2+
% Analytical Extraction of Optical Constant in TDTS
3+
% Single layered, neglecting multi-reflection
4+
% Phaseshift corrected on 2019/12/09
5+
% If the refractive index is negative, modify '-' to '+' in line 76
6+
% Revised on 2020/04/10
7+
%% Analysis Parameters
8+
clear all;
9+
Dref = fscanf(fopen('Purged-1-1024-10um-cut.txt','r'),'%f %f',[2,inf]);
10+
Dsam = fscanf(fopen('Reference-SIGaAs-1-1024-10um-cut.txt','r'),'%f %f',[2,inf]);
11+
fileID = fopen('RIX-SIGaAs.txt','w'); % Output file for complex RIX
12+
d = 600e-6; % in m
13+
%% Frequency Domain Analysis
14+
t = Dref(1,:)*1e-12;
15+
N = length(t);
16+
dt = t(2)-t(1);
17+
fs = 1/dt;
18+
f = linspace(-0.5*fs,0.5*fs,N);
19+
% f = f(f>0);
20+
w = f*2*pi;
21+
fTHz = f/1e+12;
22+
Eref = Dref(2,:);
23+
Esam = Dsam(2,:);
24+
Fref = fftshift(fft(Eref));
25+
Fsam = fftshift(fft(Esam));
26+
%% Spline Interpolation to Smoothen Spectrums
27+
% fTHz2 = fTHz + 0.5*(fTHz(2)-fTHz(1));
28+
% Fref_db = db(abs(Fref).^2);
29+
% Fsam_db = db(abs(Fsam).^2);
30+
% Fref_db_spline = spline(fTHz,Fref_db,fTHz2);
31+
% Fsam_db_spline = spline(fTHz,Fsam_db,fTHz2);
32+
% figure(5);
33+
% plot(fTHz,Fref_db,fTHz2,Fref_db_spline);
34+
% xlabel('Frequnecy(THz)');
35+
% title('Spectrums');
36+
% axis([0,6.5,-inf,inf]);
37+
%% Extract Modulus and Phase (unwrapped)
38+
H = Fsam./Fref;
39+
A = abs(H);
40+
phaseshift = unwrap(angle(H));
41+
%% Extrapolation to correct phase
42+
ReliableRange = and(fTHz>0.2,fTHz<1.2);
43+
p1 = find(fTHz>=0.5,1);
44+
p2 = find(fTHz>=1.0,1);
45+
m = (phaseshift(p2)-phaseshift(p1))/(fTHz(p2)-fTHz(p1));
46+
b = phaseshift(p1) - m*fTHz(p1);
47+
for i=1:N
48+
phaseshift(i) = phaseshift(i) - b;
49+
end
50+
%% Plots for Verification
51+
figure(1);
52+
plot(t*1e+12,Eref,t*1e+12,Esam,'Linewidth',0.9);
53+
legend('Reference Waveform','Sampled Waveform');
54+
xlabel('Time Delay(ps)');
55+
ylabel('Amplitude');
56+
grid on;
57+
title('Time Domain Plot');
58+
figure(2);
59+
subplot(2,1,1);
60+
plot(fTHz,db(abs(Fref).^2),fTHz,db(abs(Fsam).^2),'Linewidth',0.9);
61+
xlabel('Frequency(THz)');
62+
ylabel('dB');
63+
grid on;
64+
title('Reference and Sample Spectrum');
65+
legend('Reference','Sample');
66+
axis([0,7,-inf,inf]);
67+
subplot(2,1,2);
68+
plot(fTHz,phaseshift,'Linewidth',0.9);
69+
xlabel('Frequency');
70+
ylabel('Phase');
71+
grid on;
72+
axis([-0.1,4,-inf,inf]);
73+
title('Phase Shift');
74+
%% Analytical Expression of Complex Refractive Index
75+
c0 = 299792458; % in m/s
76+
n = 1-c0./w.*(phaseshift/d);
77+
alpha = log((4*n)./A./(1+n).^2)*(2/d); % in m^-1
78+
alpha = alpha/100;
79+
k = c0*alpha./(4*pi*f);
80+
%% Plots of Complex Refractive Index
81+
figure(3);
82+
sgtitle('Complex Refractive Index');
83+
subplot(2,1,1);
84+
plot(fTHz,n,'Linewidth',0.8);
85+
xlabel('Frequency(THz)');
86+
ylabel('n');
87+
title('Real Refractive Index');
88+
axis([0.2,2.2,3,4]);
89+
subplot(2,1,2);
90+
plot(fTHz,k,'Linewidth',0.8);
91+
xlabel('Frequency(THz)');
92+
ylabel('\kappa');
93+
title('Imaginary Refractive Index');
94+
axis([0.3,2,-2e-3,2e-3]);
95+
%% Write out to txt File
96+
fTHz = fTHz(ReliableRange);
97+
f = f(ReliableRange);
98+
n = n(ReliableRange);
99+
k = k(ReliableRange);
100+
O(1,:) = fTHz;
101+
O(2,:) = n;
102+
O(3,:) = k;
103+
formatSpec = '%g %g %g\n';
104+
fprintf(fileID,formatSpec,O);
105+
fclose(fileID);

‎Terahertz_Waveform_Processing_Ver2.m

+68
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
%%
2+
% To Process terahertz waveform and smoothen the tails
3+
% Version 2.0 Support pre-cutting also
4+
%%
5+
[filename,path] = uigetfile('*.txt','Select a THz Waveform with Echoes');
6+
D = fscanf(fopen(filename,'r'),'%f %f',[2,inf]);
7+
t = D(1,:);
8+
E = D(2,:);
9+
N = length(t);
10+
11+
figure(1);
12+
subplot(2,1,1);
13+
plot(t,E,'Linewidth',0.9);
14+
xlabel('Time');
15+
ylabel('Amplitude');
16+
title('Waveform before Processing');
17+
subplot(2,1,2);
18+
plot(E,'Linewidth',0.9);
19+
grid on
20+
xlabel('Point #');
21+
ylabel('Amplitude');
22+
23+
Prompt1 = 'Enter the start of the THz peak ';
24+
p1 = input(Prompt1);
25+
Prompt2 = 'Enter the end of the THz peak ';
26+
p2 = input(Prompt2);
27+
28+
for k=1:N
29+
if (k<=p1 | k>=p2)
30+
E(k) = 0;
31+
end
32+
end
33+
34+
figure(2);
35+
subplot(2,1,1);
36+
plot(t,E,'Linewidth',0.9);
37+
xlabel('Time');
38+
ylabel('Amplitude');
39+
title('Waveform after Processing');
40+
subplot(2,1,2);
41+
plot(E,'Linewidth',0.9);
42+
xlabel('Point #');
43+
ylabel('Amplitude');
44+
45+
D2 = zeros(2,N);
46+
for i=1:N
47+
D2(1,i) = t(i);
48+
D2(2,i) = E(i);
49+
end
50+
%%
51+
fs = 1/(t(2)-t(1))*1e+12;
52+
f = linspace(-0.5*fs,0.5*fs,N);
53+
F = fftshift(fft(E));
54+
fTHz = f/1e+12;
55+
56+
figure(3);
57+
plot(fTHz,db(abs(F).^2),'Linewidth',0.9);
58+
xlabel('Frequency(THz)');
59+
ylabel('Spectral Power(dB)');
60+
title('Spectrum after Cut');
61+
axis([0,5,-inf,inf]);
62+
%%
63+
filename = strsplit(filename,'.');
64+
filename = string(append(filename(1),'-cut.',filename(2)));
65+
fileID = fopen(filename,'w');
66+
formatSpec = '%g %g\n';
67+
fprintf(fileID,formatSpec,D2);
68+
fclose(fileID);

0 commit comments

Comments
 (0)