-
Notifications
You must be signed in to change notification settings - Fork 11
/
nelly_main.m
178 lines (154 loc) · 7.24 KB
/
nelly_main.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
function [freq, n_fit, freq_full, tf_full, tf_spec, tf_pred, func, spec_smp, spec_ref]...
= nelly_main(input, t_smp, A_smp, t_ref, A_ref, varargin)
% NELLY_MAIN Runs Nelly data processing code
%
% Arguments: input -- gives input geometry and other settings for the
% calculation. This can either be a filename for a
% JSON file or a struct
% t_smp -- time points for the reference time domain trace
% A_smp -- amplitude points corresponding to t_smp
% t_ref -- time points for the sample time domain trace
% A_ref -- amplitude points corresponding to t_ref
% varargin -- additional arguments
% Output: freq -- an array containing the frequencies (THz) at which
% the refractive index was calculated
% n_fit -- an array of complex values for the refractive index.
% The ith element corresponds to the ith element in
% freq. For the imaginary part, positive values
% correspond to loss.
% freq_full -- an array containing a finer mesh of frequency points
% directly from the padded fourier transform.
% tf_full -- an array containing the transfer function
% (E_smp/E_ref). The ith element corresponds to the
% ith element in freq_full
% tf_spec -- an array containing the transfer function
% (E_smp/E_ref) at a coarser spacing. The ith element
% corresponds to the ith element of freq
% tf_pred -- an array containing the predicted transfer function
% based on the extracted refractive index values, for
% use in assessing the accuracy of the extraction. The
% ith element corresponds to the ith element of freq
% func -- an anonymous function which takes two arguments --
% frequency (THz) and the value for the unknown
% refractive index--and returns the predicted transfer
% function values at that frequency assuming that the
% unknown refractive index is the value given.
% spec_smp -- the spectrum for the sample pulse (i.e. E_smp(freq))
% The ith element corresponds to the ith element of
% freq.
% spec_ref -- the spectrum for the reference pulse (i.e.
% E_ref(freq)). The ith element corresponds to the ith
% element of freq.
%% error checking
%% load and process data and input parameters
input = load_input(input);
% process additional arguments
assert(mod(numel(varargin),2) == 0, 'Extra parameters come in pairs')
extra_args = struct;
for ii = 1:2:numel(varargin)
extra_args.(upper(varargin{ii})) = varargin{ii+1};
end
% fourier transform
[freq, tf_spec, freq_full, tf_full, spec_smp, spec_ref] = exp_tf(t_smp, A_smp, t_ref, A_ref, input);
fft_sets = input.settings.fft;
% determine time cut off for etalons (relative to peak)
% (for non tree version)
t_cut_exp = t_smp(end) - t_smp(find(A_smp == max(A_smp),1));
switch fft_sets.windowing_type
case 'gauss'
t_cut_wind = 3*fft_sets.windowing_width;
case 'square'
t_cut_wind = fft_sets.windowing_width/2;
otherwise
t_cut_wind = Inf;
end
t_cut = min([t_cut_exp t_cut_wind]);
% calculate t_cut for tree
end_time = max([t_smp(:); t_ref(:)]);
t_cut_exp_ref = end_time - t_ref(find(A_ref == max(A_ref),1));
t_cut_exp_ref = min([t_cut_exp_ref t_cut_wind]);
% estimate of time it takes to pass through reference geometry
ref_ns = arrayfun(@(x) x.n_func(1), input.reference);
ref_ds = arrayfun(@(x) x.d, input.reference);
ref_nd = dot(ref_ns, ref_ds);
t_traverse_reference = (1e12*ref_nd/3e14);
t_cut_tree = t_traverse_reference + t_cut_exp_ref;
% estimate starting refractive index
% real part
delay = t_smp(find(A_smp == max(A_smp),1)) - t_ref(find(A_ref == max(A_ref),1));
n_est = estimate_n(delay, input);
if real(n_est) < 1
warning('Estimated refractive index is less than 1')
n_est = 1;
end
% imaginary part
k_mean = mean([min(freq) max(freq)])*2*pi*1e12/3e14;
d_inds = find(strcmp({input.sample.n}, 'solve'));
d_tot = sum(arrayfun(@(ii) input.sample(ii).d, d_inds));
n_prev = [real(n_est) log(mean(abs(tf_spec)))/(d_tot*k_mean)];
a_cut = input.settings.a_cut;
% build transfer functions
func_smp = build_transfer_function_tree(input.sample, t_cut_tree, a_cut,...
varargin{:});
func_ref = build_transfer_function_tree(input.reference, t_cut_tree, a_cut,...
varargin{:});
func = @(freq, n_solve) func_smp(freq, n_solve)/func_ref(freq, n_solve);
%% perform fitting
n_fit = zeros(2, numel(freq));
for ii = 1:numel(freq)
err = @(n) n_error(func(freq(ii), complex(n(1), n(2))), tf_spec(ii));
opts = optimset();
if isfield(extra_args, 'SHOWOPTIMIZATION')
if extra_args.SHOWOPTIMIZATION
opts = optimset('PlotFcns',@optimplotfval);
end
end
if isfield(extra_args, 'PROGRESS_BAR')
extra_args.PROGRESS_BAR.Value = ii/numel(freq);
pause(0.1)
end
% Since fminsearch creates the initial simplex as a factor of the
% starting point, it can get stuck at values near zero. This can be an
% issue for the imaginary part of the refractive index, which is often
% near 0. The k_min option sets a floor value to avoid this (i.e.
% the optimization will start at (n_prev, k_min) instead of (n_prev,
% k_prev) if k_prev < k_min).
if isfield(extra_args, 'K_MIN')
if abs(n_prev(2)) < extra_args.K_MIN
n_prev(2) = sign(n_prev(2))*extra_args.K_MIN;
end
end
% By default, the initial simplex for the fminsearch Nelder-Mead
% implementation is formed by adding
% 5% to each component of the starting vector (i.e. (x,y), (1.05x, y),
% (x, 1.05y). The simplex scale options changes the size of the
% initial simplex and centers it on the optimal point for the
% previous previous frequency, i.e. the starting point is (x, y)
% such that the centroid of (x,y), ((1+s)*x, y), (x, (1+s)*y) is
% the previous optimal point.
if isfield(extra_args, 'SIMPLEX_SCALE')
n_start = 3*n_prev/(3+extra_args.SIMPLEX_SCALE);
err = @(dn) n_error(func(freq(ii),...
complex(n_start(1)+dn(1), n_start(2)+dn(2))), tf_spec(ii));
dn_opt = fminsearch(err, ones(1,2)*extra_args.SIMPLEX_SCALE,...
opts);
n_opt = n_start + dn_opt;
else
n_opt = fminsearch(err, n_prev, opts);
end
n_prev = n_opt;
n_fit(:,ii) = n_opt;
fprintf('%0.2f THz: n = %s\n', freq(ii), num2str(complex(n_opt(1), n_opt(2))))
end
if isfield(extra_args, 'PROGRESS_BAR')
close(extra_args.PROGRESS_BAR)
end
tf_pred = arrayfun(@(ii) func(freq(ii), complex(n_fit(1,ii), n_fit(2, ii))), 1:numel(freq));
n_fit = n_fit(1,:) - 1i*n_fit(2,:);
end
function [chi] = n_error(t1, t2)
chi1 = (log(abs(t1)) - log(abs(t2)))^2;
d = angle(t1) - angle(t2);
chi2 = (mod(d + pi, 2*pi) - pi)^2;
chi = chi1+chi2;
end