From 75ffb0375f74eb8966b4c0d377b3fda2dedf2d6c Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Thu, 16 Mar 2023 12:46:34 -0700 Subject: [PATCH] Tutorial for phase of mode coefficients (#2433) * tutorial for phase of mode coefficients * add text * updates to figure and text * provide more details on necessity of factor of 2 in front of L in phase correction factor --- .../Python_Tutorials/Mode_Decomposition.md | 258 ++++++++++++++++++ doc/docs/images/refl_coeff_flat_interface.png | Bin 0 -> 56457 bytes python/examples/mode_coeff_phase.py | 239 ++++++++++++++++ 3 files changed, 497 insertions(+) create mode 100644 doc/docs/images/refl_coeff_flat_interface.png create mode 100644 python/examples/mode_coeff_phase.py diff --git a/doc/docs/Python_Tutorials/Mode_Decomposition.md b/doc/docs/Python_Tutorials/Mode_Decomposition.md index e12340c1f..7f1d86172 100644 --- a/doc/docs/Python_Tutorials/Mode_Decomposition.md +++ b/doc/docs/Python_Tutorials/Mode_Decomposition.md @@ -143,6 +143,264 @@ The reflectance values computed using the two methods are nearly identical. For In the reflected-flux calculation, we apply our usual trick of first performing a reference simulation with just the incident field and then subtracting that from our taper simulation with `load_minus_flux_data`, so that what is left over is the reflected fields (from which we obtain the reflected flux). In *principle*, this trick would not be required for the mode-decomposition method, because the reflected mode is orthogonal to the forward mode and so the decomposition will separate the forward and reflected coefficients automatically. However, this is only true in the limit of infinite resolution — for a *finite* resolution, the reflected mode used for the mode coefficient calculation (calculated via MPB) is not exactly orthogonal to the forward mode propagating in Meep (whose discretization scheme is different from that of MPB). In consequence, if you did not subtract the fields of the reference simulation, the mode-coefficient could only calculate the reflected power down to a "noise floor" set by the discretization error. With the subtraction, in contrast, you can compute much smaller reflections (limited by the floating-point precision). +Phase of a Total Internal Reflected Mode +---------------------------------------- + +For certain applications, such as [physically based ray tracing](https://pbr-book.org/), the phase of the reflection/transmission coefficient of a surface interface (flat or textured) is an important parameter used for modeling coherent effects (i.e., constructive or destructive interference from a collection of intersecting rays scattered off the surface). This tutorial demonstrates how to use mode decomposition to compute the phase of the reflection coefficient of a [total internal reflected](https://en.wikipedia.org/wiki/Total_internal_reflection) (TIR) mode. The simulated results are verified using the [Fresnel equations](https://en.wikipedia.org/wiki/Total_internal_reflection#Phase_shifts). + +The example involves a 2d simulation of a flat interface of two lossless media with $n_1=1.5$ and $n_2=1.0$. The [critical angle](https://en.wikipedia.org/wiki/Total_internal_reflection#Critical_angle) for this interface is $\theta_c = 41.8^\circ$. The source is an incident planewave with linear polarization ($S$ or $P$). The 2d cell is periodic in $y$ with PML in the $x$ direction. The cell size does not affect the results and is therefore arbitrary. A schematic of the simulation layout is shown below. + +The key item to consider in these types of calculations is the *location of the mode monitor relative to the interface*. The mode monitor is a line extending the entire length of the cell in the $y$ direction. In order to compare the result with the Fresnel equations, the phase of the reflection coefficient must be computed *exactly* at the interface. However, it is problematic to measure the reflection coefficient exactly at the interface because the amplitude of the left-going wave drops discontinuously to zero across the interface. For this reason, the mode monitor must be positioned away from the interface (somewhere within the $n_1$ medium) and the measured phase must be corrected in post processing to account for this offset. + +![](../images/refl_coeff_flat_interface.png#center) + +The calculation of the reflection coefficient requires two separate runs: (1) a normalization run involving just the source medium ($n_1$) to obtain the incident fields, and (2) the main run including the interface whereby the incident fields from (1) are first subtracted from the monitor to obtain only the reflected fields. The mode coefficient in (2) divided by (1) is, by definition, the reflection coefficient. + +The phase of the reflection coefficient needs to be corrected for the offset of the mode monitor relative to the interface — labeled $L$ in the schematic above — using the formula $\exp(i 2\pi k_x 2L)$, where $k_x$ is the $x$ component of the planewave's wavevector (`k` in the script). The $2\pi$ factor is necessary because `k` does *not* include this factor (per convention in Meep). The factor 2 in front of the $L$ is necessary because the phase needs to be corrected for the monitors in the normalization and main runs separately. The correction factor is just the phase accumulated as the wave propagates in the surface-normal direction for a distance $L$ in a medium with index $n_1$. + +With this setup, we measure the phase of the reflection coefficient for two different source configurations (polarization and angle) and compare the result with the Fresnel equations. The location of the mode monitor is also varied in the two test configurations. Results are shown in the two tables below. There is good agreement between the simulation and theory. (As additional validation, we note that the magnitude of the reflection coefficient of a TIR mode must be one which is indeed the case.) + +| $S$ pol., $\theta$ = 54.3° | reflection coefficient | phase (radians) | +|:--------------------------:|:----------------------:|:---------------:| +| Meep | 0.23415-0.96597j | -1.33299 | +| Fresnel | 0.22587-0.97416j | -1.34296 | + + +| $P$ pol., $\theta$ = 48.5° | reflection coefficient | phase (radians) | +|:--------------------------:|:----------------------:|:---------------:| +| Meep | 0.11923+0.98495j | 1.45033 | +| Fresnel | 0.14645+0.98922j | 1.42382 | + +The simulation script is in [examples/mode_coeff_phase.py](https://github.com/NanoComp/meep/blob/master/python/examples/mode_coeff_phase.py). + +```py +import cmath +from enum import Enum +import math + +import meep as mp +import numpy as np + +Polarization = Enum("Polarization", "S P") + +# refractive indices of materials 1 and 2 +n1 = 1.5 +n2 = 1.0 + + +def refl_coeff_meep(pol: Polarization, theta: float, L: float) -> complex: + """Computes the reflection coefficient of a TIR mode using mode + decomposition. + + Args: + pol: polarization of the incident planewave (S or P). + theta: angle of the incident planewave (radians). + L: position of the mode monitor relative to the flat interface. 0 is + interface position. + """ + resolution = 50.0 + + # cell size is arbitrary + sx = 7.0 + sy = 3.0 + dpml = 2.0 + cell_size = mp.Vector3(sx + 2 * dpml, sy, 0) + pml_layers = [mp.PML(dpml, direction=mp.X)] + + fcen = 1.0 # center frequency + df = 0.1 * fcen + + # k (in source medium) with correct length + # plane of incidence is xy + k = mp.Vector3(n1 * fcen, 0, 0).rotate(mp.Vector3(0, 0, 1), theta) + + # planewave amplitude function (for source) + def pw_amp(k, x0): + def _pw_amp(x): + return cmath.exp(1j * 2 * math.pi * k.dot(x + x0)) + + return _pw_amp + + src_pt = mp.Vector3(-0.5 * sx, 0, 0) + + if pol.name == "S": + src_cmpt = mp.Ez + eig_parity = mp.ODD_Z + elif pol.name == "P": + src_cmpt = mp.Hz + eig_parity = mp.EVEN_Z + else: + raise ValueError("pol must be S or P, only.") + + sources = [ + mp.Source( + mp.GaussianSource(fcen, fwidth=df), + component=src_cmpt, + center=src_pt, + size=mp.Vector3(0, cell_size.y, 0), + amp_func=pw_amp(k, src_pt), + ), + ] + + sim = mp.Simulation( + resolution=resolution, + cell_size=cell_size, + default_material=mp.Medium(index=n1), + boundary_layers=pml_layers, + k_point=k, + sources=sources, + ) + + # DFT monitor for incident fields + mode_mon = sim.add_mode_monitor( + fcen, + 0, + 1, + mp.FluxRegion( + center=mp.Vector3(-L, 0, 0), + size=mp.Vector3(0, cell_size.y, 0), + ), + ) + + sim.run( + until_after_sources=mp.stop_when_fields_decayed( + 50, + src_cmpt, + mp.Vector3(-L, 0, 0), + 1e-6, + ), + ) + + res = sim.get_eigenmode_coefficients( + mode_mon, + bands=[1], + eig_parity=eig_parity, + kpoint_func=lambda *not_used: k, + direction=mp.NO_DIRECTION, + ) + + input_mode_coeff = res.alpha[0, 0, 0] + input_flux_data = sim.get_flux_data(mode_mon) + + sim.reset_meep() + + geometry = [ + mp.Block( + material=mp.Medium(index=n1), + center=mp.Vector3(-0.25 * (sx + 2 * dpml), 0, 0), + size=mp.Vector3(0.5 * (sx + 2 * dpml), mp.inf, mp.inf), + ), + mp.Block( + material=mp.Medium(index=n2), + center=mp.Vector3(0.25 * (sx + 2 * dpml), 0, 0), + size=mp.Vector3(0.5 * (sx + 2 * dpml), mp.inf, mp.inf), + ), + ] + + sim = mp.Simulation( + resolution=resolution, + cell_size=cell_size, + boundary_layers=pml_layers, + k_point=k, + sources=sources, + geometry=geometry, + ) + + # DFT monitor for reflected fields + mode_mon = sim.add_mode_monitor( + fcen, + 0, + 1, + mp.FluxRegion( + center=mp.Vector3(-L, 0, 0), + size=mp.Vector3(0, cell_size.y, 0), + ), + ) + + sim.load_minus_flux_data(mode_mon, input_flux_data) + + sim.run( + until_after_sources=mp.stop_when_fields_decayed( + 50, + mp.Ez, + mp.Vector3(-L, 0, 0), + 1e-6, + ), + ) + + res = sim.get_eigenmode_coefficients( + mode_mon, + bands=[1], + eig_parity=eig_parity, + kpoint_func=lambda *not_used: k, + direction=mp.NO_DIRECTION, + ) + + # mode coefficient of reflected planewave + refl_mode_coeff = res.alpha[0, 0, 1] + + # reflection coefficient + refl_coeff = refl_mode_coeff / input_mode_coeff + + # apply phase correction factor + refl_coeff /= cmath.exp(1j * k.x * 2 * math.pi * 2 * L) + + return refl_coeff + + +def refl_coeff_Fresnel(pol: Polarization, theta: float) -> complex: + """Computes the reflection coefficient of a TIR mode using the Fresnel + equations. + + Args: + pol: polarization of the incident planewave (S or P). + theta: angle of the incident planewave (degrees). + """ + if pol.name == "S": + refl_coeff = ( + math.cos(theta) - ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5 + ) / (math.cos(theta) + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5) + else: + refl_coeff = ( + -((n2 / n1) ** 2) * math.cos(theta) + + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5 + ) / ( + (n2 / n1) ** 2 * math.cos(theta) + + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5 + ) + + return refl_coeff + + +if __name__ == "__main__": + thetas = [54.3, 48.5] # angle of incident planewave (degrees) + Ls = [0.4, 1.2] # position of mode monitor relative to flat interface + pols = [Polarization.S, Polarization.P] # polarization of incident planewave + + for pol, theta, L in zip(pols, thetas, Ls): + theta_rad = np.radians(theta) + rc_m = refl_coeff_meep(pol, theta_rad, L) + rc_f = refl_coeff_Fresnel(pol, theta_rad) + + rc_m_str = f"{rc_m.real:.5f}{rc_m.imag:+.5f}j" + rc_f_str = f"{rc_f.real:.5f}{rc_f.imag:+.5f}j" + print(f"refl-coeff:, {pol.name}, {theta}, {rc_m_str} (Meep), " + f"{rc_f_str} (Fresnel)") + + mag_m = abs(rc_m) + mag_f = abs(rc_f) + err_mag = abs(mag_m - mag_f) / mag_f + print(f"magnitude:, {mag_m:.5f} (Meep), {mag_f:.5f} (Fresnel), " + f"{err_mag:.5f} (error)") + + phase_m = cmath.phase(rc_m) + phase_f = cmath.phase(rc_f) + err_phase = abs(phase_m - phase_f) / abs(phase_f) + print(f"phase:, {phase_m:.5f} (Meep), {phase_f:.5f} (Fresnel), " + f"{err_phase:.5f} (error)") +``` + + Diffraction Spectrum of a Binary Grating ---------------------------------------- diff --git a/doc/docs/images/refl_coeff_flat_interface.png b/doc/docs/images/refl_coeff_flat_interface.png new file mode 100644 index 0000000000000000000000000000000000000000..6a248f1971b6f77bad049ae05f151de921b538b2 GIT binary patch literal 56457 zcmb6Bby!s0_XiBqjf8ZAbaxCXDIz6ENlHoA&@J5{NHjDRdV+{KZ@8uyE170MB%IQEg94w&jCeG$??(XiKR`xb7 zrY4T&oDR;G83$ryaBwtmin7l&Ju>&^-F-CQ-1Hm|2b)=p($Y&YH@rBc)qEarnACt# zg7yKb7Rj7UC};MD0Qvc^A6<*Ajjumnjx+|-A-zD9#^!*RE=f$QZH9%8V&prx>o28|LOK!+-BM4vxbAw4c}O zdHL`2jQ@Kz1X4BgDv?q7j6)tw-qC*(vj%CLAKIBHgMUB$jTRmr9->VTe#nFK-vm0( zZwzKADk&kFZ5M(&Crte}DAyl3U%TY zM?tJHopH+5mA?GM>a^UBtdcEO>dpcShk(7(|6gb(sxiSH6AT~a-zxt557sMb5J<2W zrAbo%zr?)fu0a5p_GD8(W%Az_zkoVfgIcZ`Vd z4K(%u9UI$xtVo%jiD?YP8HHron=-9eX}lqFQA9k>YrPbUmqVh5BaQA3-LaJYXWL^g zK8pise5L9IC|M9!Izd4}g-Z`lPdsw+{)sZ}i16^<<0|W^ib3$nnV(iqD_d@?V-I&1 zTP2?<#I9vHw3z%(Vi9RM{9A1!;`bFiao_|bCbb$_sU=qXi_V(6QTD~cLjT{g~%(XD(3ziUZ3Eb(@brIYg9_$y#*h9Xw zS8&2vL$_X933PIBbw1Z8zVq%Ed%AG1U%x)sZTZNcv~g7pCbu(($11&O+W<_m$43l8 z*7)q~p=6G?&OWN|-!nfvJb?Idygpf{Sk%s!MQ>_u7I9iZa%&VN1L(9t2$uF&tp;uFEc$CoxXrU65HCN3`Ca~TpAW-6E=;(20b;c?H$m>AJR z<&C90ChT?gJ~qd^T8{;6tf>lv{^MvSDMv?cWK1G@4h{}J_D$K{$**8%uCDbZinJ(z zV5lR9?Wcr(R?U*OLZuX8r)VeuV3r)gu`o!ihIgV6tCwp|&_5L{y4S=(j* zD+7rP313Kod`!MnAOZoWzRsaD2=%8Fdo&zh0~X9-L!~1I!DsqiN2l7FkdRf=W9c$M zEoInNkeqmw(=lr|b3nM%U1UWy6|Rc%rBScnJ+ARs0nKuEQ-CAQQIXAN00GUV#iskd zv-3rHM2>&xq8}eE+f{U5e}7b4n^b#8$2q4Qn2hwI1bZbW5pmta*PV#YMIy&Sx2ulA zR6D-({{9$5EWM>aHQuN2JHpB9z}diC!%x`zP$BNYRqSg0>G%#O5P&v5ymA|8 z*TylM+Z4c*`Ryh-r;i*?acv`Rr~T?KfJdq^qe2;30<~h}uf{vmythIz8-vSbmHdA9 zgb%0GSNE6t-od@7T4kR=;yB%}U0i`ej!!TDt$$#cd6E9)*fWO4ax1dLvnm zT%Pu5AhSANsBxCPnDH4?k>H%`cjx5OQp2AR2}22p*29vD0i-LS3o%?Zf@wf9dN|!* zK5NC5o0xmo;ID7vyH@GMNo8hP+NNe|&Sb@seoMwg^hyMwhm)GXJ}YUL2G1{52nkZk zZK%_Tckue{B!OpIM7a>w&w7E)B?0L}cB?#jEsBwICqzH4l;?(!kP`397?!%v=c*ix zCVOaUxjnw!5E2)+&vUi5eK)PYd69tnit#CUx^LgO`ITXzs)RV=18?onocaWtr?Gyr0;Px8@Lk$^)l zN%1SwpOc7WBs~OdA*z7Sr;4c=tmGkN9hMyC8Q8Vu(Wu#6C|;o^L^;&j8g*S}$!!Ls%?yz;HT zaz8KYh-X^#IgV+1k<-nnmEqePvhArc^zXr9XWnGH?ow-9m}=D3vVTygk!^eG;Cm6q zpASOx)B_LPm&Cd$;ymc*ImYWfJfC-j9^GF$ZeR^$qI2j4{CapG0*NtTPes&;#i^<0 zOhAx&@l}ZtvqnzV%dt^qY|<@#c-N7Y=eEK>0zXWCunrx{4IC;J5pg<3x+gxjyiYhc z6XtVX_l-Gkq&bO0_P{6H8Yq#fMO;Op(atSHzsgUJmoK={A#oC0NpL-)(H~wm)G{*K z8b!{}U)oZ?+T&x^=5*KvM^s)8mIBJ{(4ML0II|$t{0sHR{8zpftbSg8Gsek?Lut7h@HhI0P_>*%sokl5s86m4CN44Bj_@^tHv-k??DQ&zrQ_f`2!* z{FxX{Cf8eE#WLQnN2kYphPf`&5vmcJ&2J9BL#zby{PJE|i$~O~|HXeVxIhJI92?3z z4fF6`pJkV<%cN1lr53s;@vri)Xlf0HBxKq2t$z(1DuLhHkg`N}D?l({9c_3`GMU5i zt0@H;(|EjOt_)TNcnfE`r49@Fc7w-bA=UiYDs(Gq*I(T=)idP?&kTmwQgusGWsjvE zme0Im6s(%z|FZF9AFDoi{NSex%{d%(`F%hkYVzw5z4&8}?y(jhZ+2w%`V%axB)T7S zO<_Z02w7aLBm1%lS=>#i*sh689pQLONtOPfG7@Idey!u}!4Zcn+3g*(P+)NFc6EKV zZ7}KUGAJCo$e&G!s;$PJOHx1{Q>4+Qa+MY%k{vvbBr$tM8ugUu=?Zx}YK-Jk{qJFu z-(62;KoRHclA@%F4z2EMS$mC;b%c9KzEfco_Hh43;T6~g*?s+MU8~5{mYPB8(DD(U zspRu@ri_m6>@n}uI&-@IC+x8KSI1U3p8D-Cu;#%_LE2f=d9t3&r5pb(K&@46!E=YQ z;mvkvYn-GuZ^hc*DsZS12JGbjy=!dHC#z?PBo!c<5`k8TL-DtQ!|HO>+)wDC2kc>q znBZ_>Bh~qs=Nfi0^|BjFN-!*DzSxcL`m&E*uV0AF` z|Q5ZLP_snA{HAXE5Vkdn$!N$!j9$Q4rBVdXAcH3bM+cs2!&Nam;ja(0uBu-0v4_AQiCVq5Jv8!#D%jO>0Sm_NJSqS+8a#VnWa{ zSqSwb1FKO-oAul~P#RLlcwWuNt^FVKd|J8QC<`sJsgYKa{QJ6gMp$-#qY5OQ@kQf5c zP?&|w^t7ur@uJ;hEb=iXGq^Q;HOfeV7mf^S!=by+aOSi~Fly*ntCqi)!r9QN6e_r> z_=A+d4}^I`WAScmXp7XFC9Co1RljkM=!Rg;aEbP{{fd+P5+{J& z9T!IJ*Skmtd%Cyv4f=lM_jVf}`@vJKU1wJ{tw6LL{$5z!)kw09oR`ipEh139xa?O5 zLMPl}El6he%0n1=x%4I0RUAt)mkt~th{sX{j=DYR!MzF!u~J?Tc0DBDo3;s2A1)kH zgg9}OvEq(l$xZ(lSjjpFcHW*yb68)is)c&Kb>zD>hLj3Z%p!SJbCo@cd2pu^Z~i$l za6a^I1y-4`?LPu>(>Wd{QOvXSj~14WFv30 zk{kisavuX3NcAnhnPlXf3G;@FllQMJGVyWQQ^_Glj<#w&`WQHS^6C0*Z;Oe5p+~&P zeHvHpjoi@4@>YvuDmwDR$jDr!B1oQfpXN<0A7}2PL9*(Kh+?5q%beP2v6z+ul}?0e z%x2kFBlBz-p^^(4!>^^Pxrk|{h*D$xpy*5CGTQk4iA}%J4LhtOFeQc9qAD#nm*zuI z(D4bc4vO1;UZcPWsQ}zO+13*{C7uJ1;9g6Px9!oSY9;YEcm9t-HFm%3r)oK2X|10y zpxch5Co+f)4F}Cq_8t#f_IL1GASaJN!f{4-6}KO)$%9Fg)}?^guZ9Vq7T&K2sqvZaEQ#n)pxJ7rt0H^>=ZrAkRlBY+Y) zS0?Pqc)4yOK&7XBOgSUX-tMu~k`?ZyBCAo5nzo5`-X}VNPKfcJ_6#_mCP*=LB3_+o zTpjvX8F}vukO-^gc!8YY!3j#SqmYj1NU|>244q<{LJR&tzrGug+8ocMP!|1tTefIv z(k?CSBSSJ_ob(n65}yWv14sx4{;F;-| z)Q@duyxaQ6;yUW6QHwsdfeDTJdTM;J%4&g}R5LPQM5Fu@ELE>aOSHLmZ`S(aG!Exkk@bwm9^Z^k>4GB=ob5hyXQq%hd%>+UMu#T{8AVo^bsistw43^V1obnU zadd8V7i?=4X3QdKqPt3cVtrHcy70H%)Sv|)AS{;V-51xph>gs)M+?o&%!~yyTuXRS z*`?@5aL<_+cL0zi>>mZg_X!RRWvP3&-?9_gIZrE5wtgu z=bw!PvbWx8w6wG&aT!qs>`qm>9L$Miu{B_VQq==rK&BSe^;=y6q5-{iH>ecjzJsK? zEU%Is{DgApxALMYE$uf)S|rq6?G!LMt4b=WMLly=wqgwlxQn%=+XtOz1jtHDo1r#tB*>W@}2(~1FNXQZlF@8kCC@A*5iR!M$pAdu&joV|F( zH!sXizYhjqg>u*Lo9U(4f!e**+rg+IN@Lr)+i7rF0Gu9I`CY`ImPn9t$(@Z+NH%(i z7$4UJgiXkV^?GS>cDr@;##W1}#C0&>8|K~a0jKb}=(DuFFQtx}Nd=L$yp7N|WKghT zTC=~T;Y#%8geIxE6Khu*nxAeAMpH{tIXF0=;84AqZ*X<- zc|2O~*q(2^w3|xdH2A$~bg=y{i4hhT`@i$7?e5s*U#|?e>`}9iG5gsiS{9Jc2(UrNe7%H6|u-LGWh zmiO1lpxkxcg`LO;TQ}#sw(||2zz$CZ3)$d)sM6|pe{?>0jWP`QQMG*8h}nvs&d#O% z)A-hXPfRe;G(#kvx;42?u6_9EP?H@ML*1mvY@e2 z2VYo##m!r7j?~LU>+pOySD4i3)%zdR>ehGO2IKWR7Ov#-&6-dwL$5(0e-}E`)m%3A z9FT*5{`~Q|*)FDlK*OIcEn&OoCMxQ-x75$oO3%m`3h+Ij^V*A-DV$&3t$t@A%yWu8 z8_`t{v}Xe4??_)V$4DhLy8)xfmb%4o_TUF(%%h8oBFE)+un^8$BYC!!d>Mmcw|iAr zx0|czClxMRbbSj1UA;=2HU-@-ryrF76sXGzfV)~J-kMhTeHiv{4@CTu7YT|xg|I;r z)63HNSAQW9Fd#?w4aq^!BXN^eGOSgE+MjcXfixEDB5S_q^VdY3)Nk?Jt=tRNM?~{q zpwIrH6&zV<%J|I8frI%*5$_B3z`#J;zds`u{U7gZrjU`5kH1jGX{X@(!RZk+SSOMq6HvBvAL>)i2jd3|TX~-vLvV;BWhcyiw)g-R~Nf ze(P~VOoHV6clZR7- z+y2YJbb-Tbh>?-e-QC?uG4oWkmH{Za-J-p11KCC<7}!SSOy?LrON${qp?porqg{bi zrE9>ZjzGj2EFZJcudbIJQXkTg3tb_th|Kno(TB`b)UC z-1fOu`Q4`k)}i)X)?#3ba!32$W42r~-d1xOfCA#Y6wt4-3(m+5R)aWYp~S2_%n|)m53& z*M)E|UcQuWv>JR5i16pnYaa`*IIRJONY8?sCsr3ey!K2a|!{$_!q(i-swPElegQeR|CH8REB^T=Ft0@3-pW5o*qP8 zT->(x{uDkVKKnUg6f6?Qo?F!%iNR=6J`z#4awIX@;2sysx7_Bl1QP@dbp(jk*xY+q zS-ZJsrGGT|QuMv7a}m*9Z@^-{G(&C>HfHF?fhw!oo>{-Yz}q3`)<}q1ektsPXFPLnW?f|bQ@l1_I9~@ zK*gal1LQxPhqt#z1JsI+j*bKZWw!hKGtYj$;d5c3_#8~~QZ@(RN`$flmN}fITbY4 z%p|#em&^(m4(Ymt5@wHGj#OcuUA-U_KqrorX)lPKs!&D(8KEEK<=C<9&$}}A2XDW> zLk6Y#klR<~3_*7Ny7-gzem>}yngrO3VP1g3uX{EXJsCXsRL9Ghg(2uA@8==h4!8ON zM>E_E{C+UrRO-rD4idH5s`UK12>zFUqPS{@M?7EWDj_8U_CuTMz0;~Gt*l@o<1fjF zg!{j2l~Vo5H{*%KKwIesiv;(6p{ago0&J)QaM;3VBv?v!Hb(D41m#SpT_+47D!x!Q zZr?zM%>&*e9abq{=PDvN=f+f>G*dvrKInO=mCJ`%e9obO<}TpDU2ae&;U>J4Wg|~#BY}CB1}mk>>i*h=+l08bsc;i z#sf=}rN@yysE~NIVrmC)V%J+xWGsmURv_pGa-vuFsG~V$L}Ph8cuu!a{{0zvyn0^2 zeEEq*0|n?+W?4Y#4WC;D6$8ZCq3^Jq6_0ZxRzjP2@@N?hsqChpA0O;zqc0BoV5*RfL z6p^PMJ?r3JWm75Z2HacO5rl(vW+MQ^SHwkG$WQp;#q2FQu~#W8qPM zLA!!2I)i4}OMwqTBs$Wvj*2(K196RAKdpW;cpngNI#31sQ07&wf`H@sneU$cWn1EkqQQ>xe|y z#;&8*AHM$bnY`l@Slu#hBJBxd@ZRaX=Rr&;VdPq9lQ0(O#+I_kMH#v3?r|z{Igf@= z(LyIHvM`{WzBeVItOxcNdo$w|1{pv`1A?oD#NCk$r%}rQ(4k2LZBrN4(P=wnlv+{; z`dy4H0LUV^U>iuP(Ri;>gY75bm->nk#umVEo}~%(B!!{OH-S_xJ-_EY7?eQJB6k+U zQL@;P9nuZZlM)}t0t6aje1)Yt&M$%0X*huz9b*nhb-*-n*FD_l>5G&YN$&ogGAfhl z*Zw+MSt02HR8qIS>4fWRFu0i-J8~enFGq2;7=uy)RCfwJ7$LYTzv>a@lULuu=I3Pw z%TMO9_HXc8VBT&xxB45N?V{Bn@#U=}WwShD1tv?$#4-u>zRKz)1+b(f3PXmcO+K1D zK%X!6T-|gQryOzgjQu6aYxrhNN3N-g=bJ;0dpeT>4i2dz8X@A(bZcx$$BI>ZdV0c! zB_G9xvsc?p2AU_N(Q1pW$ohyAw|Qz6v~q3A5Dao>QsK8WZ%PNGh|SR7NtdC6WlIJ! z1Z_{^Tl(<`$xRGswVKI$oNt7L0nt}3;iLk|kdkMl7w{PYm?SP@bdSNg4FxU81E{eG zk9C`rai2Acb4?ZMe8zY;;&FBAa50F+D=ScsA8rp#t7hKMHwb_lcPPnYGYbMd@U#=CKeW$^=aVUUUdH+}4YiPG`V>P1+!NAHC1){ByB!PKrh66R9#ka-o z?l)`FaoIQ_%}nRhUg-BIAinPZke(IQbbL%@Ozuc!`cJC~!GL(xDV;E-trBW&7ML45V9L45bpXwndX#ZwxWVQi!zs1Gy1S;%X^K%rg`jGvxl_ z&IV`#=`BIBCcphH9InTzX%6}0U?7ExKaR^r+w}LD%4%!#CYj$oG2}4+gy-por{LSE zI6neeWmCgtj7{O5Hx(aWvc4&ONsF=G=K{vqbMP}UKD{m6{F#ownK*3=?JIk}v}AKV zGUw1CF#M}KlTc2|-#ect{t*xB#9|s3q%tUw;thHc@iZ7K!)iq+sBn;ZYg)dDoZzn? z$hp$p(-=reB04f5jjSm-j%C1)ztx%37^>%-CVf#!h@Bb)$(Eg)P@w(=5=84oH#vX= z>T_amFPFopH~RR}|6o0x@gMnNoF~3dI=Cg~S%lkuI-TYAZn*p8;ETJiuLn$&@|oYT z0rZG+fJmET;1E!Sj^>`rf2EKc&trx`IJthV_u4I#&nIa*kRf4XzrO!8jn?P+2%yG> zbpe^_y-N6e$Tqgg?bkzT3XuRkdXOPS5S`7)6)mX$KwKff?w_ixkJlGAl>EU9bZ1b zfho~fM5~@PcK}t9RQ4CzbUU5>{R4w-FjDIMF)}49l%^LjiSAeJCJS%xsesUdXZ&QT zqgax@jmiCqqR`}FXs`D_P)fK3wq%G6P8QopS1!UABuIjWdB9-FVN|sl<+2HC&6Z~3 zFa2hWEXc2jr@B)*{V!;2plTW$!JJR2Ue4>kc9_Jjd_CFZFe>PC{1 zrUPjaP%>VLnEG^m-QLc`3aXQ8IGUIhU52!B{v9x8*TdDvB8kMxjEUj4Th9b5bm*1MB?-rIr>p?>#gWj&P5L%ZzeEdI+} zpsvNx8p%4aHNjLpp?B)sHpmB+ZrUrxmnfgs22mVNeAYLCcn!d-_wb(i%47t4v=K~G zZkkx)Ff2JifZv$YElSOSA^pPcW1udkJrxZD0VOX$(FFvzp>y$A6l->llYzm_8c-?E z#yWwD`3X^_@}vm#N0`c55;{o;hCOk>rkVoToGpL%-GU>c0@?k}XY3nL}Y z>k{9J8aGtLm|#7w80oLzf3+GMzAKjj<{x(=D%_#dPdP&W;e4P30Z@}hU5qdv-ypN4 z-mH*~l-snSCw*ggRy(jDnE_}L!?IJdW}6GfUi>7WR;pbxSLI#|i$~m;2%@+$rX=7H zRR;b=T~E)ur9_9)SxN`2OW}-PKKG^3nu-)lK21JUh(dn2cK{{t&eT*f05er12_J}@Vg z$o2N9FR3^)7A2HMnYsDlH8hL}5Egi!_G_fQ-l=oCb9c3DcAKn6S zFo4gC=xhrOJA{xeSDWjbp2>lJHNVE3t_<0a5RRk*idE8Ku(1Wqu*-1UQ}qK68Yvff zfNlkr6s9d2MtCr&BB~iWAvaLF-L(ze^TEnu7>XvJNujY>h*glDFX$x(BSDlXA^nP7 z;!+z;62LC?W*ld&LaG)X2G*~46~`svG-X}u_E1K}B_1`)Q+?ekba)8&?(WXY`sP0s zvS-Y^dPc^Rm^MnPfoAfef+1zhZ7wwrlUb#;*Y7NH%lt}8$aE2_{jPmWqYCGq?iSA$ zd+2misl!44cpf;W1@6D#T z4w}6-;MpDjQ^{EXc0%$G0wW*K5N(nUlI>`aKiGiZsV+&lFTBiU{Qm1uy}1X z2!K`Qoi&}2!)?0mttH+L@yo~FhX$ahscl~CD{*El=>}FHvyM0mfx`j4dbXA&!NhDk%Ipe{E4 zhn$zHlykn*_`%EFMICv~Pk9!GIw8&nx}aV;_ou*;0?4Q+lTXbJ6<+Zfs9{5;+zIZ1R zqq}DJzMQPL%fJUL*fZC{h)UzRj3102V0*p%9OzwhwV zZZ`RK;mc)UU+oCRKIi1v900N?8c<7cT%TmABgV7k3evLWa(2|M;$0Wu?(D8@JW08A z{rW;IHyX>;a6S-m=hS8@SG%%&SZ?U7$lbznAo~m!h~Z-HV0bV<2r6=dpG0d>$wnik z@0jjA@FQaNdEgB<{m$q~UO6b63JduF0wDI4_B>SJ92}tm_Xba#2s> z;85HTD1JKNwB8pX7(jn65Rz)Rq*uk*iS<|cbSiY>^{Czj6P7I$qSJ~E^qnO(H7vcE zq}Ek^Zu6Y+XL?)JdN!~d1@xH#3is9x82eB-sANISjxYf4GkW**LGQQ{_z~v;RlH-5 zhK`|HhnNDck(U7P!V|AV|S7|{P9?|_uDrc`x4!aRh=U^qW`whU(g zH8K~XA1g(c?-5=nT%v_-Z-D*S`P6YQHT}<1{|bur>2@N>p9s={CE=d59x-Y3iAH|l zx2+QyL4Nz;`o~1d|1j=bo!3U5zafH9I=iF+y=)IisGBOwjRGIY0UFF-vONk zC|q1ztIm=#7@iVOxesfQ6CpWhdB z8yly$tcTCJ+vM8J;9Cr+Sxz_Ke31|!RX4$?`uORrNa%6ZU#u^Y-GXPGrPAF>8V%6} z7&Jsid?$9;{wCS2;sp*Z?hO|YN)-BrT!2vM>|gT}ot~?t%hLI(g+a{u3b@pd&d&>h zIoK=ZY||ba+Nkaw&S^XVx27kpNA61BVK{E`XY*eLO@cBY1&-cuqy`14$A@noJOr#r z&E7%kL~rUaF^P)!vT&~B5h-%vv2}*Mn)~skNs=^NK4xaW%0c%1f!r`@{`RREwsSum zl&c~D^8Sxo`&>&G$h*4gP)q(W1aB4mD9+#(BeRLs2*bakm&E65cM;ka4fCDT~`_{$O-U^ z^^U=TQ7Atyj{XiXq}ctgp`9k%tAXFL&rH{F-PadO{iMiNWwG_U(LEPY*DuY|?`j#< z@A3X}wt|U)0dx$6gn7Gw;nW=;Q@5$WXp=xr{3Bp91G>y{c#}5T74~}!TbZoq?s`*^ zvJ!p00j&>&HFiNgR2oxF%_bMNu2R>2+fXkf!}T{*n@-uAo?atCviR5WI6Njzg|DT< z@77^#6)7(rR#qS*|H>3$oz9Z1bnulQ%iORvd$#|+N5Q8PFZKxF((ZCtED4!| z^X!%KT95!Wt}8R^>GTJ;^}qE7+P#^MgTfDX=Zm~?99KWDi=pMyf9x(t=3j!rZ7(#7 z0%E&)|0Wch{AyWJ((&M={?=bZ6$@HjCq7aLET?m6&^iL4y#;J>uM$#7CAt5UR!MGp zvDhir^#0AHCI6ceY5qUZB0kCErxElo!HLsJ31jP3V=U)xoY+O{qTBsKuIsjYfV%4^Qh$C z08jv7-m37MXC8SwnEyKtVT-2}u?m_{EPzIbaZ6|coM27*_vQ#3fnF(>QqkXEoT;y4 zNV)bwbbLg+pFjoFt73)R>yxB_g&j@1Qo2~h@y=sXfSNj00P}b9L?!T+9U%dqrg-t? z(O`{(4oPeLQGjFUB>8Wp2}52pIA;#ptRHq5F7GyRcZh$FG{|i078KL=tkUHT8Ps27 zPAG)riI+_7%RhnQl^oYuY}L>jaE9R4udFF_h zcHACj^7Fg%i9{u`e33vz)EZy1s>@+U;2Zl@jXEZpF9djr=6-MHA9p9A zpj8H+ntB*?0KtbxMEu%@I)i2;V8lTVsI3)BNJxN*u%2(D2}|JlW$)4V;cMWmU?5B( zrUzjUvF6=e04vo?(9i-TH&hx4iEO15&TajKBLY!(=~lm^ey*jz80w|?(WHpaFq1858gN+U4y z=QxmIXT2@B77wqoL6xiy#1jCJ{uEnKe<%lwLX*cimd0n9D(>xO5zy#$&iW)Q=qZ!R z(Z#F-;?m*aa1XVAW^LXOaR~}Pf)6NY1>uQ#EZ}~?y$S3tajc~l$UX&3?@Hc8sYV!Z z9Kny&IhhLDPUpJFwI|+PwyG6kG+k!j`1IieH;=E)Zw?^r6&SVp!BI&1Wy1RZ^b;ib zl}%SVBVMKRM*!oSxuHFenoE9*yp
Vz*-hkt{RRTl;_?r5vOKj^apq&sT9)lPWO+YbFegrlof$E|n?SJLNR zUJ%K<+$VwM+LhU}Aq}8&j>n`O2`*j0rW+W*PA{n6hIn+={j=0^00iYNn>!0kF^l_G zoU>S2#mErQ-E_kT{{Vzrgc3M{kne~F$ho8dG!_gI^6lWizUJ(E1DgE0&3H|C(?CaK zFJCQ6Lz@&ZPXMy&b>Ug*w#oVaEHGb<@b4!Wc-31DXAk(+0W;+%DJEG?=fU^NGlC3N zEHHEy`ULu&<@>9c?B6iZc#3&CEroir`laPn2`Lbc2S@4=Mm~AZZ4LpBN5s0%obDR| zw>1|QrkegC=mf+mctr}tb_g^eUk*6z{a3RR+0^=F3?d2HicA@#J&1t9+WGZq8nt8& zNE`r+y`L6`0_=E|B zYtT{(Of?%@e+mGXytKY9*PqN`HC}>eVqyZ~H}HIS3N#nK&UByYCBVo3IFn}&`ZvKF z@M=oAY+B_nO3BBVf$8w)`#&o+x1fzIYTmwe-wX8Cd0xy|UiTY4eCT>d%fJuxEL~lb|BQ>9s+L|P~OdNsb+zw1;+kNM1&N~#m?kc zB_P#**ZW%A{0L(3g^J41ozLULeYt);@2m~5Zz-qonjh#aUhL26);lxw_4NS@sam-X zPddMqH0Yq@ssEd>K*zxmv$gde4k-O-*w_r>;#ok@zQ!)W-Tq2_6>-u?nOOR*w7Me3aCUXR2HUa_yFt1{KeOud7IWq*-uIgPE z#ye9}@U(Ox7=-y|eILc{&&D*%wDz+viwXq|K zsCq?}V2U|ZUYuV(|5%VL>(7x4V$j2h+_RN#kA>$2QNMajX>6sW$5HMixL#wd60Qlnou4i{(SGygLjQ42E3W%Tli`t~!?X$>l&{N@KI)=LHB34)fwlD;-mE8vl-r3^c%go{MNLxp6AU z6xRH^Xi+LmwI}Z*{yt4jx{SlZ!C_)$jcgCWtnxU1_k>EU7uJysmK6;HW2GQYf{uqL zL8DmZyICJTXfE#tl!c;-N|AI32EWZDR#aCz@Syf)i+lgLx2frWsw(!YJBFMuh23tk zB}>3&a>-E`*60rQe!fZ;q}-@AyGS|X1XKre)eE3A5%k~70gDspQpX1C0vb|NVLkk? z;sF3<9@B0tkO$m`(^@kqg&hMxV#OmRRRl?(4)lOy60z6qyONTUrctA!p!~3(hkOCY z*fkYQC12@^l9iPmsj+*T?XiMR#Gxk(@?4sb!-7GM1L$n6bEK`dnc7(T6$CuAgvDcu zU-hbwPCkCxk-waLtJWg+(HgR%qvMrNp1LMiMlfi2z+8tN`|v00)aF-v(N~V&Wyz34|7br z*%}7WQVrmLx}e=lGc$V7pegQs@eNEoNF0iwF?YJjlMnP1)Y;6dtE+>~Z)9LP29qEU z4D%u`8(5$dTC+l*l$x3vG-AVAy5R;hgqY~p9%n=kk4hpBaat97?Y(HB1bENA(DN=C zb<}W?J_kVyS+;Q0WS>i^rTc6nu;!j|!f%Jef5uQmQ$)+6|BQgnkTHLK?=+LazD-YC zZ?dpC!z6k+Fz0PiHM6&wo#uj3wwkE0@bGZ!0QOS}NXOueDF)dO7$7S2>)Z2W!s=W$ zpMbqX3+%R_8@Xs_7)%4qTGa3wmL;36H!_;QS3QXNT*ZR|_=m%y1Zc@mdRLD7;o<)F z6r>4+->?7FF8XDp7xLKf^&x-E=;;}QD{*`W zvzNpPT-kqb#z<^-SW{c`N?Y=i*5iAJl=`ch&lcJd)Nl=y+!XLP2J`o?g%2c|kqNxk zre=_y!rNE)DW*;bUB&asJ-*(;(^NPj{eX^^frBTo?;`oyk~ggr#D*+5lsWdy`;!Sl zg$x|IA3l8O1kc>{-Lyycc%_TuM!>Ick^BLuR}YAlTI_4aG4V!FMmbt+ju%cnEA9N=~sq>0CE@KuHK}^B`o+ zu5LICn{>pdfM%*Or|dBLIbz2J3jN8fo-f-_sr2q$KCrKqeEad^$JQxS^gJwjeEl-4 zif^SIx&Qrq5Nf#Bd3GyTJ!xsxk01DMZf@^DvL$j$?4xhFTERT0*l4X@)30}ytL8Fl z5i{@yu`CC|3Yo^x@7mny=4_k57wl1nNFokaICWqrC%9NkfUvQ#sTas&JAwLuz%OH_ z#kcYABk;^BkJVnEZbpH{EV7P@iz5a>tU2rgYFy3R=HqU1>l`psUvj4Uz*jrqgJVX; z$!DR-b7$J5BeamjhMS}loEPH~s*e@)9|um#6CeeB%$y2FM2D5WM#3tt1eF7xa_o{- z`g8t1k}2Xtl+la@b@(oF(8Xy3VPzM5u3?Z_g}>7*Dp0uazGjk}cabM_&1U{E#B@FtUw1?!lXbjwT^HElg+q@it#2{4b|yA|@IG+GD|OZf%gWZQI`gMDME&z-D~k4;3+TnQ_Pa^ z^}_I0iD&>Nd1OxTM%C=hc6%bA5%K`9q3xYSVaZtg>V!5P z9$s?$!G)dd_UQ0b+#G>xKq{Z5%9qU-FJ7og%7Ak#r#+icUC3&b9v+RHoZM++AdUYP zTtQKMbF9BLF$#i;is}T8OlW1Au6hby9JM47IF#|?Z$=M&Z!b!x6dy=Gk^Kmsak#E}uaPtr#UmeQW23wh+MI=DFX=myc+jzlS~!PaNp`zP2KpLyfYqJXkmPV7Y#X zsao98;-h5&pQYtMNlr!vxsd&w+Uy3G6Zo5ko11&kB{(e11<(p&;^GrLrA|BJrTlie zW|3@fYuN4W?YA%bp8z<4D24XxR}kvho+Uq-h~xjm)ms2n)kgcnl(dwTga}F~4bn(= z2?j_?C?FwH(kVzODBT=Nx}_ThB&1tFkQV8d`qqB$y)*y$X5Mk;74_`1_p_g8tzWG< zG$;W-1OFf_CZ_+qm;V|aotmqwYw-ztm*m;Gkmnii*_q=IFE=;yE%WMEmUnNNTRkq) zqww%Jy|=;LxKD-#$9KQy2k_`Q0J4Av6?x>VXRvqN)hX_U(hC>z|{*y5Krtj)}1-(p`K+R$H8*3)LhfEnXU@1EUk^1y(mjjb(`VFFoA z&_giQ)YlhWm!#s<`vy*(NB4tHIv%@Ud=5PM{He0i`5lIs2>0(e&91#$wJ68YdcK?) zx~Uku99{J1lukIsbnxnYJ1S~knF-zcZMBl$b+U-{a739e=!X2jsC0{va96QL1B#M!Psb=h>ALZRMoJG-fDxOy*kr+*Us@!x^W21vogI zA%h{PPXFEOpI=@EAEwj00<<9->x-AO`mCQ8DP(Zlez)Nlq8vv=cD>fOUS7I;O#huO zt#Y>?4s595k9&Tj0<|YY>eSR?|ApPRS5;{pBOJ~(oM4%xy)LZfr6BNAc9$;eA>~QO zqE7FHh^XjpH+xCe&vEgbAYz8@7qMgn8JRIt&dXbygM;0CYIjzLa^TGWW6=e&_`L~Q z)x)h{#8L&?MME>0v9lMyA!-LLOTUDl(+JF{Qo2@F8Z^wG7O$|F{O*{cC7h7(y)_=5 zMdGA)P0HL))MkYTmX#gj>h#8pfDnUk@NqC*yU>0Sm9)w2QOz6aJ>`cX4rXiC%r;%WiweJ7#EM8Q=3hd6uYkhSLxKUD@2! zSZ7^d2lMPJtp)=o^lUV@>Q~+6JsltBLs(oH2iFh&?2Q~IChTl)Yw{btf@9IJIG%w= zv2%K6X813>Cf!hV41GzGPs@JSp2X2P zVylRW)stsmk6hT;C{mcY*t)!kV3H5bnzGwn?E2Fb4+ZmbX?3K3+cO{O4>a3Ov6rmU zaa7pE#Kg`%22dH)c`1INxx(gC0dm|pXLt7@1KdB~$cb;s_9b#5pz^3D zhHlBZu%Lil@@yLux-Fb&ZdaePoleK(5IJk{6vA?SApo}&BKl~3b}>>ka$*Timg1g0 z{bET9C|9zCS0bpE^`kY`1~YN&8V@MZ2&lL~QLJ$%mOW981!|3x&H632XHA}U%vT*k zg?@Qn^MwpJez=igUiQ^?zdnNj6h4h%E!N=9qtAo@zG2KkuYB+2Pyi+l(bH(T;IOa- zs33BHccwHeDlA02>h6uHr$<|eqMQMkuPD+$|&o zhjn*+shQ2@Sv-!&wD8wush*JiiFD}m{ zo&CDFG&{}Sn4l`}FZ&gu7R{vGHSE7B{u+Bj;mgBoPC2e7SSbkr{yb{(ber`$GMx3? zQQp`T?8=g8MaJUWk=69@2n!0TQL4`rB&QBZA*Zbgy6+ z=!jr=;OZ(2495zzP!P&TBP1jQk#F*kA3v^h=A^x636y8j>SP5goZQHlM-m13UUbhb zF3JH&2@t{KqoFSvss2xx+!P2|FsW{F@vR)w9b#>z=bC9Xy}YI5?&R{FYA@u&2WA!)mKkPccr;KJL6bdlZl3q)<44q?s<+T&FF=oL zBYD?g@vlHX@>Ew>tJpvUg2#}0G(A22J*`kU6tmg7#N7NB{th@|@Q;l$%VWQlKHwMl z#2hdO7?c^LCnBEExjAF7SZ9X2xy$aro&lnyo^y?|vWdOK*x;92krX z0wdFwp+Eoj?ORs8kJ7H`{3emt>0N)KmJV#~l6NXGyJOsqMiKS(Vw;vf=9htP@gF z!rj{13YM~k)m0_HfpsgbViLc6M=jcom$5#!TFyOOWNwoDcmF45^`nUt|@Q8^!`}?;<6jqY`Vx(4{cY0&7pr5dY1l)FfoQKgF?iZ#yR0d6x z(f#1}y%#6+^z?3k=AS?s8i~yF_!w<=IRynWI=WtH4%i;&k^tfL8E#NQv02ff8-W!m{O2y^`N7 z(-W0uafxn6x*vPEnems~elLqUbwy=ZYSqTC54?~apm2h68-UJkYaxedXr$n|e`t6k z`w(sWC`iJ~c+Wih4Lw`sx}dGU&@7s6frN7?`JG7fp^OHU2~Cql&?bPqaqsj-06JK^XZ1{4}8%}T%YAju{2I}eYQWE(OL(JFSQ4*zsL5^ zczgV-`O>WM(%oRg;JTsDy>qLBuevL}n6iB{`<&F(*vyk(=?3|zejmy!GiAJx*17gV zLW9kDRfFkDiX2lH+0UxhZ#h+>_pGri%$Gh$l7+}c-e68q;VbOpaeus(u8Y3t(iw4u zm8$gQw|S@txtVL@{C?i3V#OZbd;I6ET~`VZ48x#%xb^3QvG0nnfv(J5`I)eAk&Mxg zS#7^r2Z~V@SFF&kD&6Jq_NU^y6AB-B{qut$0c(zks|-Y6-~I1LQf4Qs?wC{(3*7M~ zdyjWKyYpf&Pu^z#mlAKMDPa+RG>AtOTeLVF>hH*<~clhQ67qw$=SuzlO^ zaf91H{*&3`+nh#TRWB^X3p~4ER7L^l&;)_-!LLJ86*tnN)3v+k?3w ze$|tiBNzOi6y};6(u6xS{sk%4#fCbebyO5$(S$cSHNT}uz9tWrPV(2ld3+@?orCcz z+J!W}DBM%q8~6Rg>F1Y!WSLb}m6~8ib3;F;e48v}jF-FCz`-4l_Do1~x=a}JK8}^q z0H;{RdjYC3P1pk%+`PY)ihgx?GT=w<%*}ErMS)@trLM!J!!{lKj2gy_@zZ8=B)|$s zaWIbl_PkA|+;oeBaMlAlVvS#?MTo5&+6MV zF0q|cm2cB!_c|CtdImZ`79@sA$z>`?Lpz5eo%di*jp=NQ_GZeHzv*O zT=rG*#>H|JULAW2Q5Gsn>2`li5OKQlH6^t#;;$@lHaAIV5`~ZjrR@F9l8=mxzHL8q zc=F9lT6wxqp7L3xe`9|Nx749LMA&Tq@=cokkna4BZ^HU@syyA&Zjb+a zdSzWk0)N!-t^YO$A#u)CmDm8;pq8Js$;rFlnBggOxlpSs6t^2K{QUZ>T)>zeU=+`> zIsSj6e{z(iEy!4q18%=Lr=rh}zE9NR-8gM(jz3XWE_=;t8HYZ~Svp5nS7wS*zEDM) z72ZZt@R;r_Xf`+PdXC;LjlDEsChH;ySB|f7-lAL}esAG1_XgS4#-G!9Ssk|OmNd$& zL1TH?_|wzzm{#!E4r$&ku1X?H!L5l^^P=ZY@Tw+$Shf4EIEgKRsD5;?*=`ow!vSSj zQb6HoepY)}Z>AV8F0yq7t2MGg0ayB|0d?BuNi(edg7KGAKz?cLTM3l5-nCf*SAw|x5^zC&1Ub}6K3Qir6H|6(`s1qhYNwJhda6d%ri2J`=POs}HEAt_nS9UuzFJnlN+QbjW76F4 zRb|6?Qm_aA{1s?U<{OW0GBK^ir=&bQoXd{e?c;{k;g%XOLu~fFkCCp%0nMV3pI@&Z zCEHGPqJs7!aJ%NUf_ua_u?boneH@nRc(u5@h5&B<(1D%1_$CUc!Posk>M4c~rVH}U zqA;z}aUbJEGiL?)K9yQb(MR4M(*M+3A3uI{D*QiuwwDw6Sgu_E%|GbLP8z+I{BI8_ zey)h2N9*?bzv^GQ4-JlmLl06~Uw83vDJ%v;JZerW^DKdR(cujX7Up7eZMu8!% z_xTgu$UDe=fuQrc#UYMW^Qr{pQ~V7m>w8=<*!XW#IlT$|8$aZLzop}z)hZ#Y0>F$Y!p#N;H9CNGc_6zFo)EWHa@tSOx;wPM+f z&@|8nEgzvJ#C;q?r_+l4_zz9X(oupeWHIo?l1jYfuiXh82DYw8u&rg?G$CV=S&+H@ z4*zR|kbQ~yMIyzHT(IGyS zor5?5k)=neNJ;~weYb)#!Ez`I56J%oaD$tU73$uzor;F)j?a3mu;C3q-hCCPVnLOF zO^nxM8~8W>2aNLmzjebcp5~2T|3>#l5}*7XpPR@XL(lB`{-Iuoy388OL-sJ=V|A)| z)+h=ARRTH#`K^xF-r#I?G<|Gj!J>cuXr#8Lh6KcKK$=hsy`qD1Yq1q`W2sLlK@HJl zK#z=+tw_WYRAET08K-{L^X<7-$S5xbF7Axc=!btXWhhFFI9aCtKd2#!JY5V=4&$@< zzABtw<$i;>GxcwqeTI11*&#t@dq0}t6)e1;$e!0>w6xclR`7!0|KJSBOm$4(@4xu> z9yX8BEH*t7r;qAd<5|krSTvHAU_}7}%Ldqi1;`?jC4D4v&OErF9R9tv^#D3W2pvxc zQP0H01R835J@0WI8K_HG#Kn_AzqC|T(RC#S6E5@v)9ROFdPYDlWkLp0u$TaNgr|D@ z`qAc{5C|oX3ZMK5ae|T^-6~-&^Y9ZaGW8HFMI8Q}Ik|?njrN1}JdFoVvnBgYELV`m zU+=F#llV6CPj1^zt(EW9oK4w|F@1?7PH2JU{58O{B61S|^h!NGySH@zpBCWgRh7+T z8~E$-aB-2Sx`xY()0(=v<>-=H1!$Fl%I=(;OaMe+9)c50s;jGWfF!Y*ZAgG-pVzp> zx50h{;!&fraA&UL$D+j4FQ}a0yt)q@n#7mOG1IMqN{F?!hKbt>$# znLOnUhz^kZp;NQKQ-cvy>nA+n9ALz)JlaDkSDUyXnK-LWrOg!+R58eIH0 z&^0~LH!z?BZc+Q^jT<-iw|+gd{ng+W1-xwj23+SL5Hyt=(0jnOP95&=XF0ErS^ezn z!~s^ka;E%+YC}WANEOH_fr7QF2X!DMy&2h-7}m!EJ(0fwm-q0CrB-DlknX>KL)LgW zjR3yBzP{8Xm6$~z-Ke>pT|aA@XGT?KW~OeLY2Xx>ZuxD5-~0*v3)~l7Tb=yJA<%CX zYyj$84n$(5VY9av5c>I5#Pt8-9_5U%Y)nZX{Bicv^7RA8pOlc&@FEjt)fnWmSD`YkS)0UODBx6r-~` zludw!_XaPyRHQXTmNlA>O$t;MmOrbC{v%|UK;~&x4;a+Q18vWZlIDx@nuH@jSss>| z%2zg;e8mPN?sFA@Jj$T889fArikc1ID3VubY-d+e1%R63L9IrP+2Z$4R;R#&ap+g6 z*#agX4YwfWw1MC~!&wUKSuib$)_cs$+Qg zk?q!WZIrJsTICEvr^T@6RE-zve$7@*ciPD3XyjTVC+W~C*QZc6tY8E&KSp?GL#`6vBD;-AzP?%DS zXbPEH;IIyXcQ6z#I216(5D;C23}z~Z!6p(6dK~^AiU7xJ)O_RwPd*cD(5!d03)w(r z7u$A&fq@%v;@hzLZo_Gkl;;Tev94OGxJWxPPN1zex1xnYnlT{8p+ImU=!Ucz)-4iG zzd3a*dD5pArZ(Q5t0ecNw!5NE8s4JhE|May1}EN>+f{52{q0OL_TV2Z;9tL8Bs~so z!W==-(cUHjacBgmz+iIx1?(z3a0WsZ+Kg%)UJ!UZeHtHVSTbH9=;7f(@Ds*jaCU%3 zX0**40<{QE>$jS$O`ZWT66gibi4s{$V{H;Hn+60PznZ)_0W7tvE)}vK>uvJ}`j6lg zgsRU?A|fJmZ1M2$D64><(tY3Pv=$kN=%py%fV70bBjUz=j=z9-d@1z-+?4XL%&eKb zRLd)g`oj|q#6Q6|- zLXUsy+H2n=p?gqHd_#P~^b9nw937CcbwqgxB<#}(*lz4h>|oUfoVo0jE;^XysSW28SsJSr3~V(18uXgf(TCG4|4uB>v&&Q>T@uvLwN`<%F_vt z{Z$i_qr=0hQX!iO&B_lnfgm*dOv3cIGD(fa@y9?PcAxPiA|V7tHpotO1vhuy?GsoFb3k%+^*NudD8Y%Vbj5}vq%y00~ZBBj2CWpTvmJ#$zKdGRp%(c&d;NdK8l;D#i@ zX`S;r$QK_x^MP*$NG+54rA2SNpGFQQ$I4YW-0kM|<#s~{*tEJ94*{Mfs4pog2|S1C z8xiN2>G>?TTHTmhvu zhr=ETW~Qbbf2&{3l$Qnt1mqspbGevcX$p6FgS#tirz_J54Tp%X)U^?W1ZF0`zPyZT zU%yU9W@-YEeyzL{huh+HW2ke4Z?pel3D8WFbl?BpR=r3np=s2)c$iZ{PeNjBA{isz zfDhX8>9z|%^7SVg8yiEEDf5mF4k}9#nIW(UdbX?);{8cDf`UKTv`k5#6kr?Q^Oib2 zWI)Owy%y-38Rq*NuI)|Gl*T}2-$~)T1(*UbSsNd9uzt6=O>+(DK{rEck;biHHiq|V z=KyHw)fMQ9QJs8HB1`W8vhZo=s1D%?kbGXN;m7BzaA5_@bve5*J^o_=k7=3=hj~1 zzk!nqcE3ga)Z-?Xd)&&hqNj}b3+vc+^_nv4PAQP^5ztEnFYrt4 z*BMy4H%aLXTQ8Haa6D8+6}CSGK@leh`>pHIjLX_UIwoL=-FT9YAhGbEqo;>aE@1bN z2%NFJ=Kw?sw1L~(t{T^>XgSh1 zgV*9Ijv)R{gW!EQ>HX%V3SAi!^`+9Gru2&O)fTv0G@sQw?F_ub=D8j60Xcd z4y!(0ucicO3ES2>?`v6d|r2Ga!Lw9tuc6RqeDG1)ayV7 z1FJwo)(vkmz|Rj%gs2G(8r>;ihXabaa@s`+0r@h586{hAg5{zAzSEuBf?ToHH!Y?+ zmhtk2CIB%h*w{ob#lY3B^nJ8ZHO^PkLhuGLybdhSuU@(UxFM@3;1M3A8hlf%lo-lF zIWM5bITu%PuXVTB;OTLZ3c8KNl>K74 zq;8Nm#@rhH8n4oY_%$;GVp%s8^B;(y_uB31!I*!Md*IC`2*xm zmRj|SsvyPzc$os2ZlF7m$yoqCW57MVA?d9K=sK(o2!Bon0mVw{twX67kytbGb2md# z!%XfdvVp|g(k{3g3Utho{-GkGIYp=TZVEfU2Qq`07E=ZK$J@ikKr;atONhvw?%$6> zQS6&f4fg!GZt}Ir3>=h2j0g!8Fuh5=ImBx<4IM;E z52ctxf)%dJWap{^;45e7FbSq5iQ-dxKeGS2IG!TJQ$TMyzp(HXnqSaS zVgW!f-RzSFBD3r;Zq+4F?I#cdD@eK zIPj^5y8>n`kkAp*3dww8(q-6MiQAJv5W6?L+*-*#rV{!XRR{zOu-}kP0M}r&lICL) zf(KmBi?;N{>7tqAu5ij&uqdj~kYm}{*jSDgUI#Ie1<>Ni4i04T-oN#d;RZKH5I)O) zvN+jFAw6Wn$fZ#B@1xxys1sy6UE{*^$(RYwBzmE&XK5tH^l6c4S#Rp$$cx0ULFwLJ zpm)oE(lhlGlum#m-h$6Jetgq0{D;YUb`iNFN5_DuYz#*n4TV#yR%`F2i^Td_bomv) zt0{BuJ3Y81Z^#E%g}5Dl0a3j9?5}cwMgv#$eWyUSr~Vqm0cc`_#Q;~AU{w13+MA#X zxh*7u<$1E0b9*=-fQ^{|=s+TRK1q#G1T0VEAwVTu_C7~sxFFGGb(TzjR^*+?p>KTe zI$j8#$>GiCFiX(wFv8vbhRH?# zq{IB29!>M6Bqtu8&@?C=M&Y3X4Hcp5{cwa@T29%BsCGf zz^fOuaGvq!U$DWx^my(Hy^>a(%S1*nP7#49y~A3p5Y>uk!_wQqAWmh0a^32h-|Ys=#1Mwyk{!%oQdg5VHkC~Xn;@ac5JA|JIa zL7fbns_9RF!3tWmznq?aQsV1zjwURLK~EAtls9{vk7YgXuU$ui!gnZurA@(m@ZwF|~nTw}>%H#B9An^hJyTo}C<~w=cUUdxX~@ zNg%=3aQeCBSSJoHqFFZ_8X;|!8!(9Ia&PYefIAOw{q?4P>q2C1_uU9U<%B_V?34nh zjzR4tWvpHdXNv(1zz@}4i6q<@^_a^V{w(FC{oj3%qC7NSWOtHMM_ms^z z0}x8TVW&NIMl?wxBuo#U`kuE9yc`YHvAv%DocS3Fx3!g1oP|D8$p|9KIDFXZ33&>U$ z_dbhITJgo2$YnQla%e7p1r%q{u<&&761Tn4BU=Xpv)EBnE z_)1}T=wfuzZoh;SCv_}(UBAP%eFwOsbost_M(7=LO24!GO1Hh#l&S}TLdS8RSmYT$N3fd@%>2|-nN%7}Il;L5xH#Y-+uxg+TRAE!dp zK?r*8Gf<2+gf8HV%|=>>M(#%wD7bR#)AKmo=w@r3=p_2|EJ4U z;hA`OfyvEXlh)1W%;}vTsdRfv3KwAlus5imGp^3&Z&Mb(DNZDpdO%Q@$C)z(vQGS; zRAY;t4xv;o%aOgAm<_(7Y#x1Qz6&5yL(pZ+WKFCT;wj)7A|J)a8CQL;SqJwaVxxNs zKr4O{35LG@Pr`m)7>5T+>poCKf?k_gL!-Yc&kMrnkOc)2vg~cl@6p@|^P@T@!z=@9 zk>P>wWP*s4d)S9>iHa&te}uhC&7}CJpzDe;k->f6ido$c1GSU2FE!Bc7;z&AENYHk zU(}sZ9R~hAf&1N{vvVI`%<|EI{5QqL2{O;xvxkuq=F&Q}%rYxD$#V4=NZ2+nt#`ULzwT7>4*x=b0BpqY(nWf;9#Vjebz1IcynC(bC$&fjZkD0#+CVUPPkC zLB0WwoF32$B`pgBU;*Yrf)>EVgEdz4UnU9KB9jRD0Fx{$jR6E`_W_>@owzGW=AaU+ zGf<)P1Kiv`IQR(!4-RL?_F?u!L`29XfGXw@l#5Vw^}?Gb;&UPHakRMvDAQVBEnshm zycKE=Xga=vTV%4@jsfbjm*F1iB7~mr%ds3^{WFB^L&?;X32;rg()PCHKYzdm^t^5X znIZt76e7P^19cLR2*!gM^03haf{_nKm>RwBP%A`7z&9jr*=@b3nj<{vm*e%5$m<3M z(Q0w>p||Nhj!8t{1n5_nd$}WLmKLJq+cL8F`^;>Ts+?FyOvl?Y6IYZ|VO6c|+s0}2 ziA}%7U!vI$nMJoA>8oeUGTSPE_EAne*!b*XXZk@+M19|VMUA^|;f03v+R=A)ROUia zByHW@DyZ+7XGa5!mH~6uu%6v?rIqHz8arvJwn%w!AM0ys@5+5F=AkqM1?$eg8`img zvoI~{7kpt;a6_*(l@MaTwVbZ8&tCw^3wvLhWZLosJjl@U7d@XvtH5hw)hsmYdP|)F zq;wmfGFW#&Qe_Icuo=j7*l-WPUWX1~x*}j-aYf_p9oS`#tTMq_L0tim5jWYt zwz|p&uQxNeCO|%Py-ylBtG<4;B(d%*A@zJ~+0ant`f*zd*n9gszav-cnhIXW|G61@ zscCTLzNrIfk~X#a)MO}D9AcGfavh3R4iXAcTGLG~`unrf`oa$! zDWj#MOP)%6y)3xSc*{9*gjaD@T)v!TBv*h_($YH!?{|wjJyt>fJErgd%FP?IxOS%t}Z-AET}}Gxx)cGeAad0 zDp&_(!zekASJHee!C{O{#yk*FPZdvq=R`_R{{!~lkdTo17#*Xs;BF+o+T;PBmbm+V zF!Zzw>+4~l4g$mzWO#}2cpnx%EhAsqnkc&kPhK`>81!vO=mE_5=zPwHn4=JlBJAp+ zphg7uB_jt1$M{O&^-a{7?W!+X%f?9k+JQSj$(`_T9XDT|N)qV9Tl2L1IUYb?4`R9K zK>!Q}k^~qR3+#*uM7AI!D4aH@6LgH+n19#rY%9HD3;vnJM;VKfx1^ zIy=sM=)oldL*E#Dp09sF$jSvoPymLqJBpu6*3GbSpeA_9J%l;NWvr#z75Rc(Oia!FmM)0*j?$k zSE;H>ZXZVkFEEf@ulJVKL7vEAgV+xF_;wRwYuw%4hcXq3B_KN42!m8c{DhgcDM5BigG}yX_?P|Z8c-+4%mHkMfK-&7-CwOnWq@Iy=?ri| z?;by4lcD$2z}g@BsP9XwJ)S=M9#z;j?ebe~<3F0lZ0JSzA}0_!V9@6lP_xI)9haTe zs?9$>?OChNsyHwp#3g>6Co+C$WCe&MjmiiTf0%ijHPHzWn-DA_oEF$1I6xx{uF{>( zsMoUpW1`p6s7+-7^<}-ow(9O}hOzy1yPu2js zD7P6@wICWqkgE!4mA|pl96*z3)2O(2!K`gl~@uB39hEh+xRccz?Y2P20Z4UOLp2`V?7OtWE7+T$36@I zDrciG23aT22Cmz4e(fH2If0P@haoYDL*X~E0Z302vg<1-Y3hc0Kprf&pA7_#V0k6F z_l}OnK>e}@XmuI9=HzD3(f((gl<#`xCLi*BE$Y1FDsj-!?AzmGwOSTfc>OU@%w7=w zi!u#wryhPy_@5Q)(ndFMkL>NBl+R+X%@V5GRYDsZ?}N+p$Y zo-O#bAE$x4(qgRv{A`i`JOSytV-=e%564J_Qwo*;Rkn;`_H7ZBrj8s$TmdR%wIso- zfWQHrCzBPlUw3V35=! zGlBm8{R_!@I`|BQMJ5MUuA)2_JDsl=m8xc9pr$s+%yLx5_PF}l<;N@({fK&wluFmh zopk+tY+-(`RkI|WaGMT#OXiZ99zOUnYV~GpvR+%_o;7HOc)s*$zH`0WI@-6~=>PNo zv;YE~C195g3+GgkE8sN_IwT1Di81P{0L_`||6s<^OQ9gWU#jSrSV(c1ZJn8+2mL%0 zol-$W^zk5q;J4^ogcjmSGiqvTDqppS?;1Mqkuv97bB!pYA^g3A4+>_nkn+xL;03`Q z_u^b7e!%(m-~FqQBGE*DfoFbqunvq$|C?iTtQPkjPuhY7GO@v>pY^lz>sM%rK>9uj zSHTO>4&>)ery6%yQj@Rq8=veZ?+P;O>8pU|bC)$cMRgd!kRSLrL0cP}Y}Dr7d_Zo` z)@EGNsv_5;k@TO!bu!^aVN)rMkeCp}x_n23nDYM#F$F~dizFY&_>Z8TUk)6>&@ z;?1=W+rw%R*eBVRmYow<4q7}!2H4_=2H0Y^0D5#1p!WS1 z-sYnxPiiJ>2mcrOEbdceQE}56%)~{-BK>OgBCZ~#Q;Z)ao30q9U-O_-0G9wH(ulxj z0IHT8=x|`NnxZ-jpVbIknApX#P5kwZ>%;2@k8?HAZ^~DFho_X^e_vdDNIXcS7QB~) zgpbq4^PjM*nniJs@9*#TlKL1~3W}dh#@wqgwt~_Tj5}7R*>OhYHS_aJSDMc-y55>o zakm9A$dr?X!mhy7B^^mnSLtni2PDwwshAG7A*QePzMved9dn-j@Fpw~q7$lrZ&z*s zZcQV(znWN@?Bw$SF^h24Ol3fPF|E$PA!x8aRzbB&F z#>B|~gzu$!>j)WazAO-h0YEXI={H<3P3J1SQ;iL<{n^sWnCJ62I&y9I`rGmUz6#pj zjXLI>xLp6@W{XRVTdQE>Th*`NhF$Nwvoy}bkw(`Eco1X%sbmSJ<&!X!h#QE$$S3IY zlrjJkpVl@ius|%}h_wMb)i>6j$_uC@a{Dd%k`w+k(JLw{vPkA-B!B$a@}Eiov87x~ z^4xC4$xF?IVOoe`bfN3HKU6mPZP2-Uf3cb58#W6h=}i(M)(o&Q_jGvux@Bzaa5CA8 zbmDl?;!^y8Nf8cIbdbs{0pel|Z7#%Z1u$)Fn5=6n-LMg)j?RwP?Hxq~kCAGoeGav& zBTUX+O9(A7Vx7x!#SoPC$Oe<5+hMojQ$g-@_`h0X0uOkVu$k*9uv{^Z{*(UUbVuD@ z<-+tBc=e5!0!bbSzj41?k-_V~AZ=}J50AG`x0*3kLIU9Jncv(bVmqDHtA6zZ!U#TOW}1cn za07=6*hkEu`9@Mp_lF;F#&)eai9){5H>hWpIIg z^kb5CI5LSqvH0`icHCMGZP^hq5URf6j`f7p^PMq<{|($^hYL}`n^yiL5AGY`n4Hrs zi#Y#z9Cq};(lQq@>&R(pesn$FPIo^8QA7sd2zNT_YimP*U(Y`WOErX`l7banX{YH+ zWu^LP5LEbKAY;k%0$EBu@)}^?^7g#!Ac&<#CE&;+UjacL^w_Uoe}6p-%qS^X zN+`!c-3vvbL8b9%+m7MAy5fetsMf&~m(s$DV! zfl4FGk8$wkSn{)5$Jm3eZe}|-=cSG6qq(p%rH_B@P`#}pCVTa(K4e!ihD!2m?`^Ob zb<)n9XK`)7^KEhIBTVhsOw(Ja}9Xo5J!$viHrW5Sbp>yfOjEfY7T**sZj3H5n0d zlD(wCA9!p4zCf7u9zhnVX(n1UWGnTpe%(;`;BZ5~`NpifY6?uZZJW`x)g~LhmT2_Dxe4PcB598o#>7sXU6N-6{Rm$IhW;gW$f4^`DZJ#V|$YOb%*_sA^q$ z<)p8XVs6r|s#r=GvE@kotv46PTH0+Kf)q`1ayd8vc8PH%z5d;URfH`FieJl%wncMMv5IBi%^78JylP<4;ox*CYP#Jvo(2wCOD=R}ayjj(s zUZ9x{vT5!1{(kPu3)Tmz9DXHYl5>iv?bv~k7O;HrO6uHUk8rpboCIO`ocgHk;`xcE zmQ&=fwo~OoAQS{KRrV$d5Y*3beei;b*xyVmc_iH>O1Ru(%#Q@kyv7gW514(HZL~x% z27;P>zWND?XW$=d?%MyCO~^wIbr0 zm*^N9eua8XcbjQXQa(EVnNA?h^&28@r^iS^U*cI5j8;K83P4B_V2mX6^j?Zb)jv1B zwYQ7V(`r7xbLUQ0R@U^+`Bp%?DHN1%T|6&gn5gcDdi@nGQR6r0?=sfRWHxJYVmTd& zoOE@&yA5Anfw`P)sxW3Sr(A>Ae3&}6gNwO^yd_EfSC+n|kot}2R~0{n9d`&9V3A(T zpB&S8UZk!te8VYv<;LjlpR}{)()$?iMQ`0oj6k_OwTpXdAFa9aicrdHZl@V0CBED%8AZ!d z)U`jO>c%$|@1}exWE091pqg|w86y#6<8kp4QNy|5{+}iB!-WW*myx0eFsHj-wNJym zB_F*o0I+3jh4ow&#pnZ+q3xaaMK{QQ5)wB$x2T7Io=k?nob zT+yI|?89#S8Pe#tHI?o}DQ4|u%<#Wj(Y}eBL>I1SP1W|%rNJ-w_T!18M+ zgII*^BnFJd)_}32^P7Do9L@6owPeM9duivilVv@$GV6z6Cj|)IgDJWh7~0%$R?n^(BZA*ag)`{TWG*h}YsW%DJ@iMgdi}(liI_f(xV$M8xNn)=V2}+; zq-R?(>|0u1&>IcT;77aar{F83f;zpXq?m;My)1B_G+K%Z)8@($2)VPwjRwguPnGzg#5IS!*jfrS}}1eL436GCjks(4#c_8pn=B#$_;PrB4z)uN?<Te{<62h98U~V?3v0-~bp2d_-s{7CpLCurtW6&y_~(W)$VuXI}j&IZAD@PDz{2 z(66-z6(WV>a4d{Q85W3`&_R+)Mn(pO$s#}{STMcycNq6BOmNLJ%&2Jd>%(~w_)4fd zG?wz&h!}<;^%XxZd_1RWKCGqbQq>lM7JTO((ia0fbCocQ7Dsv_D5KyQC35>^{ZMAHJDz%^IINqt`ke1NkRT3?862PuAjguI#6mN zTiN6I#ug=p`&YUiJ@++gPZkfxT9bbfNsIi9*@W>A)56!nuQ6Zyc6a_Tn;qL<`Ip^` z7&0=(v*H63|GM5R6FU8cG~|3_4tWi2*bf;~1tx~`3(}t7vGdm3OK-Iq2GY%O&eSis z7<`LgRO%{stNi}@gQ1+wTnkQGdTSS_uFq^jOzA|WlQt7Jf}n*~g7i6HFaURgEU=fA zs*cO64%Zc01)TJVoze{$RwD;tEhjQ`$5vbTp0_wkMdnk4ly6oyd;8!m#lMXhVQ>G% zK*{NkLiv~V#LWj%wTw2!KrU9St_@}QXZodK+qCj~Z|?YgmYTJ_;7YM}>F8h0sWl@g z6lCw@kKi|1$D)+uZqeWRamqQnvQ)=yZPl(NLq7+)9S9Ya1{S+^CDmhB=a0eN5s^9r z%2%P~1h#gygzYCZqIgb1I#-^C?~3TZ-DT1fHNe-|M<^;@M?0wq1*xTe{-7B2%Py z#q5pP{wlhbR`#<9xMv08);H{5+IyCr6Hqefs(ri@Y!X3~soKZ|XO;lfDjj?Hv$K;c zd&K>n<8wXR=3m$JnhTbRRK^wdIZPmW0Pu1Ihm_^BH^5`gQRpA zR-~u`vFeiW6i<$bw}_M0cL67%;qY}a1HC6jO2O;u7_k}2hDZE!Puw5O zWa1FgEQ&)o=k|P1w$59)7xi#22J;3zX%9FFGN~TU){A4(5@&f>x}3fDwu;kv9!}9| z=i=QnE$fKKgU!R{NdlGHg5B=80Sa@{Th|m~@?Aht_h#l8eywp@6C@^cEVL9rS7Y;g zh11>oTaQr69mUn{Nrw~GQ+7p&m7Xh5d&2KViv)^U#iDnI`SQz;q_2xzeP$j{9gUp%sca3wZKsKYZ%)W z_$T;R)%T9n#_nc(Kpk1R=ZORGcg{1%VCVnd*LQL1Qht>UJy%4fv*#;YkXlaW569Jg zDSFiP?oczrD>z}QD97c%^4*C%u;vUH3Fy4WP>Hdqo;ZeN$H^=B7~#kYNjsN0R% zsKs1}QhjbXIH+v81(IZN*kje@)#va*O-JE%7YYxTi_qAyMcut%^{}O&`;t7Tb2e~= z`&I*CNX;uc;71!!2N@t`(ELJvh9SasGX=#E-_m!iH7=gAPIuUB@R z8q*WES8>wZ>-=3{{qhP|h_khQ%~f;<7cn58U@R8AE8O zfhyl~UV!+uTYrqnCtnKTIkTWx256WG+ffagb-sd_Dd@TXUvuvrmg67454V>FO${ZH z&@O53L>VPY(ojlTG*#MzG-*=WrF>{>8-+w^ik7ytX-GxtcV2v-@ALfe{QDfo({X%` z4|liwzOMK69#lF`Ul(16ct)uNHAIut7oR=pXuxBv;07>_AnF``5 zI}NrDXw*vW^dZ5fQr>)kJ(Wm3qqY-M_^VeA#FI`MdiIQO`#%=XoD z11=#|_DO*p`tto1YSqe~n*%Bi-G$-~cb9smKZ@IR4X9MW1|1Rf^8u_=MzM^bmwC#4%O(cn(*=o!ih!uVm5^ z7hC(Y=<3a851~iy>Q%vULYsL-Srj$3>+*JetPhU6@Rx@eM;`t_=hiOny++?OwxBSV{hpeX6MJnJ=Dx5fZ# z+e7V~0d)S?A7w3(Jtuum(MlQv2kfeJ$2zC@!$*JX)eH>=`fsqXi&A8lSO|+S^i_oJ z`uP^F>dxfsDWR)RJe7wwy)#!Ayb#iqmcuSCl5%o#a5a&Xl9EFWirPoc`U*F0v#sn4d zkdRRk0~-FOUC7KzBCU?~$Okq^M9|7mONsmJ{F}7(kl|Mg@5g)=COFnzRyhzu6FXF; zY{jtQ?yf)9Kx%yer(CqLA=L7qMS?Aj#(l7YDTw z)x3`+U#jp($V|@+Uk<;Fq+k_?z7|9ID9rcq{|1dBytuAQW zcD?oT8CLCi@~|p`RXwb<(WcCK^d+tG`X7zN*TVcgJn!2DHKus`Oyj5{`e$WgOJnaM z#U#$Vwn^>HH-8TL7q8CU)>}zJkumn}-6Mc`H}1Z>xJybZ?*02?Zt^?cyu9xJ&5(_* zH8VoH-$ayw)vem?rfQxxty^2h`-#&MN>0A<0$Xxks|Ts-QmR5HP}=-Pmk^!+#LEX9Ak1KKq3IAu|Sp|h2LQnEN#h@fj9Re zF0Ec&x+xrH7En~5lo>;xi>#}XiLfFVBH%#$ch38yt*Y@rxCX9gS;$mj%OR}}my8dY zh=mX7Xm?_oD_@(Tv|wy{v9{IXl!W4gt>-MM2UllSGKH>Ob;JL;Jr}gVFP#8$VtB4n zxAX|IU2f!_?|R=oE&nEB$ZLA*MMwIbn5ZuaHAS|P?Y(mzr(8~XUKzVC9Bv#^IhZoi za>L$Z%Ty|)@$KGe$=!4k>0@*&7m~dx$~#V(oJej;)0C(yb(QSw%Yn_`^rU`TUK0I` zv~F3mE)=!5godQpNp=*yFb<{?cYRR$r;}gC(73$p?KIZ?=koabtW<@(Z~x1FG%vNf zaN!Y?SM|4M1S=m!=c04I<}7}ysXWVCMpM@M9uZU?4)5bsT?%^+xZ52iiw@L|ruGow z{rrG>I>XJ;jYs%7m&>0N)b$LICo*$3u7T_;l4r{`C!u`H7Bx247=$d`Y~ua8)v#Do`m(Li({lAx1X#bk4cua>k(CBwP$)EjgvQAGy{Xhn)VNo4DW0w45e9r z^0xDhUmx}lIBKzYqak8s@cI@fo02+RvHFhe>G*viTfcolI%8(OCeO`$Xc&=koWJkO zEeFk!Hy6DliH}YFo8^1C?Ga49B@q*oiFzLrrt%n;TqJE?M& zyafC1$MEJeNNXEAMo{~aOA7V*m}y^TXW1gGiLlcda_C+x*8HBHub2#=?l6%;&q>C8 zCEMTu_H7f_$L)^CK8MBB-#mHF|0R*r+0woKpmBXhabNE~dSiu#BN;5ChM!1}TfBCzsxbQaV+t}HcC4`f`!-&qx+1oxo$_9=GsCLK<|uu`Czb&!`@Bh9*r<@kGOl9%4A9-rn_QYwuQ!Taf&MjR zj&Aza%JH?E!>>wOt^C(X|NMSrKt$DBw2j_w9QZw9uCem&h4M)C&GH?SHwt$&C1=}< zPe(>N$|WMBSM~b72sKLQeYH{?m2Cl~i(NE(+INsvxL$(QlPsC3RUO}+c1QA4&Q5MF zwZAfWCiKppj9&sbt@89E_V>w@X`y~xI~1tq`Zb&iz1>k){?t1wbK+9!-^$7Hvk@_a zH=PR_FUw$g{#<_`^eDRUWxj#oWZy?yZOU}I}p4deK?bON&$=4iwzZAD3NfD&WX^*~?vmT?fKV@FKP1C0X4QI=xRNHzh$8i}h-BZOX|#7V=HYp}zxO2~Av~<#(D&xsqXx@x1WLqJi#A6N+b* z7A3~Acf5j5R`|Ru?Xjk~OSXM!b-R4mbjG{C`QBd9`@>@@$~GN$v`}E+Te^QR*vvj+5t9FM4i6qrl{Z^%J(n;N#b!3;oNO|}D zMw`}u7XHNZnJfi~~R7CSKFpb~6^Tx9y=sFmIFEtlWqP zc6z?BE<>AY+cvJi^p2BrrQyPD$=W{NibJ#OH_RC~)dTh=-#FU$e{lf<>PLd7dNcNJ z9Nv84Pr_k^XO`xC|4P>L5hMX41uHgst2=*5&jz<@+@vo3bBcAZQWJZFwh)cmxo<8H zagI>PH_{E9*xgQ>89Quo_Bjt3`Ia}OWeF#e_U)VgNbEQ{LDwjM9h}x)Y5$f?uj=(l zakO)kh?!lPTK2TF9(i&2&*Y9z-2QPtzv8=(ygzbrvGCpdK9wze&zSyjEob|WZ6MXF zj=%P>Crwv96<>T@pRqMMbC^dSLA;V4x3zLKTg{Horb)W+1p6>$vEk!6nHeNeyTf}q zh@85MaBKSfPU-nd`y(|m?}L-My{a53?9E*R25W2WZ4~!{5y{?s&hVSrf=;jY=`L3X z{?v*~Tv97-ue*nJzZVbE)9~r-LS@zHlbqp`WOAE$)k=mV8Y!eln7Y*O#fTBrzBwm; zwwDY(%WUg2=BS}*?^ioXgW_TSsJWJL=7-^MLnI14OC`tprK04hC$=K7LjL;0zUyR4 zV5=R!vOwC@^#@XLpk1f;`$02TKRjQLvuAUI^jknkP3ijp7f%n}x)A2ciVTXHQMu30 z-KGM#>a2uJj*JBs*k0T6aH|uncPEF2yc|Z8ef_NkS<4<4cz>GqT2rV8?A`b!%icqO zF{4E_kD8b2T($#$l=i(CLnEHhtK=f*>aGpR#mx1R57g}I7YjvR<`^xK6iRNo@@rsf zsLKCVfo;O>h`A!R8s86V&PB|VHJ@gh5Mt6wQYTEh>gFGah!i57JQ~kdBWodWmDG_+ zo?LH-ZE$I*H0dhe?PxITL%2QX)y_Vlp}!GSQg&%%ZDQoqf+>e6ZYrAx#<|FF?gzUr zk`DyTlcFF=J4JhJ$e>ugyWGrZ(01|gEo%5V8U4H)a8r*gbsMrhk=wa8`q%17PLrn{ zrgghjC%=;g0Ym=lHl{?8;KlJRZwk)`a+n1`u7Q{-A@;>R@gm`I*Nz(|XWb2&+}5Y@ z;egw%satq>$E$A?#*^`6E5TES>-K$P_J1 zuHC0uBjw~xArfw8ij47mu$9(KO6rjKO8fbs_x=a9E`L5*7WnP#8;J+AX)h!H0Q3s; z-j-WBUMf!OZXGJnUan@GsX!&HT^ra+qe!C|%shF5tT`aUj(tza>a3LKPveP_80Awz zX^1OITkOhQgvYj3?R(rBW0dI(9>~Vb4KbP41|64gN1^?m4e+8=_vDz()w`h`TNVzjO3C_&qr7B~ z;Cfc3^}OP%*N!ard@W=E?4ZKaj12ofNB#u8iB@`{Mt0fctle32B1zuZhDL(Ey)ZYv zBgM7lO>6qipL6`$kA|*b@~{GWYweA=#94lsH6mL0B?`F&Ty-fpEm=`$z^O9se(L*1 zf1mfLZb?iXqo?D``E@W5*Pp;UZn3ae&v25t-2*v<+nf*AI8V)39W$}zNt&AjHE0N@ zwDR&ov+i!ezlB(u?%T<6f1#nH@0p)24?gW*H~qWgt4)zih^1cYx{`hOnkO%n82J?$ z-w%5#-~sEZ?VTahPG=qZO{kb_j!IsMk~&u*TI6gXF`2%OkbcP1c;?+Y_UHGjuPxJE zyMd*;=UmJacbn5Qd|gn7d&zG+ykIwAgy;K^fTmiG5+K?9%pG%y>By%aW_Nf?PmVB? zo>5P&?I3p#WM@U}Pdeuj>h@T+Len@?SG9j)ns(%{@-G@Q?{{9aWr*1Mx`@(z$jp&< zCoa5A=6@V+9%_8W=KYw}!`;wenbgTM{MCg{jo$Ka{Gv4g~!~o6{cBJ z8WM;bk2zl^uO7VXZu8~|57e1<&`8%FCfWaA*ZtEs3i#e2mRJk?mLOlWb#_N&-p#ze zaZ)^|o22~NsJN_?GIVb7W|ih z(U+bk4r-;g1{ZMEgHRUH!|;#YIzsV`orpgrL<6lmaUCHE)QtSAp%>|Oc# z!8tj*6NQb{otz4F62eR0JAfdZp_33)@>J3KfQ`>nA~e*z+}zSaK>4#Uu0g-Kg|N?q zEd&b({qup~*n<6#YdOMR@AMYS0yFt`xklj@xAAAD_F6L@1Jdnlr(b{gu#D$>qM2T` z$NYs)fUU(|9g@z`3ooKoyJ_aNs|Bx8-wu+QaNPpYXK4sJgxt8n%+7uoWJ*Bo7JW}H zpU5l(nhxq@6w03CI8VxMg6#lL3$}wg9MF3~lSugESU0m1 z?L25ux)mKQC?*zz@;jO?68s5ypN_IiCqc!CSCN6P3ZVpmqI{0}_@1`E(X2WFA?+{? zHNIry&|LBt$x2;HYDI>Ux6-LVDKBf`9h~|a0%3c^O$0#n+pc!J~VP0Aj zU)$Ta79f7)nz>#s?t6T{RlcQ|0f5%gOP0r)C6Gk?4@4Fl7wcb`)@smL+cM<7MpVJ_ zj{i#{`b;}qC~9u7K5o|%$$7|c{D=^}02D2bVaD$337Y<@y<@y z+k)f+#J*g1InsqOHR;YZlc%~syW>TYM3sqkU|>i#AO}2b*RMvEjXRqQP#z)VtD!Pr z3Fk+}Uz6Er(j+4(>%Gr6Ndgw`sp<0Y(tXmbe!pnx%$0B7&41iyEc`bVwSy1lyUi-- z-w<+VKYqBSK}nr3_|0Rh0S9+I25ffy|J)O$BU)ZAM2QQPp*6M==i){yv9TzJj0#wP}UBj`>` zfE6t7F?$zRINe>8zUl1bB=$DCJwV_LSSyn1Og=0;OaB%eAF{*9t0}??z*r z@^4{}+Em6J^7EfeuE9khq*)khQbgWi|Ev%Z z^TBHb`0qTL&T{X%ttiSGl_}4TYyP1*ULdQJPmWdp`d!>!NKTJ~3X=m5-12i!6v)GB zf^5d!{h|HO9@aZMs0q*|XW1&Xf6UEjw&tXF8yp-Y7_cOdJ~lzDWfhtzoJ!u(5MH@j z<7jJZyU~aqNatopo|om@>Baq(KVLJizx$D{=qcGvrbh%snSUA5=WBe>VOk=Za!}+? z0&Ple;2ZR6&TjOz0_W3!+mmj2Ypi9-IoV{t_@fGf62MeQ-;msp;nZ278o&v6S+t$V zZ7=a&>Wn_R?lXUhIAc4$FnG&M9?dvvLY=FdPo9;PwK|oT9)^8FKmhHOVh%Fg)z#Gm zM%-f+(~+yo!6j{*aUw8UBz%OGLSE#MmNSnq?+G=7cuDq2l5-Uf_lC&4(4h2mbrgLM zIp?qCkPiQNl`+r6%CFN@IrP#hMlR;rn&9K}4++vlX2s*j!uU#(@zOAm@ud7We9$%Q zG2pl*hHb;E$ogN_vno2R@O4L;HexfU3{UvwsNT_l5_r}(SPv*AzMW+fY-ieD#mwxRswotuhD-Yp3-vYX+f}p z+vC=PtIO(>lIc%`FWP-ROSRShl(<-IAy-+ z;Ep?bH(j1}-saaJjNzCR=Sk#K7cW<1-@#4N^C{m)-^*gBuGj(yoJnREqS*?x_tCid zXQ4jWK$1gkm;JtxBTPe1$n=rPk?*HQt}FtK%7aEfkoAR{{?O3S@fY`Q-@pGJ!sWUJ zhxLCdF0=Y9-dgHO6q)}uGu?&TD$JTcaSxh@olu_7vme~CxVQ+{HP?;wE#~%lT6kiX zt&QUDZKuETChLN!HhAILRDW@iAo=YEP7 zkHL-569itl<{uqS_-e8pdg?f0vmf*Dl~d*B+!!}KRb}rBF`bi9T;f-FS?!)QEsV{N zMqFL~&(0bJFkK;@X}JAxg0IMMrP6ic1dv4%&*4C7kPEA)tor6W2n|4x*3ZJobK}o> zq!`DMl%_Ln9K#)&F?e?r>8l*4+Szz8P-XT{;Jx3MMz8$H0}4nL%CxJUHQ&SQnwI*& z(EK8_cg43Np8nYBea4E!jvQ zv^gJo<~X9aa?M0+*;@%NFob;jHi5&RP48qKqh-7?-prrT-ZnYwfQUMiaS#O!LD0G5 zAv?dWq7A`CgQ|m9`9l>h4%*FF1J@zCd_e(z`m@iTaNoqv{p#|}gJb}j$g%ZMg5k|n zn(RK?AL#UJaO!=}A#zaW2BxoBy%o@)k+xT84-Ba>DLq#}Fyf#GkNAn`-rUSNoIy2D zVWDo^@r&WDYHi@R3&Tg_1;uFV1e@G4$R}o|HrAKgxXmsyG==s=&&&c7;FfWynzEvz zr0mE$$;a!rdf`7jzo(R2Dj#&e3Zm3oHM!qHRrH*Wn&}@$BmAWrY$8*jQZ)FZB2=Td zrNXSwK?Y4vwou8($|PES(cj~O?gOHYt~E``#HV41wDF*u!OF1~CO1HxA6_{7A1|~u zcIc+ADO|3tIi5jY{ovaJT+6+C_h1iu3kCs8^RJke&>=YV^Uv7W*il{G+7VADr*Bc; zFOu=9)2behXDe7%W$&YRqsow;KIJczT=Gtb4to93bUc%8Nrk4Qkj0^w>PP zP=F@?M8+ZXA&@Zxs+92uj2QJqAD9(e)E~dn|Jb-}J?0QXhfQbnj?*wPO~4nk=9{vN z38xtuo0g1i`|@!;gdKv%^!|+mHrI7nVXyfU%8gFm{UAqtx%QSGhP2{(v}z!iK3!^6 zK))Py02-V>i$8HklPY@6WAEJ2zi&|Gh)wT%*7Ig2H=v1~Mb0WMR5h=EXrcm=?9X6| zqx^=Vm9aU=i3TQz`_q3p>XyK8S88nV#U;CQAxjy0Ue+|$gaTucAmoV%^Y!3;7 z7e`{uAD-HS={hp=23B>)&z<06uXWkEhZo$4s~Nf!h9q)Y-a6eQK&??V*{cVJTldDG zGE?}$F!n+77hX#?rwk0Pc=rn)PyXUR+e5e=AAi2SsGIjung82smXmc8=?XD(gAA2D zXOOH~w_*-0ao=72^=oiZzWU|yWuByy1P2N)OM!R3{-XA0u?uzM0u9<-pU%`oGmo+t#-|vuWXgNG>n4PqV@7Lz-uD|&lzjXltr-A ziRKd}gd`mZ7HSXFgz=D0E~%(o^84XKP(=^`lr^%SMVy7Q5+EpVjQx2)#}r zPu?2%ZFe7&+Y0GO-RDuSBM?u{%ce>cGxz(@@9pVH*aK}&V9gT7Z|TauIxC(GBOdF2 z5@yDmlEk*x`f!?glAcI6B;by#TA>;PL!kvPCtA?ik&Nz*CCo7?+xCnXzPvDm#cH{W zi$KVD&Ftzn0&XyqeYaXLX1x8`ZQ%9fTud2%`6luJbss)9Ac9yx`P3-;79I473XGzz zam*l}^1Ih&QwA@IMv%?hpovEaJfH(wsMFm|Yjl!sYSLwjF|bgc(nNHCzg>w$F~X(L zOaktL0x;5sc9x9O=gBVKwR5m;1p$aXKBxNw7MKZ@8vjOg636FkEgFH+GV_kqf0 z;MQM9Io~+(Hm|>00zX*t$vZu>E!U5uyf@oW7C9uDhRDw9(>>01PZIK112}%TA9t}) z;|PUU)m=p|NoX2mU=a~bi+!Qy%dC~x^;;cXrL{U154w5Bt4`R4XDy2)$q%ayZ|Z$D zh%^mFGtl{*Lv%GBb4{1hecwPHLJr!i>=C)E;oGb`{Wr)rSNvP8tQx_t&HszZp;k4t zYuVBB&vG|3Dj0B9goK2A{vljzKL@XCDELKTzaN-i%!sw9>_6Ii;_E=v_tM=Sw94Lm zHx6;1K>2=^JHn&)&Vk~8kd@8S_OCmVZz7?d6wsLJKayf{Btl1hht5(g^i9Ma`=kkF z7NXw=eKk5ik20|Gm=uF0nz1$Bo~)(L*lFw307sjJutrXrlhIt ze5OFpIjUG;yUk}laQ`@pP$v7QjZ3g!62pUF3X-0n^vvXb&0($A}!LtbkVbH9Lm#?I2&k7WxF*DiJFw}*CdTQ}%8n^=Yl<#@T? zRqXY0_-+;!rdMH32bzE{==4PRiXoOW)n>X1>rQTu%VLb64$vD*g`Ckg8nK$cFgs>7 zbvlm)n9C`XeZnUm$~~kvQHk6mp({RhEYemv2!ub1?T9GT75Fn5-@bdtf8f9tMEy{s zc!hwrs5S~!nzXdE29D_yq;cu;6! z`tmX3wXb`ehRB)fie9LX913VlHIO&``Zk3fsY^I`9GKQkIY;9b+sa%C@pP4(6muYW{cCU}7Xh1WL~Ji0_RPW-YI3C4c$)33_2?6vd&%vR1lZ9L#^&i|GNHjAJ4 z(YmWbh<2){WK87B*vLp=y19g?_C(fWlFQ<4nbgxHE!buChY#&!nQP4xY$wazr?ab= z^sQ>#5+Nc-0-xX0)yyMwKE6PaLN^u8SCc_V`Bd!R(v@*yT-+{&ipe!91o@sXMUS#f-%@%GK6-vyDs*!|U6-bzK4cK81VZn+GO zV|;5dz7TGG`fW!n^Cj)LKvHs@4EylXEO*w`9UyTpIX{+_$==!ihbkwtdTl#B6_%S0 zggDO@Tac0vhEvFQFF-P2ztP^V;VD|zuc`wztz+cWFMrPctpC}kZ(~vii1Xm+Fp_*_e;(^!QB7{!COq}d!=W3bh*N*`n({5 z0I(aFJV8LX-~V;XYf_VrizH<6sz%(lWUJZdSyrj9@*Z#9d+1OS+TNvnf4PVzp$hHz|%h! z1T}4b&|&cD$sK}TjPT)bnB!1FjjP5=kZQ56NZn62HmJzF0fbLLfWp)m8~d~E1$$9(8g z($8%5Fy>*|vp$o3nU57~BsaO_3#jC^Rd=8e-E5gOex-U=?-4n8n4{xUxgNK)58RY6V-@D%DF=8a*8{g(aDbo>%xtGlMg&=F7V*^#Wt@q^flUcp z-YhO=SfRU0uSnF{-umHq#n!AZsCJyVnD!SecYbe)YDc7(4?W1+yOEPst6nNMmP4-% z-}(kAC~_8`?s-)<)VdW8+*?#0hdpD;z+t-+`>hVejCTvXh_!!r#w>eq53Tzs)aEl3 zf&sX_KB%ZT{9-foO;lU{I^b4`Zg75`!%!;H{$lj2u}9m zPN&B3v5*250Q7zv^`oZK9{fEf*G&~1J-E@IW4-=qfjG)j^BBXY(lt%eY4zm+-EoZ! z2A(jS4B_Tf;*!7S!TWo$d8Vj<1uZ9OsADZmEv!$pZ5E?>2Yo$*PCsR~EeilV$FB%t z_(Iz)0-Oh}Nugyh)8vm6h95OW;P!W`9CxYAZPKRZMFYxS_~gaWJZ3*r#;-G2ZYpZk z;g|6!STfz7JXJb!t5cTm_3a0nyd3^`ZoSVWOc(kV)s-V_82iv~6b7(UPEgG46Mj-Cmb{pY_1yvQPZWNR?t%H3-w&sDZvX5Teb#LE8f^Ial>)zC^KF|YFC$PpWaj@%Z4Q5LS#}kmUfi{l zUmKsPSpJagL{12*Zp@%TjJCwwi+H`gtjHb&fw44_DYXAAp)5xTFMHF8x^(Hb8oR2ly8#IY9v@W{70;!qBr#>tAV8K#lM;ENtU4DUu`nm6@lTpt(rm zO+SCQH>-+GHU>xvBMs@f3sz9|XZ__!me0ZQ1C)jGGar5Q$$G>}+y6FnlXWgP= zGXMqJA7UhZ+uUrvr5^WD+z3k z)PepW1t>8($p}Uy=cN z;=$a#M1O{XN6Ti`;FO@JAo`p=0Y+D~;9JkK`t7$4s!_8fHGw6s2Xl4{|ya|K{UqgL24^?B?dolw!UQrEJh ztw#-O<72JFJX5m{vA$CHDf%YXslI=>1=Z&KGmk#)P>KG%Ui2vp(Uv+dWBrE@aj=mr zZ{BAJ(JCHgpY3J6$81ye@Phv4>DmCxic4+0kXg||qP(!Y{O7kqSpR94O@mXXh5A--LbfYS#U7_33H*J~;_X}gm5lfLn=U9#7I3ZGP5pg@NZ zB+wq2MmNX$k^|nxfkTsdfq9zvKk=7WySuPX{{H30(R>_*fB(V|r?>Zi-oYg}Af{D9 zv?2pc(3RD)%=UmPw zhM?8=?*+Dt_xgcI`}Z2T1T#P;Ps9QPR0hX=N+$VPJam%gbkKHJ!J>6Jhr?PoFBHe=rQ%br9ttP_VLIqF#?uaLCP@2B?q`9x|dv>rr?kifC|Q zIf#D+F{BUBOTxt*x=?-&#Z@bOFxsi_BqW5M$P<96J#^4-!#su=)*9u_Zk){1nzu^( z=g>KCSa|*&su+YKt4s>R$&)7^8y@(F@9Dvzm)2cW7@#a1ZWbFj%+bN3iVuDWU?xPK z!xfk@fmH$u)I<}Zd@3V~2dMDOFD;D$yv>D50IvCwsog!?gsk7XEMa;I;WfcMdnj@H zKsWGCE6DNCb4iC`-O_r9+8{~ufzCX{30%>lh~2^PQ36&7q+=<00rr^rtyyF6Z<-Yu z)t-0Wf>!$X_vPZ?A`Pn}7}3VzsuEICbdMdQM@nRgH0YwW^#$l!^*VeRABWwIrl)6B zNR2v@7iMc@_ciVZ(tQ?bA3l6YXlO#dYHVUc`@{)GMafT&QVqEr3(=Eil*-@PN5tf?ygX{hf= zA_)o#4xak>hyx01V(;F)gZajNvFEEngl`bc5g4U}LWmE+@_*A7z5}Lb;9~;zfs}+P4~FCI1QMWtf8qELb046nB-z*tf>rz%-|V==#YY$ zN;nH)k~BhI2aJAOL~gpg-3 zuJsVXqR&Dd6z)ueQ`Bi_Es>So9``QCGI!L*#nqK~2NW1koDhV%&8Z5G>SD~CC79iB zWo8~_?SW@BF+d|DBWH!#3De&e8-tS0(sA?v0Xe#gMZ@bo$q(tO+b_5-p}+&;2JT}I zKR;o~bSEh(&!ypIJ+`~~`K$8o)35)m41NFMK8|)%XkY1~wkqd5N&y>cL_F*&(?wCrPo_G1M82BZgL0; z=M5Hdrp?OTT&rvSCingOL0Q?JPkTBHj>Lw9htu)MyU{>-4E`X#yV^!`cEzWi{@T%pWBAgDYJm$FQ=;-uJO`|@2GJIR=Fmq~ZAKjKM&X8D^RZtjP z?ly14M>v`?D+%bcj!aL}4eK+>E$@J-b8IYg{zp4JL}rwPOrsxtw8I{wv6bS>moFIw z1ucgYPMy45ey*p#Un9{P<1dr=L&K(w6Q#wsEiDZ=Sb_4@4|e%IE-z=p4`g)MF?DZ% zwbRwrwNpif4*-dfsOVt>1D5n1we?}T&mKH`wgY;+q$E)16aa9~fo}X0c{4s$DCxd} zRAbs?y+FUzc6wR2$*tJ0pg}3_JgRx|(xv_yVQ;F~-SxK8S;fU$gLdwO0L#+UuiG8A zuTjv;_-I#a>U)}1D-tGOBzmV#kwr9gcJJww+Ak`4!>s6LM8pR}`uc$|Nr4J$sK&+L)wq_=wbz(R zg|kabnbUoLa&mET5tsSsV>rywb{#y(2IZf#Wi||zEMD{)JCqVmNe zBa6$+p>X~_y`Rm_&hB<{GTqTwanL28!~648c^5_`>&X*Vo3d){?&{jwSN;9m___h! zD~^s@etyb0ti>(rNC;J&13#7}aJfH!dW2=IcBS2q_G)NpY2h)I(deo8U1dxKWjhX7=4l>3emtg7_s9`d zKvdu&VPL``X|)3a(5mpzFsc6)(Ccs*c6j08;od*=Ucj{OvDq^&NSK}EvA4H>1vI<_ zpxkLABP#^McycXE!YLb6YFl^GkXH{F2BnI(4ju*k^S z3XVjFysRuTbQ=c>H{V48oO>@@j+J)+g6h%H(fK#rS4c>pZih^gk9a&dP-Q*mcOp#W zoZfh?r+qT_)WZzJN7tpbRSOKS`E93E3n(j-f{Xkv0RfP~>O_}3Z_slFA?jddWzBf} zm>w{iMgDb03blaLRMVIV54ZQ1A@fU|9fbWll<1yRRH)~DAoYAQtYKopeeT>jJrk1% zT(SD2kBeKD<)%D`yqTvlPBp1~-GhVCSy{Ur25Q&{4dEYCAJCT}gDz#ow%DmFdXGNB z4ruyw@vkCYzK_GRv+vg${@nd&r=_pI?NgC?42V>OL@+|Y`mk*R4yIgx8o*f$`TY4j zbUC+T&VPRzb5KU+rG>a^;&|&HmydP_ku~+^cQWOFoNje2#Hfga`M0LBss9$w{mfXg9ut^W+%WEMTJcQgq{3?0^iZ9-ACJ}+-_UnyWX{I6}vmL zriKfJeqw3Cp)?AE)~e-{@4pYv!qHJgQd07EN=pA)%}2WiA zLv%edHFfz_I4#Nj`}c{tj>F^2mo7(*9H0}L`VjS+?O@IHSGXZReGvMj+g1FC$**r@ zw4(VvShyJ?6SnsD!EtdcMhr81Awi9kCmdmF>tHY=rm5}a%Yi#-@(T-Z=;Xe{b8tMp z&s{rq_vtfd_~qm{k%qkj1Q!X`k1KB~1F==KwJ{0C(P9jW&lQ*krQi``>%3@u`07>h z3HB9@QemV|V*u=5B0xm4aRVB77oekj{%A#dr(HQFNe~VKv!d@Vk8^W#UC&_IS`y~1 zAA9{?HZ)LfY;2J5K0AqxscfP!81KYD+LiP4&Trq?o8K)fE4x!tk~wJSFyane`4qj! zj#p6#Q!OO#i97U&SJl?`S6Pw;R1W!Zz8|Ih^}=Lrc35n|r3P7Q-7IPQ6u^`n9UEi4RZPHyAji3|y$fTLpX zC+j~Gx3a=u4M)%)FyKdy9h<#R-H8*v=xk5Lg!W*l-vW)qmzINz%FOeecz z7&&btTo+252IL9(?D>@yZENd2*ix}ef+i5s-mYhoq?6-`%77Xaw=su>n;{RzcNLPk zq7z|RVR!?PNsMI1&;$h_oPB(*oDLd$b+}i8yG%h%o&Dc`bQp)=q$CcHZ=Y$Qad#us z4{J+_u#tUlBZd8<7?G!_Ym@?^(ulYfGp|%FL35MW%*+gQ#lvQ1yoma@VFlB&2xe4Q zbAq;rwL-9D2xW9cZ2dK#{erK#O&<)}smiK#ELoZs5$c;mLq5SUaLmZbp*UB-4jH=r zM05nzL_+Lc?ZzaZWJ(J@ha`Y7YePB2dz2?Zhccoe2oj(YUCT(56_8NvAe7&s%1Gq0 zAf?nn`W37sXazPB3c@Nme4R7Psn=7GXz^wVd0c;iKl!X5wpy1K6F)NCDwf9en#aVNrwg(4i1U2|Wme9-kdeYX` z*5b;_iQq_e=|hKJVb^8V!u7idUsth`yr+k~L!}LS9<+LgTZ&`el}6(H2S*9-)rj4l zfBL^xwn0P|Lvel&tklG^J62UHxt!hI#}irUb{=IuCdh&WNSccQEXIp7vH9FPmL|Bk zxjTA#da7?nVOl>qCxYLq*!$P(X=@TcITn#sfB&FK_vU|JqIduQ;3b<&Gh{dfHd=6z zh|TrQ%otJAK)gno)l+`_OrkX$A|ZB46||GEF}&*PdVN%J*RHH@^G={)cSR&8C*Q^s zBh=!tEh6^UN*XH_eq~uQj69PdYZY5HGU8lOOSI-bb4*iH=M`}~T@gs1B;a-ne*P|m z<4DbMP^30u**iN+5F)p>(p)}RUfbo}gih~IV0a|BYuBC&`2xFkO?7PMx63$eS6=(Vg0D77($&zwXF`Og19!sxVSdH&MCK5@cT-j>rN9y0@Cn(d3QRHU#&6(c4b|9Ig*42@B8+{ z2Nno0l8|6+3!g0}!SBnQ&tR@x{5vsxH509g>;{JvN&BZyH)dxYq@|^Msw8*qIwg5S z9?wkKQksiU7rw9P#SCv^M99Q01f~kj(O3slUR)0AL>DrrWzwUuSHU?!9$V)iBOoEM z4dEWT_ zpr|Nk-ll(2;_rt`#=Q$M_gka9eEITr2}=rudBc;FGzEqcs44B*y_*7ASUoWI3o6bX z4(?EpUo){i;Q)Fy+7OXx_=r@I&K;ov4U#*^a$NmBsD#ObSgWVX4*!z*-@mk!Nx2gj z7mP_pobk|jCnUR*ToF7PKGO2>QN$^56*(`lI*>ET{FfXS7JuWBIWvZaIii6q^~6nt zhuB0c*H;&|BYj|*cvnh)H1>8}9681SGr complex: + """Computes the reflection coefficient of a TIR mode using mode + decomposition. + + Args: + pol: polarization of the incident planewave (S or P). + theta: angle of the incident planewave (radians). + L: position of the mode monitor relative to the flat interface. 0 is + interface position. + """ + resolution = 50.0 + + # cell size is arbitrary + sx = 7.0 + sy = 3.0 + dpml = 2.0 + cell_size = mp.Vector3(sx + 2 * dpml, sy, 0) + pml_layers = [mp.PML(dpml, direction=mp.X)] + + fcen = 1.0 # center frequency + df = 0.1 * fcen + + # k (in source medium) with correct length + # plane of incidence is xy + k = mp.Vector3(n1 * fcen, 0, 0).rotate(mp.Vector3(0, 0, 1), theta) + + # planewave amplitude function (for source) + def pw_amp(k, x0): + def _pw_amp(x): + return cmath.exp(1j * 2 * math.pi * k.dot(x + x0)) + + return _pw_amp + + src_pt = mp.Vector3(-0.5 * sx, 0, 0) + + if pol.name == "S": + src_cmpt = mp.Ez + eig_parity = mp.ODD_Z + elif pol.name == "P": + src_cmpt = mp.Hz + eig_parity = mp.EVEN_Z + else: + raise ValueError("pol must be S or P, only.") + + sources = [ + mp.Source( + mp.GaussianSource(fcen, fwidth=df), + component=src_cmpt, + center=src_pt, + size=mp.Vector3(0, cell_size.y, 0), + amp_func=pw_amp(k, src_pt), + ), + ] + + sim = mp.Simulation( + resolution=resolution, + cell_size=cell_size, + default_material=mp.Medium(index=n1), + boundary_layers=pml_layers, + k_point=k, + sources=sources, + ) + + # DFT monitor for incident fields + mode_mon = sim.add_mode_monitor( + fcen, + 0, + 1, + mp.FluxRegion( + center=mp.Vector3(-L, 0, 0), + size=mp.Vector3(0, cell_size.y, 0), + ), + ) + + sim.run( + until_after_sources=mp.stop_when_fields_decayed( + 50, + src_cmpt, + mp.Vector3(-L, 0, 0), + 1e-6, + ), + ) + + res = sim.get_eigenmode_coefficients( + mode_mon, + bands=[1], + eig_parity=eig_parity, + kpoint_func=lambda *not_used: k, + direction=mp.NO_DIRECTION, + ) + + input_mode_coeff = res.alpha[0, 0, 0] + input_flux_data = sim.get_flux_data(mode_mon) + + sim.reset_meep() + + geometry = [ + mp.Block( + material=mp.Medium(index=n1), + center=mp.Vector3(-0.25 * (sx + 2 * dpml), 0, 0), + size=mp.Vector3(0.5 * (sx + 2 * dpml), mp.inf, mp.inf), + ), + mp.Block( + material=mp.Medium(index=n2), + center=mp.Vector3(0.25 * (sx + 2 * dpml), 0, 0), + size=mp.Vector3(0.5 * (sx + 2 * dpml), mp.inf, mp.inf), + ), + ] + + sim = mp.Simulation( + resolution=resolution, + cell_size=cell_size, + boundary_layers=pml_layers, + k_point=k, + sources=sources, + geometry=geometry, + ) + + # DFT monitor for reflected fields + mode_mon = sim.add_mode_monitor( + fcen, + 0, + 1, + mp.FluxRegion( + center=mp.Vector3(-L, 0, 0), + size=mp.Vector3(0, cell_size.y, 0), + ), + ) + + sim.load_minus_flux_data(mode_mon, input_flux_data) + + sim.run( + until_after_sources=mp.stop_when_fields_decayed( + 50, + mp.Ez, + mp.Vector3(-L, 0, 0), + 1e-6, + ), + ) + + res = sim.get_eigenmode_coefficients( + mode_mon, + bands=[1], + eig_parity=eig_parity, + kpoint_func=lambda *not_used: k, + direction=mp.NO_DIRECTION, + ) + + # mode coefficient of reflected planewave + refl_mode_coeff = res.alpha[0, 0, 1] + + # reflection coefficient + refl_coeff = refl_mode_coeff / input_mode_coeff + + # apply phase correction factor + refl_coeff /= cmath.exp(1j * k.x * 2 * math.pi * 2 * L) + + return refl_coeff + + +def refl_coeff_Fresnel(pol: Polarization, theta: float) -> complex: + """Computes the reflection coefficient of a TIR mode using the Fresnel + equations. + + Args: + pol: polarization of the incident planewave (S or P). + theta: angle of the incident planewave (degrees). + """ + if pol.name == "S": + refl_coeff = ( + math.cos(theta) - ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5 + ) / (math.cos(theta) + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5) + else: + refl_coeff = ( + -((n2 / n1) ** 2) * math.cos(theta) + + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5 + ) / ( + (n2 / n1) ** 2 * math.cos(theta) + + ((n2 / n1) ** 2 - math.sin(theta) ** 2) ** 0.5 + ) + + return refl_coeff + + +if __name__ == "__main__": + thetas = [54.3, 48.5] # angle of incident planewave (degrees) + Ls = [0.4, 1.2] # position of mode monitor relative to flat interface + pols = [Polarization.S, Polarization.P] # polarization of incident planewave + + for pol, theta, L in zip(pols, thetas, Ls): + theta_rad = np.radians(theta) + rc_m = refl_coeff_meep(pol, theta_rad, L) + rc_f = refl_coeff_Fresnel(pol, theta_rad) + + rc_m_str = f"{rc_m.real:.5f}{rc_m.imag:+.5f}j" + rc_f_str = f"{rc_f.real:.5f}{rc_f.imag:+.5f}j" + print( + f"refl-coeff:, {pol.name}, {theta}, {rc_m_str} (Meep), " + f"{rc_f_str} (Fresnel)" + ) + + mag_m = abs(rc_m) + mag_f = abs(rc_f) + err_mag = abs(mag_m - mag_f) / mag_f + print( + f"magnitude:, {mag_m:.5f} (Meep), {mag_f:.5f} (Fresnel), " + f"{err_mag:.5f} (error)" + ) + + phase_m = cmath.phase(rc_m) + phase_f = cmath.phase(rc_f) + err_phase = abs(phase_m - phase_f) / abs(phase_f) + print( + f"phase:, {phase_m:.5f} (Meep), {phase_f:.5f} (Fresnel), " + f"{err_phase:.5f} (error)" + )