-
Notifications
You must be signed in to change notification settings - Fork 0
/
wavesight.py
2454 lines (2233 loc) · 98.2 KB
/
wavesight.py
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
#!/usr/bin/env python3
import numbers
import warnings
import numpy as np
from misc import *
from metal import *
import cmasher as cm
from fields import *
from templates import *
from convstore import *
from datapipes import *
import diffkernels as dk
from fieldgenesis import *
from wavegraphics import *
from fungenerators import *
from collections import OrderedDict
from matplotlib import pyplot as plt
from scipy.interpolate import griddata
from matplotlib.patches import Rectangle
from scipy.interpolate import RegularGridInterpolator
real_dtype = np.float64
complex_dtype = np.complex128
############################################################################################
############################### Step Index Waveguide #######################################
def equivCurrents(Efuncs, Hfuncs, coreRadius, m, parity):
'''
This function returns the necessary magnetic and electric
currents corresponding to the given field functions. The
currents are calculated in a plane perpedicular to the
z-axis, and are returned so that they are functions of
cartesian x,y coordinates and represent the different
cartesian components of the currents.
Parameters
----------
Efuncs : 6-tuple of functions
Having the structure (ECoreρ, ECoreϕ, Ecorez,
ECladdingρ, ECladdingϕ, Ecladdingz). The z-component is
not needed, if it is included, then it is neglected. The
first three functions describing the electric field
inside of the core and the following three the electric
field in the cladding.
Hfuncs : 6-tuple of functions
Having the structure (HCoreρ, HCoreϕ, HCorez,
HCladdingρ, HCladdingϕ, HCladdingz). The z-component is
not needed, if it is included, then it is neglected. The
first three functions describing the H-field inside of
the core and the following three the H-field in the
cladding.
coreRadius : float
The radius of the core.
m : int
The azimuthal mode number.
parity : str
The parity of the mode. Can be 'TETM', 'EVEN', or 'ODD'.
Returns
-------
JKfuns : tuple of functions
Provided as (Jx, Jy, Kx, Ky). With Jx and Jy being the
functions for the x-y components of the equivalent
electric current, and Kx and Ky being the functions for
the x-y components of the equivalent magnetic current.
Each of the functions is a function of vec, which needs
to be an iterable with *three* elements corresponding to
the x,y,z cartesian components; this is to account for
the fact that MEEP requires these functions to take a
single Vector3 argument.
'''
(ECoreρ, ECoreϕ, _, ECladdingρ, ECladdingϕ, _) = Efuncs
(HCoreρ, HCoreϕ, _, HCladdingρ, HCladdingϕ, _) = Hfuncs
if parity == 'TETM':
def Jx(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (-HCoreϕ(ρ)*np.cos(ϕ)
- HCoreρ(ρ)*np.sin(ϕ))
else:
return (-HCladdingϕ(ρ)*np.cos(ϕ)
- HCladdingρ(ρ)*np.sin(ϕ))
def Jy(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (-HCoreϕ(ρ)*np.sin(ϕ)
+ HCoreρ(ρ)*np.cos(ϕ))
else:
return (-HCladdingϕ(ρ)*np.sin(ϕ)
+ HCladdingρ(ρ)*np.cos(ϕ))
def Kx(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (ECoreϕ(ρ)*np.cos(ϕ)
+ ECoreρ(ρ)*np.sin(ϕ))
else:
return (ECladdingϕ(ρ)*np.cos(ϕ)
+ ECladdingρ(ρ)*np.sin(ϕ))
def Ky(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (ECoreϕ(ρ)*np.sin(ϕ)
- ECoreρ(ρ)*np.cos(ϕ))
else:
return (ECladdingϕ(ρ)*np.sin(ϕ)
- ECladdingρ(ρ)*np.cos(ϕ))
elif parity == 'EVEN':
def Jx(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (-HCoreϕ(ρ)*np.cos(m*ϕ)*np.cos(ϕ)
- HCoreρ(ρ)*np.sin(m*φ)*(-1)*np.sin(ϕ))
else:
return (-HCladdingϕ(ρ)*np.cos(m*φ)*np.cos(ϕ)
- HCladdingρ(ρ)*np.sin(m*φ)*(-1)*np.sin(ϕ))
def Jy(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (-HCoreϕ(ρ)*np.cos(m*ϕ)*np.sin(ϕ)
+ HCoreρ(ρ)*np.sin(m*φ)*(-1)*np.cos(ϕ))
else:
return (-HCladdingϕ(ρ)*np.cos(m*φ)*np.sin(ϕ)
+ HCladdingρ(ρ)*np.sin(m*φ)*(-1)*np.cos(ϕ))
def Kx(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (ECoreϕ(ρ)*np.sin(m*φ)*(-1)*np.cos(ϕ)
+ ECoreρ(ρ)*np.cos(m*φ)*np.sin(ϕ))
else:
return (ECladdingϕ(ρ)*np.sin(m*φ)*(-1)*np.cos(ϕ)
+ ECladdingρ(ρ)*np.cos(m*φ)*np.sin(ϕ))
def Ky(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (ECoreϕ(ρ)*np.sin(m*φ)*(-1)*np.sin(ϕ)
- ECoreρ(ρ)*np.cos(m*φ)*np.cos(ϕ))
else:
return (ECladdingϕ(ρ)*np.sin(m*φ)*(-1)*np.sin(ϕ)
- ECladdingρ(ρ)*np.cos(m*φ)*np.cos(ϕ))
elif parity == 'ODD':
def Jx(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (-HCoreϕ(ρ)*np.sin(m*φ)*np.cos(ϕ)
- HCoreρ(ρ)*np.cos(m*φ)*np.sin(ϕ))
else:
return (-HCladdingϕ(ρ)*np.sin(m*φ)*np.cos(ϕ)
- HCladdingρ(ρ)*np.cos(m*φ)*np.sin(ϕ))
def Jy(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (-HCoreϕ(ρ)*np.sin(m*φ)*np.sin(ϕ)
+ HCoreρ(ρ)*np.cos(m*φ)*np.cos(ϕ))
else:
return (-HCladdingϕ(ρ)*np.sin(m*φ)*np.sin(ϕ)
+ HCladdingρ(ρ)*np.cos(m*φ)*np.cos(ϕ))
def Kx(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (ECoreϕ(ρ)*np.cos(m*φ)*np.cos(ϕ)
+ ECoreρ(ρ)*np.sin(m*φ)*np.sin(ϕ))
else:
return (ECladdingϕ(ρ)*np.cos(m*φ)*np.cos(ϕ)
+ ECladdingρ(ρ)*np.sin(m*φ)*np.sin(ϕ))
def Ky(vec):
x, y, _ = vec
ϕ = np.arctan2(y,x)
ρ = np.sqrt(x**2 + y**2)
if ρ < coreRadius:
return (ECoreϕ(ρ)*np.cos(m*φ)*np.sin(ϕ)
- ECoreρ(ρ)*np.sin(m*φ)*np.cos(ϕ))
else:
return (ECladdingϕ(ρ)*np.cos(m*φ)*np.sin(ϕ)
- ECladdingρ(ρ)*np.sin(m*φ)*np.cos(ϕ))
JKfuns = (Jx, Jy, Kx, Ky)
return JKfuns
def multisolver(fiber_spec, solve_modes = 'all', drawPlots=False, verbose=False, tm_te_funcs=False):
'''
This function solves the propagation constants of a step-index fiber with the given specifications. This assuming that the cladding is so big that it is effectively infinite.
The propagation constants determine the z-dependence of the fields along the z-direction. The propagation constants are always bounded by what would be the plane wave wavenumbers in the cladding or in the core.
Parameters
----------
fiber_spec : dict with the following keys:
nCore : float
The refractive index of the core.
nCladding : float
The refractive index of the cladding. If NA is given then
nCladding is computed to be consistent with the given nCore
and given NA.
coreRadius : float
The radius of the core in μm.
λFree : float
The wavelength of the light in free space in μm.
grid_divider: int, not necessary here but when later
used in the layout generator, this is used to determine
the fineness of the grid by making it equal to
λfree / max(nCore, nCladding, nBetween) / grid_divider
solve_modes: str
either 'all' or 'transvserse'
drawPlots : bool, optional
Whether to draw plots of the mode profiles. The default is False.
verbose : bool, optional
Whether to print out extra information. The default is False.
tm_te_funcs : bool, optional
If True then the return dictionary includes the keys tmfun and tefun. The default is False.
Returns
-------
sol : dict with all the keys included in fiber_spec plus these others:
kzmax : float
2π/λfree * nCladding (no kz larger than this)
kzmin : float
2π/λfree * nCore (no kz smaller than this)
Vnum : float
The V number of the fiber.
numModesFromVnum: float
The number of modes according to the V number.
totalNumModesTE : int
The total number of TE modes that were found.
totalNumModesTM : int
The total number of TM modes that were found.
totalNumModesHE : int
The total number of HE modes that were found, excluding parity.
totalNumModes : int
The total number of modes that were found, including parity.
tmfun : function
The eigenvalue function for the TM modes.
tefun : function
The eigenvalue function for the TE modes.
hefuns : dict
The eigenvalue functions for the HE modes. The keys are values of m (m ≥ 1).
TMkz : array
The propagation constants of the TM(0,n) modes.
TEkz : array
The propagation constants of the TE(0,n) modes.
HEkz : dict
The propagation constants of the HE(m,n) modes. The keys are values of m (m ≥ 1). The values are arrays with the kz(m,n) propagation constants.
'''
assert solve_modes in ['transverse', 'all'], "solve_modes must be either all or transverse."
warnings.filterwarnings('ignore', 'invalid value encountered in sqrt')
nCore = fiber_spec['nCore']
if 'NA' in fiber_spec:
if verbose:
print("Estimating nCladding from nCore and NA ...")
NA = fiber_spec['NA']
nCladding = np.sqrt(nCore**2 - NA**2)
else:
nCladding = fiber_spec['nCladding']
NA = np.sqrt(nCore**2 - nCladding**2)
separator = "="*40
coreRadius = fiber_spec['coreRadius']
wavelength = fiber_spec['λFree']
kFree = 2 * np.pi / wavelength
kzmax = nCore * 2 * np.pi / wavelength
kzmin = nCladding * 2 *np.pi / wavelength
kzspan = kzmax - kzmin
kzmax = kzmax - 1e-7 * kzspan
kzmin = kzmin + 1e-7 * kzspan
# split the solution domain in at least 300 parts
dkz = (kzmax - kzmin) / 1000
sol = fiber_spec
Vnum = 2 * np.pi * coreRadius * NA / wavelength
sol['Vnum'] = Vnum
sol['kFree'] = kFree
numModes = int(Vnum**2/2)
numModesTE = int(Vnum/2)
numModesTM = int(Vnum/2)
numModesHE = int(Vnum*(Vnum-2)/4)
sol['numModesFromVnum'] = numModes
sol['nCladding'] = nCladding
maxHEmodes = int(np.sqrt(2*numModesHE))
if verbose:
print(separator)
print("Approx number of complex HE modes: ", numModesHE)
print("Approx number of TE modes: ", numModesTE)
print("Approx number of TM modes: ", numModesTM)
print("Approx number of total modes: ", numModes)
print("Approx Max n for HE modes: ", maxHEmodes)
print(separator)
sol['kzmax'] = kzmax
sol['kzmin'] = kzmin
# secular equation for TM modes
tmfun = tmfungen(λfree = wavelength,
n2 = nCladding,
n1 = nCore,
a = coreRadius)
# secular equation function for TE modes
tefun = tefungen(λfree = wavelength,
n2 = nCladding,
n1 = nCore,
a = coreRadius)
if tm_te_funcs:
sol['tmfun'] = tmfun
sol['tefun'] = tefun
if solve_modes in ['transverse', 'all']:
dkzprime = dkz/numModesTE/2.
if verbose:
print("Calculating TM(0,n) propagation constants ...")
tmmodes = findallroots(tmfun, kzmin, kzmax, dkzprime,
dtype=real_dtype, method='brentq', num_sigfigs=10, reverse=True)
if verbose:
print("Calculating TE(0,n) propagation constants ...")
temodes = findallroots(tefun, kzmin, kzmax, dkzprime,
dtype=real_dtype, method='brentq', num_sigfigs=10, reverse=True)
kzrange = np.linspace(kzmin, kzmax, 1000, dtype=real_dtype)
if drawPlots:
tmvals = tmfun(kzrange)
tevals = tefun(kzrange)
tmzerocheck = tmfun(tmmodes)
tezerocheck = tefun(temodes)
plt.figure(figsize=(10,3))
plt.plot(kzrange, tmvals, 'r')
plt.scatter(tmmodes,tmzerocheck, c='y')
plt.plot([kzmin, kzmax], [0,0], "w:")
plt.ylim(-1,1)
plt.title('TM modes (%d found)' % len(tmmodes))
plt.show()
plt.figure(figsize=(10,3))
plt.plot(kzrange, tevals, 'r')
plt.scatter(temodes,tezerocheck, c='y')
plt.plot([kzmin, kzmax], [0,0], "w:")
plt.ylim(-1,1)
plt.title('TE modes (%d found)' % len(temodes))
plt.show()
hemodes = {}
if solve_modes == 'all':
m = 1
if verbose:
print("Calculating HE(m,n) propagation constants ...")
while True:
approxModes = maxHEmodes - m
approxModes = max(2, approxModes)
dkzprime = dkz / approxModes
if verbose:
print(f'm={m}',end='\r')
hefun = hefungen(λfree=wavelength,
m=m,
n2=nCladding,
n1=nCore,
a=coreRadius)
# sol['hefuns'][m] = hefun
hevals = hefun(kzrange)
hezeros = findallroots(hefun, kzmin, kzmax, dkzprime,
dtype=real_dtype, method='brentq', num_sigfigs=10, verbose=False, reverse=True)
if len(hezeros) == 0:
break
hemodes[m] = hezeros
hezerocheck = hefun(hezeros)
if drawPlots:
plt.figure(figsize=(15,3))
plt.plot(kzrange, hevals, 'r')
plt.plot([kzmin, kzmax], [0,0], "w:")
plt.scatter(hezeros, hezerocheck, c='y')
plt.ylim(-0.01,0.04)
plt.title('HE(%d, n) roots (%d found)' % (m, len(hezeros)))
plt.show()
m = m + 1
numCalcModes = (2 * sum(list(map(len,hemodes.values()))), len(temodes), len(tmmodes))
if verbose:
print("")
print(separator)
print("HE modes = %s\nTE modes = %d\nTM modes = %d\nTOTAL modes = %d\nFROM_Vnum = %d" % (*numCalcModes, sum(numCalcModes), numModes))
print(separator)
# put the modes in the solution dictionary
sol['TEkz'] = {0: temodes}
sol['TMkz'] = {0: tmmodes}
sol['HEkz'] = hemodes
totalModesTE = len(sol['TEkz'][0])
totalModesTM = len(sol['TMkz'][0])
totalModesHE = sum(list(map(len, sol['HEkz'].values())))
sol['totalModesTE'] = totalModesTE
sol['totalModesTM'] = totalModesTM
sol['totalModesHE'] = totalModesHE
sol['totalModes'] = totalModesTE + totalModesTM + 2*totalModesHE
return sol
############################### Step Index Waveguide #######################################
############################################################################################
def boundary_test(Efuncs, Hfuncs, fiber_spec, modeType, tolerance=1e-5):
'''
This function checks that a given solution for the fields
satisfies the boundary conditions at the core/cladding
interface within tolerance. If the mode is HE the test is on
the relative difference of all field components being within
tolerance. If the mode is TM or TE then the test is that all
the components that should be zero are exactly zero, and
that the relative difference for the non-zero components is
within tolerance.
Parameters
----------
Efuncs : tuple
a tuple of six functions which describe the three
components of the H field in the core and cladding
regions.
Hfuncs : tuple
a tuple of six functions which describe the three
components of the H field in the core and cladding
regions.
fiber_spec : dict
a dictionary containing the fiber specifications.
modeType : str
must be one of 'HE', 'TE', or 'TM'.
tolerance :float
the tolerance for the boundary condition test.
Returns
-------
boundaryTest : bool
True if test OK, False otherwise.
'''
(ECoreρ, ECoreϕ, ECorez, ECladdingρ, ECladdingϕ, ECladdingz) = Efuncs
(HCoreρ, HCoreϕ, HCorez, HCladdingρ, HCladdingϕ, HCladdingz) = Hfuncs
nCore = fiber_spec['nCore']
nCladding = fiber_spec['nCladding']
a = fiber_spec['coreRadius']
coreBoundary = OrderedDict()
claddingBoundary = OrderedDict()
# the radial component of D is continuous,
coreBoundary['DCoreρ'] = ECoreρ(a) * nCore**2
claddingBoundary['DCladdingρ'] = ECladdingρ(a) * nCladding**2
# the radial component of B is continuous,
coreBoundary['BCoreρ'] = HCoreρ(a)
claddingBoundary['BCladdingρ'] = HCladdingρ(a)
# the azimuthal component of E is continuous,
coreBoundary['ECoreϕ'] = ECoreϕ(a)
claddingBoundary['ECladdingϕ'] = ECladdingϕ(a)
# the longitudinal component of E is continuous,
coreBoundary['ECorez'] = ECorez(a)
claddingBoundary['ECladdingz'] = ECladdingz(a)
# the azimuthal component of H is continuous,
coreBoundary['HCoreϕ'] = HCoreϕ(a)
claddingBoundary['HCladdingϕ'] = HCladdingϕ(a)
# the longitudinal component of H is continuous.
coreBoundary['HCorez'] = HCorez(a)
claddingBoundary['HCladdingz'] = HCladdingz(a)
coreArray = np.array(list(coreBoundary.values()))
claddingArray = np.array(list(claddingBoundary.values()))
if modeType == 'HE':
nonZeroDiff = np.max(np.abs((coreArray / claddingArray)) - 1)
nonZeroTest = nonZeroDiff < tolerance
zeroCheck = True
elif modeType == 'TE':
coreZeros = np.zeros(3,dtype=np.complex128)
claddingZeros = np.zeros(3,dtype=np.complex128)
# the z-component of E should be zero
coreZeros[0] = coreBoundary['ECorez']
claddingZeros[0] = claddingBoundary['ECladdingz']
# the ρ-component of E should be zero
coreZeros[1] = coreBoundary['DCoreρ']
claddingZeros[1] = claddingBoundary['DCladdingρ']
# the ϕ-component of H should be zero
coreZeros[2] = coreBoundary['HCoreϕ']
claddingZeros[2] = claddingBoundary['HCladdingϕ']
zeroCheck = np.all(coreZeros == np.zeros(3)) and np.all(claddingZeros == np.zeros(3))
# Now for all the ones that should be non-zero
coreNonZeros = np.zeros(3,dtype=np.complex128)
claddingNonZeros = np.zeros(3,dtype=np.complex128)
coreNonZeros[0] = coreBoundary['ECoreϕ']
claddingNonZeros[0] = claddingBoundary['ECladdingϕ']
coreNonZeros[1] = coreBoundary['BCoreρ']
claddingNonZeros[1] = claddingBoundary['BCladdingρ']
coreNonZeros[2] = coreBoundary['HCorez']
claddingNonZeros[2] = claddingBoundary['HCladdingz']
nonZeroDiff = np.max(np.abs((coreNonZeros)/claddingNonZeros) - 1)
nonZeroTest = nonZeroDiff < tolerance
elif modeType == 'TM':
coreZeros = np.zeros(3, dtype=np.complex128)
claddingZeros = np.zeros(3, dtype=np.complex128)
# the ϕ-component of E should be zero
coreZeros[0] = coreBoundary['ECoreϕ']
claddingZeros[0] = claddingBoundary['ECladdingϕ']
# the z-component of H should be zero
coreZeros[1] = coreBoundary['HCorez']
claddingZeros[1] = claddingBoundary['HCladdingz']
# the ρ-component of B should be zero
coreZeros[2] = coreBoundary['BCoreρ']
claddingZeros[2] = claddingBoundary['BCladdingρ']
zeroCheck = np.all(coreZeros == np.zeros(3)) and np.all(claddingZeros == np.zeros(3))
# Now for all the ones that should be non-zero
coreNonZeros = np.zeros(3,dtype=np.complex128)
claddingNonZeros = np.zeros(3,dtype=np.complex128)
coreNonZeros[0] = coreBoundary['DCoreρ']
claddingNonZeros[0] = claddingBoundary['DCladdingρ']
coreNonZeros[1] = coreBoundary['ECorez']
claddingNonZeros[1] = claddingBoundary['ECladdingz']
coreNonZeros[2] = coreBoundary['HCoreϕ']
claddingNonZeros[2] = claddingBoundary['HCladdingϕ']
nonZeroDiff = np.max(np.abs((coreNonZeros)/claddingNonZeros) - 1)
nonZeroTest = nonZeroDiff < tolerance
boundaryTest = (nonZeroTest and zeroCheck)
if not boundaryTest:
print(claddingBoundary, coreBoundary)
return (boundaryTest, nonZeroTest, zeroCheck, nonZeroDiff)
def calculate_size_of_grid(fiber_sol):
'''
Given a solution for the modes of a multimode fiber, determine
the half side of the computational domain that would capture
most of the energy contained in the solved modes.
Parameters
----------
fiber_sol : dict
a dictionary containing the solution for the fiber.
Returns
-------
b : float
the width of the computational domain.
'''
goal_fraction = 0.99
a = fiber_sol['coreRadius']
nCore = fiber_sol['nCore']
nCladding = fiber_sol['nCladding']
nBetween = fiber_sol['nBetween']
λfree = fiber_sol['λFree']
kFree = 2*np.pi/λfree
grid_divider = fiber_sol['grid_divider']
maxIndex = max(nCore, nCladding, nBetween)
Δs = λfree / maxIndex / grid_divider
# calculate the side of the computational domain: START
allnums = []
for modetype in ['TE', 'TM', 'HE']:
allkzs = fiber_sol[modetype + 'kz']
for m, kzs in allkzs.items():
for kz in kzs:
wnums = (modetype, m, kz)
allnums.append(wnums)
(modetype, m, kz) = sorted(allnums, key=lambda x: x[-1])[0]
(Efuncs, Hfuncs) = fieldGenerator(a, kFree, kz, m, nCladding, nCore, modetype)
ρrange = np.linspace(0, a + 20 * λfree / nCladding, 1000)
fluxCore = Efuncs[0](ρrange) * np.conjugate(Hfuncs[1](ρrange)) - Efuncs[1](ρrange) * np.conjugate(Hfuncs[0](ρrange))
fluxCladding = Efuncs[3](ρrange) * np.conjugate(Hfuncs[4](ρrange)) - Efuncs[4](ρrange) * np.conjugate(Hfuncs[3](ρrange))
flux = fluxCore
flux[ρrange>a] = fluxCladding[ρrange>a]
integrand = flux * ρrange
total = np.sum(integrand)
insideEnergy = np.cumsum(integrand)/total
b = 2*np.interp(goal_fraction, insideEnergy, ρrange) + λfree
numSigFigsina = sig_figs_in(a)
b = rounder(b, numSigFigsina)
return b
def coordinate_layout(fiber_sol):
'''
Given a fiber solution, return the coordinate arrays for plotting
and coordinating the numerical analysis.
Parameters
----------
fiber_sol : dict with at least the following keys:
'coreRadius' : float
The radius of the core.
'claddingIndex' : float
The refractive index of the cladding.
'coreIndex' : float
The refractive index of the core.
'free_space_wavelength' : float
The wavelength of the light in free space.
Returns
-------
a, b, Δs, xrange, yrange, ρrange, φrange, Xg, Yg, ρg, φg, nxy, crossMask : tuple
a : float
The radius of the core.
b : float
The side of the computational domain.
Δs : float
The sampling resolution in the x-y direction. Assumed to be half the
wavelength in the core.
xrange, yrange : 1D arrays
The coordinate arrays for the x-y directions.
ρrange, φrange : 1D arrays
The coordinate arrays for the ρ-φ directions in the cylindrical system.
Xg, Yg : 2D arrays
The coordinate arrays for x-y.
ρg, φg : 2D arrays
ρg has the values for the radial coordinate, φg has the values for the
azimuthal coordinate.
nxy : 2D array
The refractive index profile of the waveguide.
crossMask : 2D array
A mask that is True where the core is and False where the cladding is.
'''
a = fiber_sol['coreRadius']
nCore = fiber_sol['nCore']
nCladding = fiber_sol['nCladding']
nBetween = fiber_sol['nBetween']
λfree = fiber_sol['λFree']
kFree = 2*np.pi/λfree
grid_divider = fiber_sol['grid_divider']
maxIndex = max(nCore, nCladding, nBetween)
Δs = λfree / maxIndex / grid_divider
b = calculate_size_of_grid(fiber_sol)
# calculate the side of the computational domain: END
# calculate the coordinates arrays
numSamples = int(2 * b / Δs)
xrange = np.linspace(-b, b, numSamples, dtype=real_dtype)
yrange = np.linspace(-b, b, numSamples, dtype=real_dtype)
ρrange = np.linspace(0 , np.sqrt(2)*b, numSamples, dtype=real_dtype)
φrange = np.linspace(-np.pi, np.pi, numSamples, dtype=real_dtype)
Xg, Yg = np.meshgrid(xrange, yrange)
ρg = np.sqrt(Xg**2 + Yg**2)
φg = np.arctan2(Yg, Xg)
crossMask = np.zeros((numSamples, numSamples)).astype(np.bool8)
crossMask[ρg <= a] = True
crossMask[ρg > a] = False
# #Coords-Calc
nxy = np.zeros((numSamples, numSamples))
nxy[crossMask] = nCore
nxy[~crossMask] = nCladding
return a, b, Δs, xrange, yrange, ρrange, φrange, Xg, Yg, ρg, φg, nxy, crossMask, numSamples
def calculate_numerical_basis(fiber_sol, verbose=True):
'''
Given a solution for the propagation modes of an optical waveguide, calculate a numerical basis.
Parameters
----------
fiber_sol : dict with at least the following keys:
'coreRadius' : float
The radius of the core.
'nCladding' : float
The refractive index of the cladding.
'nCore' : float
The refractive index of the core.
'free_space_wavelength' : float
The wavelength of the light in free space.
'totalModes' : int
The total number of calculated modes.
'TEkz' : 1D dict
A single key equal to m=0, the values are array with
the TE modes propagation constants.
'TMkz' : 1D dict
A single key equal to m=0, the values are array with
the TM modes propagation constants.
'HEkz' : 1D dict
Keys are m values, values are 1D arrays of kz values.
Returns
-------
fiber_sol : dict the same dictionary as the input, but with two new keys:
'coord_layout' : tuple
The tuple returned by the coordinate_layout function.
'eigenbasis' : 5D np.array
The numerical basis. The first dimension is the mode number,
the second dimension only has two values, 0 and 1, 0 for the
E field, and 1 for the H field. The third dimension is used
for the different components of the corresponding field in
cylindrical coordinates. The first index being the radial
component, the second index being the azimuthal component,
and the third index being the z component. Finally the
fourth and fifth dimensions are arrays that hold the values
of the corresponding field components. The modes are
enumerated such that first all TE modes are given, then all
TM modes, and finally all HE modes. This same enumeration is
the one used for the 'eigenbasis_nums' key.
'eigenbasis_nums' : list of tuples
A list of tuples, each tuple has 7 values, the first value
is a string indicating the type of mode, and can be either
'TE', 'TM', or 'HE'. The second value is a string indicating
the parity of the mode. The third value is the value of m.
The fourth value corresponds to the value of kz, the
propagation constant along the z-direction. The fifth value
is the index that the value of kz has within the
corresponding array listing the propagation constants. The
sixth value is the transverse propagation constant inside
the core. The seventh value is the propagation constant in
the cladding. (modType, parity, m, kzidx, kz, γ, β)
'''
warnings.filterwarnings('ignore', 'invalid value encountered in sqrt')
warnings.filterwarnings('ignore', 'invalid value encountered in multiply')
warnings.filterwarnings('ignore', 'invalid value encountered in divide')
coord_layout = coordinate_layout(fiber_sol)
nCore = fiber_sol['nCore']
nCladding = fiber_sol['nCladding']
λfree = fiber_sol['λFree']
kFree = 2*np.pi/λfree
if 'coord_layout' in fiber_sol:
del fiber_sol['coord_layout']
if 'eigenbasis' in fiber_sol:
del fiber_sol['eigenbasis']
a, b, Δs, xrange, yrange, ρrange, φrange, Xg, Yg, ρg, φg, nxy, crossMask, numSamples = coord_layout
# determine how many modes there are in total
totalModes = fiber_sol['totalModes']
eigenbasis = np.zeros((totalModes, 2, 3, numSamples, numSamples), dtype=complex_dtype)
counter = 0
eigenbasis_nums = []
iter_fun = lambda x: x
for modType in ['TE','TM','HE']:
if modType == 'HE':
parities = ['EVEN', 'ODD']
else:
parities = ['TETM']
solkey = modType + 'kz'
for m, kzs in fiber_sol[solkey].items():
cosMesh = np.cos(m*φg)
sinMesh = np.sin(m*φg)
for kzidx, kz in enumerate(iter_fun(kzs)):
γ = np.sqrt(nCore**2*4*np.pi**2/λfree**2 - kz**2)
β = np.sqrt(kz**2 - nCladding**2*4*np.pi**2/λfree**2)
(Efuncs, Hfuncs) = fieldGenerator(a, kFree, kz, m, nCladding, nCore, modType)
(ECoreρ, ECoreϕ, ECorez, ECladdingρ, ECladdingϕ, ECladdingz) = Efuncs
(HCoreρ, HCoreϕ, HCorez, HCladdingρ, HCladdingϕ, HCladdingz) = Hfuncs
funPairs = (((ECoreρ, ECladdingρ), (HCoreρ, HCladdingρ)),
((ECoreϕ, ECladdingϕ), (HCoreϕ, HCladdingϕ)),
((ECorez, ECladdingz), (HCorez, HCladdingz)))
this_E = np.zeros((3, numSamples, numSamples), dtype=complex_dtype)
this_H = np.zeros((3, numSamples, numSamples), dtype=complex_dtype)
for idx, ((EfunCore, EfunCladding), (HfunCore, HfunCladding)) in enumerate(funPairs):
ECorevals = np.vectorize(EfunCore)(ρrange)
ECorevals = np.interp(ρg, ρrange, ECorevals)
ECorevals[~crossMask] = 0
ECladdingvals = np.vectorize(EfunCladding)(ρrange)
ECladdingvals = np.interp(ρg, ρrange, ECladdingvals)
ECladdingvals[crossMask] = 0
E_all = ECorevals + ECladdingvals
HCorevals = np.vectorize(HfunCore)(ρrange)
HCorevals = np.interp(ρg, ρrange, HCorevals)
HCorevals[~crossMask] = 0
HCladdingvals = np.vectorize(HfunCladding)(ρrange)
HCladdingvals = np.interp(ρg, ρrange, HCladdingvals)
HCladdingvals[crossMask] = 0
H_all = HCorevals + HCladdingvals
E_all[np.isnan(E_all)] = 0.
H_all[np.isnan(H_all)] = 0.
this_E[idx, :, :] = E_all
this_H[idx, :, :] = H_all
for parity in parities:
if parity == 'TETM':
for idx in range(3):
eigenbasis[counter, 0, idx, :, :] = this_E[idx]
eigenbasis[counter, 1, idx, :, :] = this_H[idx]
elif parity == 'EVEN':
for idx in range(3):
ϕPhaseE = {0: cosMesh,
1: -sinMesh,
2: cosMesh}[idx]
ϕPhaseH = {0: -sinMesh,
1: cosMesh,
2: -sinMesh}[idx]
eigenbasis[counter, 0, idx, :, :] = this_E[idx]*ϕPhaseE
eigenbasis[counter, 1, idx, :, :] = this_H[idx]*ϕPhaseH
elif parity == 'ODD':
for idx in range(3):
ϕPhaseE = {0: sinMesh,
1: cosMesh,
2: sinMesh}[idx]
ϕPhaseH = {0: cosMesh,
1: sinMesh,
2: cosMesh}[idx]
eigenbasis[counter, 0, idx, :, :] = this_E[idx]*ϕPhaseE
eigenbasis[counter, 1, idx, :, :] = this_H[idx]*ϕPhaseH
# normalize the field
E_field = eigenbasis[counter][0]
H_field = eigenbasis[counter][1]
norm_sq = field_dot(E_field, H_field, Δs)
norm = np.sqrt(norm_sq)
eigenbasis[counter][0] = eigenbasis[counter][0] / norm
eigenbasis[counter][1] = eigenbasis[counter][1] / norm
eigenbasis_nums.append((modType, parity, m, kzidx, kz, γ, β))
counter += 1
fiber_sol['eigenbasis'] = eigenbasis
fiber_sol['coord_layout'] = coord_layout
fiber_sol['eigenbasis_nums'] = eigenbasis_nums
return fiber_sol
############################################################################################
################################ Ray Refractor #############################################
def poynting_refractor(Efield, Hfield, nxy, nBetween, verbose=False):
'''
Approximate the refracted field across a planar interface
using the Poynting vector as an analog to the wavevector of
a plane-wave.
Any evanescent fields are ignored. All cases where there
would be total internal reflection the refracted field is
set to zero.
Parameters
----------
Efield : np.array (3, N, M)
The electric field incident on the interface.
Hfield : np.array (3, N, M)
The H-field incident on the interface.
nxy : np.array (N, M)
The refractive index transverse to the interface inside the incident medium.
nBetween : float
The refractive index of the homogeneous refractive medium.
verbose : bool, optional
Whether to print or not progress messages, by default False.
Returns
-------
Eref, Href : tuple of np.array (3, N, M)
The refracted electric and magnetic fields.
'''
# #EXH-Calc
# calculate the Poynting vector
if verbose:
print("Calculating the Poynting vector field...")
Sfield = 0.5*np.real(np.cross(Efield, np.conjugate(Hfield), axis=0))
# #normIncidentk-Calc
# calculate the magnitude of the Poynting field
if verbose:
print("Calculating the magnitude of the Poynting field...")
Sfieldmag = np.sqrt(np.sum(Sfield**2, axis=0))
if verbose:
print("Calculating the transverse component of the Poynting field...")
# calculate the unit vector in the direction of the Poynting vector
if verbose:
print("Calculating the unit vector in the direction of the Poynting vector...")
kfield = Sfield / Sfieldmag
# calculate the transverse component of the Poynting vector
Stransverse = np.sqrt(Sfield[0,:,:]**2 + Sfield[1,:,:]**2)
# #β-Calc
if verbose:
print("Calculating the angle of incidence field...")
# Assuming that the normal is pointing in the +z direction
βfield = np.arctan2(Stransverse, Sfield[2, :, :])
# #θ-Calc
if verbose:
print("Calculating the angle of refraction field...")
θfield = np.arcsin(nxy/nBetween * np.sin(βfield))
# #FresnelS-Calc
if verbose:
print("Calculating the Fresnel coefficients...")
fresnelS = (2 * nxy * np.cos(βfield)
/ ( nxy * np.cos(βfield)
+ np.sqrt(nBetween**2
- nxy**2 * np.sin(βfield)**2)
)
)
# #FresnelP-Calc
fresnelP = (2 * nBetween * nxy * np.cos(βfield)
/ (nBetween**2 * np.cos(βfield)
+ nxy * np.sqrt(nBetween**2
- nxy**2 * np.sin(βfield)**2)
)
)
# #ζ-Calc
if verbose:
print("Calculating the ζ of the local coord system...")
# calculate the unit vector field perpendicular to the plane of incidence
# which is basically k X z
ζfield = np.zeros(kfield.shape, dtype=complex_dtype)
ζfield[0] = kfield[1]
ζfield[1] = -kfield[0]
# normalize it
ζfield /= np.sqrt(np.sum(np.abs(ζfield)**2, axis=0))
# in the case of normal incidence, ζ is arbitrary
# so it can be set to the unit vector in the first direction
normalIncidence = (βfield == 0)
ζfield[0][normalIncidence] = 1
ζfield[1][normalIncidence] = 0
ζfield[2][normalIncidence] = 0
# #EincS-Calc
if verbose:
print("Calculating the S and P component of the incident electric field...")
# decompose the field in P and S polarizations
# first find P-pol and then use the complement to determine S-pol
# the S-pol can be obtained by the dot product of E with ζ
ESdot = Efield[0] * ζfield[0] + Efield[1] * ζfield[1]
EincS = np.zeros(ζfield.shape, dtype=complex_dtype)
EincS[0] = ESdot * ζfield[0]
EincS[1] = ESdot * ζfield[1]
# #EincP-Calc
EincP = Efield - EincS
del ESdot
# #HincS-Calc
if verbose:
print("Calculating the S and P component of the incident H field...")
# decompose the field in P and S polarizations
# first find P-pol and then use the complement to determine S-pol
# the S-pol can be obtained by the dot product of H with ζ
HSdot = Hfield[0] * ζfield[0] + Hfield[1] * ζfield[1]
HincS = np.zeros(ζfield.shape, dtype=complex_dtype)
HincS[0] = HSdot * ζfield[0]
HincS[1] = HSdot * ζfield[1]
# #HincP-Calc
HincP = Hfield - HincS
del HSdot
# #ErefS-Calc
if verbose:
print("Calculating the S and P component of the refracted electric field...")
ErefS = np.zeros(Efield.shape, dtype=complex_dtype)
ErefS[0] = fresnelS * EincS[0]
ErefS[1] = fresnelS * EincS[1]
ErefS[2] = fresnelS * EincS[2]
ErefS[np.isnan(ErefS)] = 0
# #HrefS-Calc
if verbose:
print("Calculating the S and P component of the refracted H field...")
HrefS = np.zeros(Hfield.shape, dtype=complex_dtype)
HrefS[0] = (nBetween / nxy) * fresnelP * HincS[0]
HrefS[1] = (nBetween / nxy) * fresnelP * HincS[1]
HrefS[2] = (nBetween / nxy) * fresnelP * HincS[2]
HrefS[np.isnan(HrefS)] = 0
# #ErefP-Calc
ErefP = np.zeros(Efield.shape, dtype=complex_dtype)
ErefP[0] = fresnelP * EincP[0]
ErefP[1] = fresnelP * EincP[1]
ErefP[2] = fresnelP * EincP[2]
ErefP[np.isnan(ErefP)] = 0
# #HrefP-Calc
HrefP = np.zeros(Hfield.shape, dtype=complex_dtype)
HrefP[0] = (nBetween / nxy) * fresnelS * HincP[0]
HrefP[1] = (nBetween / nxy) * fresnelS * HincP[1]
HrefP[2] = (nBetween / nxy) * fresnelS * HincP[2]
HrefP[np.isnan(HrefP)] = 0
# #Eref-Calc
if verbose:
print("Calculating the total refracted electric field...")
Eref = ErefS + ErefP
# #Href-Calc
if verbose:
print("Calculating the total refracted H field...")
Href = HrefS + HrefP
# #kref-Calc
if verbose:
print("Calculating the refracted wavevector (normalized) field...")
ξfield = np.zeros(ζfield.shape)
ξfield[0] = -ζfield[1]
ξfield[1] = ζfield[0]
kref = np.zeros(kfield.shape)
kref[0] = np.sin(θfield) * ξfield[0]
kref[1] = np.sin(θfield) * ξfield[1]
kref[2] = np.cos(θfield)