From e4fe3b6b448c1049460de255bc0242fcf1c32838 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Fri, 28 Apr 2017 13:51:43 +0200 Subject: [PATCH] Enhanced readme, changed input file to YAML, and improved functionality when interpreting the chi2 components --- README.md | 96 ++++++- input.yml | 65 +++++ paper/{ => figures}/3D-Si1x1.png | Bin paper/figures/chi2-001.png | Bin 0 -> 19402 bytes paper/figures/chi2-full.png | Bin 0 -> 19615 bytes paper/figures/chi2-shg.png | Bin 0 -> 2804 bytes paper/figures/epsilon.png | Bin 0 -> 5919 bytes paper/figures/polar.png | Bin 0 -> 8362 bytes paper/figures/refrac.png | Bin 0 -> 12775 bytes paper/paper.bib | 12 + .../{SiH1x1-chi2-yyx => SiH1x1-chi2-yxy} | 0 .../{SiH1x1-chi2-zzy => SiH1x1-chi2-zyz} | 0 sample.in | 43 --- shgyield.py | 266 ++++++++++-------- 14 files changed, 306 insertions(+), 176 deletions(-) create mode 100644 input.yml rename paper/{ => figures}/3D-Si1x1.png (100%) create mode 100644 paper/figures/chi2-001.png create mode 100644 paper/figures/chi2-full.png create mode 100644 paper/figures/chi2-shg.png create mode 100644 paper/figures/epsilon.png create mode 100644 paper/figures/polar.png create mode 100644 paper/figures/refrac.png rename sample-data/{SiH1x1-chi2-yyx => SiH1x1-chi2-yxy} (100%) rename sample-data/{SiH1x1-chi2-zzy => SiH1x1-chi2-zyz} (100%) delete mode 100644 sample.in diff --git a/README.md b/README.md index fec11bd..c007a93 100644 --- a/README.md +++ b/README.md @@ -3,21 +3,31 @@ SHG Yield for Semiconductor Surfaces [![DOI](https://zenodo.org/badge/11697217.svg)](https://zenodo.org/badge/latestdoi/11697217) -SHGYield is a python program for calculating the second-harmonic generation (SHG) yield for semiconductor surfaces. This is useful for conducting theoretical second-harmonic spectroscopy of any crystalline semiconductor or nanostructure. The physical theory behind this phenomenon is presented in the references below. We have applied this software to several silicon test cases and have achieved very accurate results that compare well with experiment. +SHGYield is a python program for calculating the surface second-harmonic generation (SSHG) yield (in reflectance) for semiconductor surfaces. -These calculations can predict effects in novel materials that are still under study, or save valuable time and resources for the experimentalist. For example, the figure below is an overview of the angular dependence of the reflected SHG Yield from the Si(111)(1x1)H surface. Experimentalists will find this very useful, as they can plan the experiment accordingly in order to optimize the output signal strength and polarization. -![An overview of the angular dependence of the SHG Yield for the Si(111)(1x1)H surface](paper/3D-Si1x1.png) +Introduction +------------------------------------ + +SSHG is an effective, nondestructive, and noninvasive probe for studying surface and interface properties, and even for characterizing buried interfaces and nanostructures. The high surface sensitivity of SSHG spectroscopy is due to the fact that within the dipole approximation, the bulk SHG response in centrosymmetric materials is identically zero. The SHG process can occur only at the surface where the inversion symmetry is broken. + +This program has several potential applications and uses: +* determining and analyzing the physical origin of SSHG spectra +* predicting and characterizing the radiated SSHG for interesting new materials +* characterizing thin films based on measured SH spectra +* allowing the experimenter to calculate and analyze the SSHG yield to optimize experiments -You must first calculate the different components of the nonlinear susceptibility tensor. The theory surrounding this problem is still being developed, and there are many ways to go about it. We leave it to the reader to find the best method for their particular problem. As an example, we use [ABINIT](http://www.abinit.org) to calculate the electron density and matrix elements and then [TINIBA](https://github.com/bemese/tiniba) to calculate the tensor components. +For example, the figure below is an overview of the angular dependence of the reflected SHG Yield from the Si(111)(1x1)H surface. Experimentalists will find this very useful, as they can plan the experiment accordingly in order to optimize the output signal strength and polarization. + +![An overview of the angular dependence of the SHG Yield for the Si(111)(1x1)H surface](paper/figures/3D-Si1x1.png) References ------------------------------------ -The complete theory is derived step-by-step in my [Ph.D. thesis](https://github.com/roguephysicist/thesis-phd). This software has been used in the following publications: +The complete theory is derived step-by-step in [Phys. Rev. B 94, 115314 (2016)](https://doi.org/10.1103/PhysRevB.94.115314). This software has been developed and used in the following publications: -* [Front. Mater. 4:12 (2017)](http://journal.frontiersin.org/article/10.3389/fmats.2017.00012/abstract) +* [Front. Mater. 4:12 (2017)](https://doi.org/10.3389/fmats.2017.00012) * [Phys. Rev. B 94, 115314 (2016)](https://doi.org/10.1103/PhysRevB.94.115314) * [Phys. Rev. B 93, 235304 (2016)](https://doi.org/10.1103/PhysRevB.93.235304) * [Phys. Rev. B 91, 075302 (2015)](https://doi.org/10.1103/PhysRevB.91.075302) @@ -28,15 +38,81 @@ The complete theory is derived step-by-step in my [Ph.D. thesis](https://github. Installation ------------------------------------ -SHGYield has been tested with Python 2 and 3, and Anaconda Python 4+ on Mac OS X and Linux. It should work on any system with the required Python packages installed. +SHGYield has been tested with Python 2 and 3, and Anaconda Python 4+ on both macOS and Linux. It should work on any system (including Windows) with the required Python packages installed. Python requirements: -`sys`, `math`, `numpy`, `scipy` +`sys`, `yaml`, `numpy`, `scipy` Usage: -`python shgyield.py ` +`python shgyield.py ` + +A sample input file `input.yml` and data set `sample-data/` are included for your convenience. Both the input file and the program itself are extensively documented. + + +Theory +------------------------------------ + +In summary, SHGYield.py produces the SHG radiation (in reflectance) produced from a semiconductor surface. As you mentioned, it requires the susceptibility tensors that are calculated with the above considerations. The theory is developed considering a reflectance model with three distinct regions that allows the user to readily simulate the SHG response of thin-films over bulk substrates, or of any crystalline surface with any symmetry considerations (see my previous reply to this thread). + +In order to calculate the SHG yield, you must first calculate the linear and nonlinear susceptibility tensor (χ(−2ω; ω, ω), χ for short) for the material of interest. The theory surrounding this problem is still being developed, and there are many ways to go about it. We leave it to the reader to find the best method for their particular problem. As an example, we use [ABINIT](http://www.abinit.org) to calculate the electron density/wavefunction/energies and then [TINIBA](https://github.com/bemese/tiniba) to calculate the relevant matrix elements and the χ tensor components. The exact method that we use is derived in full detail in [Phys. Rev. B 91, 075302 (2015)](https://doi.org/10.1103/PhysRevB.91.075302). The program does NOT care how you have produced the susceptibility tensors; you can use different frameworks (such as MBPT, DFT-LDA, TDDFT, etc.) for producing the linear and nonlinear susceptibility tensors. + + +### The nonlinear susceptibility tensor, χ(−2ω; ω, ω) + + +χ determines the nonlinear polarizability of a material and is responsible for second-harmonic generation. This relationship is expressed as + + + +where a, b, and c are crystallographic directions that depend on how you orient your crystalline structure. We can see that a material can produce a polarization response in direction a from two incident fields (E) in directions b and c, by means of χabc. + +χ is a second-order tensor, and thus has 27 possible components (unique combinations of a, b, and c; for instance, aaa, aab, and so on.). Second-harmonic generation implies that the incoming fields are identical (two photons of equal energy in, one photon of double-energy out) so it is also implied that + + + +for this particular phenomenon. This reduces 9 of the possible combinations, reducing to 18 unique components. It is very convenient to express the crystallographic directions in terms of *x*, *y*, and *z*; therefore, we can express χ with all 18 components as + + + +Symmetry relations are very important for determining χ. A given crystal symmetry can greatly reduce the complexity of the problem by eliminating many of the components. For instance, for the (001) face of cubic crystals, we have that + + + +which has only 3 independent components. There are many [articles](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.35.1129) and [books](https://books.google.com/books?id=bdFju3af2FsC) with tables and extensive discussion that describe the form that this tensor should have for a given symmetry. + + +### The linear susceptibility tensor + +The case for χab(ω) is considerable simpler. χab(ω) is directly related to the dielectric function of the material + + + +which is directly related to the index of refraction as + + + +The χab(ω) spectra should obviously have non-zero regions; otherwise, the problem is not very interesting. + + +### Calculating the electronic properties of semiconductors + +In general, there are several *ab initio* formalisms that can be used to obtain the electron density/wavefunction/energies of crystalline semiconductor materials. The most common are: + +* all-electron methods (older, declining in popularity) +* density functional theory with the local-density approximation (DFT-LDA) (solid workhorse with known problems) +* time-dependent density functional theory (TDDFT) (good for optical responses, but hard to implement) +* many-body perturbation theory (MBPT) (state-of-the-art, but extreme computational expense) + +Each of these has their pros and cons that mainly relate to accuracy vs. difficulty of the theory vs. computer resources. There are [many free and open-source codes](http://psi-k.net/software/) available for download that are under active development by thousands of researchers and groups. + + +### Calculating the susceptibility tensors and optical properties + +Once the initial electronic properties of the material are determined, we can then proceed to calculate the optical properties that include the linear and nonlinear responses. In general, nonlinear optics is now pretty well understood within the DFT-LDA (see references within [Phys. Rev. B 91, 075302 (2015)](https://doi.org/10.1103/PhysRevB.91.075302)) and the [TDDFT](http://dx.doi.org/10.1063/1.2790014) frameworks. MBPT is currently (AFAIK) still at the linear optics level, but has the highest level of accuracy available. Of course, you can combine different methods to exploit the strengths that each has to offer. + +The code that we use is called [TINIBA](https://github.com/roguephysicist/tiniba); however, we have not yet created an official release for production use. That said, we are an active research group that uses this software every day to produce high-quality scientific work, and are constantly improving and adding features to it. Publishing TINIBA is definitely a long-term goal that we have. -A sample input file `sample.in` and data set `sample-data/` are included for your convenience. +There are, however, some codes available, such as [RT-SIESTA](http://monalisa.phys.washington.edu/feffproject-rtsiesta.html) which works with TDDFT, and even [ABINIT](http://www.abinit.org/doc/helpfiles/for-v8.2/users/optic_help.html) has built-in utilities for calculating the susceptibility tensors. License diff --git a/input.yml b/input.yml new file mode 100644 index 0000000..4335b6e --- /dev/null +++ b/input.yml @@ -0,0 +1,65 @@ +# Sample input file for calculating SHG yield +# This is in YAML, so please preserve indentation below! + +# Establish reflectance model, see PRB 93, 235304 (2016). Available choices: +# 3-layer, 2-layer-fresnel, 2-layer-bulk, 2-layer-vacuum, and 3-layer-hybrid. +# Anything other than '3-layer' will exclude multiple reflections. When in +# doubt, use '3-layer'. +mode: 3-layer + +# Required calculation parameters +parameters: + theta: 65 # angle of incidence in degrees + phi: 30 # azimuthal angle in degrees + sigma: 0.10 # gaussian output broadening in eV + +# Multiple reflection framework +multiref: + enable: yes # if multiple reflections are taken into account. Can + # be 'yes' or 'no'; if 'no', next values are ignored. + thickness: 10 # thickness "d" of the thin layer, in nm + depth: average # depth "d2" of the nonlinear polarization, in nm. can + # be a float or 'average'. See Eqs (21) and 22 of + # PRB 93, 235304 (2016). + +# chi1: linear susceptibility tensor. +# Paths to files for chib (bulk) and chil (layered), and normalization of chil. +# For information about normalizing, see PRB 92, 245308 (2015). File columns are +# Energy(1w) Re[chi_xx] Im[chi_xx] Re[chi_yy] Im[chi_yy] Re[chi_zz] Im[chi_zz]. +chi1: + chib: sample-data/SiBulk-chi1-xx_yy_zz + chil: sample-data/SiH1x1-chi1-xx_yy_zz + norm: 1.2659296143 + +# chi2: nonlinear susceptibility tensor +# Components should be in separate files, and listed separately. Symmetry +# relations can be represented here. There are 4 options for each component: +# +# 1. ZERO if it is not listed or commented out (i.e. '# xxx: ...') +# 2. ZERO if you give it a value of 0, i.e (i.e. 'xxx: 0') +# 3. Loaded from file if path is given (i.e. 'xxx: ../SiH1x1-chi2-xxx') +# 4. Has some relation to another component (i.e. 'xyy: -xxx') +# +# File columns should be organized as Energy(1w) Re[1w] Im[1w] Re[2w] Im[2w]. +chi2: + xxx: sample-data/SiH1x1-chi2-xxx + xyy: sample-data/SiH1x1-chi2-xyy + xzz: sample-data/SiH1x1-chi2-xzz + xyz: sample-data/SiH1x1-chi2-xyz + xxz: sample-data/SiH1x1-chi2-xxz + xxy: sample-data/SiH1x1-chi2-xxy + yxx: sample-data/SiH1x1-chi2-yxx + yyy: sample-data/SiH1x1-chi2-yyy + yzz: sample-data/SiH1x1-chi2-yzz + yyz: sample-data/SiH1x1-chi2-yyz + yxz: sample-data/SiH1x1-chi2-yxz + yxy: sample-data/SiH1x1-chi2-yxy + zxx: sample-data/SiH1x1-chi2-zxx + zyy: sample-data/SiH1x1-chi2-zyy + zzz: sample-data/SiH1x1-chi2-zzz + zyz: sample-data/SiH1x1-chi2-zyz + zxz: sample-data/SiH1x1-chi2-zxz + zxy: sample-data/SiH1x1-chi2-zxy + +# Path to desired output file +output: sample.out diff --git a/paper/3D-Si1x1.png b/paper/figures/3D-Si1x1.png similarity index 100% rename from paper/3D-Si1x1.png rename to paper/figures/3D-Si1x1.png diff --git a/paper/figures/chi2-001.png b/paper/figures/chi2-001.png new file mode 100644 index 0000000000000000000000000000000000000000..59568842076f3f9f4394bab4bcc8ecfdc8cc191b GIT binary patch literal 19402 zcmZs@bzIY5^gllPEe1*|3Mh;a1e8>elG^AHkXC_FQqnCd(%?opkdo0bq;ud~2~l#y zXe6XYNR1fry~FqW{rP-;k01ZAz2e?;?mczS^PGihsw-WfVxodTAQ#{;d2I;f3={%6 zX7blbaAnS#h7bHZX(g*B3xO0zpWAzS3Vc8QOhri^au^)l)&RbpafPcXoS8Uz<@D7X zeaF>xArN>9TweB(=g`uqhWq5V-NPM?-rSz-GYx@)j?~_U-7BykIoT2x1#XK1 z2#iyc+s2NaV0Ro4(EIlrD!M^?TDQZw-LJB+AQVPS#Rsqt@n)NMkon?tH#b+wJn8`cB5ttrJ#=#FiLpQs3qfxCe1l-W^le&Fy%~0f*kT z&Yd|>ba;3gTsKxuUg`I~#%N-6_X#KL^kC>l60!UkxUmqwFeWZqy1T%EY<$N@->B2D zx(9>I2?>eSaw_tje7RhzF(eIt;n3{Xr%cd0ueY!at0_M&)9Lr&_v0yVNyz7*M_c_v z%dr>*niG~^pph6OUdYCWvy-XmM3G|grd(j?cE`%!=RJ3pU`M#Q@(Uoi8o5w)#$E0p5q(Ai2#}w+3 zmPm%(Ci1@AA+DvWIUaC4a{@K7rLxitrL`osi%q464K0j0aI~tH#nR17Ywy)^kw?@` zP!`R(N|m|W7lpdDk^9r$*ZPad9!3oXmPH>=qh7~onGCS&|y`JT`Ed`}I{!HM#BMw+#oBR%gQEUTuFi@~(niyepypZdW^GTayau;A|8F>&Gbo z>OWbhqeS!Cs(aE;IVZLmiJ9ov^&VW?!p7+tW)Jnb8#XmaB13ld6Ni7G8b@tCdYEd) zvn4FLbb{Y4xt;8Xsd}h_u`^MSZ%C zuX|*Mt>Vp#>njs6z>`^5oPnrRj`1;=Q7pPPT(5JH%*$@Otp2=UiA%w?e4IE?Po7rh`o)2unaeY`i63V;3Y}(qGtGSsAPlV}4mJ z^Lp+*o-SoQV#&3jg#um!gnOL;dS+p?vS1!}yGLbZA>Md4DfFRDKIZLUasBn@+Ldyt z%uhwpzq5MygF>h)FU|K?-V3$x8%8-g3WUeEuz6RBbt{~~_RrfhyK8#rWF=)A^q0MB)x^=GT(rWq*{BOqbvHGFwp>x~^ zUacZpZVv1R1A{2GyBrS%eK5^7dyiGMGcsWHn0v3+CU~PMur)f7QoGtVwY=x<4%(zng`uV(yw0bU`p0E{}g#Q+0V$egWa0=Tz zv>+%{2$dX^GEK2|5l~6JsiEdMQ45E52Ou!;g-os`q%!OIJQkSt4>Y`KD@_7uz$KgQ zaL1L;m>Okltxkf0HUIpW)VQ98Qh8N^<#J5jM>Ui@&3?id-U9tXfVGE0<*#I_3L-@MGSe* z%vnMCyzS3xY*KXcOb3(05)Td|0u)7zLyFS*X>ALOD|U0OqRIEHy#z|HWmZlH&9*{O zbF~}%#jA4*y1O#JJ<2OBmJiMsUG+P6@xQH@7&MDG)EE?K|5j(#6}L81(+0B(^6nWe zUgdAY-!^f;GL4u=$ASgZ67LvpyiAw!lc^_$4IL4_A{C zkG|n@W6yLaQp+}xU!LXWq?{~fLZBY!fd1UpUR&Ipm~4=wFk@OHzL1JeF0O9DvmP{? zajg6$692uUkfLADFy$CVvXk#Q8~#=z16;|+KwLXUrP+Z4DIbl#A!M3gw1n&# ze^{``U!~L-{0jEf%-4rmkakH}CkLjThRiZF=`9pC8^;g0GQxxi#(@|EPFYE)y=7sD z21>bEt*{h}>Cha9;;sjbD2&XLQ9x9TU{y|<)DG19$wW08 z_V$BveG3J|sjWagDrgYGG9H0=QHbCab}6J+-fd6H4oW!ZfSsrj$IpIaEQ22m{GNe) zxGOtTTu@}nY{>5qeJ1eCP{eD=Dw&9s4 zPSkP7+`Imk!x%%jnRI~G{!5+hz!GkwSHcLa*9ntc1~VtJdom9gH-(zRAZu zM=Bp?X2*2H(iLd-<8UUP&X#-*d+h5xd3*I6s8o}Hzl>l0tC=-|`kV?SNvd33O_o#7 z6GLZE%Roa-OXFPQgKm1W>zMwDh0x*jVvEOJQwp2pAwVZ zM`zSM`%R<$YYmoB6oq!Zm!ikK69;rE&LMild*df{ZVInY2RVLGg&!X!j`HQVU0+K8yg`1p$TaUp_XV9N7($w;?;jwX1{ z#*>w$q{!zMshgP_{6E@PgI{?zj+(~XJE8S_XhWgi-!0C(B;`xeKp;QA;_B(NF&j8! zMxH$ifr-OpKeKvFPSc00tUKr-EDN%R7CIa_$TZY3rGhI6%l;voEZA0yPJ=ixOzD{) z%)_AE$S2#EvEOPqSQwu^;_-1(P6h_;)*m$aoE_+nNyHAW{j`cc>GtShqvWMtt+U9t zw6e@`zn5?yMj3*$TO&&i=A$YsrhC@n8_)P*tjzZd#3Tm<$PKn;n{bVxxzc8;jjS`| zc~|KoCL)uHH_zCPHS0!BBnEa(a=vFZ(E9;xuU2cEkHN0-{TxxYJFl105>(%Ng>_uh z&jry~^j?Z>+=mnmjwz(Zc$BovR|nVWtC6Dk_1s_=Vf_TZxch5Kj4MoMYaEdCUAqFh zSE@o7;2S!m@4dRO&3$f>=L}Da%t`-iqJ|?>nJhcTdq`^OqNK4++lb&jxglag(zy?!TUt?cO=s4xY*XuyRL@Ep+T|Mw^U;l>2jV z5V^-SO27JzI|-48#`x6cYTrD){H3E1nAp@sLOnIJ(Fp?Y%lK1ReF)!K<06Pt`R7e4 zneQVUSL59rM9lf1k}tGIRb+_)cYXSOOTL!L1Cz1Kv2IJzA_e_=0a(LN)4!tSk3QX) zWC_9gn!WP}hEU%17x~r3xq6&i$*`(^8@)!8w;UGz%xVx0Q zQ*9k|UPE!OePD0j{;Q4dTDvcHiU!Tl6azMtt-|w0XL+mS`H41kEfn||vHiJ$T~j$` zb;w-lf!7{mYC_4GET=A5I8KmEQ@4?s5(nf9;Ga+}z8zbQMjN z7wMI4^t;rjc?S9Rn$`uk-%Ew)$k=-FUB8PPf3}`!k_}uRDNCRiuS}`CB+C|RdUsG5 zPlzhqc`dVl0y^%b=Ek1?0F2aMrbO&5`Kd0@W+EX@+}A~%up z0`G716$osmNou1O^*UR*J|vGGC!0Jd>3o&&pQI4p=OR$cvGp&d>MKNDWw}TtduUgi zVGy|P-W(6!E6TT>zfO+bGv6kliQ&Qd`qZx~|9dy#eK7-yWl?edn8)16!u>6d4hPJ+hgo9da)S=1~egnu(?)4h2(+C; z!^+g?6+}VyBgU~pDK3r3unh_3;N4(~ZF6{ic9S1oL+mW1Hu+phpfANgHX?I~iTN3h z|8Pma`Jyrb?Y4~m%WRFBz>buses3}(TZhlCmn!`JhKXpJ{_s>-S)Ye}*`%35z>xr| zu_K>f93La0GdaUHbPk?#xjfzP6V!8XtQ!2$#&-4l?ANK5{i<9J7<8g`dgF;eNvd%X zEhJ3^IB!-27OiX;3>E5L`~97XU`N$}&SgMf4GO1?7@|7-rNl{&CenkXhV?77E$4jw zdqtkD&!-g~nj+V%dV?J4LZN#J0dJ=l@(e+KL%dzJ?A?auSZ=qXw;ImZ6<6`wDE(0J z$*aZlH1KW5SRb`Y4hW=60)ffIaFM-Tk*@lU@PjheD$V*)35pvy+liKW>cC)PbPOWX z3AWpHrhA>FrnT_v(PzH3NTF)emr&I2(avgumYEX8vB_B;r$p*0#7~JV`Y|1aRqbZg zgz8cUFdP4PaU2N-W)mByv(}kFvdniB-aN7VlXJ(%*~My z-a*Lfp$oz8Vk$jrf`y=1EE4v`V~3%&0h%q(W)fk$hHmHoBBSo_be6*tHSPRd_R&8h_xQ}wDHJA* z{Bx2~ik{$MX*5M01;iV6_(UWH2aAnnLELltdW>Tl2v;RhD+5f`rZ+)?&ij2RMfDVN zQ~lZ_Mr-bk5D;03-1rr&G-F+H2d=|4juZK2d`)iY4uThCBu~-W-kNW%H{DeyzqfHS zMr(qTTBD84^lr@G*Mf|i@n4tT!M-@|Fu<4JfehyUC75?bC!wxyjCgB4n|@(*v8?A2 za#_G`wDB{2lC}2i*grIH>XrkVT_Z*xo2O#^w@UXhN;unVT(KZRy>~ zBIKVj?u?)+1?fTvZdv0^!JS z1OFUmXYxHrn)TuM-!G~!L-<$qmqqwcbD_TBP_DlqFJE-hB=13fKIEd<$Tyu=O6eh* zp?)RJQZ=UMAiq^t7Zc@aa5CYP%w$*W4NwDdE+@qn59 zv^&9n=5#n;@!IrO)F6<~BbxsVtl4;Z;en#$$C z9j6RN{(R~EOThdj%Ml|Nw4QcY0djb_+a~_GA{$H#r%m)Ib;2=R@T-V3Ncy~K=()VU z1k&dqu+1NM(lg2c7wR`mkH-=+ZmL9o&cLv3hk3n+Y(J2*PL0A=diSd?Jl|`#GqcOeyi; zZ>hvo5&U=G;^CRf=ZDNbk9fd9FRxKpX)%Ph!+7CFWKZ|j;)DELOeD)9myZ?eVzP*L zmtG)ZW`9pcw36^Hj~CDw-t+C}E9yO6RO4UR1oJ-<${fHr^M+XkI;CaE+QEltHoWz2 zG=oucKqJv? z|8MF3iE&p)9PIzju3irFhuUk*zE!*spk?|ua)a^p6fDXwyT)bzClYe8|de~Aa_ z?x+aqg?olR{O=37bCGR z7Deh`>cxrg6 zkJvzBlN}F-qvyO4GAYwr)E9a+8=^^_x3=JVr0D+?0Y?|tz9`Rj)kAa?y7aPTY4}Vh zU-jX>0C-xW)-g0$^sv&K%?U@bBA?&6u_O1c;Aa-ug;jI*{od&Ie7@l`JFz4pJZ=6t zrG?=Dc(4-|GYBh&O2P@qJ1DK|4Wz1~;{Md)#|(>w6(N-aleH9YTkGZ)7e2YuBe`=& zx^vnx>t%{GL^Jgw6t%lK6)fpN__)z&!INi9mbSyDUnTbsw5pY#zZ!iAEJ-ab`(}Mi9vPj(RSYpWX_PG%hI#4QPzr>)B!U zB9(cp8r0FJ^OqFOyl71~N%p5V#p=Gi>5>BD%*Z*=eL8ppoi~QOn^XNZtMN7Yigqhe z9-KtM-yHk3xOTTx*jS5nr}XX=>ZV>jt5+xB+cFgv?Reyo^*HYcl?qCjx7QjZIxl;; z0Qh9FZK$nt{I#K#$!cvD!oHvwq3eEk>~_Sjj>l=2|Bl8Io{T$}ebJQIuh75@Ul}^X zBE-I5AV5XGP|IdtUZT4!@b*u%>l%l);fD{+)%y3#)0Xz_x0E)t%d5bVW6=EqG-9Sz zs!}HbiY|_acNQ1-?|8fq%?jgrkbRwT!n3bF%E5(LWzEH#Z<1+A9TOv@qJB|ep>X0G ze#+v7`!DyuCPT}KAgnFpZ>xDxOMXOW&9Wk}Q_}lwb4cqbUajVi;2QXt;7}b`;GW1? zG+R?CZgf^z?@q4V&V*+~+ zoC-9?Ek9P98h8FKviCOH*t(!gIlVbo?gz*`xnV9jQz4gaZzAhn)Kb6||1(pEf3K{< zElL#1!?^U?32~t|^|vR-qst`fKK8Ep{qz+1X>TfMH}aiTs#ux)S|-!I>XE3a~pcelXFM={}q0X$QMB(Ik z3QYbMqyII)i!X(jT9)TrApxRW_jZQu9Ee(5-7V4kaxu;ifdCCz950+Y^Yy`S@MMkw z={`#dPuGyDaS@-xQGzsMn39m(EdL$8WnpV)hMFgT9YGMo>^@jS=M7R98GRgAmXh}| z0^VKWqNrvj_|Sk)I>Ys%q1msZ9?39MQ8wF6?u0f`Sydt*)}tp0Hu7D8@Ag?{DHqz-?4>7&G| zNkcF?q1>K@e&!^Ww)N*+o-P@ZP(MpzN-RB;Sjzr1cc(OhnjKhm-&^FY{d;cXS7t^W zkB)jZ?W%Nu2Ni#}D`ryr!(2EZk@fd31|++29N>qG@0zq53n|nH1gtUhx2)#*`9F6d zgic#xZA!_eE*E6N+vLd|RO>w|ocj5%+H_4kj}j}2v-^4G;~o3+3!?jvFI*e+v0F5m zXFt18VAs&LUhI(X%(N31?B5$SrAzi5oB;&SeaxgLWWEwoInA+>oBeiX@sKT)D6(7A zN5O2gx1z3>HgJk`Nd3HOdXxz3_Df_!Xl)t(a>lE}ubWh5UsYoyr}myg{Q}Msj!7y|j^-KU z#n(9Aq*QTTlj7E|L#jKEOEjBX~wxY{h~vxFNdtS2CUx-VWPYYyMu{GTgQEC#QnpTQ+(6IwLM7 z2t6+kC=TG%iH?g`%Knqx9XjY}w7)QK80m&C?x))Lazc8c2@sqN2rgYY{NWZpE#;gi zu08MXa^Q=o-d{OUOgX@eGiZwApRG9E$??7c=WCgIDf#C# zZj>A4LR{{(WJpx0ktW@{H1d+NHA6EdaH#ImAQwr*;GnTezI8=8#lY3)!iB+rgCy#< zx!|@fp}g&^jbpR=gj7hvpR{NB*R8Gp0|iMOn>N%pse$TYLGv9t$ z3OT~1eppaZx7D2cNfe6cSuvbSnrjNOLJpaT}Dy))_^uvcjWuK@A@I;|Hn|nJq54;w; zKDF@?t73MNTHxkr@YIz zlzP79a>gl5$Xt!mUs|_K$_A0m`&@zj2| zh<8@rb*i(cM36@I7Kq$goND_>gsRh^sLI9c0OK|OjfAs;hr9Lq{^qfXS#qIDM>5dc zbQ0k(7PDU7l#67=<^>?!n-jedGZT(`x@%FI=D7s|`oKP@M3lvEJc~Bn5DT0UAG3a0 zN?FWWdxxoyUV;XC{fUoFjLN#y^D`Qf8)-An*^V?=W!u;;KU!HGaV{0HHuloBV)RS- zy<7MI^U>sK|6mAN_r>_9?_|Zj{e5l~Lh7rdYzVKyMCNeaFozq1%ljSw{p!T{5-RN? zYOA2!!y-eD7>-$9!*8~xf(^|z+N&e2w$9H&kz3w={^HUv z2umWp3+7I#_7^I#8LY3OdIKv?GUuL9J14GvrVf#1sO;U(MTUF8jpTDIv3*6he&^I) zB+p&63X&~(#IT4_^-0}DSF=t%-YRXq!Q7_6p9+o zcjno7s$3wc&3xxGU7m8P$thRgE?6_+yg&%7I5HAd=(FYTHhw$b9CK(R^k#fT%#79> zKOR@rftyiVNu$daot{d%#XtM;UjStf3UsvO+J6O4ED&QBFefw6sYJp$=~e4 zUvrW!x8##{D%uN$Ba^tssoM_|nQi)v%?AFeiQ;n47Pq%`wK?I?pp2L0aH?fFh&|?N z-|j9_H1g~I;8a0@J68HsP>4>5tG;q-O@(c)@of3X*22fmE2*z74#Ih0oH4@z!ZCe5 zFPh`4ot4;^T&EE!02@l!o};V9KMs7AVO4-$H~g+UF8sf?wrj#5lsec>i?s9Y>r*Ah z^F&|k$Ou=Dyq>Jl#(L5GX%W|jscycmmXn#BQUIQrFC1+|1u!adP zTx=-xySv_7Ni0@o6`);6!X&Wsg;FCh0XCRoHk^$!Y;-Yy2b1fso1;0rxDr+1?B^EA zAxWQ5408M}o^|2UiQ?HzAJU?f@0prFb23n*ik)bbAMZ-Dt=No056Y~~jT2cX92pCh zoHN0XoPRzmmY#+BsB)UACXFO>XkQdK{9!H(98H)ue0I_Ov$U-E7uC3F#!%gGz z-t~)L0!-Wn3FfbvhZ0<@vLOqkudUvs zgNYsCkt!T2ZJ=&@Yz}u!xY%bB=hs6-`j`+anMV{}aykz%Gab2?3-HtBJ5E^&z_8%; z%uI?z#3LyVZ+Ir#2x^(C((`V$Sk_rJP37%CD|wFe$rYzH^z1kvheAt@Pw(GS0~-Bh z+VumK!vsdLwMEFAq@$$kJpoDU@oyil?04T5*29pWj;g-tcdD~Z&9FcCI>_PAn`|^> z`Q@Zxi}`2M7{%`nt^&7&<=`9%RM5BCilHl&`dUVO)+Z<^U}W=v+RV)xO2g20W1#X^ z*{H23<5&&vJ!i8q({)eVDZj@XqiE{oS>os&E1cU^IqtI>^KSY6yAzO?Ot-arl$uNL z=(IxtMjS_V>}SFC^Po$p zJ8tD1UVtKbHEVvdy9z$Tb;{uw<3wiesZ=24CvuOp)=3k$!}A0@^wq=7B-p z3zMMwg9<8^;sINgH~%1liYdmCER&vnQk@lE-V9UX zkoH{e7SlqA@xM5~V(EdZnw6r12&h%qKwP_ii$#(WL#EiTD~Do6<1bXLcW;)33a?E$ z4hi^44h;H`KVOo#cPlUbV!gRksBO;c5b_Q1g)K7UIc8%t8YGs->1@qcOp0zkpaD2y zOa2c4Dt1qb7WK0kjB$F68#A3fZ%n$Oxj!#64~i*R?m8XC_KETSZw&F(eFehWm?EtH zMN@Eeq%hJaU$W>>Z>9=kIxEcfFHER^><>(sVKdlK7zjmYHivmSi10W5p*xI1w{STg zyaBD0)@@*mt8m}HrOA9GU0H+NLNAB{D6}1$3HgoVdB#5;BEJHrY^iabYBKdfKSPd@ zcPdpLj$_MJKZRz7@knRX$m}BfPqGr>-al0@=>G)76#+{Tm#`<-l0x+CG`Fv9i@?^U zIHOm>s`?G`#UfGj0`%iCUfaT}zt8`f8E2cK8Z8pP@$E#56aAmKNM8dQ%@K}4G9!-f zZ@oI`t2`*@A}{XJZhbS$o&PlHga6^YG$?koKum~A&hL9G0Tlv-{-D6ABLcjBVGBih zaBk~?!X-Wa4g!&}jcJo^$}0wh>#zZ#KNw?%7l5ZkNv)d;pN;;6ct!rv7CYYx?OhWa z`#^1`q0jAk(Z0CTVn5C&Xl}}&mX8KhzlYfeVNAJ4JY;cXa~#MdAU%WT{W&h8BRYwmiVHISoVE$y2WCl zRS~PttpDdIc0;JsuSk);w&j>5CHzaWH@%B|kf?}%dfYli!9Yw);_=WnVa4k+E6h^F)^p&A%X`9{7uODJvj;R%EwB~+p=vcG?*x>oVUw4m$ zrTq>C2OdyE6!?i(PM)SL`w!8gcQ6G0TaHDbZDru`s^)9y*zLprp`&p6M;O;T?+nrq@Ri`K(zNCGr6@YxL;kj{NN}39|V@1Q=gAB|G#NN4C zW{y-(gSBC};C&BV79E@VP5kKX)}0^UHWP(OXQev^w8x=*+6|@^4T`b;Kw#y@XDqG8 zv&Y}5utz_7CgO}szx9XMI&;Sy<xw)&gCvd&!Ev{Wa6-p1ye~&>{6OG6vHTF&1u|K~5rh)VFMwTn%d^PNYk`Y)! z)KN_RBXcDbqgUcB(5BXKCe-&$g^6Opb$n>8fpKOp6{~+1^0ER=0kyvS{2JhPSv~CD96J8{aqFN5oXv|fGM3K{ z8Jl7pAb$Sa7C2)&Q}BYa9}#7-8?MuBswjrs`A|@Ym;0@BG<3}|gdaCS}W`YT3F7_zSUc+(Y08WV1t}LoeJ7#;+ zyh&J_RK6hvw@In(FPDC@-?FI$*!`VE(eyC~e$O7q+<7gB*!^lu<++99F9)crh1Y() z+G45w+sCU)2Y-9>?!I^TG3tAT+t*tIswj9*B7zg38J`$4{YfH1?nE}X6wyBUZX{9w z@)N^4d65R6X@gnR+sbVCt4x%Z_IV4jn-6wYmhea9z!jv`BUHDpw70j zsjILPQqR8(u0Qm&2U|#LAF&yqdCy>t_VVkz39{=yIRi1Xd15y$2YrNk8YU+M0fYk! z#LVIQ)gS|8SyugWEDU&khebxog-4ZrVtx~sf)W*i?gtWek3-T7+|48XYIa9ENii&& zdG>o|t{GhWjrW;e+U;AXwy}mSl6H9(!%d+={W!(l0vSrlJP$tU;O7{Ck6-d{-?KXL zl1W}#c@yJJio6RtIb^t&&e()gfXnUdJ)2VLDpS6nRO60C$02vFlZ%-g<^xYbU;ww5 z_H_p-{FmxJ@pJuWA>}*j6!5{jQ_8@c)uiYYWaA9L0qLquDNw_c8_MJZP%grHoB+UE zkW3GXgh44cAdEMWZ=wNMc7z4UoD8~u8PsXK(mH~-+l|{@1a-dh1h`)o=mVW7 z%(=tmbcF56P@g&?*H_^Wj`?=dixdyQK7flJC*%Xv8M2{ygu!`+|A7oI1g$^61~7Mh zt*7e$5X3_DqKG{dytG?;XO>mJm$D2$-=MO~0}?23r2SL%xj&ST-SP{Af(51jN1XtE zkbDGI0!Wgr^@#*{Pi8c1|Db{<6#Bq2JQQ%846eN#O9$XD=XV_FUjbipj5Ob`SIPqn z_g5+uw131Mjh=}mcE9r&&$UICXr^MR4T{336MDeIERKsN)>bUT2qK7bz$ zqd#$dSP{&$_rJ(N`xZ6KO;I*ocX3GhXh}4SZfuE5oddgA6c#$2z-vnxO_LeChq(Iu z4ESsvdyBLKL}>-g55DH;te%%(Spx;y2j43}Q8Q&V+&+RbfwEQ{*mGcaAYHsoS};Hw z?u!YE^Phc)rb}f5#T6_2y6@#?j}T9(j1i*s5uUvYg@d=!-~jYKJi^8KYHzi%At}3K ztNmi?F+LVR=&1_T1T$>zNk~@!!oeLe;mrs6EYyTZ-1L7XPa_42`1S`ksg6O!ZL?%_ z2)xH2DYrOLuJ#H7dt6(I*h0ul8ghrh+e2E2Q)q?;5(>%`#$1X((cwwR9Hw)%^YdmH z~O-pfV(DvSf zNRj@6ww&q8bLJZ5xWvzGZQmCoL{<PS%VQ+{w?T`Agjx*Y|JyX9aC0)&n{6yp+!vg?F4z-f z)<;f-`WciIk{^d`cU4q32N}pNRAR+02E9<}78d)uoR>=9`ZCmxFDNqH<{#t1Sx3!$C~3%@}jfM|)XRqOOa^B_Gc5%uH^{!%5nuociJ z;h2g3a6Ivc~g9J!;56%HeQZbCic?jBq9XYZ*AAi*Ko^!(AH?I)eSZS#Ze?aCg}J_6z$;jv`Y)aiptH&uOl+?Vsw63?1U+Go(}C(=0$G3x74DOGm2^S4kZ+9xoFk|JA9PacaJ{X==>!tcU& zu>Ly96S$U`b#|5O;*H&>>dyeG-IBP7pFnBJJf&~uSB3WWD~-_tCB9V4T-vV*)>ypG zep@{E{dX8~4$Ol3sEN331?Z&7?xkII=$kOQi7fE)Sw7rKnLGF(-fEfeV+MlwO0TEQ z+tZG>=o*xtRxwAf%y%!!b2V2NxRY0fRQ#p2vf8;fljMk63!_wJsJ@zlQ=s+cAMW*n z9gmu6K+{BBs)57!{rceXoTv0`QoZ*6%M5~)zM1Mmy0^pg1|9|5I{qJqlJ6bEYbzO= zQ*S-gJ+Kq$|Fiqc8qA}&p(qiDA74WntjG`!T3pSPyEWj}BZ;x2#Kr1_bd`TIhf~Jk zv}rUpFSDmj8z4ZXNipH}iz@Sb`wU*SKQ;Od82J?&gR`DB;&M)$_`R}icdH46{z|VZ^0X2 zZ%l?}s-^tn!K!N4PbIgUP<)T+QJRUh=|?EcZD;*^UYAA>^#hNb}#i{KRNz z9b!rGZCkvHgcCKVUjk}`hU@P1s<+a<1?9wNc6g!EX1Xe8D5GuIu0*E{sfn~kr@NI^ z+%YKGTb*>#T|Q|Hj4sIJ{idTt!eJSx2W+mI-u8=(B~gnN?RUQsn~(Al<; zx@&d?gXnUD_ix{>-(@f4zzqaX*?p|Wp!7i7Bmnsz2(5og35cs>?|sZu%Te1=vnStJ zRP;kbO1{p{?TP1$CE9V`hP%pYJ$>EmfcIn|mj`~9on>Wo%G6;D)mJl^h$)SczcTD-DKyKrktJ+G?zmWsNgPE1U%1gcu8V!&sfssC|8G-X15 z#X2Lm?Z+^v-N4(JwjV@wXZ@qZyhpXGRTA&F_QwhQzr0^-ac{2ZcM`M2}okdmW_l|tpV4gw{*NQJWvyC-{cobe_Jr`_nHF*%ZjgF5f}%uVgaMT-1+5^_aF>Y^HpmvvNPI>dV}Fr)h2?#lWWCSbK$H zZ>@m6eR>VC75Qn6jAs=112r6R#OT1pYr;57S-WgH^$^?p23@sOTjG}3-2#CeVg}hB zeom=wIoH90O2z0xe~bE3SzG@l5A3-Z2cOp5!dKfj+=MUZl=2@^UFMn`jqvnX_+1X% z#it7i=mbLYWWOsq;b9}$pr*9l#RJ9yS5J1lLzS6=noqOR?a9S{nN|qo*#O8?k1USH5X)m)fv5iS{_SQn@I~`*Uz69>-X zpkrfk4J0_f-pXnLPT-Prn$%!$UVuAJ-1_=uN3+r42R;|~-Jjp3oD=-o+pQNd9g50_ zg)KcQNt6?Cel6c?Bam6TIVp3N<`ww2L~Illlg;weF=R0Ufvg`&jf3mLLFTwl4?8W3 zSr;Xj?j!E*3F5U2m(MKVLv_r{LV*1lvp-N--b`&_fP$`sJ){?Xiug`#artBy{*Qgo2B`re^O+an?P z*7>EbbHMC{+=Dy`De-{uM+tx~Z~xY@izCxrjPscCvZbB~0L@vC=jny2gQSyOW1tQ= zq8)Y1Jb8JXF`#!^Ao;+{*8p&@Ldt`?OjUe$h;w+&U8esYsfS8gQrQ&mwmh2Q9j(bu z{f+%swV{$Ue?Enn{Yip1MGo%~+VU%P%E1ZH>;>P6tiq1M!DMyLUNLwPNMoLsjXYZT z8<9SvAPZhubZ)rrMGK|)8uuT(^nI$ejksJT`+ju!1(2uN+OM*2~5*q0krAIULS z+OIDLR1kX3#lSvVEgB)X&a1x8wU{|p2p`S{CDS3nFDJ( z_niw!VpL`I**KPr2PI{&XHF%_M%S)8v?v1R0+4OSw5v{QW5D%-U3fjRuAN^Qhbc4( zpn9Zq#GX}W;^70(ehl^I60;5{Wv|Ri+GCvzj*5HBiL&j42Q{?Tf}c6ddY^?WtF?+>kkwDj=oizKTg zQ-;}b86aTaQx^8|Kh2c4SjPcz$AFszEa6;@T1B#}TwOLUTy{+TVK1eomdb3c#0Hu{ zEPdR>2CLoT60tR3zgiuM)1@EhrFGpLVICd%U0(wEfnK>(w!4^zJzDUS106~w!)^3v zq|M@nfD4Y(dU+*{_m2p1i0nOT8dmgI&k`LGn7s*hN`H6{lrqH8Xai<+E7VF_m#wvI zYg#Wq>!t^7ts&I1v5B0MS1Me270QPd6lOetk0juE-@ zQysint2#WpJv7wi_easJD!s=vk_P83)!uzL4LXE=(*Qd~9v)E@^mF7HzoZEI*T|Bf z);=Sjc*|0XMg8E37+NX5XSEPi{!#2(gr72h;T}=-+5RkQH)d%Iyr@*FiZUP(VU-ohw2g()4;>5Fx&I$z(W4;gZ-b;-}Jm+O;Z-nsK$jj=s%;2rZl>h2iu9xYn3onU88k4QvT02BsGuK zUvX#yh5$!S!Fz%ova3wUQ|YI>WL!Egh8k7f<^^4mT*dmkgDqP0Ch zDjUT00SZqa+Z0G%%BBX?i-Lj_}wcRJ7|j zt#yL;jwk%}K9WhF$Y|T+n<9C%$}~ahuDN99GInUeMUZ5U z8K3_-APv@6&C#V)C30Z*xj>biHn33ALK-52jG%B5Yo8m+wjUdy9c_^H_NQlggkJwO z%07EaBVnQ~|F<9JnF{EKA&FGBw+pY2Y-%J(hJw|`CD_II=1Ao@!%bM>^Nh88dx(#Z z4&CgG18jR=Q}%3y7Zmw(*JMhf{&v3u5&eGN0yIJK4oQF5tM?4Mb<0IJM2ndmc(l8a z^TB8g{@Wa07!3USNve+ARy0=fm~hPY#B?U}+2^rd7WdiXHz7yc2m!FqYXyFFcJL!7 zbQ(%oyh7}tnL2syO472q)ZIX0bi^ww0SN@rwL8zbGtGL>MTbTxa{*9l?g}A$dG2vr-YJsBkZ^;@SzUQt+d^JDP8pxB! z_;#KK6a7=6(6>|HMd77*n311k17?&*vS2T0DTPv~1zaoS&TRzBHrltER1BH~uAHw` z+dRCe)_Kcq7lMu0)YG@dS}O5AGwb1BY~|r-6WO$b|kkq{>f2~ z%qpD@i9_g$(z4@wN{~w!^?_&%CvmqFQYfghbY>`YqbJwy4&q}bCRnAV zSTMZJ{8;@!6@9O2`MoW%i6i$V1wD$8mWO--A!)3FhfS!%r-U>G8aWqhwZGKe97QGu zgJW2qvI_Ch@X{)C-Z=-5E~l_T-GA22vyD4WQECjO2=>nDc`W#$8j*~QbEi=$8j>IO z8po>$Ixkptb396W)OA3`imNvfDasu63#s70Rf;odTO%7~G8G;D%Biu0 z(~z!$Xj&FK546u#h`w)XGl%h@u!fqj0S+s;MHpG>+48r8mFo;GdC#Qst{ zrcR*59LY>sf1@hUpm>$4X;k*v=evLng9mTi{69W9;*+-uwDkVPw7(f^Z||DzlQEHn z@BEaWFOz#?*7E`-0@|-+tles3Z4_+YsWM>~>(8UWjY8r;i3eqB7xd=kBt)E)kFB>^ z!2M=!q}p~xkZb?_UbA+IN4MHMzYqVJ*R43RZp*SJbD+3tiEBhjN@7W>RdP`(kYX@0 zFf`XSFwr$K2{ACUGBUL?HqtgQure_CSjr8YBS6xSo1c=IR*73fgT`i{&5TMQ8*+d~ zS1`mR<>#anSLOn#>X_t=%-XnW4Fdy1%~}OpUapv=)bz}}>Vn+FlA_FtT3#*%5HJ7{ zY9IoGsi0N^k$}p8OpL3psHmz%Qm+sb6Kz~j5=*3USk1s{97JVRMMbR|CbO!l7Tt9i zUI6+yH7^D3-z^yvR|12V2NYV#`8oMT3RKbLh*2~7YYcaR1E literal 0 HcmV?d00001 diff --git a/paper/figures/chi2-full.png b/paper/figures/chi2-full.png new file mode 100644 index 0000000000000000000000000000000000000000..9789603d1ae1b75f37e6013a43aa7790bcb3ffc3 GIT binary patch literal 19615 zcmZ^LbwE?!A2ucmD4;MJOb`J9>6Qj5r5guIO7|!g2GRo3B`qCebm$OJx<)ffy4mP@ z??CD>u$Bc}1gqU0mD6SMhb zTu=6!N9tC>*AmH2@2-1$bO+uuvpQdjzUD2{3>`ANBMT|Rke*5o@-M_RO^}9h}iut1?$Q8Y;yT6?Va1rU6tofZ`Z3DT_D`26K06H z$4q@qG0X15z_ZgkFEkUl>j$En$R+r1$Uj%$vhN!A?YMT@XSbO`hMz;u%kwIJGeTRk za*L;*=_1Z|T-f6SBKw|mI8PD=3;7FlD;{$g=u=4W<0X{iQ9QxJ`9Nyz^*x`0pAH9S z{XxRpbQFp$z1-$vj4Ub1<2NKYx0vLgAB#}%GyJ*@?e}!Q0SbA$gIJF|myUy@wUKgI zPtJcAxXlan=HS?w3@c)DZBuSmlhb!+wYYj}8}Vkwur{4^$mj(ycq+z>>y0YGEo?#H zd_Q-PH)r15@I)Db=yPSgIKI-_6`O5@2jirI@$@gZRS7GFGq$z74(4)$apj`dzI|K7 zuUdD6w&!6=xNg{=%`;B%NPp=A zm)qOTAh4$_QE9T$=NDknGA=gq~&0*Ng_MsE!rRkVT5{q`r-s)Tr^D zxrc!dE#t#knSYY@i@~L7W_rt$L;{f$*GE2=)w+B|EUmaZ>sT$9c1!8ZU7gpD-|Xce zG0};K?_|Y*{m>ec~_&?Maf z8I0OvGTeefRrT?aw2q(W_v_Gmw$nfAvB~LRw@8Fcn`d;0bRK)88w0t~g?$EcgC=OT z<~KjshR*KIAH4_p-~5Q$7H!yEl> zn#&N8R7qtjNf-N*&584dJU`0pol&P`4N<-Kxy5tnoQWo&qy6& zz>N0upy3y}roX(?ftlXS>%6;uNpH6jL5Io7+R@|4t~=>-J`enMjnVz4U-G z*c}NkZ>@?#iuVy>y5PBNEfcLm{U2U<*qAuC8l>c}x|ahvjR{7HS!+BF7|z+s4~rgG zMP8r5Z*MMit^mpA*V{O*l;`JpRxw8T+}7e78IN!c;mj{v^@;-Jvd&MTm<}@zPcm}z z;{5~c?9n)^HnfG#A7d8{;ueRF-Q8FV2bGIPkn zdK!Zi=EVzpHN0DuN~|`-`l}CAVuC^Z^V2KAcj<(0aoqWwpOeN|pWHM}NT3pwSg3oP*wOH-&b1LfVG&8V=qN z_@|WJnAs7PK!!g-Pkc+;xvPE8%4>4D>R1ecc*z!_Z7?}*?Zhe35oiZHpx0cOiZ=Wm z*W_)Xu*!c*q7F{wLOlnSYJ=0Dvq-IXdv6Zkrcv>7lMJ}%7BZGlMh$TR!>gN8tY|xbW5p8IKGzKp!un#e8gUp=6jD;U?Rv+^2JM<<*wM)cb6wlHnbunFeVHuokk@ z@G~TND!t1XP!ftTI-6u`KQcSl(hc0n)|=cS9jB^#|AGu}OLc_BYP0x5B=RXVG82Bs zUS&6pfh*_XT*D%LVH~r1SwpkDdbQ$O*#6V*9%U5Ytfr_Hzi_;IJmE%oc*f^)ycLvc zkCX)qSK-MIWcvP1A;N{bl2jU+F@_89&F7szH9dQkN4$RgP(Y7vl|gg^zg;!Y=NBJL zTDztvQxy0zXNOduct=1Og@<;VOzIcb*Iil5H&xOxae<{i9plZ8=n!m&*i$7IlCZGM zR7YBPDD!$fOQM~;Pa~&dZ7C2@K5Uw=Fyv6euXV7u9+n!+$Se=T^dFT8R~9AxC_K!} zaLU`#=l4Z9*5@nLl*yRy)g#6y8ZTDMWJ3D9R4P*8lxlJfOKJ)McDJ6eFet8TkD;_% z%r(Xi;a>}9V5uGFii=Mjud^*PBGnZc|F&oE`c*|$L|?g>Wa|;>uxX|klp`q9tor5! z8KGU7UMW%#y(T30?pltC9W$_0)0osohtCu^i|f7wTZC3I-cylt&)KFUIP& zT0_j6)z&hM`Ig42ZV;M<^*L?xxZ?))S}k3v>)JsLRhlkuRJy`AkyQyi7Al3Cs;1fb z+mLGFNAq_=kptw+Rx?ydCXk2q^#ijjdoQ*o8$)FeU1X6{NF8+DQX#agWY{8iCU(de z2>3t~326HeQ471Nwy?79E-n<|yr&?y60=xqt4RxSb^%6rI=j}4q$_pJlNToUbLzRP zHkNC-&Prj$8&w!?@nSsjDl7DMll3m$ z@~x4r#>}dA?%uM5Pil7^yx)mK)P3LLyvgNXbxx>M%S4b)k<$30^BF zNH}H|khpNZpgx4^g(X~y0qtL~w4)VM#ID)aBH+5@GFcag+8Ygk>emD{8un9F%MEoC z-U+V8D0lu$=k`U>dd@YN2Q=l(MHBSU>^}yZT+C78x%tw!ZTRBXT9|%h?^;O3WLR=z zZVl|^YYU3;MFk(&^v@#~iqOVf9;fI63LhwYwEkjjLSIwY!K?wFj6(POJ0O!#IzSe~ zOty%4r43{#w{xm)eT3)jwoc*9XOPvGDK%7Ve4g#>V#YE)Xm@M+l_01PHK2oz^ry02F(zWxXNaaE@;aeNiujE0wfL%Gp#Gt1Cv1|GX zm!uK-LToNCR?>Tp+?C6eV&=&CmsSq!(fED^4T&^V9aCy`c@L902UJy9a&qk56s8A4 zqc=2bAV$ug;GQA-Qt()w`6Dp++kfskUw@`j>Km5iF48*)2CcopH+_%p zuDd8WAO$m@%RN%hu;xMr9`Rz}Z=QY2uZcBJ^-tHRMd=>{nZ&hSx>o*zepnXokmt_j zsYkVa9G?CN>m$V!C zzVeeN&xBdKsCuA1@i#vwlFieE(7|gFZTp1xnS)4FK0g{^1$UYd+EtgFJdBN;4L>$C z@T&t6-T2svkMVU#!OUJXpSuL{Yp7d%Cbv84wM^&R!HiiSX!4v#i7k+sEN__HWx4M$ zd(EHV$zIJ9YZYd@c^*g&Fku+S_^$YWdy{yS_zDeh$6AXE6@;f!c0n zmpQw1b$tqk2~q6d_;LYj!~fd6I3@>OsH_88Du+(fvH7T)J3QV2QIOL;NG6-Vwjj39 zezmUyA9G!XBRs~LPD+#fwsB5sR8eJn=QZktmA!n$TcX!HxLZP#y+V3M>=Uu^1CNPgI>?9Q?v#U-6u6j!t~26K%m8|os9Kcl42n8(eX@_>VblB8%0TFrKI)Hh%d@=6Ob;&xKvgr5`cR zt*-D_BdO9a{3E(Ih@Rgs>;|Yw+J1P@YxL*L!!~)PaN!ZMO4h@kw<3^#f}PYKcR4Xo zm0Z9wuYM9j9v%r`K80m+c@+ERPh_xt+!A8^ox*9?%MH5>WEa}6U}NP; zVB91?_T~pUxaBCN)9!8GMeNR#yy$| zch=*fnO#!A17@LGPhm&1fo~lx5#HbRDe_^~&Y8xV3~)&b{X=G+3*W0(qH`M-QQ}ml z#bsRc$zA^wPda`E55qG7)~*Baps}=yv;i|1{F>8f=aaZhyOc`p@G##}qwn=Y+jpBl z(R_IH5tVRSmAJj{8Qaobaj=c8-+1{7L~2=NdJC5c+swU9z2N2ZkL*cKqX6|T1x<3& zd>(wlZK7peYkWV?W!?tX4sM#=w0~;zLXqWiIqK0p=Ab~oEz&5UrMZ$~ZWAYu!$c_l z-#!wauZ8+|2aQ}%;5C|^L8K(26NWF)>>6GAJ9-UW+M5Mzl0t8b*dApDAKkjy5tH2? zUN=lv>nF%Jf;LWhtNM4<@4Vr3GodhfZ$kNUlWJe;C$A`R!6 zE(rfUx_3TtR>xDOkzB<$no5%S&I~`lRoXZxz(4G6cLTO{XIJc?7aKT`C zPQMSnF*dV2q>yYD>fzv#w<43!apE_B#~f@{*FN-K@CqQlGuUIz8dS$Zol~D6(o@&p zQAfFr_)`0K)ZI+ABRV)wM~lC)^`wZ^W>2f_WEZ|)W|fk8<_ZIsaf{rS?0=Lz;N{$2 zD)wu#zn9QTn*|HP00qw}Dag2Eh_f@%EB zdSx^#1u&o{aeG8;NC9dj8+5(4X>L4WB|E$GB3#+c=dr|ylAL^@;1PqiNaeh`|%!yQrtED8~)ivoMBaelf#dZBBm;US|_ zun&WEQkHoTP2J*1>+t-l#~y3B(fKpq%^*+{=L3V#0?Sa93`n5zHq?1~^!wlR+yuoZL zmQ=E_C;qrNZ_q@b^Eh8JVl0i~V(GTgwh^&1M&O>OuC-qQ4E<;0+RTYx3I_lpciS62 z0+3P7U3-qV*X?j{c$+t@p)=3#;NT=o(*=8E83h!vf#PbRhi?I)a=Y}5an9Z6aBdlD z6Z>Wx5&<1!xVNw3pTn6rzpP=g9GzkhfXG+7LfEe0r+K5qPNCT=(gKLjRhht*FQp6# zL)Q`kq!;**6Vlj6AG(j;Vh*}6d;1iteN)6~L~uDws;M+mO!LZVD`B%HaxZ=>B2d%U z&U_NYnZ$Kjj;048RKRM{MJdse13hgy5!ruuKa&~`OkcQd;#k%I2M4^Mm-?=V?fMnK zLf~_Y*-sr03>J(@jdURASHfPxVS7BBs5C>(06ge@*(5o=L=-o!!$^lZzRMi+)=6`B zEd_>U9B^=^Q;oFc0P`-^o@pdX5;saNV!LxZ_>_+z@&0KWavC8f>}-hhr;U2T2KUi@ zFqojU0Wog94#XVz*aSEwu!mj2pZ-J&a;hOXv(|>L8fzXWN$^oa7sYAAH^g?P8G!Sm z&4$|qXF?@z-|hMzy|?|;TX6l0*e(aC7Zh3tUB80wHy4`!(TnhxXWuWV0EXbK z0d)m@z8LYADhmS_o^(H*96L*k3&AD8P4|lN(`4d&xl~&HX`UtM{N=F!y5IqJKhS0* zxO3X(4N;d)3;|Gyr)qK10P`5H;Nv!v(D-A{^gOvTT--!=6+bQAI!yMo)thQSGazOX zObbv#z{iE;#41A`FgNReUHB-@1Z6b2^Jz@yd;G6-`8t(8Cjwnp>i(awF#;7)3nhMH zefv}vfwbta!uvp&nE1>GDy~A1cv?F+->DFYlszM+88(Ng)3z_tsjF^W1*ST3{tKrU zKJfTpHC;7QH?~gl?EpMX4+7lrn{7CH1Qef9)WcG8{+l>BAGm#>8(+@@bd=5Q%m*jQ z1c|^!nTy8Qn2=L6vokbyQJlkhT}}-8w>cuyFEhOB3*tbOs_%4 z7J0?SMwv9p!W3>_JyyTs1|@4QQM6cJIn4f9PxE%_*uR*R(jzeaL*r#GOKsWQ6!gzY z4bqvlVRI&?+*sm56|?Li>{7@ zHo$+16O#`AFi%6YJQdK=T7iV3l8?uRznd3I!t28$k@_l3ER@OluU>Jj2}MgAXg(M- zLVH8iNBN9VVe^6ydL`Vn+7)M8yhBR?G(IyIWF{100Ki0+OeZT>NJ4e24FZ}6?<~mN zVWu9Iq5xUF+)`Z;ilhqcFk(O&n>1$U=Zt)|ZFhRLk=Wm1@e0mon=iD7KrGP4(=nG$ z(KPqv3;JaaSDO&b+~vBR?vpJ~24Mtr=KE zl6n4OT=`af!68OgeMEh+Wk{>9{75#GpsM-~GzRVb4~6L~-Ok-+6E}E!BMU=1QN?d!wNl*kbrv z63Le|R_27&;=Qk8Osmwz{nWhcy<#A~UK=y2r29}pnVoh}A6=aq{z8d-*M zX7ZPmBwAM)hL(0ip~3OrT{C;wr^C^976-w>m5`*5yh<;rlV_-vMheJa9spXuc$SO~ z_SWMRW$lHZd6vQ+So<{$ie*mie{7Cx3fNzRwZ9*96%6*_E+~2irR;rB9~!?Vc-%*R z>^KhRU4fygXg$D_4x#OYAngtO&h37|DyI33$-C-6YBwSGH9LbU4>MD?(?*CN2$!Oc zNkYpddq_dYlaq4y5qzF>{1=?pMGGO2aju&*WX2wgg2yU*C%Hln?)!%I&Yny2=)Hi1 zJ|A%x!4v`E7?%&q?@HIc2TVNZalei? zO3S&P1^3%>^XU=|c2)Mgn)6?YI43u><+Qgq?rqNh9@5Qe1x1XQ_)G}0dE`WJfBl@% zoqUMzrW@U{Dnu11$+g=}>ZadV!jrWbPgb(1M&BqMm#ep&2slIuZerPqrP%1$&kuMC zv*ojp{YT|NN(Sp8tiX`E;Bm1MVVvV-mSb{&rV&j6_)fb!1tt&bQ5%Kf1Ax75^nhSN zmT;$wn+L7z8HXa`CO^oGQ)MfM#-1M^f18QZX*Jp5_r|?2A-6Ml_=FBzP!wyH zD>lB8)mncOfg2Pb157-H*6;e%Xn_!~Mw9ELJcr8+6MsKZ(7ZCat895Rytijr-oFvG zVcW2axg>oS$-7@jxz9L?M@LfRbZn5+e88}CS2hN|=m0L)-ZWz3+yq@+1Bvm5OlKv# z;q6>2k8)kB8ArtnCD(oV8VywFTj3f5FghKYuC#q!0pN}Z-j^n&9A5tFr8kogOoHd$ z1CyX@pidFe`ALV5)(~Mn zX9LmC!ne^tfc$O&0Xk&4fe^dWSQ%|+T@63hJ@#5>vxO-c*nInCgFcv}TPi|ej&^I} z!sFWrAahhfL|)efOYc4tBw{ZUOS2AJx4P{!<+Aj7ZegH#^JJvS;WzC-_l=BtBg5^r z=O@BYjX`-hQlj&y{Vt4qX)ZEOeLV5;>ze5vCZ@k~HM;EqybZ zA+t*Arq9nOS>9&Kjl^S8u$Zf3yx0=~Vam&KaTX-%Djgo(a`7C`dl;CrVf0F%j(rI| z+p7^56|%aq50D%xx=<%;;nHqxx%MS5JFe_wiB9o-ncgN4DD$`fsb#Q*Y?`pFsqA^0 zoEdVBz;qF2!OK?FLU#df4(7C+lJD^fgrhshe3nU|`2pv)4 zkBXa`+2M)He%5_3F!$427RUoU3<~v+)EDuzW2w*nyi)(m{gN4SmB4gXWDAj3xIGh=1JdBYk`5VE4SCFfkh%o4g;6UD3kRy@iK0_9m4c(pDr) zXW=GUnATMb7~5{|j&pVg*3RL10HnmWI{yTu2MIhsWq&A`!dZyYUwyYEt6^Dj2@}Go zH1KUbGVy)4IEExDX@GEvbyTRJ*=iSsXmZ)+U_CTX!#8f~QMItFZf+FlJ(fa$q*n^m z02ZpC#wj)GYPkUd2)@yy+R;mV*;uxu{&AI{p3wW!IXU$sjvW$Yf|xv~J?sHtr}z;| z&&{gZ&TkGW^&wLR?Y06rkyFa0nvS9Eoi5Dkjl}|v?Z69VDaMZM7c*fb=sH{;{^$z+ zT;dN8rfbdT(=o~7lBl4*=`q=pqn*PQq$#=elpUbKOUSX(Dww5!XmPaz*c|M4 zw0KVH$kChh(ne9R)m=5grwuY$BSpo9x;mMWOKaBWB~wb0pQ)%Rh!z zFf0Cv*CPcq@CM2JKk=#(9HRv(R)yy(UBwd$1Ox#THNZz)!}l}qhA=oPEWGbD#MaYv zu_kHrY8NN)n%~pYs`nvcm>e>GNOiYWMdPC3vScn*Tq1q`x6W<9mw1-Box;P>qGVus zBA7-uu-zdU_s5UF=m{6^=4^UkXCRX4{$2Gctazi_`i+y^Z4v^qT8EcH$^mb88T8?| zFS^dts3cY0`5{@;oEK&0NBd5gliw#H-~E%Bnfz2g2N~S^tp|g!NMg>OIR@ z_~I>eG-h60Hz9b^-wHdn`$)oH_Qx%J3f*RNn7u_lpG5WG`X4|&K64sjo*ssy8ndD=9i`(p4=3je6TCoYkmL=iwYM`TbUQ# z8GT~=hlbdh(B+*8dO54TZ&P<^{0u|4(zBi(W?Uthrnk5x1ElfQMX z>+;V8EsN?2nI_Zq!?v=TFG0hFh1PZ1C+pS^4C}8a(zp=)?Fwy7V1t-rr_}y?(Zuk3 zzz)F>Gi$_p4l>PCKrP$m(?-#A+!N@FZCl9GNG*UlvQ>fe4LGqRc2IO>^36|nT;xaq&y3^Y6{1023vTh+< zfSCvo*QQRee7ZKs$2;;_$P>5%5lLV|h;b}Xelxp5Anl_ykVJhr*R}Vy`ORwO+o}I> z@svtlW=?2XVniL=dk4;s_@7A%AyAP#PCvp z%HD|$4RlC8yh5Nd%cFpKEh4~ZVM1lujQ6+Qcc&yYP4SM@3GwbDjdt&T+x%h39k+<= zzkGgVZkQOr`-vPnuaWru;oO=ft+DMh%_`BfB>g}S4&%EDH(~wqxE`-PCa#0%=jmwv zF+TVhg;@PmMBJdk?0bCMMhIXD=vGp;h|y1i$^Q{9t)yvMHYS>K&rwkTZ6JO^R1pPC zUrO8XX>lOr_`*;Un%`R`;R09PgM^Ch79%X1aipO3Kg#(u|8@ZoEuiYXDV?1 zAI(=!@)P^6;I~UP6W1JxreOMFt0?z(f_(_VWAVRPit7E|x5+QueS9*xypJ-|4-T8ZMuqYAv#&kF@E4?ve+Z^Vp}TbsCUG6Dpc_(}~y?Q}^<`j8TKLXNhbibe}Ttq186*T|mc zN(p)`Zy(KdiQ2BJ0wY>?Bcq;gq_;^I#>54phLZn1xjNtIgWwe|cIPHN4oLK=SnlWzLG!lr?%k|TakKKB+8c6#0K@~|8J&{sc&|M~ zuS%VV0R5Vod(Xjgc z$9oFF>e`E#sVx&We2?zZooN$~;zEszR)V?Oj=STXk2d%ugvs{PvxPr*i5Bilj7h+6@{ScEHXN|0MT8IU)h>QAISM~W&bCj# z0YJ%5f=&q$xPuQ=v`&UfQz8}1K_$B;r`+;-sIf|S^~)WFJ##_HPQ~CL5NQgU>uw`(i;Nc$#npZKgD`+ zUW!-WTm-OS!Iewz1OGtz8KxOsI^^f-|46wE`usi|K)s!RksSn*Of#pwPzb|$=2KNq zxPc|ZWf?qAeJc40_hIp=<{JRYeTHF^-v@e)17wKJA5{A)^D|<`m2&)8hs&qEXfLdR zq`sy^D*xd;uzj9rl3)0R*?%!IZKf>5$^gcexXIZDeJR=Qr(hy8k}jPyVkb0()M60W ztWQ9#&Dhwbb;fo)SLWO%ya8N1td87lrm8vfKOx-+iW$c;9xj#{@_$H*4;^aw01r9= za0dTfON|t>2Hz`n!;iH)v-Yr;HhT9Vu=D1Pr5_YnCeIWRBX0G)=i}U7U^%09PDlWz zLk>$Q=~fnz^NXJ#7g>99l)ZdlK+DAX}h9=y6I2E=ZqXnpRt z!@59R1T`}gPILC=OU!l*pz$-9i1E~d-3*#nA0rF^l1UU+|DgydV&HQKxc#$+b2~7F zdi0vZ!x!%PbW8D zLm(a$vF6|V&ACZ6TbQfOk$}O>)RYIzri3^+WH0^@9rvD8Pn)@6xjg7BTaU=K#M{8Y zc$b04;@bh3bMJ`&^GhN#heKF#=PiJ=;cc~j@?o+NNKitWGZNVJ>N3dn%kr!1zfXG* z3AGROUEj7D1>OUo{plQnduIu-=A6c#N#=)KH0PQgppRbj;2Hsl0@35_2HCE3l1&(>l#l*Na zzidr9`toYjIjec*f}PRsev7pJ#EsS3x?`Oi4``!!wNtuFW#TFw42)o!_cAT-PMXbA zyYB08tXN!aZJM*!T2XTM(z1FDXB1lDw5_v%4`lAYc74nWfXgWk8Yv3dJ!O3_d)%U2 zE7`Bj0o;g-yKv8Ag_dhi*H#8>*0IWpinKMPU94n3xr_!JbPFVb`~%c`+V&hwD~p~6 zWJw{M0DXIE=#X!b=sgKHWWdYm11jJ({qwd4K)_g-XGl)N;CXqoq3cj+j~*1*FZ;Qu zq(@O*y&8pHZVi6*0Lv~+qr>8eK`4()xP^u8KQH0uIy(ByMFBW`eg$2vk$ISD|I>#F zN`4@M+;mOO8&d19KbseD<}Y7Gs=6i4lOffSma1~FVY7!!S+`Y%{fdWub^zi+k7_4k z5&UDS&2={A@I?w~Xge0_pI_MLT6wz_8CvrK41W~5J{taHqk?oN!&#msJZA@O3O6;f z`2l!l#*__^vCb!%(_aP?xra@c9Dc3=yvIsD+g5sraPOvUdwi`(nk!*Wc?uS%9{|V{ z(-~NK9sQ8ojJnSTK>G*^(=7l;77_=nf5$O;d_hHP zNco_8RRg}d;6DqovuZ~zDsyP94D&$kmmJ4l2A^LYCj%)vuc%s#s;fb=ZXjWmdZ&>y zLWiXSk#pR-^Z6`taJBc?z~xnl)J~g+eH7Om&*xTkwCOIosm$X0151+h`-gRFXQO>{ z^@0nChu1ovQ;&DV$MzQ071>l-qzu~&%r(v{dLf%s$bY5gOu^;fsd4q$In%&eYh)2t zdc&AOq|p}!j0JxlWw>BBSAEZ(YYcxfmK#a1jVy@A3> zh;aI>*pxXDi0T5=#IJ;u;q?Lmav|Nx!Uo8mHSgLZ%Gvkb2@b}k;IQq?@kbjQ2X??LS4jw4#dT7O;WD09$1VH7 zPC3tZfHE5)8%}yaGgQ8dVoT3 zm&mGjOrDdGry5tgM&z-prw?5PHJvkr8Hj~%KDPcIHTG_84FV}FA{fV`sgg7dc(}}5 z+>RzV8fH*;BXvN@kayq7oUtfv6{w}=5p3&O>?NDkaW?rnn2JVE>%ORN-UawM+O$y4 zAVw2G{nk@EztXySFKl$eEi>(^HZb32mu^c4G_E77$qb+&Wi zgdN|ioC_D`0kH-<=A@lT(sz-{rAnc{24prlIHjD4&&6ZrF6-dk!oqr>GO|AH)YMZq>r@vTNicL|+)U;w~BkoVBd)OiIg+kE4fDl_-auWGn zlREc;{h4*FL62wT9bQu%_}9C<&Kx#RojV@s&THQ{bZ$O>)TyYX0^>l7x{G>OI~FFH zZV8SYXdZj5W>_0|P7nyMVl;u$l2E}GZo;>Nty@$FY~)nbndlFu2pq#VB|h&EPKz{_ueEyIy2pEVL% zHSCiWXI~a}kU7Kq^n;_p@Xy?HG~+QCJ(lD;nu{hUHOmanp=G^v0y8xLtE*52yen<*7o>$e@lWu+d}8WKFoE(CoY6!5Pe z!B6R5Ofc2Bx#GC*lz?Sd4n{gi>J(FunZpSf-EuYQDro~o9lRFi6mn%bH3vsDQa6+TtGD}!tz17)g0@(=qsDx)S;OP#od;K0~5ZaiAD zjEjZYJStyjDFU&_VdwOU&cIU{R!S1u2pn!GD`Ly<(NUQQfzWr}NC@b^ zjh%`Lr&6>-8p>pH^oP?No%pZmOB7vy&b@UpG-<$m^t_^UB{H%-pya*dB_r$l7iqD0 zr1>)(UeXYQsJo5S=r+O1Ux-%zh?R}BH1HXQTB_lOY8DcFB zlowppq0m)8SYV_+`b4?dIhha>5tWhV!VK^o0PX+dqyAjAGm>L3YQmv(djy!=21yB# zcGM-3JuWbtar?e5$%{)wuX<$GL7ApKp8OMSc`qGp8&cg3d%|FwjB-n6(p(pMyA|tW z4xBKG@tZ>}t1;*tD)%e}SE>~y=g|X9WM)R1+uea?^O8)ZnGM9oPTgVp%j^e0wUt^r zfWdcwQyxoHOR))ZHczz=EPy95yZ}FQzU#<4X0=*_9IB4%ciuU+loYN&+^g4@haOS( zx|@44v?F;w*{9^AGVki^vpq1$apaTk%k*^1?NZiR*>Y0Rd;Fm*t#Cs#p^B$YpjRck z6X0R8-XZ-o1Rrh8FjxjEZD_iuTN>d59Lo&Av7}<0+;=h|zSt6|14(3-#WG~5Hi__|Yc@!pE{NhoLFuSei^<`DC>saSD7oG@X zSmqRtat>C)tY=Rjr4?>t!vM(^?^NX)8}MFGEJArz6u(~RyU=)auXKDcu<^jaxAO5s zsHP*23XEF^ws`nDb??>Q)|`ppw;bLtu7)hYw)Q3RAR(>wGmZ2T`2;J z#YwBE{oLZW42`CF)b^?P2Oiw?%pQTyof(v93cWQV1di=VLG4p2)YUq@^1v(QZ&{wrg;W)prg?szQ}W&h50}Nv&UZ#>x=og6 zS(msEvmq`~Z^8jV#`i1|iXrz8QfE9o_9)^Tvv{4#Gxva(k4WH@h5|=qdZgwN4Jn4x&y=sl1D2fkalt^r&8y2l|^>H z+l7ZY&X%mSu4}EGyUg5;#nq(k^vpqp|A{FB#5BFW&rltExn5PXkW=#R?RDiPIGHll zbO?ax4S?<6Tmx?`=&Qd_inGYc2B|^IY72&bSye1R+Rpg5Rp%`5Os?c2V@(a$@u7)C zS+wK0bu^UVd4NCp4!Yhr)vp&g+&TxG=g@hmXyDoOiv3~8Sca)0umJ~qo`deR9(S&( zs;FACMmg=98sTcyUwHJj&AfG(ZrSUi@OYEf>gN9(>DYS1_AR^m*O89l=*goXe_m2Q z9rl|m_hjkpA57Xt0g)&Amd)j3bQnPRTqPzls$-S#VDIOoxOC^>4{|Sqk!ykn%FeiT zz>#?L@S)2oEJqhT)2Ab|(Eg)~pGSP-mIfRp?K%rBO}>mDJ9(pSV7u<3@W|if2`1ys zi$4fE9nkIBsK9K#^XB)j$GjebVDP9Qa(fo6WDFeXDE*%!9YjRe{~*8@v;piL_Y*DO zA!K!3gdov?7wW;q8dtX#zgNk7qUdScu{|UsC!6f5Vu<-i70;HfZ?3%`$Bg{X7|HzR z^&w9@NZEAXr8W2BacjFvSKnee0;4;Hj3(IA7yZB&Pj(f{#4&usE5MFb#c8LQ+FO`pd%c0;{`?o;3T-4c>pq7U55h zR-Y=AK!BAC#sI!S@;UV4bs?(kYUQb)#;hMG{{XJIr5)Euc6_sycpIE+C zT9_;1c#I5?n6MWWJmq{ezpspm`2(k8OcSXdtm4BQ8S-nEwR?18&Tkzfy7U%tO$64% z$%%+@O{m7~CnjE}@uoL^9m>Rn!~=&s`uS{Z%AOw@x21zfgbhG1Zi49FctkN7ca`I} zYeQiN$>JwYE(HHr)D$VFQX{Gvb9r@H6sk44JVyWk^=tAC9qxFxfAC)SCrbd46DxbmdR)GUA}#-Ec|KBX3RFa}v==frM-DIU*VxnbdZ zymD4MbaZ=d5>FS}&%!Eu*Lc^ERXaf}wVaJ&S#&%>uhW_=&fhZvvywkCNap-KC*^`y z=|APeaWeCo?hzs0ReTi$jnqRJW^RyV_%ratDZIT4Ff323A0xW{=^|hiw7z3e_x)Cn z51W!AD{yewkEdQoVe7RBA0vk)YG?in9RoJL!it;Z|5+EP)gP6%OXyZ{Sd9T?`UpgP zD1Y^^4LIddRi37r{5|UWA^W@Ebkh*4t+yK!_P9Grxv*LxIzYlttumA@= z9R4wk$N0HU+*Ho(bRTH;WQzbhI{xS*z=C-$yF{%A!aR0dZ$0`$Ob}vgix$TjXCcx9 zIB4MyCWH$kPO65&(#`)Lqzn37KEPoespvEvzbK^DkkO1G-;U?7T*D+$HMRj$WpcGJ z9G8E+`RZ?$=B$JF#5O5vvn00h1_w(*Vtc!OW=V;560f292#fFWSPY@~&L}}+vo&y{ zp*O{T>L~g-1jb(4(f&;EyIP`vK(&q_unC%ayBbwJVp6_EhCPc?KE5AX+_`F(y2-i5 zNbrxk3zK!aPQBkg`al_{oH7dwPO2D7)l%rTb+Zn4lI+F8{|=WRfpMl89Bt?|I^|u& zCb|qmKU~H5?#Xzmq>vr#E!!$}Z2ORTM&SNk&R`pXfMch4d<25CGOd2ya`M^H%wgEJ z8f$qn?$r5z+OpYa6|Pr%RJP*-4;_OVCT;8AblIo&D{SVWOtBS8=KQllTq-IA9xM|# z{z**G32c~ubXz~=p79-pfzn*=9WpNyfJT}9ziE^oP25--g-(142&T8*PzDTEF-BMm zUDy=5Vsc^0 z+}_JmGTM$SVh?HJnzN#qgK2y<*XYaT&TY>`&IGgChr(oKb>^OKT;A26q5jEzha<5OI|SXAe8uX*hj4Jz4kvpjnWZSjEqde8U~&vFM+n#STV~jx3UXn05%e_cWZk19{D?XoS{qhv zx=ZGBCI^r#9t5J49eSLRI=W{!V37&zL8r#kCnFxz{MM@%kFi)9^{mCBq?943Tf$RJzK%CrlJrK8Q@7e?Sist9Rq+x;tw<70R zsByh=lWugjFVjk6m5krpzyp35NWZ3DZMRYSnT++M3>$}KOZ%6bou{z%mk%+o;O~!M zM?aexYN2UggqL(G`Ph}j=gz=;#0z@0XD%Pg$-V=UuXQMsCNHCZ|ANs;On&A}X0z7k z1r}<;oF`3}&7cw5q-wY@TS zeC7;8Hk{_Yu#mf^cng9rEd}?T6g_F6{qws=Ut{0)wQ57T0A9o>4RUp>`slHP{+a!_ znO$|G1Z-b_()k{{PEZ|m=jD<;KgI1|5V&{crh6xe?EVnYnlX z7l1qZZZ0k^9`wKdb#?voF9hOoi{737*6%xzKX+XI=U-QkTe`YB|7XSpC9yE$DiOw2 zB8;mQk{+LDZr6XTXs@<=NKTNHSuQ6L9 z-g_hD70rvM_Do#;sbmcPo|=`041c-TzTKIPxe*HpHFbhVd8MxU3#7}xkM~cpxa@7{ zka)>ojuCg(a93aRmY&J?=E%E(pL$&F_iE}u;3qNF#Y-xrt}px!WKTo&MCLgyhPx{T zFqNE1SSX7y8xt=NBLt2xT&2IgidG~-8}4}FHhCP81w2naIF4P$7mTEsk_8<%mEaG} zxGsy|Z9k|>&XpnvzW8qt%4|B%1KW}30fhlaNvdSA8$z0&#b&e(AOSmX)DA>-s^3&m z?B)A)s$XK-ec&lQnuzbzoW`=Oz4%;I4OISf!+H6-GkGr#sdyE>FsglqEtnR;u9tk& zq*XnDbWrWz@VlnHVldyHnGT2BL|Vda_57*2rIj@-aRfYps$!uEHR0E{Sy828%;z^u z%O$l>GgtZ%t@(Z@jqIBbm9ijr3?ChEQ@dvi{*m6yr?YiF$@*CH3ph7N<3z@b?SKALlL(l%_(iRh~If z5*f97&8awR6}o{F%sxQIy!%JBb;+VOOMS3w0bAPwX5HvYhj~XWwDIw$gF$CvB4ZOu*BtvPk(03@l{dniJrf#IIY}V_3`narbujXOw&G^m{c#=2N~2~FuT+Z zDs{W)#2Ggmnbs*?owXDtiYD2%h6SHIbfWonmVQ|#>BooCspx|zVw4W%ram{unl*_@ zL{CrTK^Sjg*uuc)dz_EgOAgrw+qg)!g<68Q^V!e^ZiVx#@2ob&-gmXMMvFkqNwqr7 zP4Z+$Pi{xRM{0|XVjYI}?<#Zp-!tAP@4k(7*wy>ew&#JEtto4&6PD;7#OW5$3S|Fq zwfl)PL?I@e7jpWST(*^(yNn|^CBl2-t*VD=jzh-I@dY9k@zRs@0Z%P^O1TFPg%pDqA%zb|IPEJ+22q@D8cRYww$L&wY;BJH@x9> zipeZJkv{s0r?I!2rS--C3Pd&I2_65N5_fafb!6k9Gdp4xQ=&08Mg@3Yz-c+Lk7z9G zWSDP(-AL2YPVIR{mE_h(P}I{aDkHCXsv4>cu9sC`@)SX-jlD{DSrVyNA?)GKuWBV9 z_L|h;+$&X>y$V-?x!v&C%U(w!2gki!uZ8i9Zx)7|(g3C634JzlH8<6#Bf>mUkg-&I zuo+*KBJD!N_f{A>j((PC0?e8BAA0$+KYmIpu`=tM%KZb8Wv=*Ie;?t#7Tw`s&6%T9 z=Y7C5y#il);Z*iBsle0!d}53f>|^f2cnoeu8VSs0_WNGGiig>DcJ@^HXT>7|NY=H$DPqo|XG>J3J%2^?Qr&3M&b3YvZ?dZQr2EIc z+Ux98x20z+%O?Bg9qBNpr8cAc0=`Vzyj3HRJ!Os6-JkZ7vRy5q3|tv>E!nT=(<;UG zFe+(R+gd($R(pADFjD5vd^id$jF8Hip92-V`zpogw&_VEE1hFmw)WLisNU^!NDC0> z6Z1i#;3xjIVr24Kl$x{hb*DUnVH)%PZ1L~ZABa3n`^15T|La34%}S`tJxAkMv@es zJ4o7VEE}PO1vIIn*0=A%S~7@Fr>uwiB79^%hijAO;ezLOU*opBVg=&L^$=;7zrBEN zXwUp|ej>dVxR_QYHp7P!dHhcFsJ(q$??oS(>#aP1kw$k;61g_AM<|1*BN~UL$ThZa z1QL_b{Dsls{ohQI#1)K_X&e#Afn`}#+5=Iei&brfr#h_1j{NwLgxjnr?5YCw!;AJLZt*lyBgotZ(*# zB(f5^-N|iHPs4#UBvn8pLj|*ONt>k>Y_{6!5_i?ycOwW+&d|New%|jw-(; z2F)z-)#MOABeRa(;$nKE%E>{FN4Cza&FKyKD*1dRJOqatP{k|6$LETwR9}s<={{Mn z%d1}FiqW?YUCRyY^D}v{LE9Xm25qXW&R9#n0cAo->Gj|{TKFoB#H`F^IFFvrq-T^CxO<9A zODX*My_xsn&Ac~nKkVIm?w-4McV>5IPNX(Sm5i8?7z+!FOzn-bE*91!o`>2U~ugfG!?M0YU4?6t?(Xv9BU0#Wvu&%_>mtE9sxv6Q-xq1mmZ(`rA&}qBo-FU zh?=s3zR$uh#9`v%_sMX!w^CSWrmkGkoOa3z?nk3;xU;r7e5xMmWAdAXZ5gCuacd9o zezc{KUEnc6S2-``6OPloPddoqS|3N8wMF?y%ax@5fy5Jj^5;SE4staT<7UH8DF|`9 zK}SAwNGvRVOVgWTi~o&jdOTn_LGtYMk3QytM#gwo&Liilo>3Cw2FrvZmVS~zaV$B{ z`<3Q&32|{07`j1jd$H2k0Ls=X@I%gEB|_j#zR>~h=N4Ldb_zm%N~nQ6~1=YK&3 z#p%fZP#O9EuO-5ADue8bjI^XY5H84jtax;IO5|PN3GbT_r#Mwhrzk$aL_O9#$CP0P zonKao2q!w@Se(sAO{0Mz!aro1y(s|;OQmN=9It-LWYVpQ|0K05W;L!eBsbIGq==#E z2JLpGBVq=+Lt|x^YV7(X{yFya?mtFdYJ>zO%yjX=Q*c?INzl*AIp>Cnd69-2&wY54 zDAK`+nbu__fSuicht@i#)s9OBJBQ@D+Ek`f<^2pDnS5TO3{UW3^e(Dm;^+BE6=Jo5 zWY#T+DJKN_CX}QgZe{r!2%@7_XRmS^ohI>3>-%b?rY}Y?U5&qZI}h&KI)c2}yMXS06 zVlOJ1YZ+>4lAAJQP@p^3dJwr0%@5Erzto&14|n$d^3|@P5yzU#1q8wk>PY?qu`-JM4H;LOuzq??k&`(okNRjHavkBifG&51^Qk_NxDCF zlVOvQ#jyPisO%0mdXDWQIxe@m{9+{TU5mk7(gh8_wjNS9JdBO_J6_323QcuQ@(X5r zo@4VDCOaPm{OqJN7{x%fQsG`Ne5~;ukYe#>(UWu~eAj;R)Jr?` zt>TRFUk6C6iVa`(gXn7Y#y$Y~hPJmvQp%ws_Es+6ZNaH71eQ{uRB~{Wt}Yk&#g;5d zaQBs_F;L1Je%_t^=r_)DU6O9lr-RxggiJ)BevnDMJqq`;&w3r`=biwJJ z^6|#;3(QvQYXLNB5u=Rf@>$iE@_MK1_v|%z?_$v3?*^%xe?>POgc1lWu)lJGs*V!Q zU#@TlG+)!f&f>P~~z+F)4#a>v%V2FN!1Io=v{;Ydim%fRW3FX0PZ0 z-j*;>=PfR~J;y+VnQM}>=hY-)<&I$94miL4J`JfC>ixUlnN|7Hs1@BZ<}!=Rt?pVK z027L~{H%S`NM=R|G$W0`Y^fq_P5gE};>c~cZu(B6dPRK6hD)`cE+3{4DR{^kb`BfN zlH3qm6%Xh3g0lXB;_%wEt!|>LfOLT;)5?1L5Xqocm6z!l;DQJ|#P$qG>6MzpH~obu zQM5FcrR63&3^DdP5^-CSTrnA^NkK=$6@mqNi>IYhe&fdUcXWocq8oB64EUFQm|%TG zB^EsCN6U(FK`8rR!TiMhIGjM}iT)XhC_KaYz4 z+gvYs(YKYRJtUrGHz|5CB)~e_y~8zMZ?zAsMb;$un6=d;FFvET8!@hH1VhstGfK%U z>;Qu4R8%6!xZhK{$i6~OQ>5!-1U01ZOF|%kuEljxOh}V|kpU87r{yd%{*1zvep=WY zdot^EEs;vvboHLc*ofJWuNx%pecB6x=CQ_ilnRVh zOT%^O*+%Cr(RCGZN`2g`31KgRYR2Nn4$b3xF`n`M#fD@l{hN^*^|cz@?FChpWo$Y; z2D*Y$pY7mT(h9q}n~efAP-|!Bps_-SxZw@V^$H14yd^>s<}@351WINg0)8yn8fBke)MvKn;wagasl5_WyD~KZ=7AoAKf^1Z|2t$j2+ff{71HQi zCTq6R%|P9a%r*qAYn_DFOsNr(R-#L9Li70N>rPxp+uYPtDee6>W{Vb+1~tg+7|ju& z{MIUz9IXxyjfhM6L0+>o3M_lA^SFc4@QWd-uBqR9Mx6%nHqw~#{$3ExJiWo-IArdU zNsG%PNL`UtY(+BH4p)9!+HQ8W?`Z?SLqi0Ylo5^sF>_7|uH0l#S8CV7XQ)0g6n{*!W(aZB=zYi&}Zeky1W{cR*h{w~B8M`=-14K6oWz!(zJBue9 zEe`EotI7wtCyp%x3}Np@vck#yanz~b>3$dFggsiidudscz`&Zv%R8zC|J_e+epXph zRDn*9Czx{~^mWP^JX&8{G}&sr3g?G-;Wd<1qw4<#Q`Sz^My6Vck^E>5D+-T>qmIb5 zK5W6XzR5CNiX}GvRj4x9No{MUZ&e>(Wi3Co>woUti#EEHGSts_2k?1wWw@6K3TJNZ z7VSj#dDiN6#R@~auzl!b|LRFVS0vg*;tDqTl*4@0%s%CXr`m1OTI^Jf8>OY%QV@)k zD+m$#h_0Hf+K=J-$R@Gbs;kIT&}6zUjH{fZTjK`DehfQH4wMufQ5Y^k^BM3sb;|&)5DI!aOin&Lh?sTN@S`x2aYn|xpyR?KB4uCTHMmHu>`MD%l}PVi=+1dS@A{<_NB$ds5%GB&y_`^aOF zC$Grqh+a+(=S&!choe=4&(G%T$r?0ZjpgFkrXj66wpE()Ok(0*c&$Ao z-?|uxMYBa}7K5H#?v?)Y1{Dq`b?S@4vp^iR;5k^`{-i5s>YAv<>F`An-iK-NlpMUs zv4&aGZDhKN$i{(>d`$eKGK2|^+S|OY>9K+B8ezjHP2J@muR- zE~*kp1Ut5Gr}5PPE@@7fE*Fx^`aaciq?Q=h;{|HsIcx+#`I(jlZA zZqK4pBFtmpAKBHEM=h!hB7GzxqEB5`Id=J`sCW0hfCIMI>h!CDxYc~;KfApCLP~VJ znco}p$O9TLX_{9ni+1vHh<3=%lwhWUHu@Nd>d#awX86+u>tF#-9_X%wJv`^vP~?x= zY)ip+*yBbf32_G7!r93hZN|aV;d!&6${|owmw49tDHauY8$tH$3Wt+l-|=@U?@7y7 z-^Q|}vW>&9%mspV`pmT%GEf)t`C_mc$Z|4IM(`k>uCR1kI<>s4bYqR1m|(e zVE&k2fnDQ({^DFR(7D<@)!QzHD^fO25O=x+z0Z=a#pRPBBd)I7$?#CFPJKg*gta;xbbA7qkY|qUQB_I=RJ*abhf7jM-l4}`y|&-OGA|H@ikd*z$1)4J-(Ds zZ0h@u+J`DX6d)s%8urg$<|*6h(H?H#biTSRFO3!RI2z<D!UAlem5W3IM-v{!3;_AFfDoj6r-wGS}9;x-3MEm^P`r3H1J$_^O4@7(OX~j1} zoS+xHh!#r2%odTA6B!_1#TN&0#V`?WSc`zB`HS_KX;~_&JB-gK+|EcrY$tQTngdj9 zES(rbK#xo-;b9q7-Ovxpk*dSmovBYI*&@-T+>R^WT`3r#rv&!x*Vlzjdg4Ck<@$r$ zu=z*rQ9$*bvC!m>dOs-o098ypk3Mj-xzs$8ea)lSMUfB#Ws#x&Vl-q_-fgt_v|~*h z!1?fMsHLcTZA_7sO1|OQgVL%-7jzXJSo*4mIVJbMze{lrA6>UE*l}^;nMW=pChgt# zbByB!GmnXZi%mo0tG+$48)fBNZ4t5~R2|f^vE~Ws>Re)_rwy;AfqZ-kl@3VyYcu?f zB9_~>U;28u-_kPkp3V-}L5z9_79QTT$&2cz!=)in&hWowLGH6NS$&gmC z?!k(|n29nFoZ=2zo=r~8eI>De%!?qoyGiW z*09h!gM>;A8s}ack`fVZ{Noe1?-#(wUPr>+0_=FQR_wd7T>bBr>icUP7tf9LPQTLL z(PuY1KgAS`h|hRLXlJy`-ij4?ep-+}j*xoMCR@Uw!7g&#nqvD|TcAk!ws>iAhF?;^ zSD1YCXu&%s?BRQ--7^8`9B@gqk1tIB*;DxvX?a`Y7Hv=pKXC5HO2 zQS@iM+EdilOX*aCT8zi8a&XoCv}DFWlnz2&22gW`-p#B@6O;F@P`WNCC3cc}4u}3A zv5$0Oi{vY%>^J~yju5?Xm$LEe(9Pv`INlB`yLBNyQDwPHq4Lj2YT~xo4_^IF{L%CM z(4H&!=TqC)seq5aF%Z*w>fhG)uL7UuIx^2;5AqN8EM;EjyKW$dd+ySR}_lxBLVkKl8?mK$Fe6bgovadHyb(p1r%L|yNOfwft zis?@GiN9*-2&1%}6=AnJ0Y&qP*Qj^X7Q!z{zaiW}h%^KuYZRX4Eq(VcD`IQ3Icymbvo@W*X`f1X5(|T^FooW~@mJ^&D*dcVaE6k+@;*QgVsk zrAbl?(Xaqux8Ieqta@}}AYS#=l-Qr>##gVtwRIf5O)hHEQNHn^2l|p4A!HsBZVPHaubogovKA>g#2xv5$=U zl~y^BqQ|GZK&NP8%V%%?PHBGeK2ABEK70ryx$B)2vAn&^b#Fpb7gUC|j1$C8F20$S zK<>Z^yIu2rzdAxOzh5+?#yGo{%O*D>G2J`nyIAm`c?yLOfSQ`?1;YOHlGGe1J<%#O z&iADTKd6`CgCNN&P^ZiZ=`AK`aUwj>h-*lhodT=Jra04x$CA}|jgCb!jrezGa`SPE z#Dg>%01$-T1gLZP#ZIl7Cq78M&mb!AAhuQzu(Ztw@Pok;0SHS90mOtv#Pk6Y(!x^G z;t~P?fHVLgADvV1e-K>VZ5{1=|K9}RS$g^hg4%y8xICox!Fmt&@N#rqq`fH1ORaMX_@;Wf~BScQm#_84E`VC%Ur<# literal 0 HcmV?d00001 diff --git a/paper/figures/polar.png b/paper/figures/polar.png new file mode 100644 index 0000000000000000000000000000000000000000..0f74690ab56542b098b57876ae4df770059e5091 GIT binary patch literal 8362 zcma)hWmFu^v-c7RzK{R`76}l7`@#}D5Q4h|2<~Lj;EN_m@Wp~V!9B=gi@QT`3j`0k zIOKuAyFBNA&%N)5d(N#7J=4`yUB9ZX>6z&YS67t-;#1)R005wZ{98=`01Jq|H^s$5 zx62yRC+G{7jg*QM08kx6aBKD$-G5}RB=;6@9~v{#j_%>ODX7TctYOjM3p{xjuE+!c z&~YfdmC}YS9vU4YHVS9`Q``wmRoNTUh}^)>dWtYJ^8;Zb$K~3wTe}@+I9*56hpokfQH8xIe!nKgOr zb27S-d=qH*uMh*8q~2%(AF!e6<^P>2)@Un<|0a2?8>>K{C-VUsbPdy!sgO@bNYlr3 z=CkEyn7EG8PPbpvcSckbrlBpC1&Jr(5uQ*HAV-N_e4xljn*dZz^?5`C4 zTIwXTmCML@&XnTxu_(=eov)P=CBao|WO?icj(UR71`m=Y?ooS~7zRbgvB{O6ziMM~ z;+Tz3myt7gRNH}1sD!VIL`Uefn~`XKgY(XK$^X{ptezGVpuYJ;(oMPKciq^6_Tzo3=97!Vl2pYRw`vXbDY2C9> z(alLzDoh7JM>~g@yx~SQa*!9M)^LTeLb{ZdUy!{HT*o7O{n1s4l<5vzeLb*U4{H6J zUN{=?+38^a=PHqJ$aYUwKb7_XJ`RBC&LOS)d&m9(0Xk_^Ns>a-S7f3GFaW55*tKBY zSzNBSv-21LZkQRSS?JzXy;xC47M;JqZc*^Mni-)7oMha=%Ef@AwcsosX~;kN0-TY< z6>-y@kd_ou-KF}Wk9jw&qWlwdmydQab)rQs%`)}`Q*K*4o3yHyKB|UIC4Kpv`^)2c z+PTej-Z^ibebIu!=&s_7GKI2g1W!|W;!2Izu#sb`wLeI8io82%u><^8^Ve|LvqWJR zd)#AM1E7&GsT@z^oMoHKsDq5l4HHA~IskBOUWg=($UfQG4KVLYB^9Y{{xRppH``Y@o>Kf}89+}!LE&Gbb zdiEHO7)H)O{#ClEWyN`C64E+At1#6dL<=E9|DojKw}^oisBhklax$sFS!KW6=H3Ky zk?+^wIPW!H+HdZ1>BF5c&Xb*tUQs2MSFK9t!TTLcM7KE>9*N_vV7IZ=2$m=+hX6sL zgVNVa#1R)i1G{Eh8cP!ckX(&+Z347Y%3R&hCG$qj{nCk_FNGH_|2|Ei78*PJR;8;O zQ;usrQ3ozDqpj{Ix)GIG+dCtn$uf{55YobjH0>HS`>BKR8q!R^gx8;0;74FmRthp3 zZVBUTO*a<3NWQtLLrnTq@p+DB7Fwi1t%T5N`QFxn!l~dZ%N ztH;sWkz=3ptCv%sDY?{W0V_=>bd;_qllRrxsF&=qt4!pt6Ukaej11P?2sW@M8pYi$U6pCM%su`kH-#byj+{F6>Yau?MRv=<4%zxY&8wlZwdX$Auy{DhhaDBU;G{XkCCBIBg9iY!r0b_&X! z0$M|{v!RYl=2)5=-qH@g_oj!)4o6Mdeq}CX=Gy(NX`%b&UlDYnk7)+}F#;K{PI|8_ zhEz08i+S7XExZ*s$?4D`g#uI3%~}82?&2Ao6m$-!_H*s|>KHX)cIc+T5oCpHiCJS; zn5u(*!g?w5FEIcE_|&-^ZlcJ93q#B+;6ZRQk`DC=WANi!3;^p%mvJi6dBmwgU*lDM z1owMP-OpO+RBHo|Mvg{pFn&$T{HeV^Olh3GUA_NpsLZMZSg8Xlmx5#C z4A0s?ohwZ}Iz3~<25%-(>O2v7S^oBuRXwz%4YKY47s|gaC#6!J6>gnj$fQu$|G0yj ze<_vr3}a{|^(i}T`dBjrQREZ))#v1vqQ6#xr8clSCFrrx@s6QuAS#+=yMoneL0x zGO;Y~)0^jwz(5KFFHWQmDsl?`sYKbO9SlfG&z0fdgl=2di9`C6Xdp+|O2Ql?d)Ca_ z9_f%=B^xF)W93C;v-(-@($Ces>gmYY+5<(EcQp+$ErCnli``;N+Vr78Y@UArq00P0 zghd}BHRKO52VseRR^zQA(TXIWI>}1a)168^We!BMKRQ9LorTUxv!4SL3+F)TOrwL| z$I)8jf3rR(cEK*sCcP9=KwmE^x#bi5Z8c=~*czoEY2||}r=%N>Ena1rKc?_75tAOR zY@S|#B4z$o!pw$Vz5L4nqvOX_2aGp@AFFOVIk*U!?hl2(X0QCTI*@3Z2{mee4^cRE zKaTD@k0{zWEzKoLW?|q~o4TsJ|)ihB1ru}T`$cFXj3v+Z}U%VSFd8w ziasdt3Z*BDpbpew{=59~-E`zXbl2-NKOU@_fTL(& zcEI&-?+h&VM>IyxJ!1$p72Vdd=%Z*%$}tt@STaJO$+vhCf9-6nA+}k1T#bzd+w&w0 z&7E88-9M`i6=6tOKjY+05IP2BgtUlWSbm|;*Kh6bEpE*b)L@D4A zbGDqC`A*WJIAHtgejWa>gtW4(qmFJt_^klPB=2whNk2D3B6C$eGdR#g=(w(dc-NdY zx$JynK|-GM@`}V`3VF>E4aR-3kS7^U*hE&RAc~4O+vcy=vw(BmABr%KwT8X)Qo;-Z z+4Gra>@@3KYjJoj%zvQRdO%l1X!IX+V%+s%Oa+z=Mb&kE!AnD)ra_F9WAnfz?w1Dr z`Efw1B5_7c<7Ve(ff+voH9w&kmFNC@aCvzaUz9^#+hZ%q%3>TH3+sakh%;}q=Nhhp zmHI%PqR&+Q7eXZR_|OJoLSCcoKheCncrC6`p&d%TjNlw^G)wccVm%cv;bSwn{TR*i{WbF`LTAtqo!-chx!n(2*NYhCc+g5r^ZYy7_O0D~|QpsfZaqT3&} z?GK|FI93sKTrLb2_x#hrpO?z44w8Kdkv4kNMK8*O*^u1Bvx31|kU(C0+RLyg<(&!`O&fDAgPNc3*0lhQE_}bgN=^`+rpFuzVO% z@179nQSaY2ou%Gyl&{d4@RA)$lUlMo*7myKBpX7-LoDD@C@=K)dxWR>Yo8y|7cHNf zsZ}J_R>roHQ6fG^r$IJ-{3>DKC?#9j<|jc5BOZwDQ;~+=oLo^Affv)OtX+0pTVW0Q zgN_1aOnA$=Opk9C4zO)>XQTa58A?U6!&HXipYxzfVZwMg(hbvt;sd@$J{H49ZO0?9 zLOhKSX9bnl8x5B09~;4{G|BH@u^Z~B)|(qFOnOnaOrE>|NhF$W?^NGyiq+Lx7O2Fy z^S!C=d%ly7B%xWp25U2mVv7h}b;>1IzIz?N=6}5uYBigwMP>DUliLJl%{dj!_7E^q zV+Fy(4{7$eOg<8c$9-8pxoM~RipcrHrDGEXez=(Pi0iJrfFX@HJ-&YVffZ9X(ChGI zWr9Am6U0FAW{^ntK>ZgUMN4iEz8%d#qWxhuKH%C&td#HIg_P5h}#fLa&{}wWb%TLQ!m7?vv^bN^%2}oXrwr zJX3nsQJc9X;m~J6=`IBtbnAHq5W5gA&gM+@WpYpA8a<>s#k^iUkfAJEm8+{p?veXTYD0qW_0`e>7~a%VWQ1E9%Cc8K?X$li#tkkZHE?QCC2Ww487wKCPa- zvjZra!xh(CT4!ROZZ0?WMPum1b${m2+d7BV4KDdE8c|>U3pf2(w5z7J&X?TwS<{a# z+W05+bSeJl!O^^Gq30O^a=s;C4>jc{;#b{p^C%%B=jzlxo zpMQVvZzZL)!lGwK&nBBjmD<*)&q1rNTV)Dqzx-_n)LyLi?NdbgJDavDtDj^rlzbcP z{b8DEH{{8Bz+Gvvxg&X`eMU46vnHGhW=B70f_1Z!S;;>!^a9rF`PKso?jlKjauytq zqR9cqb;GTJmB=C6-jEd$xNcOs5|{w_lz7vfuin7N7=o&r#c$D9vgL80y(`dnkR%$= zAhJ*KlcyZl%xkCn3S{M$`j}BAW-%f_%p@ICFk40ytIVuif;6M;DdxGmdIbLyWgGB55-f0Dz!oE7!!tB%A^N>}sS zH^jSsRao-?Sf+RclyxAH-=IwWOVz|-OeK2WQm;_7#g*Z%L^HPN;%(aJN);xM*z;FB zV$V`?iK0^B?q`Gjwm}gXj`~M(yJA6=6emmXZq^6}sX+9pX&1v9J3R|5N$R&wZBYEv zL^$ryh@zE5d7bxwP2|mM1B0I{FD-s^+4d&`Zq_+JjA_QCyD&TkbMj}Jb{iR=M2Syf zadj8x&+_Ny;0}(<(EKe<-yC?#wP1?*OQ4)`x%Bdgz4RpdppeUsD|z8jqaK`QKO;q! z4AmD?TWLb+!8a45DN?T$bIyXJ80kZ0njG>d-jK$<`;ha3yyNB#=_`+AnNXh8aGQfv zW=r($i)ce&pSlGm64mh9N@?7L0`fb^JMdmZKYSNXv$JJ8mmu_2t;mHmZ0(TK45+>! z)dS9puGknfeFTX=qI?}^ewIo{x}zxw%(A{3rXZPav!HZWeY`Lq(@N8kfn@bv2l*7@ zy?=OtS%l_wN~pBQ3r(Er&(}(8`^2!;k2aiq!d}*v=L2D%vfm!RS9p6`wF93UN?qzv z6{#-89Aq2#oGo2aFNamD`7K4P_!a(Wnx@)r7>kwCT)oi5Cih(o%R>&D4f8UWxa$Mu z8K4rV0LO$<7YA|c4gbB>)JnSBXKC1V>CAtWR#hw+&B1`zSH`I)3x!TeB-IhD&!25B zh#u-~|BPS1}n782kGQ$ez)}xX>I*qikPi1xk^}R$?Dp{fqJN;$%Laf~z zs;~QES8pekP^PW%{v17w>|eKINe59EtbY ztn}6tM=jwXN(W#+Cf*e;Z&8fTgH?WmqZDk7k_Y%!Wq&&kGqcHl>n#UIn>L!mzp+$! zEt={=_EpzbW*A*cIP7JFo`)o_|~W*dmCAj7)fT zh&q}OUaTuu9v2X3{snHXS?_qHdFDMK6V~#4FFOAj5yIJ%D{_#9tT7^sgNt!uw>8|V zSyQzFS2M*DMBXr#;*zW;W-%?)(n&plOXW*efTd7;QVQ0~{v^544DP3~BLo6;R@P*1t&fl}Sq;kVGEQkAK2KVx6c?=|fjtX)cDy zF%cr$T9ehn)~!>SnS*aiBSVd$6Z&&LWQi|hEf^OihX-WR9#1z=Q7%HSk0_yJIgPg7 z*KA+n{M>}RzPAl!Aa8$xj0(fGLQdo>4oSS)n6G2P=v@xqQ(85JI`5I)UmW$@-4s2R z8AjOn%@o!e&L~e7L&}|*)a?yy)7R=JLhDzw(g8TiS(n4K&TpkN)ssJP896Dj9Ejv^ zzcA__ofqLOKpz-%@HZ=ymJwRS=`S%|m|^9S^0^q?28BJs$>{vp3%a>uve(=A%+c<# z^X52y=&FO7Le!SIpWKJWxep|FP`yj-EMNaQV2p?DsVs-bM@c20cd8k6VF|aY9?J-U zaWRppo)?Itz{hQ`q!L(Cf0U&a6aAhP#uHV*?qG=7*10LG4%JKaz4A$Xp}wd^Wseu& zk}Q;?wDVG7)Y$flkF2zSaJ6m&8wS+*b}Ds0T#WOWNb@+TL2|pu zTi2*ILhlU|qt{gzYh*@|SdPA)61VW=>zk_q#S+S?``VazS2ZghX|&ZiDa}cDo#RHL zEgyU18r#O0K>+?+_x|VuUr8zUPX@&A|6%KBb45x?i6p&A77&Z+iTUT`WOmdg^SDcbEAuN%(uD>CP;{9oW?c*43~Sj} z7;*=VP3{ioi7NXO26r2ULv4e(5pva3(-!6G)}Jz@eg{Jr(xktWR*h%WT|-!WDk)oR z*|2do8-KCD8;t2+J0ab;UXBDxHgfR611AdV&kMsr4K~`2cH+k|5_U~hrnnMRAu$S# zW|igv7ILbWDH7(XiL!7uVmXBdEZt zWwOr?rB%EB^mRAWrRGMn38D*uEoms^BYwp4(#Cr6S=JljcZ=XlUZSl8P254xhdfE2 zCLdoOd8g4^m8Z00S$MjlyhIGz?z!)Iefs6ZPjpK=g`--Q3L%qzyIt8MP_5u_GBRom zYH(gO+w~pMoITt~fV2b?5++%)FQVBTpj;QkXUV$%6N&$a zA6R}8NTZg+>fA|pon=WYdb@#!(T%HNOI$*Q8!C`y=G2AhiS&U2kJ4jfpK#u;qldk5 z)of!nPUSIt`lvsP4p`{C7gB5sf>bI5IFCEfXfIJ;O^iwfKDKK z;CSp`fhW%&1RR2y5upzPkE^b6h6#rFdL;%Me~Us_%2qiPWJ9@P$x%#ATRH!|4(N02tYwb^=;)F)4=}+Vc*sM literal 0 HcmV?d00001 diff --git a/paper/figures/refrac.png b/paper/figures/refrac.png new file mode 100644 index 0000000000000000000000000000000000000000..acfd2737aa6fc497e6706864a1a12c3addf1be37 GIT binary patch literal 12775 zcmbVybzD?U`!^~fBBFq(AfYZRA;+nSpB8di;xf>oXj#{;=G zh)6KcB{@G>23zo4ChSUf>}S9KP=mCsOQfT*9Mvp9o%cjyH-W_CQ9tr&^t#~}YdJ0+mCao>B;+gG*P*q_DVybaXO^~5!Nqk9<^fL+yYleTH`;Sm^Y=~}(H zBH?@hyz>~(qFjqq|Iu|VZhyKPAA?liFATJN;JNmUQQduz|AQ|}4~A*4tt5I2;Pja9 zlkL}+k=J}O2?3#Y5|c<}em}UTa0k=fxSOHshwFVJT>b2t_Wc7udl~}>_0jg4GJGG< z%|)?2AO}YHujs>og@Zal$uij@*CP2-&`INbU*;~q;rX9XVHxhD3rI>~s3VRTrV!-^ zqx7P&x>j=USVT1!bYHw^;396m*&5Um0L~JT7yHyqM4C8hhfq}QzNM38>C7Hb=C6Vr zyt7sN)gn{J-@dK(ZMT%=v7b2AEy)GAlY6H0>4IC3!Oq3)@ou#4_!q zBjF9ESqfiXlm1uZBdFMqgyz4b=4GdgyGCue`5*Luv-<==*@a_>m!~GHnDzE3Wud8` zDWQ?!7FAW~;F=c*3r6M=Nh?7H_n~2Dfd_bqo8u}BY2-e(N0NKhud2sXQoe<_d7)H@ zSdiv18Z}STY)MKVmz41c_?>4V9A*foZpvtkw3mUrG+SG$gOj>T~ zilEu?MJR7_s>H~v=D@VR@i=*BZR&=y_CDDX9-eCI&yocR!*6SK zWi2gNUd^e$V}1u^i+J84%RMR%AFg?4lm5oXc1O~&Y+Q$bB`XLQ-M#+a(`8gvLxX@` zLr3qdCO&TbRQO*pa??tA3P_mIJ?jB}W^B5BnVN)#(qeuW6Xf_1+TMVv9pR-@kPZk7M)9?D|L8577jLf!6l z0s3b04(mP37*`A{P2Dcm_U*_&X^`FO+(+!ozmu<9_Z5L*xS zqHO9%d*BgdkQNZ!~#b@mAy<8ozZd+ z55nI#{YT}!;^4QkKXhZs#Dz8Pl2IlDRlmDbpGmXV1v)6}@G1Qa#1(v)uOT?EmWIk1 z+-F*=DOYLy)A-;K%kB!w^OVi5FI7RvT{W^`+{)>>q^+=(dPUrwpu8>Ag83q7Y_fz- z%nfeqvMYetajEl|T2<)bCS#4n89%7wGn@i@;eOBqv}3kpeD_6k&wJXXWlzw(ygK6u zs+?_=;u?qkH2q!)Z5TYV z*N&LF?;2-E*h=N9oLdqwPs6-wpmY;DuDa%v(<76B#cS=sN&aM;m`uqg7U;*oJrteG zdqRr-ax&Shh$;hS+Gw-*aF9$Z9@mZWn=ZFn@ZPMXd&OY4Z8IBa0_H!ZM4}s&q|zd6aKt={~%2SX5x7z|H?R zTUr;gE|cj=Y%NTtYU;;Hr~$}3t-;;&Vw2<3&@1_;==dU23TvA72_x|JRA^KRno@Ii zt<_Elo>z!0z;36^&YxTT!ZOOCn^c8sGkd~JNzbIh47MU6Rqu3)pP5$K+@%DqilOU7 zRy2|7R&T&VYw}GHdFo=~V!Ty+AELMNF+Vm~H@IOp89ocABJ8g! znm8;h$h6?r4l*v;wK23nGa)0+o|b-W;vXwo1(`7&kiwKXZ@{GKlABCdE;frcE4kV|nP|S!`|9)oI`#6QSB$0nGu=RRv z5(QMH$wG1eR$tMRK29@YZ3QwlQ;}qfQD0eFWnz2sm zs~IXC%{_zN-J4jR5KLlNx$9rGQd+ZKkFS7w7R_y@ZjoqN27=Yhe;Z$5w}yPUOZp1g z>W5pMb@M>IfJc){hbH7(|0t1BT`sK}(J=85#m5A6*aY=N$o%+E)76||+=8gK<-zTd zocZbU0vvAbP_dHTk7509iKzFBFa>SP5IqiS2tM!w)Ass%H1{^8!Ts0M48jTXxiNl8 z?7@L|;y}*TFU)A{)ZzG0Yph+ju3@K=Wrhw8|9)gqa8ocHu=rxhO-7;$rJC@ z=YNbpW{p{}JKb)7RXFJB(;pHhp28{F|EifWr19>5*!8CxCZ01LkAN&Rv%XjI0RGGG z;_uECtWeHCwwlkH)Tj?ZtemnFC3HKy#tXc|Ic(vWnh#Xq@=j5)n`C*l^`?am-(0^+ zr7a$agw~ixr|L}Nj^4cviO!EP9{Q&Oh+PMn85{khS#ti3MFW~}hs%K$b4x+-A@JXo zp5Iz&WIG}nfR`A)gv`7W@q~Qn1@YDXS^b+B!<%ZpOJDtL^gF11?>6w$=UvccCk0l$ zA0dAgb6CSBE?Bi*pyZjnZbRYa;yaCusz*)^*~`hOgd~spfytPLu;=h^T#m}~Sgc6~hSjl`ly7e1^;b{g z5P5@7Vu85ep64HLMGW$>R7f*=Y05|YiV*UY-}ui^`OcXVK1fzjaq`l>k&>H>rMlKo zu)GB4)INI?**dbN1n$9CWZss^Ec4pk$m;1aehwv<2_NVz9YFT)o$aPfmR763S~Kjh zu`zgKTvo{O4&441Jn>efb$dXYDdPjI;b;XWjs#ifWgQJaNSlFWP@t7^e?_n_;E>P3N8*QQC>F zUJ#+gvZ=&4M&V8{XQfR4%&UH6X+PR-w|OfcVc*{=baJR-FchB-IEwE`FvUy2uYyx< zpnAy5Ysk|D9e!Z@(_;>bQorlVE3lv8HY}_}AwTT|e41*TfcXsvMX^ zDo(D@c=OKbzT3E1P5phO$9<~eR^tR6lh0Pqf>n4GQ!?^63+w;tKh~v%vtagpl%$fP zD-L;7KWOvIbQYjnVm}?~4#l&!*V1u=kz&!wHjBk5Qu6|u^Ch1%&LN%ELM1l`3TukC zl~LQKFlQ9E0G-EP@7g~NN2wgm>!RNC5G#T9E~50)%4n9IB?edm7Q%YHI# zu_O})3!+~as6x5^G&$7Tn?C@;iaWKIt&_Sd0KLo z_2zkto~R))=8n(hgZp_C$P~y#$U{sOrJHqHSoWAf$#+kpEi|&{Z&xjYFyoXqG?%!5 zN6Eb4b<5~-=|Dg8kY+Qm=IZTT-_hEGqFhyZeKqHIvDQatwR#@p%!bZo#ZO6^87#wk8linfeKc%vf)%Vk z8|~v?^|y+?ROJ4u$CREkGFP`1SIL=oizSs zY~Ia_5+rTdT02qEw|0>cZkgD|Pqv{8;X7vfDlZ4EcCGPbTA^#lIDK$h@g5CfDa-Ko zDcLiXFM4mFLqw{kUPD?o`hFg~?HOJ@ra_QqP$bk;HPBFQ)483jLv!j{ih=M7Mhq1N z-a}g%o5!7L(Ob=_2wPQUQQ{ZM4s2UX5JAmdhc)`7Y^)4|LMfOkgk2zW_=NX*h2IKr zBz7}>3bTu7wdC*ff>rd!S3Fz7{%pv)WkBtr>$v6YQX02S7K>bjQCKhbVaD@p3*uUD z<)CG|_}yAeGx~IT#vj#6S1xO^^$QL}lyikc`dD?>lK5Z_t)}12H-@wsGIye6!npPp zB49PA^hS;=5BgWe^^`F_n5z9M0bmiB$HJ24>=Pje66cONq=a zDtHVIza4Cd#js4i7Es5ZP9QA#!z7*g*A+5oZIBx=e-ZEc7hHFMKKh$yq`T)&o$RZu z;j({Ib&E5eI>kccolxZA-)WaocsidYOn3iVhz;{J)rv)5bcbTzZ<9^Egqq0674`?bB?sEKIwwo)IR{wOKeKh|0r)Ap^IuFm^e-&tGw%01XOjwcXzAa-^; za31ni=x9TWx3?b?GB*%9T}Gh5Qw7T{zFWgfrm;vbz4N_D6@LjI6@}$UIW|j=TOFUH zeftv`7ihR0@n^uRCACw(jJ?;K56g_MW}DhraUT08NGNtS>d3q38OzT%jhra)E}XOy63spE$Pm zbg17H!r_Qk@-(j@pvZ=Z zb)QHidt9lHgHGuNwjio;2`xR|;E{z+Mm4l|5_Gj)Kn@&73sM#_x*O)Z^=8urD#z%b zdXw9%8-+FnTU$@mvp?NKQ@Ht0R`g<;S@{_U)E^YOyoiKZ~YBa(KU_FTB;SLq&#Ao8TuttHpOMP0ht5 z^hp=Ti6ZuoW|PkUwmT(sDw?l_g*o5#;gA39=-r?H%R3^mo(qG4%f5G^NZGZ}sw6#g zRP9PRvlE$As!iS+i@&p=tJzokTi3LK_HK+u?5x91qZsmC(nWuMIflvtIAl9GUZdXV z2?7)JGE^iK{*gF-$dYj53srh)6Ej=Eg=FN>k312Kjp?y{ahI~}hLEeA@*ieIOfXVK zpvNQp7dBCYCJL(vBhZE%BjZlMnkRxx#bwqdDcu(fve8t0f5dpADO$R?5@g<`3wb^< zHr}ZGsQJi$ur^+BdI-W$m|J=PEpm1HQ#PtI^T11c-IBCAXVdr5X?avfK+6=5y2UYU z2Q0|qy*z&myRa)sFF}l$aitIx_xi?LmyGm)*turiB9`1#jw+CR-}AU$Wpww@XDdHr zGI>6rc1C3Kq#`=kPgA!ijHF{6^;FQJT;EHcg1AIryjB{nmYjNs;L>lUQPU;zmM@FT z{VnF(b2^BH7=q%=@n))uzHU-e4fONdOPwq$kD5QJs5m?rw=6oW>>%@zBXtjcfc#u_ zLhQ>ck<7H`E+YkVY~4@pE*u68mw zJ|K?Hd z%!QaAQk~kr=hndeIz$=VW)dH;t z;d!qbER5_-lZjN*C?M!BG^q;F%(l~M)oeIn&0?b=AJ_3kzSw}l6y@TCLB%b- zkH?SRv|sQGlX^5|*2|rMlFd3btl-=+)rZl4w7>d>a;XRwb;op;R>xuodWSp*5N(~SXa4onE-!15ST5w3pR$Pe8rIn!{~Wnmg^LBur7 zXb3Y}ahiq~Z$2qhR&o5YKpL{^cwr+qu?PpX#TTbQn}W2-D;Q6 zvCYqEnbhpT45aMc;cDHRos4vOoOaBTTCUR|HO0J|CR_xIY9c6(15If<_t$7+k1s%L zBt0!jFt5(+Z?9rs%2cr2V^$6Ekvot3W9|I+&!7v(dDce^lo4;fAyDv}^zHIZUoOSg zlaj%2!ISmC3FdvJvTHdu6*n~TM7-8zhJB0`RwAu+Fk`))uX;cq{EUk4UwHch3 z8pV8+&R_CA(%d}XUL9#(sFbRS4FYyCgk^vjvaRrn<^5TMo`cs#-65I%CU1*m)cD_Z z7Ivm!?ko?=TPHXUN=K(toAAGI^LE2`T=(s^?J0f1UXf?%#piYs1R+`b(HJR(ZgrYf zmMK?7pu1jfapx#cu>CD3r^3(|L!GcRF2eBLP4Zu>W zDcPtWoiBoI-jzuO%02CIXu=5e|2|3EI#n28GuW3@fYYL_wuIi>>pZE4-5z^-d!!%9 zpgX&|On6ANaY8~JXv4Lc;tVQ(_$RlNOT4i)wLmYlaAm25z{YNvO=&xbYDx z{1q4!NBYw5EFM^KzA8kT$J#Uzld7$n(rw9Z20J8>clGj|R~)jYlbmiF$gwpgG=Yl| zDDLi!OateDQ5}Uw{;<@wd`47PuWhEzqZRpviVHAO9uu<5>DKNlcamYd?N0SH2wI;} zBQeO`(m`!3Br)hhMeKasSGh|G^g9EiF5!%fn-GmLc06>mD;Cl)Fa<=oL#nnBMK@8A zx6P_|FUx^Ir{IQ$&(R-jR{>lQ{%&tG)HDYn6(+p=4cD=wT?>>1B36SunGgl{4xrX{ zzvGSbe5bA(hISmlb9oEA`{YVWlDzQ;*edP;pD)9K<%siS?QlyUApNT zG*A3b7-0W<;ZY@|5ZQ*Ux?xxz^S1BLUBr;z$FYdYs1?dezzm9c^M@-HHggXoN(d$6 zx+KbMonttZVOmQny(J{U9yVt0?GUtM`IhU63XX8R8hV9vQYCFe}LB6-OX$_cBtns%yrhtm}hdh)v3F<0v;H>`WKZKYjmCN6`vtM3jFKI{z~~XGBfiGRR_`$AnWIa7 z*I8oZJaI%0m#~I;c2CwNC0}-oVU0=pU87Xs*+A{DI145sQi=u&9j45J>}?XmT7g< zpje;tHqZJmpNMmPGMM8|rXHuU9L39Rv_Q&F+A|;P_HM#<8jU&dF2&%mP|toVGm(UQ z;gWx$Z{qXo@%uEX`^@!^=pDZiifC|*f#_RBEDwAbwVlfIa*r#h6AD!#5^kJK`F>s2 z?~UBvQK-R){$ogls)%SzP=*)3gA&yhsYmQR%*|Yk?{iX8;khY3+h?X*=k`{U`atLh z$O3w*GW@}msAUNWtcEvz*`MmDxwW=$<9Fi5hJgvX{!tI3RQer+YmYXb3JRK*`nw-= zKKgw>a?Nl)bo^}Kz@A-h#e08xr2z+sFWiJXto$eRit+71Y zH#hAb6qkiH+rQt5xO2i;ncz@ZDT90gCz-8Jo?Wp!tN&JfQBWU!r;nta!l#eE&Oyw! zx@6PWS^nae`x_(5UZH%yk>Bi(YndX~=7E)=faJ8DljVM$Rf2*&kzdNJ{$5kun<<|1 z_(%+EWz%uw(k%VrnEep+5vsD>;(|LO-5Cy--H>3JTFAP+-2l$h89<;EsRne%f1s0O zTMmZDFQ^$33PV9g55*twi6M2_5SgcL<}I zi)X9AdFbwyFZwgRE9BVc)IlzxmLJab8pH_PAssbEn!8jnddOl^NddHq(<#^%Mq%~ ztL_E1cY=&yS?U3(`Ncv7P%oxU};ZU3+g5Bmb2Bmxc9&`d*r9HmZGD zi$Oe)n+}t5u5=hh{lsvYI;MTfrX<8fGMsyq^U&*;>$hz!vRdUna#>81cA0Ga>yqpa zwh4YqT}Y`lXL0i#%$S#CtFC*?gI;~z^5`m?V8&FLNZowr)V^~41bVrLajY;Aa+PPV zrV%WV&+IL(RUZV!mE@Y~#8Up_XyLYvNpfIhP*Dr3O<~d%v_5r+1u_Yf>hbvP-^ogo zrCo=AU;Nv>UfKPV9Mz2{fQfr6-r{2kWorAAkV+N%wBlgzHqwjDT9_VUd+rfdRLJ9p zG8_Lz%yzhSsvMKcz)Ph*r)Yl`zOw!k<>dX&EseP+P}*V5z||EhR|XGr2on_dlbiil z9vcVY?*GVRN2^}bY#E~txOIVZ%~7<lcRd5Z}P0x7!@sP|7AO2<0S z+%2rU8at5G+%)HK2ZCrSU>E?3nW1;aDo1sObgVrc`hX+J-oHO;fw57|O-NducQrN+ zP(jU~Kj{Nwqy0eHeA*kL4WYmzoYB11`rWbIZ26ysg}wJ1cBooI%x{d+bjoM41os5! zY1N{Z;)!a;;u8OT6L(Ks-efIXUyYBXEL;y0+*f&@j4wM8-rGt~D9uulSwHEp*iK$C zh{$^uT0+BxiEHG)t;5F~R2HAisFAC+C%+f6!ouqm9GYg}Z7>=iKgiJH$^`6|HInfN z(aM*{5}Mv8Hw?Qv%jCJI+HR^%9lTc1_#`P7*!Ko*DR^6Qs_)0Mn76;{G5P44QbEv< zyr51X&&y zE!AwUbre(+`yLBy!76>`>BPKZ-Q0cAGmY%-4?u;sa<-P_mwx_^g-f;Ab>NX79nFb+%0tZ$rCSsEt0oReJsaO?6R~>6EQ8qCIRGZS}Gw$W2GY_ao{# z+^%4y#y|uvb=ofZL;PVzzv~#{O8v~>x>Pp5WVO?~$bY2YzMJgvt;x%?VD+&=B3^u!ic9^p~F%pKABg#qmqg# z2Q+%D&-{X?Zg`I!y6h&fQ;fCukA`gv-9I)lM7Ucd;66n*=}u)otDNS8<3g=4UzwqLQ1g0EO(H_V`Ig8|BPCBi;+}1LomdQ_R5UJ) zT@r%tU{ANUyLB|-&!hGexMO$tWRvuqZO~=YmxB||9clCsE#Z)#kItYM$K7F-AE*Tf zCasAi%N>_19R$Iu=C@LW(akW#{N_to4`-YOI_Vj#> zz-68ONw|uozq)u@XHzncyeEuUgEsEOXa|Lm z>J4;r%go6gQ&MQ!RNx@g);>x8rV^loh8c&OzmnCzd1>Fzsi$)8;OcC2&L>#)J?RiyJa#)CQwkbui{4sZoC?{1~^pLlRzZzxA&lCxkZ~}4pA>W-< z+ok5De{2+TFBq!8fXFV1;SLITaP?o#ACWEAgVo4`82OuiJA$iJR;0(@$HpcTL$ETtb)NvPal`FTzVyn-QZGUWBWHw;JAc@B6AD#xn%d2gYkVGJ(dfKA>s zQlg10mikD~*RHoXf)5yA9WB2tT9f$(UU(F#>RN%b(1_gm<8&Q7Q&Z9w$ivhqyK3&p z-CO5gU%^ecj;OK3+=OSw#Wr_oKwps0zC^dUX_}Fyy*l|ur24|?)-wwj`Q3`j_OPfD zlo0bsN?zrSCXpA1Wkp86N-nhK3=fv2kqw7u?c13SnKG*+F>0nHfw*oNHR&Ea%2v8} z_&YweaO(KH*`1G4B>6J`!ONPu)@DE)xcOIIP4j2uw4YDD5A2nuzlw#pX4(F0dG7Ja zAi4lE!%@X&UWuK?8sYTKy9x+K&0EjXAKTCbbe8IHjly#?K2Wv=WN{hHCzIE{kiCF8 zjHHB?-}v}e`1QYh$j2umF8N@eTU&VXfE$x+8P^5?{9q45oqfzP=jO!lNHGsG#n;aC zKBvFeM^r61cy=T4B=Rh1Zl**yv z`x!pu6;sW*HyoeJ4qT#FPh3zP9AXT9E|dDyzK~{=yt<7iramcF@)}o`YL5R6?%=f* zr~6l{>Pxs!fSJfzV=h8QT=#jBG?k>+=u6g{TTxWE9YJq8ukjrpKEvR^$T(^DWhEqB zi2=2S)_9};)rPUCQd8_^kMBrCt_U+WaONQZaX{$dqL-h7A4}4-6JFssP!fqt%Wp|V zX;}Y_5&@e{?JuHLO2^r)Z%F*jOkiW00-!a2)u5ooI+|O`S15(Wu6}ZU_=`JqB#-_b z0uLdSRBj~3PKjE-p54d60rl9W$_JY2vM?*;27#;{{O&e~0#M2CZFa6UvVVSt4(GRx z(TM_-m-5lJb~_;;0{}3(gDtwK55;~R�N}N;7hRP95Lr5#&=j8q0maD~R*!6ZsEvs6tL@KDj*zl8c0_R>p>c7*_6 zJDABF7&5``V9HLO91XtY>~afu-wr-2ubAqWzi>>xY#>T^1t%7z+vIIF=ZdNczSKwj z8t44&&v?@3YCeK3#s|zN|L97; zVJ0vqGaRD3w>d@lb6#`LMIXTE-v3A|ehHe!G5rQY=f{LmEcdj;mOC{iM~ZI)8n(Fr z_+8*Z{~nr2b>IP*;b3oGbECBIjUm`kYe4Q(o**yXn)%{dI`5ffVo=-!|A#q1egCy< zoN+_05Zl7y_$j!$2dFw4a*f&Uf=U~m00@m)Zr-C8SVf@ZUs-7!Cie1js807RKAEnj z0R@SySBhprHeV;m?G&@+mH4zIU$9_bUp9JjEfRPrTJdLq(gB=hJIW|M#TzlN3kUHW zh)wa#1GU=Tv!_TMcF0M6zviT1;1bI{ngU}soa`^}tU&&fI}cW;@s)iMDW+v5&PhczzQTt_Yr*QibuXiaisN18vH<$MOS+S}vxA+hi?xG24i7u~uJ7#SbbuoxsVGq{ZWQ>x E0RLE5wEzGB literal 0 HcmV?d00001 diff --git a/paper/paper.bib b/paper/paper.bib index 3dfd196..1014c0e 100644 --- a/paper/paper.bib +++ b/paper/paper.bib @@ -6,6 +6,18 @@ @Misc{shgyieldrepo url = {https://github.com/roguephysicist/SHGYield}, } +@article{andersonFMATS17, + author = {Anderson, S. M. and Mendoza, B. S.}, + title = {Depth-Dependent Three-Layer Model for the Surface Second-Harmonic Generation Yield}, + journal = {Frontiers in Materials}, + volume = {4}, + pages = {12}, + year = {2017}, + url = {http://journal.frontiersin.org/article/10.3389/fmats.2017.00012}, + issn = {2296-8016}, + doi = {10.3389/fmats.2017.00012} +} + @article{andersonPRB16b, doi = {10.1103/PhysRevB.94.115314}, year = {2016}, diff --git a/sample-data/SiH1x1-chi2-yyx b/sample-data/SiH1x1-chi2-yxy similarity index 100% rename from sample-data/SiH1x1-chi2-yyx rename to sample-data/SiH1x1-chi2-yxy diff --git a/sample-data/SiH1x1-chi2-zzy b/sample-data/SiH1x1-chi2-zyz similarity index 100% rename from sample-data/SiH1x1-chi2-zzy rename to sample-data/SiH1x1-chi2-zyz diff --git a/sample.in b/sample.in deleted file mode 100644 index 8a77581..0000000 --- a/sample.in +++ /dev/null @@ -1,43 +0,0 @@ -# Sample input file for calculating SHG yield - -# Establish reflectance model, see PRB 93, 235304 (2016). -# This is mostly deprecated now, so when in doubt just use '3-layer'. -mode: 3-layer - -# Required calculation parameters -theta: 65 # angle of incidence in degrees -phi: 30 # azimuthal angle in degrees -norm: 1.2659296143 # chi1 normalization factor, see PRB 92, 245308 (2015) -sigma: 0.10 # gaussian output broadening, sigma = 0.10 eV - -# Multiple reflection framework -multiref: yes # if multiple reflections are enabled ('yes' or 'no') -thickness: 10 # thickness "d" of the thin layer, in nm -depth: average # depth "d2" of the nonlinear polarization, in nm. can - # be a float or 'average'. See Eqs (21) and 22 of - # PRB 93, 235304 (2016). - -## chi1 & chi2 components. Can be path to files, or number (like 0). -chib: sample-data/SiBulk-chi1-xx_yy_zz -chil: sample-data/SiH1x1-chi1-xx_yy_zz -xxx: sample-data/SiH1x1-chi2-xxx -xxy: sample-data/SiH1x1-chi2-xxy -xxz: sample-data/SiH1x1-chi2-xxz -xyy: sample-data/SiH1x1-chi2-xyy -xyz: sample-data/SiH1x1-chi2-xyz -xzz: sample-data/SiH1x1-chi2-xzz -yxx: sample-data/SiH1x1-chi2-yxx -yxz: sample-data/SiH1x1-chi2-yxz -yyx: sample-data/SiH1x1-chi2-yyx -yyy: sample-data/SiH1x1-chi2-yyy -yyz: sample-data/SiH1x1-chi2-yyz -yzz: sample-data/SiH1x1-chi2-yzz -zxx: sample-data/SiH1x1-chi2-zxx -zxy: sample-data/SiH1x1-chi2-zxy -zxz: sample-data/SiH1x1-chi2-zxz -zyy: sample-data/SiH1x1-chi2-zyy -zzy: sample-data/SiH1x1-chi2-zzy -zzz: sample-data/SiH1x1-chi2-zzz - -# Path to desired output file -output: sample.out diff --git a/shgyield.py b/shgyield.py index e8f13da..e59c77d 100755 --- a/shgyield.py +++ b/shgyield.py @@ -16,34 +16,25 @@ """ import sys -import math +import yaml import numpy as np from scipy import constants, ndimage def parse_input(infile): ''' - Parses the input file specified on the command line. Returns dictionary + Parses the YAML input file specified on the command line. Returns dictionary with variables and values from the input file. ''' - params = {} - targetfile = open(infile) - data = targetfile.readlines() - for line in data: - if line.strip().startswith('#') or line == '\n': - pass - else: - key, value = line.partition('#')[0].split(':') - params[key.strip()] = value.strip() - targetfile.close() + with open(infile) as data_file: + params = yaml.load(data_file) return params def broad(target, sigma): ''' A function for applying Gaussian broadening on the final output data. ''' - value = sigma * 42.666666666 - data = ndimage.filters.gaussian_filter(target, value) + data = ndimage.filters.gaussian_filter(target, sigma) return data def savefile(file, freq, val1, val2, val3, val4): @@ -52,10 +43,10 @@ def savefile(file, freq, val1, val2, val3, val4): Energy(1w) R_{pP} R_{pS} R_{sP} R_{sS} ''' data = np.column_stack((freq, val1, val2, val3, val4)) - np.savetxt(file, data, fmt=('%05.2f', '%.14e', '%.14e', '%.14e', '%.14e'), + np.savetxt(file, data, fmt=('%05.2f', '%.8e', '%.8e', '%.8e', '%.8e'), delimiter=' ', - header='RiF 1e-20 (cm^2/W)\n1w(eV) RpP'+21*" "+\ - 'RpS'+21*" "+'RsP'+21*" "+'RsS') + header='RiF 1e-20 (cm^2/W)\n1w(eV) RpP'+15*" "+\ + 'RpS'+15*" "+'RsP'+15*" "+'RsS') def epsilon(in_file, norm): ''' @@ -73,31 +64,15 @@ def epsilon(in_file, norm): epsa = 1 + (4 * constants.pi * coma * norm) # epsilon with normalization return epsa -def shgcomp(): +def shgload(infile): ''' - Reads each chi^{abc} listed in the input file, with the following columns: - Energy(1w) Re[1w] Im[1w] Re[2w] Im[2w]. - Checks for all 18 components possible for SHG; those not listed in the input - file are assumed to be 0. Sums 1w and 2w real and imaginary parts, and - returns a dictionary with the component name as the 'key', and a numpy array - with the component data as the 'value'. + Loads chi^{abc} listed in the input file, with the following columns: + Energy(1w) Re[1w] Im[1w] Re[2w] Im[2w]. Sums 1w and 2w real and imaginary + parts, and handles scale and units appropriately. Returns numpy array. ''' - chi2 = {} - components = ['xxx', 'xxy', 'xxz', 'xyy', 'xyz', 'xzz', - 'yxx', 'yyx', 'yxz', 'yyy', 'yyz', 'yzz', - 'zxx', 'zxy', 'zxz', 'zyy', 'zzy', 'zzz'] - for response in components: - if response in PARAM: - try: - factor = float(PARAM[response]) - shg = TINIBASCALE * PM2TOM2 * factor - except ValueError: - data = np.loadtxt(PARAM[response], unpack=True, skiprows=1) - comp = (data[1] + data[3]) + 1j * (data[2] + data[4]) - shg = TINIBASCALE * PM2TOM2 * comp[:MAXE] - elif response not in PARAM: - shg = 0 - chi2[response] = shg + data = np.loadtxt(infile, unpack=True, skiprows=1) + comp = (data[1] + data[3]) + 1j * (data[2] + data[4]) + chi2 = TINIBASCALE * PM2TOM2 * comp[:MAXE] return chi2 def fresnel(kind, i, j, pol, freq): @@ -133,16 +108,16 @@ def shgyield(gamma, riF): ## Initialization: Parse input file, establish relevant modes. PARAM = parse_input(sys.argv[1]) # Parses input file -MODE = str(PARAM['mode']) # Establishes the 'layer model' to be used; - # see PRB 93, 235304 (2016). Mostly deprecated. -MULTIREF = str(PARAM['multiref']) # Whether or not multiple reflections are - # considered. +MODE = PARAM['mode'] # Establishes the 'layer model' to be used; + # see PRB 93, 235304 (2016). + ## Energy: establishes energy range ## Assumes chi1/2 components are taken over 0-20 eV with 2001 (0.01 eV) steps MAXE = 1000 # Upper energy limit, usually 1000 = 10 eV ONEE = np.linspace(0.01, float(MAXE)/100, MAXE) # 1w energy array, 0.01-10 eV + ## Constants and parameters HBAR = constants.value("Planck constant over 2 pi in eV s") PLANCK = constants.value("Planck constant in eV s") @@ -152,19 +127,17 @@ def shgyield(gamma, riF): M2TOCM2 = 1e4 # Convert from m^2 to cm^2 TINIBASCALE = 1e6 # Scaling chi2 in 1e6 (pm^2/V) SCALE = 1e20 # Final yield in 1e-20 (cm^2/W) -THETA0 = math.radians(float(PARAM['theta'])) # Converts theta to radians -PHI = math.radians(float(PARAM['phi'])) # Converts phi to radians -THICKNESS = float(PARAM['thickness']) # Thickness d of the thin layer \ell -DEPTH = PARAM['depth'] # Depth d2 of the polarization sheet -SIGMA = float(PARAM['sigma']) # Std. dev. for gaussian broadening +THETA0 = np.radians(PARAM['parameters']['theta']) # Converts theta to radians +PHI = np.radians(PARAM['parameters']['phi']) # Converts phi to radians +SIGMA = PARAM['parameters']['sigma'] * 128/3 # Std. dev. for gaussian broadening LAMBDA0 = (PLANCK * LSPEED * 1e9)/ONEE # Energy range expressed in nm -CHI1NORM = float(PARAM['norm']) # DEBUG: normalization for chi1 -PREFACTOR = 1/(2 * EPS0 * HBAR**2 * LSPEED**3 * math.cos(THETA0)**2) +PREFACTOR = 1/(2 * EPS0 * HBAR**2 * LSPEED**3 * np.cos(THETA0)**2) -## Components and pre-calculation -CHI2 = shgcomp() # Dictionary with all chi2 components -epsl = epsilon(PARAM['chil'], CHI1NORM) # Epsilon from chi1, layered, normalized -epsb = epsilon(PARAM['chib'], 1) # Epsilon from chi1, bulk + +## Linear responses: chi1 and epsilons +CHI1NORM = PARAM['chi1']['norm'] # Normalization for layered chi1 +epsl = epsilon(PARAM['chi1']['chil'], CHI1NORM) # Epsilon from chi1, layered, normalized +epsb = epsilon(PARAM['chi1']['chib'], 1) # Epsilon from chi1, bulk epsv1w = 1 # Epsilon for vacuum = 1 epsv2w = 1 # Epsilon for vacuum = 1 epsb1w = epsb[:MAXE] # Epsilon for bulk, 1w @@ -172,7 +145,8 @@ def shgyield(gamma, riF): epsl1w = epsl[:MAXE] # Epsilon for layer, 1w epsl2w = epsl[1::2][:MAXE] # Epsilon for layer, 2w -## Reflection model, see PRB 93, 235304 (2016). This is mostly deprecated. + +## Reflection model, see PRB 93, 235304 (2016). if MODE == "3-layer": ell1w = "l" ell2w = "l" @@ -191,6 +165,44 @@ def shgyield(gamma, riF): epsl1w = eval("eps" + ell1w + "1w") epsl2w = eval("eps" + ell2w + "2w") + +## Nonlinear responses: chi2 components to be used in the formulas below. +## There are 18 possible components for SHG. We have 4 scenarios regarding the +## input file and how it is interpreted: +## 1. Component not listed (or commented) (chi2^{abc} = 0) +## 2. Component specified as 0 (chi2^{abc} = 0) +## 3. Component specified with path to file (chi2^{abc} loaded from file) +## 4. Component has symmetry relation, such as 'yyx: -xxx' (chi2^{yyx} = -chi2^{xxx}) +## +## The end result is a dictionary (CHI2) with the component name as the 'key', +## and a numpy array with the component data as the 'value'. +CHI2 = {} +CHI2_equiv = {} +all_components = ['xxx', 'xyy', 'xzz', 'xyz', 'xxz', 'xxy', + 'yxx', 'yyy', 'yzz', 'yyz', 'yxz', 'yxy', + 'zxx', 'zyy', 'zzz', 'zyz', 'zxz', 'zxy'] +for component in all_components: + if component in PARAM['chi2']: + value = PARAM['chi2'][component] + try: + shg = shgload(value) + except (ValueError, FileNotFoundError, OSError, IOError): + if value == 0: + shg = 0 + else: + CHI2_equiv[component] = value + continue + elif component not in PARAM['chi2']: + shg = 0 + CHI2[component] = shg +for key, value in CHI2_equiv.items(): + equivalence = value.split('-') + if equivalence[0]: + CHI2[key] = CHI2[equivalence[0]] + elif not equivalence[0]: + CHI2[key] = -1 * CHI2[equivalence[1]] + + ## Indices of refraction for 1w (n) and 2w (N), where n = sqrt{epsilon} nv = np.sqrt(epsv1w) # Index of refraction for vacuum, 1w Nv = np.sqrt(epsv2w) # Index of refraction for vacuum, 2w @@ -199,13 +211,15 @@ def shgyield(gamma, riF): nl = np.sqrt(epsl1w) # Index of refraction for layer, 1w Nl = np.sqrt(epsl2w) # Index of refraction for layer, 2w + ## Wave vectors for 1w and 2w, see Eq. (6) of PRB 94, 115314 (2016). -wv1w = np.sqrt(epsv1w - (math.sin(THETA0)**2)) # Wave vector for vacuum, 1w -wv2w = np.sqrt(epsv2w - (math.sin(THETA0)**2)) # Wave vector for vacuum, 2w -wb1w = np.sqrt(epsb1w - (math.sin(THETA0)**2)) # Wave vector for bulk, 1w -wb2w = np.sqrt(epsb2w - (math.sin(THETA0)**2)) # Wave vector for bulk, 2w -wl1w = np.sqrt(epsl1w - (math.sin(THETA0)**2)) # Wave vector for layer, 1w -wl2w = np.sqrt(epsl2w - (math.sin(THETA0)**2)) # Wave vector for layer, 2w +wv1w = np.sqrt(epsv1w - (np.sin(THETA0)**2)) # Wave vector for vacuum, 1w +wv2w = np.sqrt(epsv2w - (np.sin(THETA0)**2)) # Wave vector for vacuum, 2w +wb1w = np.sqrt(epsb1w - (np.sin(THETA0)**2)) # Wave vector for bulk, 1w +wb2w = np.sqrt(epsb2w - (np.sin(THETA0)**2)) # Wave vector for bulk, 2w +wl1w = np.sqrt(epsl1w - (np.sin(THETA0)**2)) # Wave vector for layer, 1w +wl2w = np.sqrt(epsl2w - (np.sin(THETA0)**2)) # Wave vector for layer, 2w + ## Fresnel factors for 1w and 2w, s and p polarizations. See Eqs. (13) and (14) ## of PRB 94, 115314 (2016). @@ -226,12 +240,15 @@ def shgyield(gamma, riF): Rlbs = fresnel("r", ell2w, "b", "s", "2w") # Reflection, 2w, layer-bulk, s Rlbp = fresnel("r", ell2w, "b", "p", "2w") # Reflection, 2w, layer-bulk, p + ## Multiple reflections framework. See Eqs. (16), (17), (21), (22), (26), ## and (30) of PRB 94, 115314 (2016). -if MULTIREF == "yes": - varphi = 4 * math.pi * ((ONEE * THICKNESS * 1e-9)/(PLANCK * LSPEED)) * wl1w - delta = 8 * math.pi * ((ONEE * THICKNESS * 1e-9)/(PLANCK * LSPEED)) * wl2w - if DEPTH == "average": +if PARAM['multiref']['enable']: + thickness = PARAM['multiref']['thickness'] # Thickness d of the thin layer \ell + depth = PARAM['multiref']['depth'] # Depth d2 of the polarization sheet + varphi = 4 * np.pi * ((ONEE * thickness * 1e-9)/(PLANCK * LSPEED)) * wl1w + delta = 8 * np.pi * ((ONEE * thickness * 1e-9)/(PLANCK * LSPEED)) * wl2w + if depth == "average": RMpav = (Rlbp * np.exp(1j * delta/2))/\ (1 + (Rvlp * Rlbp * np.exp(1j * delta))) * np.sin(delta/2)/(delta/2) RMsav = (Rlbs * np.exp(1j * delta/2))/\ @@ -241,8 +258,8 @@ def shgyield(gamma, riF): RMminusp = 1 - RMpav RMminuss = 1 - RMsav else: - D2 = float(DEPTH) - delta0 = 8 * math.pi * ((ONEE * D2 * 1e-9)/(PLANCK * LSPEED)) * wl2w + D2 = float(depth) + delta0 = 8 * np.pi * ((ONEE * D2 * 1e-9)/(PLANCK * LSPEED)) * wl2w RMp = (Rlbp * np.exp(1j * delta0))/(1 + (Rvlp * Rlbp * np.exp(1j * delta))) RMs = (Rlbs * np.exp(1j * delta0))/(1 + (Rvls * Rlbs * np.exp(1j * delta))) RMplusp = 1 + RMp @@ -255,7 +272,7 @@ def shgyield(gamma, riF): rMpluss = 1 + rMs rMminusp = 1 - rMp rMminuss = 1 - rMs -elif MULTIREF == "no": +elif not PARAM['multiref']['enable']: RMplusp = 1 + Rlbp RMpluss = 1 + Rlbs RMminusp = 1 - Rlbp @@ -265,72 +282,75 @@ def shgyield(gamma, riF): rMminusp = 1 - rlbp rMminuss = 1 - rlbs + ## r factors for different input and output polarizations. See Eqs. (50), (55), ## (60), and (65) of PRB 94, 115314 (2016). ### r_{pP} -rMRpP = - (RMminusp * rMminusp**2 * wl1w**2 * wl2w * math.cos(PHI)**3 * CHI2['xxx']) \ - - (2 * RMminusp * rMminusp**2 * wl1w**2 * wl2w * math.sin(PHI) * math.cos(PHI)**2 * CHI2['xxy']) \ - - (2 * RMminusp * rMplusp * rMminusp * wl1w * wl2w * math.sin(THETA0) * math.cos(PHI)**2 * CHI2['xxz']) \ - - (RMminusp * rMminusp**2 * wl1w**2 * wl2w * math.sin(PHI)**2 * math.cos(PHI) * CHI2['xyy']) \ - - (2 * RMminusp * rMplusp * rMminusp * wl1w * wl2w * math.sin(THETA0) * math.sin(PHI) * math.cos(PHI) * CHI2['xyz']) \ - - (RMminusp * rMplusp**2 * wl2w * math.sin(THETA0)**2 * math.cos(PHI) * CHI2['xzz']) \ - - (RMminusp * rMminusp**2 * wl1w**2 * wl2w * math.sin(PHI) * math.cos(PHI)**2 * CHI2['yxx']) \ - - (2 * RMminusp * rMminusp**2 * wl1w**2 * wl2w * math.sin(PHI)**2 * math.cos(PHI) * CHI2['yyx']) \ - - (2 * RMminusp * rMplusp * rMminusp * wl1w * wl2w * math.sin(THETA0) * math.sin(PHI) * math.cos(PHI) * CHI2['yxz']) \ - - (RMminusp * rMminusp**2 * wl1w**2 * wl2w * math.sin(PHI)**3 * CHI2['yyy']) \ - - (2 * RMminusp * rMplusp * rMminusp * wl1w * wl2w * math.sin(THETA0) * math.sin(PHI)**2 * CHI2['yyz']) \ - - (RMminusp * rMplusp**2 * wl2w * math.sin(THETA0)**2 * math.sin(PHI) * CHI2['yzz']) \ - + (RMplusp * rMminusp**2 * wl1w**2 * math.sin(THETA0) * math.cos(PHI)**2 * CHI2['zxx']) \ - + (2 * RMplusp * rMplusp * rMminusp * wl1w * math.sin(THETA0)**2 * math.cos(PHI) * CHI2['zxz']) \ - + (2 * RMplusp * rMminusp**2 * wl1w**2 * math.sin(THETA0) * math.sin(PHI) * math.cos(PHI) * CHI2['zxy']) \ - + (RMplusp * rMminusp**2 * wl1w**2 * math.sin(THETA0) * math.sin(PHI)**2 * CHI2['zyy']) \ - + (2 * RMplusp * rMplusp * rMminusp * wl1w * math.sin(THETA0)**2 * math.sin(PHI) * CHI2['zzy']) \ - + (RMplusp * rMplusp**2 * math.sin(PHI)**3 * CHI2['zzz']) +rpP = - (RMminusp * rMminusp**2 * wl1w**2 * wl2w * np.cos(PHI)**3 * CHI2['xxx']) \ + - (2 * RMminusp * rMminusp**2 * wl1w**2 * wl2w * np.sin(PHI) * np.cos(PHI)**2 * CHI2['xxy']) \ + - (2 * RMminusp * rMplusp * rMminusp * wl1w * wl2w * np.sin(THETA0) * np.cos(PHI)**2 * CHI2['xxz']) \ + - (RMminusp * rMminusp**2 * wl1w**2 * wl2w * np.sin(PHI)**2 * np.cos(PHI) * CHI2['xyy']) \ + - (2 * RMminusp * rMplusp * rMminusp * wl1w * wl2w * np.sin(THETA0) * np.sin(PHI) * np.cos(PHI) * CHI2['xyz']) \ + - (RMminusp * rMplusp**2 * wl2w * np.sin(THETA0)**2 * np.cos(PHI) * CHI2['xzz']) \ + - (RMminusp * rMminusp**2 * wl1w**2 * wl2w * np.sin(PHI) * np.cos(PHI)**2 * CHI2['yxx']) \ + - (2 * RMminusp * rMminusp**2 * wl1w**2 * wl2w * np.sin(PHI)**2 * np.cos(PHI) * CHI2['yxy']) \ + - (2 * RMminusp * rMplusp * rMminusp * wl1w * wl2w * np.sin(THETA0) * np.sin(PHI) * np.cos(PHI) * CHI2['yxz']) \ + - (RMminusp * rMminusp**2 * wl1w**2 * wl2w * np.sin(PHI)**3 * CHI2['yyy']) \ + - (2 * RMminusp * rMplusp * rMminusp * wl1w * wl2w * np.sin(THETA0) * np.sin(PHI)**2 * CHI2['yyz']) \ + - (RMminusp * rMplusp**2 * wl2w * np.sin(THETA0)**2 * np.sin(PHI) * CHI2['yzz']) \ + + (RMplusp * rMminusp**2 * wl1w**2 * np.sin(THETA0) * np.cos(PHI)**2 * CHI2['zxx']) \ + + (2 * RMplusp * rMplusp * rMminusp * wl1w * np.sin(THETA0)**2 * np.cos(PHI) * CHI2['zxz']) \ + + (2 * RMplusp * rMminusp**2 * wl1w**2 * np.sin(THETA0) * np.sin(PHI) * np.cos(PHI) * CHI2['zxy']) \ + + (RMplusp * rMminusp**2 * wl1w**2 * np.sin(THETA0) * np.sin(PHI)**2 * CHI2['zyy']) \ + + (2 * RMplusp * rMplusp * rMminusp * wl1w * np.sin(THETA0)**2 * np.sin(PHI) * CHI2['zyz']) \ + + (RMplusp * rMplusp**2 * np.sin(PHI)**3 * CHI2['zzz']) ### r_{pS} -rMRpS = - (rMminusp**2 * wl1w**2 * math.sin(PHI) * math.cos(PHI)**2 * CHI2['xxx']) \ - - (2 * rMminusp**2 * wl1w**2 * math.sin(PHI)**2 * math.cos(PHI) * CHI2['xxy']) \ - - (2 * rMplusp * rMminusp * wl1w * math.sin(THETA0) * math.sin(PHI) * math.cos(PHI) * CHI2['xxz']) \ - - (rMminusp**2 * wl1w**2 * math.sin(PHI)**3 * CHI2['xyy']) \ - - (2 * rMplusp * rMminusp * wl1w * math.sin(THETA0) * math.sin(PHI)**2 * CHI2['xyz']) \ - - (rMplusp**2 * math.sin(THETA0)**2 * math.sin(PHI) * CHI2['xzz']) \ - + (rMminusp**2 * wl1w**2 * math.cos(PHI)**3 * CHI2['yxx']) \ - + (2 * rMminusp**2 * wl1w**2 * math.sin(PHI) * math.cos(PHI)**2 * CHI2['yyx']) \ - + (2 * rMplusp * rMminusp * wl1w * math.sin(THETA0) * math.cos(PHI)**2 * CHI2['yxz']) \ - + (rMminusp**2 * wl1w**2 * math.sin(PHI)**2 * math.cos(PHI) * CHI2['yyy']) \ - + (2 * rMplusp * rMminusp * wl1w * math.sin(THETA0) * math.sin(PHI) * math.cos(PHI) * CHI2['yyz']) \ - + (rMplusp**2 * math.sin(THETA0)**2 * math.cos(PHI) * CHI2['yzz']) +rpS = - (rMminusp**2 * wl1w**2 * np.sin(PHI) * np.cos(PHI)**2 * CHI2['xxx']) \ + - (2 * rMminusp**2 * wl1w**2 * np.sin(PHI)**2 * np.cos(PHI) * CHI2['xxy']) \ + - (2 * rMplusp * rMminusp * wl1w * np.sin(THETA0) * np.sin(PHI) * np.cos(PHI) * CHI2['xxz']) \ + - (rMminusp**2 * wl1w**2 * np.sin(PHI)**3 * CHI2['xyy']) \ + - (2 * rMplusp * rMminusp * wl1w * np.sin(THETA0) * np.sin(PHI)**2 * CHI2['xyz']) \ + - (rMplusp**2 * np.sin(THETA0)**2 * np.sin(PHI) * CHI2['xzz']) \ + + (rMminusp**2 * wl1w**2 * np.cos(PHI)**3 * CHI2['yxx']) \ + + (2 * rMminusp**2 * wl1w**2 * np.sin(PHI) * np.cos(PHI)**2 * CHI2['yxy']) \ + + (2 * rMplusp * rMminusp * wl1w * np.sin(THETA0) * np.cos(PHI)**2 * CHI2['yxz']) \ + + (rMminusp**2 * wl1w**2 * np.sin(PHI)**2 * np.cos(PHI) * CHI2['yyy']) \ + + (2 * rMplusp * rMminusp * wl1w * np.sin(THETA0) * np.sin(PHI) * np.cos(PHI) * CHI2['yyz']) \ + + (rMplusp**2 * np.sin(THETA0)**2 * np.cos(PHI) * CHI2['yzz']) ### r_{sP} -rMRsP = - (RMminusp * wl2w * math.sin(PHI)**2 * math.cos(PHI) * CHI2['xxx']) \ - + (RMminusp * wl2w * 2 * math.sin(PHI) * math.cos(PHI)**2 * CHI2['xxy']) \ - - (RMminusp * wl2w * math.cos(PHI)**3 * CHI2['xyy']) \ - - (RMminusp * wl2w * math.sin(PHI)**3 * CHI2['yxx']) \ - + (RMminusp * wl2w * 2 * math.sin(PHI)**2 * math.cos(PHI) * CHI2['yyx']) \ - - (RMminusp * wl2w * math.sin(PHI) * math.cos(PHI)**2 * CHI2['yyy']) \ - + (RMplusp * math.sin(THETA0) * math.sin(PHI)**2 * CHI2['zxx']) \ - - (RMplusp * math.sin(THETA0) * 2 * math.sin(PHI) * math.cos(PHI) * CHI2['zxy']) \ - + (RMplusp * math.sin(THETA0) * math.cos(PHI)**2 * CHI2['zyy']) +rsP = - (RMminusp * wl2w * np.sin(PHI)**2 * np.cos(PHI) * CHI2['xxx']) \ + + (RMminusp * wl2w * 2 * np.sin(PHI) * np.cos(PHI)**2 * CHI2['xxy']) \ + - (RMminusp * wl2w * np.cos(PHI)**3 * CHI2['xyy']) \ + - (RMminusp * wl2w * np.sin(PHI)**3 * CHI2['yxx']) \ + + (RMminusp * wl2w * 2 * np.sin(PHI)**2 * np.cos(PHI) * CHI2['yxy']) \ + - (RMminusp * wl2w * np.sin(PHI) * np.cos(PHI)**2 * CHI2['yyy']) \ + + (RMplusp * np.sin(THETA0) * np.sin(PHI)**2 * CHI2['zxx']) \ + - (RMplusp * np.sin(THETA0) * 2 * np.sin(PHI) * np.cos(PHI) * CHI2['zxy']) \ + + (RMplusp * np.sin(THETA0) * np.cos(PHI)**2 * CHI2['zyy']) ### r_{sS} -rMRsS = - (math.sin(PHI)**3 * CHI2['xxx']) \ - + (2 * math.sin(PHI)**2 * math.cos(PHI) * CHI2['xxy']) \ - - (math.sin(PHI) * math.cos(PHI)**2 * CHI2['xyy']) \ - + (math.sin(PHI)**2 * math.cos(PHI) * CHI2['yxx']) \ - + (math.cos(PHI)**3 * CHI2['yyy']) \ - - (2 * math.sin(PHI) * math.cos(PHI)**2 * CHI2['yyx']) +rsS = - (np.sin(PHI)**3 * CHI2['xxx']) \ + + (2 * np.sin(PHI)**2 * np.cos(PHI) * CHI2['xxy']) \ + - (np.sin(PHI) * np.cos(PHI)**2 * CHI2['xyy']) \ + + (np.sin(PHI)**2 * np.cos(PHI) * CHI2['yxx']) \ + + (np.cos(PHI)**3 * CHI2['yyy']) \ + - (2 * np.sin(PHI) * np.cos(PHI)**2 * CHI2['yxy']) + ## Gamma prefactor for different polarizations. See Eqs. (49), (54), (59), (64) ## of PRB 94, 115314 (2016). -GammaMRpP = (Tvlp/Nl) * (tvlp/nl)**2 # p-in, P-out -GammaMRpS = Tvls * RMpluss * (tvlp/nl)**2 # p-in, S-out -GammaMRsP = (Tvlp/Nl) * (tvls * rMpluss)**2 # s-in, P-out -GammaMRsS = Tvls * RMpluss * (tvls * rMpluss)**2 # s-in, S-out +GammapP = (Tvlp/Nl) * (tvlp/nl)**2 # p-in, P-out +GammapS = Tvls * RMpluss * (tvlp/nl)**2 # p-in, S-out +GammasP = (Tvlp/Nl) * (tvls * rMpluss)**2 # s-in, P-out +GammasS = Tvls * RMpluss * (tvls * rMpluss)**2 # s-in, S-out + ## Final SHG yield for different input and output polarizations (in cm^2/W). ## See Eqs. (44) and (38) of PRB 94, 115314 (2016). -RMRpP = shgyield(GammaMRpP, rMRpP) # p-in, P-out -RMRpS = shgyield(GammaMRpS, rMRpS) # p-in, S-out -RMRsP = shgyield(GammaMRsP, rMRsP) # s-in, P-out -RMRsS = shgyield(GammaMRsS, rMRsS) # s-in, S-out +RpP = shgyield(GammapP, rpP) # p-in, P-out +RpS = shgyield(GammapS, rpS) # p-in, S-out +RsP = shgyield(GammasP, rsP) # s-in, P-out +RsS = shgyield(GammasS, rsS) # s-in, S-out ## Write output to file. -savefile(PARAM['output'], ONEE, RMRpP, RMRpS, RMRsP, RMRsS) +savefile(PARAM['output'], ONEE, RpP, RpS, RsP, RsS)