forked from complexphoton/MESTI.m
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmesti2s.m
2032 lines (1903 loc) · 104 KB
/
mesti2s.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
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function [S, channels, info] = mesti2s(syst, in, out, opts)
%MESTI2S Solves frequency-domain scattering problems in a two-sided geometry.
% [field_profiles, channels, info] = MESTI2S(syst, in) returns the spatial
% field profiles of Ez(x,y) for scattering problems of 2D transverse-magnetic
% (TM) fields satisfying
% [- (d/dx)^2 - (d/dy)^2 - (omega/c)^2*epsilon(x,y)] Ez(x,y) = 0,
% or of Hz(x,y) for 2D transverse-electric (TE) fields satisfying
% [- (d/dx)*(1/epsilon(x,y))_yy*(d/dx) - (d/dy)*(1/epsilon(x,y))_xx*(d/dy)
% + (d/dy)*(1/epsilon(x,y))_xy*(d/dx) + (d/dx)*(1/epsilon(x,y))_yx*(d/dy)
% - (omega/c)^2] Hz(x,y) = 0,
% where Ez or Hz is each a superposition of incident and scattered waves.
% The polarization (TM or TE), relative permittivity profile epsilon(x,y),
% frequency omega, and boundary conditions are specified by structure 'syst';
% epsilon(x,y) must have homogeneous spaces on the left (-x) and right (+x)
% sides, with an outgoing boundary in x for the scattered waves and a closed
% boundary in y.
% The incident wavefronts from left and/or right are specified by variable
% 'in'.
% The returned 'field_profiles' is a 3D array, with field_profiles(:,:,i)
% being the total (incident + scattered) field profile of Ez or Hz given the
% i-th incident wavefront. The returned 'channels' is a structure containing
% properties of the propagating and evanescent channels in the homogeneous
% spaces on the left and right. The information of the computation is returned
% in structure 'info'.
%
% [S, channels, info] = MESTI2S(syst, in, out) returns the scattering matrix
% S, where 'in' and 'out' specify either the list of input/output channels or
% the input/output wavefronts. When the MUMPS function zmumps() is available,
% this is typically done by computing the Schur complement of an augmented
% matrix K through a partial factorization.
%
% [field_profiles, channels, info] = MESTI2S(syst, in, [], opts) and
% [S, channels, info] = MESTI2S(syst, in, out, opts) allow detailed options to
% be specified with structure 'opts'.
%
% In mesti2s(), the boundary condition in y must be closed (e.g., periodic or
% PEC). Given the closed boundary, the set of transverse modes forms a
% complete and orthonormal basis of propagating and evanescent channels.
% The inputs and outputs are specified in the basis of these propagating
% channels, with coefficients normalized with respect to the flux in the
% longitudinal (x) direction. Properties of those channels are given by
% mesti_build_channels().
%
% When an open boundary in y is of interest, the appropriate input/output
% basis depends on the problem, so the user should use the more general
% function mesti() and will need to specify the input and output matrices B
% and C.
%
% MESTI considers nonmagnetic materials with isotropic (i.e., scalar)
% permittivity. Even though epsilon(x,y) is originally a scalar in the
% continuous system, subpixel smoothing (which is needed to ensure smooth
% variation with respect to parameters and to reach second-order accuracy for
% TE polarization) creates anisotropy where epsilon(x,y) becomes a symmetric
% tensor. TM polarization uses the zz component of epsilon(x,y). TE
% polarization uses the xx, yy, and xy (equals yx) components of epsilon(x,y).
% The other components of epsilon(x,y) are zero.
%
% This file builds the input and output channels using mesti_build_channels(),
% builds the nonzero blocks of matrices B and C, and then calls mesti() to
% solve the scattering problems. For scattering-matrix computations in TM
% polarization, mesti2s() can also use the recursive Green's function method
% instead, which is efficient in 1D and for 2D systems with a small width in y.
%
% === Input Arguments ===
% syst (scalar structure; required):
% A structure that specifies the system, used to build the FDFD matrix A.
% It contains the following fields:
% syst.polarization (character vector; optional):
% Polarization. Possible choices are:
% 'TM' - Transverse-magnetic field (Hx, Hy, Ez)
% 'TE' - Transverse-electric field (Ex, Ey, Hz)
% TM field uses syst.epsilon, and TE field uses syst.inv_epsilon. If
% only one of syst.epsilon and syst.inv_epsilon is given,
% syst.polarization is optional and will be automatically picked based
% on which one is given. If syst.epsilon and syst.inv_epsilon are both
% given, then syst.polarization must be specified.
% syst.epsilon (numeric matrix, real or complex; required for TM):
% An ny_Ez-by-nx_Ez matrix discretizing the relative permittivity
% profile epsilon(x,y) of the scattering region. Specifically,
% syst.epsilon(m,n) is the scalar epsilon(x,y) averaged over a square
% with area (syst.dx)^2 centered at the point (x_n, y_m) where Ez(x,y)
% is located on the Yee lattice. It is the zz component of the
% discretized epsilon(x,y) tensor from subpixel smoothing, used by TM
% fields. We choose
% (x_n, y_m) = (n-0.5, m-0.5)*syst.dx,
% with n = 1, ..., nx_Ez, m = 1, ..., ny_Ez.
% such that the lower corner of the first pixel syst.epsilon(m=1,n=1) is
% at (x,y) = (0,0), and the back side of the last slice
% syst.epsilon(m,n=nx_Ez) is at the edge of the scattering region, x = L
% = nx_Ez*syst.dx;
% Outside the scattering region (with x < 0 and x > L), epsilon(x,y)
% is given by scalars syst.epsilon_L and syst.epsilon_R for the
% homogeneous spaces on the two sides. Note that nx_Ez = 0 corresponds
% to L = 0 (ie, no scattering region) and is allowed.
% Note that y corresponds to the first index m, and x corresponds to
% the second index n.
% The positive imaginary part of syst.epsilon describes absorption,
% and the negative imaginary part describes linear gain.
% One can use syst.epsilon with ny_Ez = 1 and with a periodic or
% Bloch periodic boundary in y to simulate 1D systems where the relative
% permittivity profile is translationally invariant in y.
% syst.inv_epsilon (cell array; required for TE):
% The xx, yy, and xy components of the discretized inverse relative
% permittivity 1/epsilon(x,y) tensor of the scattering region from
% subpixel smoothing, used by TE fields. It has three elements
% inv_epsilon{1}: (1/epsilon(x,y))_xx, size [ny_Ez, nx_Hz]
% inv_epsilon{2}: (1/epsilon(x,y))_yy, size [ny_Hz, nx_Ez]
% inv_epsilon{3}: (1/epsilon(x,y))_xy, size [ny_Ez, nx_Ez]
% The third element, inv_epsilon{3}, is optional and is treated as zero
% when syst.inv_epsilon only has two elements. The yx component is not
% specified since we only consider symmetric 1/epsilon tensors where
% (1/epsilon)_yx = (1/epsilon)_xy.
% The different components are located at different points:
% - Hz(x,y) at (x_{n-0.5}, y_{m-0.5}); size [ny_Hz, nx_Hz].
% - (1/epsilon(x,y))_xx and Dx ~ dHz/dy at (x_{n+0.5}, y_m).
% - (1/epsilon(x,y))_yy and Dy ~ dHz/dx at (x_n, y_{m+0.5}).
% - (1/epsilon(x,y))_xy at (x_n, y_m), same as Ez.
% Here, (x_n, y_m) is the location of Ez and syst.epsilon above.
% inv_epsilon{1}, inv_epsilon{2}, and inv_epsilon{3} should each be
% (1/epsilon(x,y))_xx, (1/epsilon(x,y))_yy, and (1/epsilon(x,y))_xy from
% subpixel smoothing, averaged over a square with area (syst.dx)^2
% centered at the points where each of them is located. The
% 1/epsilon(x,y) tensor from subpixel smoothing is given by Eq. (1) of
% Farjadpour et al, Optics Letters 31, 2972 (2006).
% In x direction for a two-sided geometry, nx_Hz = nx_Ez + 1 =
% L/syst.dx + 1, and the sites on x_{n+0.5} start from x = 0 (i.e., half
% a pixel before the first site of x_n) and end on x = L (i.e., half a
% pixel after the last site of x_n). That means the first and last
% slices of Hz sit exactly on the front and back surfaces of the
% scattering region, and the first and last slices of inv_epsilon{1}
% extend half a pixel outside the scattering region. For a one-sided
% geometry (where syst.epsilon_R is not given), nx_Hz = nx_Ez =
% L/syst.dx, and the last slice of Hz sits on x = L - syst.dx (since Hz
% = 0 at x = L).
% In y direction, where these points start and end depend on the
% boundary condition.
% For periodic, Bloch periodic, and PMCPEC boundary conditions in y,
% ny_Hz = ny_Ez, and all of the sites on y_{n+0.5} are half a pixel
% after the corresponding sites on y_n.
% For PEC boundary condition in y, ny_Hz = ny_Ez + 1, and the sites
% on y_{n+0.5} start from half a pixel before the first site of y_n and
% end on half a pixel after the last site of y_n.
% For PMC boundary condition in y, ny_Hz = ny_Ez - 1, and the sites
% on y_{n+0.5} start from half a pixel after the first site of y_n and
% end on half a pixel before the last site of y_n.
% For PECPMC boundary condition in y, ny_Hz = ny_Ez, and all of the
% sites on y_{n+0.5} are half a pixel before the corresponding sites on
% y_n.
% For a 1D problem with ny_Hz = 1 and periodic boundary condition in
% y, the y derivative vanishes, so inv_epsilon{1} is not used and does
% not need to be specified.
% syst.epsilon_L (real scalar; required):
% Relative permittivity of the homogeneous space on the left.
% syst.epsilon_R (real scalar or []; optional):
% Relative permittivity of the homogeneous space on the right. If
% syst.epsilon_R is not given or is empty, the system will be one-sided,
% terminated on the right with a PEC boundary with Ez = 0 at x = L +
% 0.5*syst.dx, or wth a PMC boundary with Hz = 0 at x = L.
% syst.length_unit (anything; optional):
% Length unit, such as micron, nm, or some reference wavelength. This
% code only uses dimensionless quantities, so syst.length_unit is never
% used. This syst.length_unit is meant to help the user interpret the
% units of (x,y), dx, wavelength, ky_B, etc.
% syst.wavelength (numeric scalar, real or complex; required):
% Vacuum wavelength 2*pi*c/omega, in units of syst.length_unit.
% syst.dx (positive scalar; required):
% Discretization grid size, in units of syst.length_unit.
% syst.yBC (character vector; required unless syst.ky_B is specified):
% Boundary condition (BC) at the two ends in y direction, effectively
% specifying Ez(m,n) or Hz(m,n) at m=0 and m=ny_Ez+1 or ny_Hz+1, one
% pixel beyond the computation domain. Available choices are:
% 'Bloch' - Ez(m+ny_Ez,n) = Ez(m,n)*exp(1i*syst.ky_B*ny_Ez*syst.dx)
% Hz(m+ny_Hz,n) = Hz(m,n)*exp(1i*syst.ky_B*ny_Hz*syst.dx)
% 'periodic' - equivalent to 'Bloch' with syst.ky_B = 0
% 'PEC' - Ez(0,n) = Ez(ny_Ez+1,n) = 0
% Hz(0,n) = Hz(1,n); Hz(ny_Hz+1,n) = Hz(ny_Hz,n)
% 'PMC' - Ez(0,n) = Ez(1,n); Ez(ny_Ez+1,n) = Ez(ny_Ez,n)
% Hz(0,n) = Hz(ny_Hz+1,n) = 0
% 'PECPMC' - Ez(0,n) = 0; Ez(ny_Ez+1,n) = Ez(ny_Ez,n)
% Hz(0,n) = Hz(1,n); Hz(ny_Hz+1,n) = 0
% 'PMCPEC' - Ez(0,n) = Ez(1,n); Ez(ny_Ez+1,n) = 0
% Hz(0,n) = 0; Hz(ny_Hz+1,n) = Hz(ny_Hz,n)
% where PEC stands for perfect electric conductor (for which Ez = 0 and
% Ex ~ dHz/dy = 0 at the boundary) and PMC stands for perfect magnetic
% conductor (for which Hz = 0 and Hx ~ dEz/dy = 0 at the boundary).
% Note that this yBC also defines a complete and orthonormal set of
% transverse modes, upon which the input and output channels in input
% arguments 'in' and 'out' are defined; mesti2s() does not support PML
% in y direction because a closed boundary is necessary for defining
% such a transverse basis.
% Here, syst.yBC is required, with no default choice (except when
% syst.ky_B is given, in which case syst.yBC = 'Bloch' is automatically
% used).
% syst.ky_B (real scalar; optional):
% Bloch wave number in y direction, in units of 1/syst.length_unit.
% syst.ky_B is only used when syst.yBC = 'Bloch'. It is allowed to
% specify a complex-valued syst.ky_B, but a warning will be displayed.
% syst.xBC (character vector or scalar structure or cell array; optional):
% In mesti(), outgoing boundary condition is always used in x direction.
% But there are different options for implementing such outgoing
% boundary. By default,
% syst.xBC = 'outgoing';
% in which case self-energy based on the retarded Green's function of a
% semi-infinite discrete homogeneous space is used to enforce an exact
% outgoing boundary for all propagating and evanescent waves on both
% sides. Doing so doesn't increase the simulation domain; For TM, only a
% single-pixel slice of syst.epsilon_L and a single-pixel slice of
% syst.epsilon_R are added on the left and on the right to place the
% self-energy and the input and output, with a total simulation domain
% size of [ny_Ez, nx_Ez+2]. For TE, the self-energy and input/output are
% directly placed on the surface where the first and last pixels of Hz
% are, so the total simulation domain size remains [ny_Hz, nx_Hz].
% A drawback of self-energy is that it introduces coupling between
% every pair of pixels in y direction on that slice, making matrix A
% less sparse. When the system is wide (roughly when ny > 800), it is
% more efficient to implement the outgoing boundary using perfectly
% matched layer (PML), which attenuates outgoing waves with minimal
% reflection.
% Note that in mesti2s(), the PML is placed in the homogeneous spaces
% specified by syst.epsilon_L and syst.epsilon_R, outside of the
% scattering region specified by syst.epsilon or syst.inv_epsilon. (This
% is different from the more general function mesti(), where
% syst.epsilon specifies the entire simulation domain, so PML is placed
% inside syst.epsilon.)
% To use the same PML on both sides or to use PML on the left of a
% one-sided geometry, set syst.xBC to be a scalar structure with the
% following fields:
% npixels_PML (positive integer scalar; required): Number of PML
% pixels. This number of pixels is added in addition to the
% scattering region.
% npixels_spacer (non-negative integer scalar; optional): Number of
% homogeneous-space pixels to be added between the PML pixels and
% syst.epsilon/syst.inv_epsilon, used to attenuate evanescent
% waves. Defaults to 0.
% power_sigma (non-negative scalar; optional): Power of the
% polynomial grading for the conductivity sigma; defaults to 3.
% sigma_max_over_omega (non-negative scalar; optional):
% Conductivity at the end of the PML; defaults to
% 0.8*(power_sigma+1)/((2*pi/wavelength)*dx*sqrt(epsilon_bg)).
% where epsilon_bg is either syst.epsilon_L or syst.epsilon_R.
% This is used to attenuate propagating waves.
% power_kappa (non-negative scalar; optional): Power of the
% polynomial grading for the real-coordinate-stretching factor
% kappa; defaults to 3.
% kappa_max (real scalar no smaller than 1; optional):
% Real-coordinate-stretching factor at the end of the PML;
% defaults to 15. This is used to accelerate the attenuation of
% evanescent waves. kappa_max = 1 means no real-coordinate
% stretching.
% power_alpha (non-negative scalar; optional): Power of the
% polynomial grading for the CFS alpha factor; defaults to 1.
% alpha_max_over_omega (non-negative scalar; optional): Complex-
% frequency-shifting (CFS) factor at the beginning of the PML.
% This is typically used in time-domain simulations to suppress
% late-time (low-frequency) reflections. We don't use it by
% default (alpha_max_over_omega = 0) since we are in frequency
% domain.
% We use the following PML coordinate-stretching factor:
% s(p) = kappa(p) + sigma(p)./(alpha(p) - i*omega)
% with
% sigma(p)/omega = sigma_max_over_omega*(p.^power_sigma),
% kappa(p) = 1 + (kappa_max-1)*(p.^power_kappa),
% alpha(p)/omega = alpha_max_over_omega*((1-p).^power_alpha),
% where omega is frequency, and p goes linearly from 0 at the beginning
% of the PML to 1 at the end of the PML.
% The syntax above (setting syst.xBC to 'outgoing' or a scalar
% structure) uses the same parameters on the two sides. For a two-sided
% geometry (i.e., when syst.epsilon_L and syst.epsilon_R are both
% provided), one can also use different parameters on the two sides. To
% do so, set syst.xBC to a two-element cell array, with the first
% element specifying parameters on the left, the second element
% specifying parameters on the right. For example,
% syst.xBC = {'outgoing', struct('npixels_PML', 40)};
% specifies exact outgoing boundary on the left, 40 pixels of PML on the
% right.
% With real-coordinate stretching, PML can attenuate evanescent waves
% more efficiently than free space, so npixels_spacer defaults to 0.
% The PML thickness should be chosen based on the acceptable level of
% reflectivity given the discretization resolution and the range of wave
% numbers (i.e., angles) involved; more PML pixels gives lower
% reflectivity. Typically 10-40 pixels are sufficient.
% PML is a layer, not technically a boundary condition. In mesti2s(),
% the boundary condition behind the PML is always set to Dirichlet (ie,
% PEC for TM, PMC for TE).
% Note that PML cannot be used when opts.method = 'RGF'.
% in (cell array or scalar structure; required):
% The set of input channels or input wavefronts.
% To specify all propagating channels on one side or on both sides, use
% in = {'left'}, or
% in = {'right'}, or
% in = {'left', 'right'}.
% The character vectors 'left' and 'right' can be abbreviated as 'L' and
% 'R' respectively.
% To specify a subset of the propagating channels, let 'in' be a scalar
% structure with the following fields:
% in.ind_L (integer vector): Vector containing the indices of
% propagating channels incident on the left side.
% in.ind_R (integer vector): Vector containing the indices of
% propagating channels incident on the right side.
% One can provide only in.ind_L or only in.ind_R or both of them.
% The above generates flux-normalized single-channel inputs of the form
% f_a^L(x,y) = 1/sqrt(nu_L(a))*u_a(y)*exp(i*kx_L(a)*x)
% for input channel 'a' from the left, and/or
% f_a^R(x,y) = 1/sqrt(nu_R(a))*u_a(y)*exp(-i*kx_R(a)*(x-L))
% for input channel 'a' from the right, where nu normalizes flux in the x
% direction, kx is the longitudinal wave number, and u_a(y) is the
% transverse profile of the a-th propagating channel; nu = sin(kx(a)*dx)
% for TM polarization, nu = sin(kx(a)*dx)/epsilon_bg for TE polarization.
% The user can first use mesti_build_channels() to get the indices, wave
% numbers, and transverse profiles of the propagating channels; based on
% that, the user can specify the list of channels of interest through
% in.ind_L and in.ind_R above, or a list of customized wavefronts through
% in.v_L and in.v_R below.
% When the input(s) of interests are not a single channel but a
% superposition of multiple propagating channels, such custom input
% wavefronts can be specified by letting 'in' be a scalar structure with
% the following fields:
% in.v_L (numeric matrix): Matrix where each column specifies the
% coefficients of propagating channels on the left for one input
% wavefront from the left; the wavefront is a superposition of all
% propagating channels of the form f_a^L(x,y) above, with the
% superposition coefficients given by that column of in.v_L.
% size(in.v_L, 1) must equal N_prop_L, the total number of
% propagating channels on the left; size(in.v_L, 2) is the number
% of input wavefronts.
% in.v_R (numeric matrix): Analogous to to in.v_L, but specifying
% input wavefronts from the right instead.
% Note that the input wavefronts from the left and the input wavefronts
% from the right are treated as separate inputs. In other words, each input
% either comes from the left or comes from the right, but not from both
% sides. If an input with incidence from both sides is of interest, the
% user can manually superimpose results from the separate-side
% computations.
% out (cell array or scalar structure or []; optional):
% The set of output channels or output wavefronts.
% When out=[] or when out is omitted as in input argument (i.e., when
% nargin <= 2), no output projection is used, and the spatial field
% profiles Ez(x,y) or Hz(x,y) corresponding to the set of inputs are
% returned.
% When out is given, the scattering matrix is returned, with the output
% basis of the scattering matrix specified by out. In this case, out
% follows the same syntax as the input argument 'in'. Specifically, one can
% specify all propagating channels on one side or on both sides with
% out = {'left'}, or
% out = {'right'}, or
% out = {'left', 'right'}.
% The character vectors 'left' and 'right' can be abbreviated as 'L' and
% 'R' respectively.
% One can alternatively specify a subset of propagating
% channels with
% out.ind_L (integer vector), and/or
% out.ind_R (integer vector).
% For example, if r_full is the full reflection matrix computed with in =
% {'L'} and out = {'L'}, then the reflection matrix computed using in.ind_L
% and out.ind_L is equivalent to r_full(out.ind_L, in.ind_L).
% One can also let the output basis of the scattering matrix be a
% superposition of multiple propagating channels, with such custom output
% wavefronts specified by
% out.v_L (numeric matrix), and/or
% out.v_R (numeric matrix).
% Here, each row of the scattering matrix corresponds to projection onto an
% output wavefront specified by a column of out.v_L or of out.v_R. For
% example, if r_full is the full reflection matrix computed with in = {'L'}
% and out = {'L'}, then the reflection matrix computed using in.v_L and
% out.v_L is equivalent to (out.v_L)'*r_full*in.v_L where ' is the
% conjugate transpose.
% opts (scalar structure; optional, defaults to an empty struct):
% A structure that specifies the options of computation; defaults to an
% empty structure. It can contain the following fields (all optional):
% opts.verbal (logical scalar; optional, defaults to true):
% Whether to print system information and timing to the standard output.
% opts.nx_L (non-negative integer scalar; optional, defaults to 0):
% Number of pixels of homogeneous space on the left (syst.epsilon_L) to
% include when returning the spatial field profile; not used for
% scattering matrix computations.
% opts.nx_R (non-negative integer scalar; optional, defaults to 0):
% Number of pixels of homogeneous space on the right (syst.epsilon_R) to
% include when returning the spatial field profile; not used for
% scattering matrix computations. Note that opts.nx_R can still be used
% in one-sided geometries where syst.epsilon_R is not given; the field
% profile on the right is simply zero in such case.
% opts.solver (character vector; optional):
% The solver used for sparse matrix factorization. Available choices are
% (case-insensitive):
% 'MUMPS' - (default when MUMPS is available) Use MUMPS. Its MATLAB
% interface zmumps.m must be in MATLAB's search path.
% 'MATLAB' - (default when MUMPS is not available) Use the built-in
% lu() function in MATLAB, which uses UMFPACK with AMD
% ordering.
% MUMPS is faster and uses less memory than lu(), and is required for
% the APF method.
% opts.solver is not used when opts.method = 'RGF'.
% opts.method (character vector; optional):
% The solution method. Available choices are (case-insensitive):
% 'APF' - Augmented partial factorization. C*inv(A)*B is obtained
% through the Schur complement of an augmented matrix
% K = [A,B;C,0] using a partial factorization. Must have
% opts.solver = 'MUMPS'. This is the most efficient method,
% but it cannot be used for computing the full field profile
% X=inv(A)*B or with iterative refinement.
% 'FG' - Factorize and group. Factorize A=L*U, and obtain C*inv(A)*B
% through C*inv(U)*inv(L)*B with optimized grouping. Must
% have opts.solver = 'MATLAB'. This is slightly better than
% 'FS' when MUMPS is not available, but it cannot be used for
% computing the full field profile X=inv(A)*B.
% 'FS' - Factorize and solve. Factorize A=L*U, solve for X=inv(A)*B
% with forward and backward substitutions, and project with
% C as C*inv(A)*B = C*X. Here, opts.solver can be either
% 'MUMPS' or 'MATLAB', and it can be used for computing
% the full field profile X=inv(A)*B or with iterative
% refinement.
% 'RGF' - Recursive Green's function method. Cannot be used for
% computing the full field profile or with iterative
% refinement, and cannot be used with PML or with TE
% polarization.
% 'C*inv(U)*inv(L)*B' - Same as 'FG'.
% 'factorize_and_solve' - Same as 'FS'.
% By default, if the input argument 'out' is not given or if
% opts.iterative_refinement = true, then 'factorize_and_solve' is used.
% Otherwise, if ny is large or if TE polarization is used or if PML has
% been specified, then 'APF' is used when opts.solver = 'MUMPS', and
% 'C*inv(U)*inv(L)*B' is used when opts.solver = 'MATLAB'. Otherwise,
% 'RGF' is used.
% opts.symmetrize_K (logical scalar; optional):
% Whether or not to pad input and/or output channels and perform
% permutations to make matrix K = [A,B;C,0] symmetric when computing its
% Schur complement, which lowers computing time and memory usage. Such
% channel padding and permutation is reversed afterwards and does not
% affect what mesti2s() returns.
% opts.symmetrize_K can only be used when all of the following are
% met: (1) input argument 'out' is given, (2) 'in' and 'out' are not
% specified as wavefronts, (3) opts.method = 'APF', and (4) the boundary
% condition in y is not Bloch periodic. When all of these conditions are
% met, opts.symmetrize_K defaults to true; otherwise it is not used.
% opts.clear_syst (logical scalar; optional, defaults to false):
% When opts.clear_syst = true and opts.method is not 'RGF', variable
% 'syst' will be cleared in the caller's workspace to reduce peak memory
% usage. This can be used when syst.epsilon or syst.inv_epsilon takes up
% significant memory and is not needed after calling mesti2s().
% opts.clear_memory (logical scalar; optional, defaults to true):
% Whether or not to clear variables inside mesti2s() to reduce peak
% memory usage.
% opts.verbal_solver (logical scalar; optional, defaults to false):
% Whether to have the solver print detailed information to the standard
% output. Note the behavior of output from MUMPS depends on compiler.
% opts.use_METIS (logical scalar; optional, defaults to false):
% Whether to use METIS (instead of the default AMD) to compute the
% ordering in MUMPS. Using METIS can sometimes reduce memory usage
% and/or factorization and solve time, but it typically takes longer at
% the analysis (i.e., ordering) stage.
% opts.nrhs (positive integer scalar; optional):
% The number of right-hand sides (number of columns of the input matrix
% B) to consider simultaneously, used only when opts.method =
% 'factorize_and_solve' and input argument 'out' is given. Defaults to 1
% if opts.iterative_refinement = true, 10 if opts.solver = 'MUMPS' with
% opts.iterative_refinement = false, 4 otherwise.
% opts.store_ordering (logical scalar; optional, defaults to false):
% Whether to store the ordering sequence (permutation) for matrix A or
% matrix K; only possible when opts.solver = 'MUMPS'. If
% opts.store_ordering = true, the ordering will be returned in
% info.ordering.
% opts.ordering (positive integer vector; optional):
% A user-specified ordering sequence for matrix A or matrix K, used only
% when opts.solver = 'MUMPS'. Using the ordering from a previous
% computation can speed up (but does not eliminate) the analysis stage.
% The matrix size must be the same, and the sparsity structure should be
% similar among the previous and the current computation.
% opts.analysis_only (logical scalar; optional, defaults to false):
% When opts.analysis_only = true, the factorization and solution steps
% will be skipped, and S = [] will be returned. The user can use
% opts.analysis_only = true with opts.store_ordering = true to return
% the ordering for A or K; only possible when opts.solver = 'MUMPS'.
% opts.nthreads_OMP (positive integer scalar; optional):
% Number of OpenMP threads used in MUMPS; overwrites the OMP_NUM_THREADS
% environment variable.
% opts.parallel_dependency_graph (logical scalar; optional):
% If MUMPS is multithread, whether to use parallel dependency graph in MUMPS.
% This typically improve the time performance, but marginally increase
% the memory usage.
% opts.iterative_refinement (logical scalar; optional, defaults to false):
% Whether to use iterative refinement in MUMPS to lower round-off
% errors. Iterative refinement can only be used when opts.solver =
% 'MUMPS' and opts.method = 'factorize_and_solve' and input argument
% 'out' is given, in which case opts.nrhs must equal 1. When iterative
% refinement is used, the relevant information will be returned in
% info.itr_ref_nsteps, info.itr_ref_omega_1, and info.itr_ref_omega_2.
% opts.use_continuous_dispersion (logical scalar; optional):
% Whether to use the dispersion equation of the continuous wave equation
% when building the input/output channels. Defaults to false, in which
% case the finite-difference dispersion is used.
% opts.m0 (real numeric scalar; optional, defaults to 0):
% Center of the transverse mode profile with periodic or Bloch periodic
% boundary condition, u(m,a) = exp(i*ky(a)*syst.dx*(m-m0))/sqrt(ny),
% where ky(a) = syst.ky_B + a*(2*pi/ny*syst.dx) and ny = ny_Ez for TM,
% ny_Hz for TE.
%
% === Output Arguments ===
% field_profiles (3D array):
% For field-profile computations (i.e., when 'out' is not given), the
% returned field_profiles is a 3D array containing the field profiles, such
% that field_profiles(:,:,a) is the total (incident + scattered) field of
% Ez or Hz corresponding to the a-th input wavefront, with
% size(field_profiles, 1) = ny
% size(field_profiles, 2) = opts.nx_L + nx + opts.nx_R
% where [ny, nx] = [ny_Ez, nx_Ez] for TM, [ny_Hz, nx_Hz] for TE. Recall
% that nx_Hz = nx_Ez + 1 = L/syst.dx + 1 for a two-sided geometry; nx_Hz =
% nx_Ez = L/syst.dx for a one-sided geometry.
% By default, opts.nx_L = opts.nx_R = 0, and field_profiles only contain
% Ez or Hz in the scattering region. When opts.nx_L and opts.nx_R are
% specified, field_profiles will also contain opts.nx_L pixels of Ez or Hz
% in the homogeneous space on the left, opts.nx_R pixels of Ez or Hz
% in the homogeneous space on the right.
% S (numeric matrix):
% For scattering-matrix computations (i.e., when 'out' is given), S is the
% scattering matrix, such that S(b,a) is the flux-normalized field
% coefficient in the b-th propagating output channel (or the b-th output
% wavefront) given incident wave in the a-th propagating input channel (or
% the a-th input wavefront).
% When all propagating channels on one side or on both sides are
% requested, e.g. with in = {'left'} or out = {'left', 'right'}, matrix S
% includes channels.L.N_prop propagating channels on the left and/or
% channels.R.N_prop on the right, with transverse wave numbers given by
% vectors channels.L.kydx_prop and channels.R.kydx_prop, longitudinal wave
% numbers given by vectors channels.L.kxdx_prop and channels.R.kxdx_prop.
% The transverse profiles of these channels are channels.fun_u(kydx).
% Matrix S is then one block of S_full = [[r_L; t_L], [t_R; r_R]] depending
% on which side(s) are specified for the input and output, with r_L and r_R
% being reflection matrices on the left and right sides, t_L and t_R being
% transmission matrices with input from the left and right respectively.
% Note that channels on the left always come before channels on the right
% in matrix S, so the first channels.L.N_prop elements of the matrix are
% always channels on the left (if requested), with the last
% channels.R.N_prop elements being channels on the right (if requested);
% the ordering of elements within the cell arrays 'in' and/or 'out' is not
% considered.
% The phases of the elements of the scattering matrix depend on the
% reference plane. For channels on the left, the reference plane is at
% x = 0. For channels on the right, the reference plane is at x = L =
% nx_Ez*syst.dx.
% When a subset of the propagating channels are requested, with
% in.ind_L, in.ind_R, out.ind_L, and/or out.ind_R, matrix S includes such
% subset of the propagating channels. For example, if r_full is the full
% reflection matrix computed with in = {'L'} and out = {'L'}, then the
% reflection matrix computed using in.ind_L and out.ind_L is equivalent to
% r_full(out.ind_L, in.ind_L).
% When the input wavefronts and/or output basis is specified by in.v_L,
% in.v_R, out.v_L, and/or out.v_R, matrix S is the full scattering matrix
% with the input and/or output basis changed. For example, if r_full is the
% full reflection matrix computed with in = {'L'} and out = {'L'}, then the
% reflection matrix computed using in.v_L and out.v_L is equivalent to
% (out.v_L)'*r_full*in.v_L where ' is the conjugate transpose.
% channels (structure):
% A structure returned by function mesti_build_channels() that contains
% properties of the propagating and evanescent channels of the homogeneous
% spaces on the left and right. Type "help mesti_build_channels" for more
% information.
% info (scalar structure):
% A structure that contains the following fields:
% info.opts (scalar structure):
% The final 'opts' used, excluding the user-specified matrix ordering.
% info.timing (scalar structure):
% A structure containing timing of the various stages, in seconds, in
% fields 'total', 'init', 'build', 'analyze', 'factorize', 'solve'.
% info.nnz (scalar structure):
% A structure containing the number of nonzero elements for the various
% matrices, in fields 'A', 'B', 'C', 'S', 'X'.
% info.xPML (two-element cell array; optional);
% PML parameters on the two sides in x direction, if used.
% info.ordering_method (character vector; optional):
% Ordering method used in MUMPS.
% info.ordering (positive integer vector; optional):
% Ordering sequence returned by MUMPS when opts.store_ordering = true.
% info.itr_ref_nsteps (integer vector; optional):
% Number of steps of iterative refinement for each input, if
% opts.iterative_refinement = true; 0 means no iterative refinement.
% info.itr_ref_omega_1 (real vector; optional):
% Scaled residual omega_1 at the end of iterative refinement for each
% input; see MUMPS user guide section 3.3.2 for definition.
% info.itr_ref_omega_2 (real vector; optional):
% Scaled residual omega_2 at the end of iterative refinement for each
% input; see MUMPS user guide section 3.3.2 for definition.
%
% See also: mesti_build_channels, mesti
%% Part 1.1: Check validity of syst, assign default values to its fields, and parse BC and PML specifications
t0 = clock;
if nargin < 2; error('Not enough input arguments.'); end
if ~(isstruct(syst) && isscalar(syst)); error('Input argument syst must be a scalar structure.'); end
if ~isfield(syst, 'epsilon_L'); error('Input argument syst must have field ''epsilon_L''.'); end
if ~isfield(syst, 'wavelength'); error('Input argument syst must have field ''wavelength''.'); end
if ~isfield(syst, 'dx'); error('Input argument syst must have field ''dx''.'); end
if ~(isreal(syst.epsilon_L) && isscalar(syst.epsilon_L)); error('syst.epsilon_L must be a real scalar.'); end
if ~(isnumeric(syst.wavelength) && isscalar(syst.wavelength)); error('syst.wavelength must be a numeric scalar.'); end
if ~(isreal(syst.dx) && isscalar(syst.dx) && syst.dx > 0); error('syst.dx must be a positive scalar.'); end
% Pick the polarization to use; assign use_TM
if isfield(syst, 'polarization')
if strcmpi(syst.polarization, 'TM')
use_TM = true;
elseif strcmpi(syst.polarization, 'TE')
use_TM = false;
else
error('syst.polarization, if given, must be ''TM'' or ''TE''.');
end
else
% When syst.polarization is not given, we automatically pick based on whether syst.epsilon or syst.inv_epsilon is given.
if isfield(syst, 'epsilon') && ~isfield(syst, 'inv_epsilon')
use_TM = true;
elseif ~isfield(syst, 'epsilon') && isfield(syst, 'inv_epsilon')
use_TM = false;
elseif isfield(syst, 'epsilon') && isfield(syst, 'inv_epsilon')
error('syst.polarization must be given when syst.epsilon and syst.inv_epsilon both exist.');
else % neither syst.epsilon nor syst.inv_epsilon exists
error('Input argument syst must have field ''epsilon'' or ''inv_epsilon''.');
end
end
% Check syst.epsilon (for TM) and syst.inv_epsilon (for TE)
if use_TM
if ~isfield(syst, 'epsilon')
error('syst.epsilon must be given when syst.polarization = ''TM''.');
elseif ~(isnumeric(syst.epsilon) && ismatrix(syst.epsilon))
error('syst.epsilon must be a numeric matrix, if given.');
end
syst.polarization = 'TM';
str_pol = 'Ez'; % for printing system info
else
if ~isfield(syst, 'inv_epsilon')
error('syst.inv_epsilon must be given when syst.polarization = ''TE''.');
elseif ~iscell(syst.inv_epsilon) || (numel(syst.inv_epsilon) ~= 2 && numel(syst.inv_epsilon) ~= 3)
error('syst.inv_epsilon must be a two-element or three-element cell array, if given.');
elseif ~(ismatrix(syst.inv_epsilon{1}) && isnumeric(syst.inv_epsilon{1}))
error('syst.inv_epsilon{1} must be a numeric matrix.');
elseif ~(ismatrix(syst.inv_epsilon{2}) && isnumeric(syst.inv_epsilon{2}))
error('syst.inv_epsilon{2} must be a numeric matrix.');
elseif numel(syst.inv_epsilon) == 3 && ~(ismatrix(syst.inv_epsilon{3}) && isnumeric(syst.inv_epsilon{3}))
error('syst.inv_epsilon{3} must be a numeric matrix, if given.');
end
syst.polarization = 'TE';
str_pol = 'Hz'; % for printing system info
end
% Check that the user did not accidentally use options only in mesti()
if isfield(syst, 'PML') && ~isempty(syst.PML)
error('syst.PML is not used in mesti2s(); use syst.xBC instead.');
elseif isfield(syst, 'self_energy') && ~isempty(syst.self_energy)
error('syst.self_energy is not used in mesti2s(); use syst.xBC = ''outgoing'' instead.');
elseif isfield(syst, 'kx_B') && ~isempty(syst.kx_B)
error('syst.kx_B is not supported in mesti2s(); use mesti() if Bloch periodic BC in x is needed.');
elseif isfield(syst, 'PML_type') && ~isempty(syst.PML_type)
warning('syst.PML_type is not supported in mesti2s(); UPML will be used. Use mesti() if SC-PML is needed.')
syst = rmfield(syst, 'PML_type');
end
% syst.epsilon_R is an optional argument; determines whether the system is one-sided
if ~isfield(syst, 'epsilon_R') || isempty(syst.epsilon_R)
two_sided = false;
syst.epsilon_R = [];
elseif ~(isreal(syst.epsilon_R) && isscalar(syst.epsilon_R))
error('syst.epsilon_R must be a real scalar, if given.');
else
two_sided = true;
end
% Number of sites in y and x
if use_TM
[ny_Ez, nx_Ez] = size(syst.epsilon);
if ny_Ez==0; error('ny_Ez = size(syst.epsilon,1) cannot be zero.'); end
nx_s = nx_Ez; % number of pixels of Ez in the scattering region
dn = 0.5; % the source/detection is half a pixel away from x=0 and x=L
nx = nx_Ez;
ny = ny_Ez;
else
[ny_Hz, nx_Ez] = size(syst.inv_epsilon{2}); % (1/epsilon)_yy
if ny_Hz==0; error('ny_Hz = size(syst.inv_epsilon{2},1) cannot be zero.'); end
if two_sided
% nx_Hz = nx_Ez + 1 = L/syst.dx + 1 for a two-sided geometry
nx_Hz = nx_Ez + 1;
else
% nx_Hz = nx_Ez = L/syst.dx for a one-sided geometry
nx_Hz = nx_Ez;
end
nx_s = nx_Ez - 1; % number of pixels of Hz in the scattering region, excluding the two surfaces
dn = 0; % the source/detection is already at x=0 and x=L
nx = nx_Hz;
ny = ny_Hz;
end
% Check boundary condition in y
if isfield(syst, 'ky_B') && ~isempty(syst.ky_B)
if ~(isnumeric(syst.ky_B) && isscalar(syst.ky_B))
error('syst.ky_B must be a numeric scalar, if given.');
elseif (isfield(syst, 'yBC') && ~isempty(syst.yBC)) && (iscell(syst.yBC) || ~strcmpi(syst.yBC, 'Bloch'))
error('When syst.ky_B is given, syst.yBC must be ''Bloch'' if specified.');
end
syst.yBC = 'Bloch';
% mesti_build_channels() uses ky_B*periodicity as the input arguments yBC for Bloch BC
yBC = (syst.ky_B)*(ny*syst.dx); % dimensionless
else
if ~isfield(syst, 'yBC') || isempty(syst.yBC)
error('Input argument syst must have non-empty field ''yBC'' when syst.ky_B is not given.');
elseif ~((ischar(syst.yBC) && isrow(syst.yBC)) || (isstring(syst.yBC) && isscalar(syst.yBC)))
error('syst.yBC must be a character vector or string, if given.');
elseif ~ismember(lower(syst.yBC), lower({'Bloch', 'periodic', 'PEC', 'PMC', 'PECPMC', 'PMCPEC'}))
error('syst.yBC = ''%s'' is not a supported option; type ''help mesti2s'' for supported options.', syst.yBC);
elseif strcmpi(syst.yBC, 'Bloch')
error('syst.yBC = ''Bloch'' but syst.ky_B is not given.');
end
yBC = syst.yBC;
end
if ~use_TM
if strcmpi(syst.yBC, 'PMC')
ny_Ez = ny_Hz + 1;
elseif strcmpi(syst.yBC, 'PEC')
ny_Ez = ny_Hz - 1;
else
ny_Ez = ny_Hz;
end
% Check that size(syst.inv_epsilon{1}) = [ny_Ez, nx_Hz] since we will append to it later
if ny_Hz == 1 && (strcmpi(syst.yBC, 'periodic') || isequal(yBC, 0))
syst.inv_epsilon{1} = sparse(ny_Ez, nx_Hz); % (1/epsilon)_xx is not used when d/dy = 0
else
if isequal(size(syst.inv_epsilon{1}), [0, 0]) && ~(ny_Ez==0 && nx_Hz==0)
error('syst.inv_epsilon{1} must be given (unless ny_Hz = 1 and yBC = ''periodic'').');
elseif size(syst.inv_epsilon{1}, 1) ~= ny_Ez
error('size(syst.inv_epsilon{1}, 1) must equal ny_Ez = %d.', ny_Ez);
elseif size(syst.inv_epsilon{1}, 2) ~= nx_Hz
if two_sided
error('size(syst.inv_epsilon{1}, 2) must equal nx_Hz = size(syst.inv_epsilon{2}, 2) + 1 = %d for a two-sided geometry.', nx_Hz);
else
error('size(syst.inv_epsilon{1}, 2) must equal nx_Hz = size(syst.inv_epsilon{2}, 2) = %d for a one-sided geometry.', nx_Hz);
end
end
end
% Check that size(syst.inv_epsilon{3}) = [ny_Ez, nx_Ez] since we will append to it later
if numel(syst.inv_epsilon) == 3
if size(syst.inv_epsilon{3}, 1) ~= ny_Ez
error('size(syst.inv_epsilon{3}, 1) must equal ny_Ez = %d, if given.', ny_Ez);
elseif size(syst.inv_epsilon{3}, 2) ~= nx_Ez
error('size(syst.inv_epsilon{3}, 2) must equal nx_Ez = size(syst.inv_epsilon{2}, 2) = %d, if given.', nx_Ez);
end
end
end
% Defaults to using an exact outgoing boundary in x
if ~isfield(syst, 'xBC') || isempty(syst.xBC)
syst.xBC = 'outgoing';
elseif ~((~iscell(syst.xBC) && strcmpi(syst.xBC, 'outgoing')) || (isstruct(syst.xBC) && isscalar(syst.xBC)) || (iscell(syst.xBC) && numel(syst.xBC)==2))
error('syst.xBC must be ''outgoing'' or a scalar structure or a two-element cell array, if given.');
end
% Apply the same xBC on both sides; the second element will be ignored if two_sided = false
if ~(iscell(syst.xBC) && numel(syst.xBC)==2)
syst.xBC = {syst.xBC, syst.xBC};
elseif ~two_sided
error('For a one-sided geometry, boundary condition on the right is PEC for TM, PMC for TE; syst.xBC only specifies boundary condition on the left and must be ''outgoing'' or a scalar structure, if given.');
end
% Start with default setting: exact outgoing boundary using self-energy, with no PML and no spacer on both sides
% nx_extra is the number of homogeneous-space pixels to be added to syst.epsilon or syst.inv_epsilon{2} in x direction.
% The two elements of nx_extra and use_self_energy correspond to the left and right sides.
if two_sided
n_sides = 2;
nx_extra = [1, 1];
use_self_energy = [true, true];
str_xBC = {'outgoing', 'outgoing'}; % to be used for printing later
else
n_sides = 1;
nx_extra = [1, 0];
use_self_energy = [true, false];
if use_TM
str_xBC = {'outgoing', 'PEC'};
else
str_xBC = {'outgoing', 'PMC'};
end
end
syst.PML = {}; % to be used in mesti()
n_PML = 0; % number of sides where PML is used
% Loop over syst.xBC and handle PML parameters, if specified
str_sides = {'-', '+'};
for ii = 1:n_sides
if ~strcmpi(syst.xBC{ii}, 'outgoing')
use_self_energy(ii) = false;
n_PML = n_PML + 1;
str_xBC{ii} = 'PML';
PML_ii = syst.xBC{ii};
if ~(isstruct(PML_ii) && isscalar(PML_ii))
error('syst.xBC{%d} must be ''outgoing'' or a scalar structure.', ii)
end
% Check that the user did not accidentally use options only in mesti()
if isfield(PML_ii, 'npixels') && ~isempty(PML_ii.npixels)
warning('syst.xBC{%d}.npixels is ignored in mesti2s(); use field ''npixels_PML''.', ii);
end
if isfield(PML_ii, 'direction') && ~isempty(PML_ii.direction)
warning('syst.xBC{%d}.direction is ignored in mesti2s().', ii);
end
if isfield(PML_ii, 'side') && ~isempty(PML_ii.side)
warning('syst.xBC{%d}.side is ignored in mesti2s().', ii);
end
% Number of PML pixels must be given
% Other fields are optional and will be checked in mesti_build_fdfd_matrix()
if ~isfield(PML_ii, 'npixels_PML') || isempty(PML_ii.npixels_PML)
error('When syst.xBC{%d} is a structure, it must contain field ''npixels_PML''.', ii);
end
npix_PML = PML_ii.npixels_PML;
if ~(isreal(npix_PML) && isscalar(npix_PML) && npix_PML >= 0 && round(npix_PML) == npix_PML)
error('syst.xBC{%d}.npixels_PML must be a non-negative integer scalar.', ii);
end
% Defaults to no spacer
if ~isfield(PML_ii, 'npixels_spacer') || isempty(PML_ii.npixels_spacer)
PML_ii.npixels_spacer = 0;
end
npix_spacer = PML_ii.npixels_spacer;
if ~(isreal(npix_spacer) && isscalar(npix_spacer) && npix_spacer >= 0 && round(npix_spacer) == npix_spacer)
error('syst.xBC{%d}.npixels_spacer must be a non-negative integer scalar, if given.', ii);
end
% number of pixels in x to be added outside of syst.epsilon or syst.inv_epsilon
nx_extra(ii) = 1 + npix_PML + npix_spacer;
PML_ii = rmfield(PML_ii, {'npixels_PML', 'npixels_spacer'}); % these won't be used in mesti()
PML_ii.npixels = npix_PML; % to be used in mesti()
PML_ii.direction = 'x'; % to be used in mesti()
PML_ii.side = str_sides{ii}; % to be used in mesti()
syst.PML{n_PML} = PML_ii; % to be used in mesti()
end
end
% Total number of pixels in x
nx_tot = nx_s + sum(nx_extra);
%% Part 1.2: Check the input argument 'opts' and assign default values
if ~((isstruct(in) && isscalar(in)) || (iscell(in) && numel(in) <= 2))
error('Input argument in must be a scalar structure or a cell array with no more than two elements.');
end
% out is an optional argument
if nargin < 3
out = [];
end
if ~((isstruct(out) && isscalar(out)) || (iscell(out) && numel(out) <= 2) || isempty(out))
error('Input argument out must be a scalar structure or a cell array with no more than two elements or [], if given.');
end
% opts is an optional argument
if nargin < 4 || isempty(opts)
opts = struct();
end
if ~(isstruct(opts) && isscalar(opts))
error('Input argument opts must be a scalar structure or [], if given.');
end
% Check that the user did not accidentally use options only in mesti()
if isfield(opts, 'prefactor') && ~isempty(opts.prefactor)
error('opts.prefactor is not used in mesti2s(); the -2i prefactor is automatically included.');
end
% We return the scattering matrix S=C*inv(A)*B if out is given from input argument; else we return the spatial field profiles
% This opts.return_field_profile is not specified by the user; it will be returned as info.opts.return_field_profile to help debugging
opts.return_field_profile = isempty(out);
% By default, for field-profile computation, we only return result within the scattering region, setting nx_L = nx_R = 0.
if opts.return_field_profile
if ~isfield(opts, 'nx_L') || isempty(opts.nx_L)
opts.nx_L = 0;
elseif ~(isreal(opts.nx_L) && isscalar(opts.nx_L) && opts.nx_L >= 0 && round(opts.nx_L) == opts.nx_L)
error('opts.nx_L must be a non-negative integer scalar, if given.');
end
if ~isfield(opts, 'nx_R') || isempty(opts.nx_R)
opts.nx_R = 0;
elseif ~(isreal(opts.nx_R) && isscalar(opts.nx_R) && opts.nx_R >= 0 && round(opts.nx_R) == opts.nx_R)
error('opts.nx_R must be a non-negative integer scalar, if given.');
end
else
if isfield(opts, 'nx_L') && ~isempty(opts.nx_L)
warning('opts.nx_L is not used for scattering matrix computation; will be ignored.');
opts = rmfield(opts, 'nx_L');
end
if isfield(opts, 'nx_R') && ~isempty(opts.nx_R)
warning('opts.nx_R is not used for scattering matrix computation; will be ignored.');
opts = rmfield(opts, 'nx_R');
end
end
% Turn on verbal output by default
if ~isfield(opts, 'verbal') || isempty(opts.verbal)
opts.verbal = true;
elseif ~(islogical(opts.verbal) && isscalar(opts.verbal))
error('opts.verbal must be a logical scalar, if given.');
end
% By default, we don't clear syst in the caller's workspace
if ~isfield(opts, 'clear_syst') || isempty(opts.clear_syst)
opts.clear_syst = false;
elseif ~(islogical(opts.clear_syst) && isscalar(opts.clear_syst))
error('opts.clear_syst must be a logical scalar, if given.');
end
% By default, we will clear internal variables
if ~isfield(opts, 'clear_memory') || isempty(opts.clear_memory)
opts.clear_memory = true;
elseif ~(islogical(opts.clear_memory) && isscalar(opts.clear_memory))
error('opts.clear_memory must be a logical scalar, if given.');
end
% Use MUMPS for opts.solver when it is available
MUMPS_available = exist('zmumps','file');
solver_specified = true;
if ~isfield(opts, 'solver') || isempty(opts.solver)
solver_specified = false;
if MUMPS_available
opts.solver = 'MUMPS';
else
opts.solver = 'MATLAB';
end
elseif ~((ischar(opts.solver) && isrow(opts.solver)) || (isstring(opts.solver) && isscalar(opts.solver)))
error('opts.solver must be a character vector or string, if given.');
elseif ~ismember(lower(opts.solver), lower({'MUMPS', 'MATLAB'}))
error('opts.solver = ''%s'' is not a supported option; use ''MUMPS'' or ''MATLAB''.', opts.solver);
elseif strcmpi(opts.solver, 'MUMPS') && ~MUMPS_available
error('opts.solver = ''%s'' but function zmumps() is not found.', opts.solver)
end
method_specified = true;
if ~isfield(opts, 'method') || isempty(opts.method)
method_specified = false;
% By default, if the input argument 'out' is not given or if
% opts.iterative_refinement = true, then 'factorize_and_solve' is used.
% Otherwise, if ny is large or if TE polarization is used or if PML has
% been specified, then 'APF' is used when opts.solver = 'MUMPS', and
% 'C*inv(U)*inv(L)*B' is used when opts.solver = 'MATLAB'. Otherwise,
% 'RGF' is used.
if opts.return_field_profile || (isfield(opts, 'iterative_refinement') && isequal(opts.iterative_refinement, true))
opts.method = 'factorize_and_solve';
elseif n_PML > 0 || ~use_TM
if strcmpi(opts.solver, 'MUMPS')
opts.method = 'APF';
else
opts.method = 'C*inv(U)*inv(L)*B';
end
elseif isfield(opts, 'store_ordering') && isequal(opts.store_ordering, true)
opts.method = 'APF';
elseif isfield(opts, 'analysis_only') && isequal(opts.analysis_only, true)
opts.method = 'APF';
else
% APF and C*inv(U)*inv(L)*B are more efficient when ny is large. RGF is more efficient when ny is small.
if strcmpi(opts.solver, 'MUMPS')
if ny > 80
opts.method = 'APF';
else
opts.method = 'RGF';
end
else
if ny > 200
opts.method = 'C*inv(U)*inv(L)*B';
else
opts.method = 'RGF';
end
end
end
elseif ~((ischar(opts.method) && isrow(opts.method)) || (isstring(opts.method) && isscalar(opts.method)))
error('opts.method must be a character vector or string, if given.');
elseif ~ismember(lower(opts.method), lower({'APF', 'C*inv(U)*inv(L)*B', 'factorize_and_solve', 'FG', 'FS', 'RGF'}))
error('opts.method = ''%s'' is not a supported option; use ''APF'' or ''C*inv(U)*inv(L)*B'' or ''factorize_and_solve'' or ''RGF''.', opts.method);
elseif opts.return_field_profile && ~ismember(lower(opts.method), lower({'factorize_and_solve', 'FS'}))
error('opts.method = ''%s'' cannot be used for field-profile computations where out = []; use opts.method = ''factorize_and_solve'' instead.', opts.method)
elseif (isfield(opts, 'iterative_refinement') && isequal(opts.iterative_refinement, true)) && ~ismember(lower(opts.method), lower({'factorize_and_solve', 'FS'}))
error('opts.method = ''%s'' cannot be used when opts.iterative_refinement = true; use opts.method = ''factorize_and_solve'' instead.', opts.method)
elseif strcmpi(opts.method, 'FG')
opts.method = 'C*inv(U)*inv(L)*B'; % opts.method = 'FG' is short for opts.method = 'C*inv(U)*inv(L)*B'
elseif strcmpi(opts.method, 'FS')
opts.method = 'factorize_and_solve'; % opts.method = 'FS' is short for opts.method = 'factorize_and_solve'
end
if strcmpi(opts.method, 'APF') && strcmpi(opts.solver, 'MATLAB')
error('opts.method = ''APF'' requires opts.solver = ''MUMPS''.');
end
if strcmpi(opts.method, 'C*inv(U)*inv(L)*B') && strcmpi(opts.solver, 'MUMPS')
error('opts.method = ''C*inv(U)*inv(L)*B'' requires opts.solver = ''MATLAB''.');
end
if strcmpi(opts.method, 'RGF')
% check options not used in RGF
opts = rmfield(opts, 'solver'); % not used in RGF