-
Notifications
You must be signed in to change notification settings - Fork 29
/
calc_dispersion_stepindex2.m
144 lines (118 loc) · 4.96 KB
/
calc_dispersion_stepindex2.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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This script uses the propagation constants obtained from the mode
% calculations to approximate the dispersion parameters of each mode
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set parameters
Nf = 40;
lambda0_build = 1030e-9; % m, the center wavelength used to solve for the modes
lambda0_disp = 1030e-9; % m, the center wavelength at which the dispersion parameters will be evaluated
lrange = 180e-9; % wavelength range, in m. Should be the same as the range used to calculate the modes
modes_list=[1:10]; % usually just 1:num_modes, but it can be different
polynomial_fit_order = 9;
num_disp_orders = 4; % i.e. if this is 3, 4 coefficients will be calculated, including the 0th order
radius = 7; % Only needed for the file name
folder_name = 'Fibers/STEP_1030S';
%% Load in the calculated effective n values
if ispc
sep_char = '/';
else
sep_char = '\';
end
% Set the range in frequency space, which is more objective
c = 2.99792458e-4; % speed of ligth m/ps
if lrange == 0
error('Cannot calculate dispersion with only one frequency point');
else
f0 = c/lambda0_build; % center frequency in THz
frange = c/lambda0_build^2*lrange;
df = frange/Nf;
f = (f0 + (-Nf/2:Nf/2-1)*df)'; % frequencies in THz, as a column vector
l = c./f*10^6; % lambda in um
end
num_modes = length(modes_list);
field = 'scalar';
boundary = '0000';
lambda0_build = lambda0_build*10^6; % um
n_calc=zeros(Nf, num_modes);
for kk = 1:Nf
lambda=l(kk);
for ii=1:num_modes
fname=[folder_name sep_char 'radius' num2str(round(radius)) 'boundary' boundary 'field' field 'mode' num2str(modes_list(ii)) 'wavelength' num2str(round(lambda*1000))];
load([fname '.mat'])
n_calc(kk, ii)=neff;
end
fprintf('Loading lambda = %d um\n', round(lambda*1000));
end
%% Calculate the propagation constants
beta_calc = zeros(Nf, num_modes);
w=2*pi*f; % angular frequencies in 1/ps
for midx=1:num_modes
beta_calc(:, midx) = n_calc(:, midx).*w/c; % beta in 1/m
end
%% Fit the propagation constants to a polynomial and save the appropriate derivatives
dw = 2*pi*df;
w_disp = 2*pi*c/lambda0_disp; % angular frequency at which dispersion is calculated, in 1/ps
b_coefficients = zeros(num_modes, num_disp_orders+1); % The dispersion coefficients
for midx = 1:num_modes
beta_calc_i = beta_calc(:, midx);
beta_fit_last = polyfit(w, beta_calc_i, polynomial_fit_order); % the fit coefficients
b_coefficients(midx, 1) = polyval(beta_fit_last, w_disp)/1000; % Save beta_0 in 1/mm
for disp_order = 1:num_disp_orders
% The derivatives can be calculated exactly from the coefficients
beta_fit_last = ((polynomial_fit_order-(disp_order-1)):-1:1).*beta_fit_last(1:(polynomial_fit_order-(disp_order-1)));
b_coefficients(midx, disp_order+1) = polyval(beta_fit_last, w_disp)*(10^3)^disp_order/1000; % beta_n in fs^n/mm
end
end
% beta0 and beta1 should be relative to the fundamental mode.
b_coefficients(:, 1) = b_coefficients(:, 1) - ones(num_modes, 1)*b_coefficients(1, 1);
b_coefficients(:, 2) = b_coefficients(:, 2) - ones(num_modes, 1)*b_coefficients(1, 2);
betas = b_coefficients';
save([folder_name sep_char 'betas'], 'betas');
%% Display the results
% We need to use cell arrays because the higher orders are calculated from
% finite differences. This means that each order has one less data point
% than the previous one.
w_vectors = cell(num_disp_orders+1, 1); % omegas, in 1/ps
l_vectors = cell(num_disp_orders+1, 1); % lambdas, in um
w_vectors{1} = w;
l_vectors{1} = 2*pi*c./w_vectors{1}*10^6;
for disp_order = 1:num_disp_orders
w_prev = w_vectors{disp_order};
w_vectors{disp_order+1} = dw/2 + w_prev(1:length(w_prev)-1); % in 1/ps
l_vectors{disp_order+1} = 2*pi*c./w_vectors{disp_order+1}*10^6; % in um
end
% beta_full will have all of the orders, for each mode, as a function of
% wavlength
beta_full = cell(num_disp_orders+1, 1);
beta_full{1} = beta_calc/1000;
for disp_order = 1:num_disp_orders
beta_full{disp_order+1} = zeros(Nf-disp_order, num_modes);
end
% Take the differences to calculate the higher orders
for midx = 1:num_modes
for disp_order = 1:num_disp_orders
beta_full{disp_order+1}(:, midx) = diff(beta_full{disp_order}(:, midx))/dw*1000;
end
end
ggg=figure;
coo=hsv(num_modes);
ylabels = cell(num_disp_orders+1, 1);
ylabels{1} = '1/mm';
ylabels{2} = 'fs/mm';
for disp_order = 2:num_disp_orders
ylabels{disp_order+1} = ['fs^' num2str(disp_order) '/mm'];
end
% Plot the results
for disp_order = 1:num_disp_orders+1
subplot(1,num_disp_orders+1,disp_order)
hold on
for midx = 1:num_modes
plot(l_vectors{disp_order}, beta_full{disp_order}(:, midx), 'Color', coo(midx,:))
end
hold off
ylabel(ylabels{disp_order})
xlabel('\mum')
axis tight
end