-
Notifications
You must be signed in to change notification settings - Fork 18
/
CELES_MAIN.m
152 lines (133 loc) · 7.39 KB
/
CELES_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
% =========================================================================
%> @brief Run this script to start the simulation
% =========================================================================
%% add all folders to the MATLAB search path
addpath(genpath('./src'))
%% import example file with sphere positions, radii and refractive indices
data = dlmread('examples/sphere_parameters.txt');
%% initialize the CELES class instances
% initialize particle class instance
% - positionArray: Nx3 array (float) in [x,y,z] format
% - refractiveIndexArray: Nx1 array (complex) of complex refractive indices
% - radiusArray: Nx1 array (float) of sphere radii
particles = celes_particles('positionArray', data(:,2:4), ...
'refractiveIndexArray', data(:,5)+1i*data(:,6), ...
'radiusArray', data(:,1) ...
);
% initialize initial field class instance
% - polarAngle: scalar (float) polar angle of incoming beam/wave,
% in radians. for Gaussian beams, only 0 or pi are
% currently possible
% - azimuthalAngle: scalar (float) azimuthal angle of incoming
% beam/wave, in radians
% - polarization: string (char) polarization of incoming beam/wave
% ('TE' or 'TM')
% - beamWidth: scalar (float) width of beam waist. use 0 or inf
% for plane wave
% - focalPoint: 1x3 array (float) focal point
initialField = celes_initialField('polarAngle', 0, ...
'azimuthalAngle', 0, ...
'polarization', 'TE', ...
'beamWidth', 2000, ...
'focalPoint', [0,0,0] ...
);
% initialize input class instance
% - wavelength: scalar (float) vacuum wavelength, same unit as
% particle positions and radii
% - mediumRefractiveIndex: scalar (complex) refractive index of environment
% - particles: valid instance of celes_particles class
% - initialField: valid instance of celes_initialField class
input = celes_input('wavelength', 550, ...
'mediumRefractiveIndex', 1, ...
'particles', particles, ...
'initialField', initialField ...
);
% initialize preconditioner class instance
% - type: string (char) type of preconditioner (currently
% only 'blockdiagonal' and 'none' available)
% - partitionEdgeSizes 1x3 array (float) edge size of partitioning cuboids
% (applies to 'blockdiagonal' type only)
precnd = celes_preconditioner('type', 'blockdiagonal', ...
'partitionEdgeSizes', [1200,1200,1200] ...
);
% initialize solver class instance
% - type: string (char) solver type (currently 'BiCGStab' or
% 'GMRES' are implemented)
% - tolerance: scalar (float) target relative accuracy of solution
% - maxIter: scalar (int) maximum number of iterations allowed
% - restart: scalar (int) restart parameter (applies only to
% GMRES solver)
% - preconditioner: valid instance of celes_preconditioner class
solver = celes_solver('type', 'GMRES', ...
'tolerance', 5e-4, ...
'maxIter', 1000, ...
'restart', 200, ...
'preconditioner', precnd);
% initialize numerics class instance
% - lmax: scalar (int) maximal expansion order of scattered
% fields around particle center
% - polarAnglesArray: 1xN array (float) sampling of polar angles in the
% plane wave patterns, in radians
% - azimuthalAnglesArray: sampling of azimuthal angles in the plane wave
% patterns, in radians
% - gpuFlag: scalar (bool) set to false if you experience GPU
% memory problems at evaluation time (translation
% operator always runs on GPU, though)
% - particleDistanceResolution: scalar (float) resolution of lookup table for
% spherical Hankel function (same unit as wavelength)
% - solver: valid instance of celes_solver class
numerics = celes_numerics('lmax', 3, ...
'polarAnglesArray', 0:pi/5e3:pi, ...
'azimuthalAnglesArray', 0:pi/1e2:2*pi, ...
'gpuFlag', true, ...
'particleDistanceResolution', 1, ...
'solver', solver);
% define a grid of points where the field will be evaulated
[x,z] = meshgrid(-4000:40:4000, -3000:40:5000); y = zeros(size(x));
% initialize output class instance
% - fieldPoints: Nx3 array (float) points where to evaluate the
% electric near field
% - fieldPointsArrayDims: 1x2 array (int) dimensions of the array, needed to
% recompose the computed field as a n-by-m image
output = celes_output('fieldPoints', [x(:),y(:),z(:)], ...
'fieldPointsArrayDims', size(x));
% initialize simulation class instance
% - input: valid instance of celes_input class
% - numerics: valid instance of celes_input class
% - output: valid instance of celes_output class
simul = celes_simulation('input', input, ...
'numerics', numerics, ...
'output', output);
%% run simulation
simul.run;
% evaluate transmitted and reflected power
simul.evaluatePower;
fprintf('transmitted power: \t%.4f %%\n', ...
simul.output.totalFieldForwardPower/simul.output.initialFieldPower*100)
fprintf('reflected power: \t%.4f %%\n', ...
simul.output.totalFieldBackwardPower/simul.output.initialFieldPower*100)
% evaluate field at output.fieldPoints
simul.evaluateFields;
%% plot results
% display particles
figure('Name','Particle positions','NumberTitle','off');
plot_spheres(gca,simul.input.particles.positionArray, ...
simul.input.particles.radiusArray, ...
simul.input.particles.refractiveIndexArray)
% plot near field
comp = ["real Ex", "real Ey", "real Ez", "abs E", "real Hx", "real Hy", "real Hz", "abs H"];
fieldtype = 'total field';
figure('Name','Near-field cross-cut','NumberTitle','off');
for sp = 1:numel(comp)
subplot(2,4,sp)
plot_field(gca, simul, comp(sp), fieldtype);
caxis([-2*contains(comp(sp),'real'), 2])
end
linkaxes(findall(gcf,'type','axes'))
% plot Poynting vector
figure('Name','Pynting vector cross-cut','NumberTitle','off');
S = plot_poynting(gca,simul,'Total field');
ylim([min(z(:)), max(z(:))])
% % export animated gif
% figure('Name','Animated near-field cross-cut','NumberTitle','off');
% plot_field(gca,simul,'real Ey','Total field','Ey_total.gif')