From 53d98ce71acddd028a8361bf966aa8d6be204cfd Mon Sep 17 00:00:00 2001 From: Marco De Pietri Date: Tue, 3 Mar 2026 10:48:56 -0500 Subject: [PATCH] Add method on Material for computing photon contact dose rate (#3700) Co-authored-by: Paul Romano --- docs/source/pythonapi/data.rst | 2 + openmc/data/__init__.py | 4 +- .../data/{effective_dose => dose}/__init__.py | 0 openmc/data/{effective_dose => dose}/dose.py | 0 .../icrp116/electrons.txt | 0 .../icrp116/helium_ions.txt | 0 .../icrp116/negative_muons.txt | 0 .../icrp116/negative_pions.txt | 0 .../icrp116/neutrons.txt | 0 .../icrp116/photons.txt | 0 .../icrp116/photons_kerma.txt | 0 .../icrp116/positive_muons.txt | 0 .../icrp116/positive_pions.txt | 0 .../icrp116/positrons.txt | 0 .../icrp116/protons.txt | 0 .../icrp74/generate_photon_effective_dose.py | 0 .../icrp74/neutrons.txt | 0 .../icrp74/photons.txt | 0 openmc/data/dose/mass_attenuation.h5 | Bin 0 -> 135056 bytes openmc/data/dose/mass_attenuation.py | 153 ++++++++++++++ openmc/material.py | 191 +++++++++++++++++- pyproject.toml | 2 +- .../unit_tests/test_data_mass_attenuation.py | 53 +++++ tests/unit_tests/test_material.py | 55 +++++ tests/unit_tests/test_mesh.py | 1 + 25 files changed, 457 insertions(+), 4 deletions(-) rename openmc/data/{effective_dose => dose}/__init__.py (100%) rename openmc/data/{effective_dose => dose}/dose.py (100%) rename openmc/data/{effective_dose => dose}/icrp116/electrons.txt (100%) rename openmc/data/{effective_dose => dose}/icrp116/helium_ions.txt (100%) rename openmc/data/{effective_dose => dose}/icrp116/negative_muons.txt (100%) rename openmc/data/{effective_dose => dose}/icrp116/negative_pions.txt (100%) rename openmc/data/{effective_dose => dose}/icrp116/neutrons.txt (100%) rename openmc/data/{effective_dose => dose}/icrp116/photons.txt (100%) rename openmc/data/{effective_dose => dose}/icrp116/photons_kerma.txt (100%) rename openmc/data/{effective_dose => dose}/icrp116/positive_muons.txt (100%) rename openmc/data/{effective_dose => dose}/icrp116/positive_pions.txt (100%) rename openmc/data/{effective_dose => dose}/icrp116/positrons.txt (100%) rename openmc/data/{effective_dose => dose}/icrp116/protons.txt (100%) rename openmc/data/{effective_dose => dose}/icrp74/generate_photon_effective_dose.py (100%) rename openmc/data/{effective_dose => dose}/icrp74/neutrons.txt (100%) rename openmc/data/{effective_dose => dose}/icrp74/photons.txt (100%) create mode 100644 openmc/data/dose/mass_attenuation.h5 create mode 100644 openmc/data/dose/mass_attenuation.py create mode 100644 tests/unit_tests/test_data_mass_attenuation.py diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 1eaf90c97..9d47430f7 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -71,6 +71,8 @@ Core Functions isotopes kalbach_slope linearize + mass_attenuation_coefficient + mass_energy_absorption_coefficient thin water_density zam diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index a45d026a0..9b38d758e 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -35,4 +35,6 @@ from .grid import * from .function import * from .vectfit import * -from .effective_dose.dose import dose_coefficients +from .dose.dose import dose_coefficients +from .dose.mass_attenuation import \ + mass_energy_absorption_coefficient, mass_attenuation_coefficient diff --git a/openmc/data/effective_dose/__init__.py b/openmc/data/dose/__init__.py similarity index 100% rename from openmc/data/effective_dose/__init__.py rename to openmc/data/dose/__init__.py diff --git a/openmc/data/effective_dose/dose.py b/openmc/data/dose/dose.py similarity index 100% rename from openmc/data/effective_dose/dose.py rename to openmc/data/dose/dose.py diff --git a/openmc/data/effective_dose/icrp116/electrons.txt b/openmc/data/dose/icrp116/electrons.txt similarity index 100% rename from openmc/data/effective_dose/icrp116/electrons.txt rename to openmc/data/dose/icrp116/electrons.txt diff --git a/openmc/data/effective_dose/icrp116/helium_ions.txt b/openmc/data/dose/icrp116/helium_ions.txt similarity index 100% rename from openmc/data/effective_dose/icrp116/helium_ions.txt rename to openmc/data/dose/icrp116/helium_ions.txt diff --git a/openmc/data/effective_dose/icrp116/negative_muons.txt b/openmc/data/dose/icrp116/negative_muons.txt similarity index 100% rename from openmc/data/effective_dose/icrp116/negative_muons.txt rename to openmc/data/dose/icrp116/negative_muons.txt diff --git a/openmc/data/effective_dose/icrp116/negative_pions.txt b/openmc/data/dose/icrp116/negative_pions.txt similarity index 100% rename from openmc/data/effective_dose/icrp116/negative_pions.txt rename to openmc/data/dose/icrp116/negative_pions.txt diff --git a/openmc/data/effective_dose/icrp116/neutrons.txt b/openmc/data/dose/icrp116/neutrons.txt similarity index 100% rename from openmc/data/effective_dose/icrp116/neutrons.txt rename to openmc/data/dose/icrp116/neutrons.txt diff --git a/openmc/data/effective_dose/icrp116/photons.txt b/openmc/data/dose/icrp116/photons.txt similarity index 100% rename from openmc/data/effective_dose/icrp116/photons.txt rename to openmc/data/dose/icrp116/photons.txt diff --git a/openmc/data/effective_dose/icrp116/photons_kerma.txt b/openmc/data/dose/icrp116/photons_kerma.txt similarity index 100% rename from openmc/data/effective_dose/icrp116/photons_kerma.txt rename to openmc/data/dose/icrp116/photons_kerma.txt diff --git a/openmc/data/effective_dose/icrp116/positive_muons.txt b/openmc/data/dose/icrp116/positive_muons.txt similarity index 100% rename from openmc/data/effective_dose/icrp116/positive_muons.txt rename to openmc/data/dose/icrp116/positive_muons.txt diff --git a/openmc/data/effective_dose/icrp116/positive_pions.txt b/openmc/data/dose/icrp116/positive_pions.txt similarity index 100% rename from openmc/data/effective_dose/icrp116/positive_pions.txt rename to openmc/data/dose/icrp116/positive_pions.txt diff --git a/openmc/data/effective_dose/icrp116/positrons.txt b/openmc/data/dose/icrp116/positrons.txt similarity index 100% rename from openmc/data/effective_dose/icrp116/positrons.txt rename to openmc/data/dose/icrp116/positrons.txt diff --git a/openmc/data/effective_dose/icrp116/protons.txt b/openmc/data/dose/icrp116/protons.txt similarity index 100% rename from openmc/data/effective_dose/icrp116/protons.txt rename to openmc/data/dose/icrp116/protons.txt diff --git a/openmc/data/effective_dose/icrp74/generate_photon_effective_dose.py b/openmc/data/dose/icrp74/generate_photon_effective_dose.py similarity index 100% rename from openmc/data/effective_dose/icrp74/generate_photon_effective_dose.py rename to openmc/data/dose/icrp74/generate_photon_effective_dose.py diff --git a/openmc/data/effective_dose/icrp74/neutrons.txt b/openmc/data/dose/icrp74/neutrons.txt similarity index 100% rename from openmc/data/effective_dose/icrp74/neutrons.txt rename to openmc/data/dose/icrp74/neutrons.txt diff --git a/openmc/data/effective_dose/icrp74/photons.txt b/openmc/data/dose/icrp74/photons.txt similarity index 100% rename from openmc/data/effective_dose/icrp74/photons.txt rename to openmc/data/dose/icrp74/photons.txt diff --git a/openmc/data/dose/mass_attenuation.h5 b/openmc/data/dose/mass_attenuation.h5 new file mode 100644 index 0000000000000000000000000000000000000000..f62785140c3ed50eac4f54a7253ffcb26b3526d4 GIT binary patch literal 135056 zcmeEv2UJwa^7o*a6Jo-!iaDX8Vh(l9sHm9ZDkj97b3#Q$MMVJxNy9Kh5|9ikDu_8^ z&WdTxiaG1o-Bk^Ld8?20ZD7B{c#bvQx985#x2n6k!mn1Dx9?aef0_I${>qzIrLQs& zZ{>eq>E&1?hv<@RdEJWM=n=Ii8i2!9dg7l-Wkmf&>HZnV_372210N;tFDZRYWvY+w zjAci4e4|j?-UQ!RQQM$BzOkit>R^1MruLlS_$Gnc#_qQrmtDw$;>5IDqX0YMUR!c8dJ`Q*0aO#NYW(|8L%*jV1pLwz=?^)S%yG z8G*m^uljv@s$44Z4fHj7sv6>}^m1T_oUXd`12VRAGM?rmiyDI}uNb4K%nGQ?zS8f{ z#d`Vm969iYc}Y@o>d|mj+7qE~m_ziRA-%izXs^ds3H)1Bs{h6xay6Ml&+U7^TCNZF zxxY68`;FlN{RQ#K=2ys7O^ z?QKZyBdI@@`u9+O5>kuf)ILM)Yt+7tbliRFe?t8)ssEjzwmovRQSsNfyn@>HFr>B- zw(ZJOe=Qc$ChqXf0OZ7pcPQ5!*1BgYBa)MJ96wsw;QwXZc*&}75uZQ5GhO|EpMoN{>9 zZtpVn^prId?tKh=o(jH|w|`ps>LwK4ZR%Fb=o;909S=Hr?lK%nY8x{4`Xvaz`M@Qj z;3bG|@^0^_po_3GZf)+=k{999iMPFP>n=dS9JN=LymTHKx7^ihQ^oVpG3MpxvvbbD zZ2mRZ&cgGN$A0a4=nR-!elnVV;WX^Hy!!L#x5;2Ky7a99rA|V}sxxzJ$$Jzwe|22+ zvH2mo&jEN5W;uS#?0vAQ9si|?@T3s`fe9di|Tc`qJ@xmuKP=&}PQhSW2> zJS_$)w!GYQM4cG0wOUtxPi#1Bx|euu;LZ>zaN426nM;A-x#Pt^uRQ)xW^7>mFE_O? z@I!L0*L8fMMd@-?x1I2TMSjCFZy07kI z7@SW0X zTz4hub&VnG_nInK-Y*&#KclQN)5%D`UtyIpEYFG3Jug3jaot|7>pkNZ@$F@}9Qo>uOC_W}b$YwGd-T)8l4er$6os|q17Cu-`xK?{Ol z{L+ArDy<(>DB8Qys-b?c{!?PY^k+WsaYjpPbptP$Q0>Nq=_9v5zuLtZ3vY(kV>rXQ z!2NNx-)z3-1dBg5pK`)}Bit@+_NvE}4X~nY>igk#j?k`O?7Z&t)`Kc^)WMF%>tNln z%8N}?9pGvBE1v=n*T96whUW(7TmuW*Tdb-wZxt+vyA#*DxqUW&vxVsIT)(%tu~)N+ z%Hy7K<(F@@RC+P;sAvnNG~@0%4azHG*?o{#S>nPM%=fTBETGSTb6HE16o{xh$yYD_ z9H1N=g%ppXx?PC};c>m4)r%C^PyI=7!v6I2lmbam+sNf?KAXKTPY6aL9IV-V&&`KzS={rp=Y!8fPCcvg+4I3nhH3T+rDFDc{NMR7u97_; z|93uA_5RKC(FA|IO4|2i`{(U@$w%>)Q!-ds|5^JU{h3*dWS`&uw)30NZ*MMDtodt} zwQ~Fzdj zn$`*|mV68j?%)QC18~%*15=Ur!Tb8Zg;u?j1as!;Oxl%8f*H{zJ|2CV2-OB$C@dPby5!q)_GBbX5sGvKod4R4oz%kU+!PG! z9^~uVQXLHC+&X=+N@4nKr`3nTR{KMH-xkb7pao^+`xiML`@;UFi;Lwis)i*ei&VDj zd|7y4dCrQf zvd?d9rTIJ0Z(?CHb5a`e*P7YL=Wj+TZ@9UlZztmOhRTC#_45V~sjL*P6XE%+dP(IS zv$DDG%b^qx?R0Zq@mDbJ`qb6y?p%kbj}n?y?s5hq8P{8mKLQt*Oz-IY<^cRMXjF?E zJ@&)vD;HSLKC-**1!G~DNCfx3-djU1?tvA9&cxp?z6(?iW0I3$JMmErEE{Xox{_ZM zh~rV|NVrz8Ny^KlFj%-DdDHZoq0q^A=!JS0Ltx3{*H72K4T27Z{D!}?4*<^1Q}XKQ ze5Qrlz4@Kh#E%|utU%cM?!UT$saKnJU-dUbH=mbfcD8qceDBp`A|E=@@7xH5QtCLW z-)(>3f*cIxr&3SY(7%}Ub@5Un_*#E^rP;Fr5;wcBXccs|y zBnXWevG`J}MBrRLY}Ou+UERUf@fy`H_ceC| zv5e2$3>E5fOYtUHX3>C|Hk?6krKQfgh%@Zg)*sZVuM^b%7@GG>>_)gUvDB7FW*cFp zfBz#5_HBT~HM^IZjobih%2*_&7To~D`}C@=UiE+EfG;T7`pMs>pNRGQdY~fcV@IUP z3z6R&jzB&?+g+K;m#Zm2>vaRAGB-rGEU%oK!&mkq3SZX!54?n~oXz*&24N=lJ_7}v z0+;sed>me!=8DZ>C_Vjx;fpp0p>6B?mx4d+1OK_)0+$3g3-E!7aQmR^&P|Scz(5%2 zcEaej^Xxm$j)T%iR?RBWD+cO27kg4m8wG1mH4oZi90ApWVspoI34>U4=ps z0^B&fa(NJNZ7S-JAC&UWVbyY(AM778zcZMZG+i=ywiASKGy4~|zVq$E zbt4G(3Ty;*(t^2qZ#F=#fWXn?7jA&k`=@o?T5&`6`pp>aZ`V(rh9h5_xgw2Qj{IGD zJg(=R`y*e!>V*A~EzsZ6wz?wpj+3R7gI$|2)9_~y%T?GTXkNL?{AKknLe<>dGMx;w zdXz6YuG>-aha3WCR5stUFx~!KML5(Fq1v5g%oDi>m^M^y(5Rg-sL-nUCJW=B!^PBB zT?WPgGXnQX42XjI!Z8s6D&g1&gT$}>B1TpYgQ~}(KG&`t0(pdkFo^t`evrevLdU%8 z{lL|%q+HB|?kpt}8zBFZ=?__g7CYxlBF5IuoA zCN~08%GPe4v=OE;?d?T{jUeM$FsetT3M)DraOJH(Ru&E>-VopA%Ge?ZOe44`f+qBgaD zJT9wW&8(bBKKc{+suj&mmg1{SJ+JKZ+uwVBGmKHhb9wBA_J@uRXlE`s4cpek6dqR5 z)7%pE+m~%ozs+Ha@fBl@6&7v5PH;IDv`4*wEweo|mA7Ao;Ynv|k1l;0ym|Oszhm%$ zhaWi|N;|(9vS^D^)sjFQ^E?wkguRvBP4R8p;q)*2P84q*3nGlIZZs6%G)KAQ97%B! z;ou(Jv&gCYq42gmHzVp)91r9Xg)5bU3Ak8zgnw=gd@FaAT*)BJ(zZouT#RFD@<=7v~It-8`;7 zvvh{}-RIqScbVz24Q_7!Fqg#)Pj_nwPJ&dz5aWs> z-h()AZ%6CZv*7?!Z; zJOBPc&~NGfo5!~KLFmr8LG32_QGAjHwg)}lS?QfORDH1JaFZ(@urN7t;Op&MU~LT% zF$%5DRcw+}DE%Ba{nfi&-4<_x$*(9vP5aI0Un@y_5_e|DXU6`WyqiC&9(J#qS$zfH&kKKLRHb0N z@;%_zxs%qHfJZeOzBmP=M4;z!c)-I#KOZ3b^*)LtO@sh$)(ziH{JaCM+H%KVEL5=J zoE8lO6ZBamXe9KRTxf!e6_W>U_GqyU-s+i8ymcoO2A|Ny1wGyhhfc4o8!$N#oa37x zw0*3DL56C-O(8n?)qlgLBUf1ba`W;HElYdDlyFV%Go~JpP;LC#c4k}P$(iS_OOCk0 zbP*W68O$z>)#_es0=9@V2g)XJdp&|#u3f--ml3ldx`3xm|B833yTFLu@v4vaouM-~ zgTylazq0(d=4;dU2aGRt{#G6?Zov5-<@_zpmWRJw9)y4Kj04I;13R35%LPc=MI$w* zyb4n<;(L;xP18^w4v)nC>HRU@v}=21{azMWU8IpxF>IR6hKZ#W?j-jvlUqSM_7mdQ z8^FDA9!1ZAA6E;eC4&f$e}4$N2}jO;;MoROE+m3r?;CqStiJ7b!koqN{X;C`VB(Py z*Y8$}0kb_@a_Ucyf&;6$10)=zU&iv_$lSGVZv+2a^*-c39Slb8nw<)5%lwcmS`-5P zz>}@&(PR8T&*-GuXR!w0_TD9qCTiHE&kf9;kPvZNTYQE)9OCBlB3shp_dPofsx`Tj zE8IQTEcd8wo1yZ{Llp|P+zd?(Hs-QO*#wGcTw?|-A3FGs?0^c_L+tJ z|1cD3s1IIO+hf~)7Pd`nkpIsNQED>x{OcN>(NEu~ImTU|s-U2rQ$P`J%7QPT+ybr< z-2!gie!1=f6e$qxH}~KvsL!%$T+SW_Va`6e9~v_4XOeFsw9w-cu?GeV2hR>#mt%oP zr!4b|h7lzuk1GE@5Y!fnZ;B zukYH*eqi@{ap9MDbr4{5Z*tp;8n}1+vSs}f-f-%E;Spo?Efbp`${i{m%Lm?BX*+V$H^`7)c~q)9%@-bXiq_*xIPU(mTh%0tsG zLp?_->2kql|3Cec@g6wOl4kS&U(5gYuE@{vw#d)Dr{ny;Ablu8N9+GO9A~ir+rE=< zoYe?Lu)m2Voevb$k3y>`)wyHla#4)OdjE#}8^1vi^X88*xC+JLz}0BGvqX;p-!`5s zI7oRYdqLzc^h|({+$|Wv*8dUb*Z6qEf$7CUS4xhKfup-?FZ~)4MSL6qXU|4N?=}wu zm*>wK*8M;EFClw*!hyNlY>XvdAQ1{qct`S-U3tH&Wx@f&itCgd$jn}$PF|+T?G{2$*;K0BVSX-win4f+vQb2c9)+zUMgEI}QRzQLyT>m+agLczn#U6D-!X z&@Fo%3(s;oUfy>&8Xm8hxTZnbC|F-Qz+tgB%cE^J@6&a^Z4mIwmo33Pw?X9bwJg&h z7#_4MbvE``AROdjtCmdvz1Y@qnV}z8-e}-|bc6=<%e8gc(n<{rMPRHa+{(RRkx@k! z$H$LrZ?>exJIyJ;Bd^?`N{3Mc_Lp@_vlp&qUZr!6uF!q)o5~$Vxq`HlXcjMBQ=*Y8 z;WD^FF_tAhxT-5O8p$)^o7GxGU~WjOBWOr+Dtq5UJEKc4>#I^*Y5^%UmZG}1p`hWte;g`+A|+58O7@1JO9 z8+(iLYA%6*>2OE}#2+~f*5bgFM0g?zuyM$jb|0SXgnJI$VG~F3)iE$|Y_I$~c16Lj ztG6E6T{{BG);jRs-SR)<`*tP|VV3WpG+YlE$BK9izHp-VFVh-cRj1j-YumMb!K_ok!+8lLD8z&-a|AqQ!4}4u-@!)1|o16;Tv!4U^dfwQ5hx~eL%SpyP~q6v;T%%3gSV%q`aKpNZ+_d z_#CHTSx=6na|8rjPP@H?zq1=QuX)$~*u3pv;(d4Sl-V($H*VaeLBpa5hbIDZJm0Y> zaZ(t3y=OXXsBH)!epLtvb{-l8!n~d@07Sau6)mJN6p@A+EeM#9K0a`or$0RO0*;!~ zG28NHX zr**Cbrti3~AKz|&`gI-^z+9}bVK_J3XOjQx2IKrnn$7=WANbSwKZ4ei8k>-xx6nCY z!FX(Mvcmb_w-M=M2h@9-*y8ovGf*yuP=3$8-k7JEs~OtGuU1t=I#1UEn6K&o4#u+R z82$Z^;G{KAkhujSP3^>4;L%`4`X}H{aUK}64@9`wtpv(f-U&xoHiNc%9F2>C+T%!om+Y zT2tx$8kp6j$IL^=)lf?uM7-0^0m6^_!~;I^Fs@!qE*|z{scr7CPKgg47Um9{A__9Z zD)%(~sQ0Fd>t;?(Z}<51y~)XC?jW|EP3|yCc{A_rL3i+a=;(ekw@2pIf64C-^qeHk z#`_a;@W0^wsrV~LB=YsXEy&Mb>`@=uK85uAp~&CCzv8%e>ye&aisO1xp7N&AiojPJ z)F0`QR+PtCOX1l*cP$Mlj_fmFK4&U8HRle_yOh6kfpFK7foHQZmk(4DfdYxpSR9n& zVP5YI`9~Jm4z)OnUO)`3=TY$3>~ycWEg~SAp*u|(6b7@VmSrizVGzZ*{kP^J5WATt z@&rK*ljpV(HT_|bSDx8pp0NBLp4K+V7uGPe!8^l!2wzYQ35(>Wh)7KG1kC&M1S@fH@ys-S z)L#$}J%T-X{>oBzSZUw4$KR$>wUcjtpE^?afDXTZUG}=t8=H=@=l7|ay61HI{pJ&kFM(nhrqyc-AV^hohF!`EU;E!I!KRktT6-0ObQw98_HsWwPP2gd;vh@m#q4 zE9tw>DewI^iVwO7Ifa8Inev79fo-h9LX#38r!ds*fSB#|TMb_o3k;0va8lc7$a!+t z%C74oDSkX0R`7IpZ7AjQgu>ciHS@+>Z3Pjw^gaMOU;+^3Gw5JaPiqD_qJgmoG!Gt+ zQUi~Miiq}txQ)M*T{4Ho!M1Kb^KPUE#Bo%JrXFN(bSHc(cPRJhR^u`W?y&7dedGDN z+~E{AvmbFMTq$>WIOmu$`ja~dSmm-F@b&@k@Swl%51#+5>6{{IHvj)Q{QrdV-K^Be z*U1#;w{HgW_4^Sh4=?**d*O2A`=La0&&SVMjHPqpV2lHQ)&@TpSx@=I;d^csQ@RVs ziyqm-UxFn=J1|^!Kg|zhC4e7JfiNSK-47i7J@!-rOb`H2JK&(_>BJ8ZOXsdA@)JdZ zrATB*$Dg{)VMfk{fPbV@V$T;rz^x?voC7EyLkDF!Tuyx*aBnGm^nn)K{MT9y9L=QC zDleGH;kSo+!dwBg|bhkDoDNiV; zBH-87z6JsJu6hK|ns(?k?fn&a##Z%^v8O@6Wv|>1Ed&sD0-VcdQ)$AS9WZ|M-g>WJ z#=uaqtVaQlwm$zM0`vugh5<`*IMO311Vngt{}ABS9NsSqf(>y8c}5WBUFaac<|&IH z&_W2yUVF9M2gcsEkA7ENot6j6vlSY*^n{14RV>9bo&R|>xkpKlG(MO7F4uRtPcSep z=E3y<2I}*&&q|LpySlWSbmeIfGk=o*M-#tGn$7=WANXnhe-e-UTu_JfofGnN`T3;R z*kHeEFxvAbQQmU}@)H>ooiG_c*Kau4y*r{lG~5jFt>Q~7SO+9G+5?l{k$>kYOzg)a zPi~OE{VZ_Q%}2Kmq~$?;<*=QWJ7E{gGCS)M3j)@?VKnH*^MVGEz@x+22#Ti(g_jjL zkYOm{cW#Aa9Trr#-W32TV}`I)9UW8_32i#y-uufjyxBTWwdE8zwx5eJEiEF&zxE_zx@%p*-U-J7JTDK%svg_Tb*z1qfyW!(io-}XW zDZb8Z6Y_n}MYx_{rFGq6s3OK&^}@W*oZYa^;spdAZON%P{_zOZ|0{Jyd-ux5K`=}xA_2faV`nz$!Seq_9UtcZ7YT9VfrEOw)M6;g&D*}5zkbZm)qO$!z z;1DgwTbceV?Y8po|MTXHIeu z@gUY%ko}_>>Al8;Z(Im*V0(N8#1(o;^~7$2D6??*0njTFfxZYF^oA2W9iychN`Gz55GU0C+xhR8-T1sD&F24~<^oY{N8ya^7EwrNY{13>rDv{!eR#Uf9^48_t5Ks=f6BH@x1rA zBE^LlQC_j?rw^*>DQZ^zv`X&RAmEBu+#z|m1XpsT$j>BFMa1eU&Aj;0M_;0j3 zL=$dyBy8{9e8AUxtS%(8E|+T+3g?d44SRSmgyP17fJZ|)sspIrqYebjTyGuJ+P|S?>dVy}+seJtN0m>E}7boA8EeA~4q* z0{Gp0ykXaE9?0lT^-t6w`~V}=@Mynf_o%Ju{Qs4|hEZIZZ3ZfFQDk0FM!#sreSze_ zwqt^LK9Z{J={aZD0Ty!dqdLGU{!(zhYzxx7UEQnjZq@AlNLPm;RjHAhtU$S0Mt-BQ zls};iMtjI!bHdB&i2I524E&t&IP8z?iE>st9X~?YRjdsaSLP*3IFeHlWm@aJqdE@v z2p{|^*$t1S#m%A}ZWq<-iG!=Jc%f5Pw}iuU_lu&sj^S`xbv;ME)!WGKv5oTHg5el< zuZFPw!a!v5Kn@-b2!6nf@T?+`FK{U+uknFbtUAf-^J@5|=@?dWQVpUEkfH`lW<4q2 zP@NV}h;h4Ye-RE@)_*G~;@N67te-L{RkKZI@3+`A8)_%yDt z!7-n2_lqHfBWJvU=8@>n%F8O$18j+A%ZrVWlg!l3R`Hi^pLhZFqNIQ1Hc!8?z@cO8!_<7ZMlCxekULm}pstV$16{CDR zeFgDU-UH^L-lw_|*Py9LFg-^8{3PJsrFpt|D70=PTiN3P^WUOD)YV!Z3D0hF19%t+ zSp3GJFn}-r*F&J0VB{e15@oFe$Zw^Ct6j_Z&g`RuTs=85odzl~bj5d{)KKE~r*1Xc zt0_+1o8*t#8%2F#Z(v5mUza<3gZH|FETh6ZZ9l48F7M$S)zrYP_qi<8z|r&G%w+q= zfrd8sz11Mfb{}T`tAIpSk-`To8{f_UY+43+aa^zHsn{zeo${ut9LbSI20Df0DPW@_ z#!33K@^x*mxQ{K-Z29_w^3|Sj(c%c7(UhpeR(w7Y5lDYUzld!qlt0x}(vy#qpUay1 zFVc7i%71;g3iqFy6wfnv9Ny2^it5UBBYVD#YnY63Sb3r%a80}P!8jq;wzSXHApEZq z3f4i+LpWid0PA%~tW~F6Og%hrNbHYPf5lgaD*taw6g{ z)D$nTPOCR3>&{`F3^kztTAdarg#KnV6cPoK)L=dP#59KrZz(5f`??d&v zGRUJX$=d{yGfB17miUm8=5@m7VNQ6FlK!MTrpaF$T(?z{X6uJy-}|rnq1_gI9+o7B zw)9-ihhTpdk}m^_FWN@(7q$bhUx-6mmHKZ`e3zdu%GI_tcs+S8Qp*XnFAhifTGoT~ z#Wr|d1-MTJk4HJ}-xtUEwnO>#s!97|Nd@bH_a7X-C!X*bV!_9U0~*cfyS)*(f0amf5dmbD9oo>-}E#Y8mAYWOpp*4@H!M*AB8$LDI`G{2+zp&FAf}8P-MQ^a$ zs=p@5%^M1Box?JEy`iWmSjqM~R(z&jJvHUMFg>z~kAdDwCZ~b%zub&Um&@%ne#eVgKsSLGlo+dppyPj$sx*=8v_xb&~w zr@zx+CcXI%Nu_?9+20vJ-S3Fbne*duDxt8k*sitzWn31&oLT>hQinpASrPg8kGyfVQ{ckNv@753c?-6!F6({a^L*P_i#dnyrs#P9Jyl#OE-B^l{gX*v=P@)Rgq`38ZI##dET1 zPb^-Kpt$n&q?eDRc3mCn+4~(Qtxb!E!}GbC>W!$Ok0^lc4UNS?AYBh{A`-;Y_4dnUIl#9XR+LSBKip0Y z>Lw*v31Kz)XPBLeMWvo7?E{s{Tx-L=6Uy^)KHGfY_6;7Gs`=Kw`;O!>nB+{-SJZCl zg3nD-DR&p?`Ahl-{U&MW;(7=1{*q?PVdms;qZ;S!5IqM2dX7;fhpUJVC7ML?H!cc4 zw}Rxe1f83IBOI_%jwok|i%_nf54+N0fge={7<{4hd}Gy}*^+!6ar z*HuKl)dFRy4sjuc758NY8}!hgH7}L&aBfq*_GA$0II07LW3UTqaCoQEaUi0lt3<;{ zVPJ{`0qeaW3`AMtAuJz3)IA9Sj)K;#auBeh*2UUi^M`J`Fa37ll@9IG?_YWPrIJ%JWKfrQB(H~T97~8PcQR!q+f-jt-W_wr`44o>u|p3Iei}x zWx<=XI#+9WVIm(;r;K1F>3v{UQcMT8LK?WSN(8WeD}QmceomqNP0|Ev8_4HNG)?~8 z@Ht*2`m^%4;}G6Q(ro$5ocxv1;Ji7|a|tF~x6fhNUxw%ks_$@__PL8Be{bn}Ly|)? z!pCC%08vM1Ea88araCA!sa}?OAHqNEi2GUHCbU0Nos=^J(f?n{9OGji)uX(Ka#%O+ zL1E=7Lth)ZBM07R%5&Hmw>Ecb$vcqTo<|6rBHW$>-Xnm<^0d;G~L z5N6oq2snM9D62Fa2Fqt}w|}-TnCcb>!^i?{3U_f0gwRnOq23>ORk}|rwZNQ^4B@1@xtQoRA|+c7_ksH?1oU96g1_;NOW%Ix!#cCvT^ zk(a^zodMU&T&(4jhO-i#(qUMeEk5w>1}})BN%Iq__)F@eQl6T+;&YSKncArt=zuLa z?xOp*^k?O-#}T}bq}lS9Ir+2I;=Ea|M|pYTK+l8ZDuLwAkm7Bgk(@cxeyG0<$FHNd z7saKubi#fU%DXF0xLY1JC@+rvFzN%958(_``G|pjtSvJZHs$30lx=#4CYopiq;Pmq~ejK){*W zr`vZcO{!B=&Bz30trwnwC2RX6aoVpcuEu$ECAnKf{?Z+^ubm}1t44Cyn66)>@x~;ltWKNQ|N1V& z``HlwZ@Y1*f0iGB{!Q~PbPjAt`=u4iVZ&b)fp5Bt>PfC_fbpP*jFnNm@O>#m+6Qwf z0xqG$E7CvjQ9YtFFowhBy*WsEe!GDcmn@dtJ`T24dI0 z4Tk2TU{Nq2K4~EA;bmh+`GZJ@YPiv4_nTuI)oK14^vkmPdXvw1e>V4l?WHHNR4nFqVNrHV_xr%!o-ZDj%c+6& zH?FXVBn`Zpz!M)eVBWA~w=vayLE78%>Z!zQ$~vUq``>`Z9di8M2<$haaThoIc$_hf zyT9?r!;W*$=HE@RnK2ewf+$XqsU@Wc%mssmalJ%WVGxpNE`!mfN%E zg|7{dinfzxvOu zJC=uVZnD<8ld0#G@%(E_dbOkh)He21#Qh|ddUhqUKT4X(c=jZzCOdpZcFr|4ZtBC#Y?Y9BowmH7>89wml4~ZG>&R^3-36`fE^sT|uo(>j|1{ z*bv*d8`JgHbRDR_BlUL|^lz`LevI_#g0h~vz4-f8rqr(xzQ7>T&mWV1tqsC)8X6~Y zz9k-YNy1sOo{i%YCZIk&eGuWSQ@w#w4biSNz5(gmLn%JE8`*cw2*<3Fg1AP-$p4ZD z@2h%4cm(%JA3jU|io?LcanGCWfyt@M8Nyf`;V4A|Gvb8XG-vUMU-e7X9v22vc{saq z2-K{URO|YoV5$!g2poNSpqD?)e|vELxg%O=DGbJ1!h2`<0~9x?`NTZW*0kGPjIoC z+SC^USK7?T6aJm^t|gscC6#hD#2e>H(s?9jmOj`o>CejBw3B!rNwejRmv*6&xiPQ=@^ zIA8$>zy2KgM~=XWAYOoX56mjVBl)%y{&Y0uy+l&}U>M=igg~d+1}%Eu4~7@4di%S# z?0n118$|ho2*1Cj1yPPYN(;OS{Ob;yG=GRZ=Sq96%maHAFm7qUN9MO|+{fW%lk_~4 zGHd)8q8h80U=-pTlurYq*;4}-O68l_Y%`1V=j9WwYCyn;l=p=@3^nJow=dPH(*C}E zNy?W(@+9eMYM;+QgGk<1Yw)>9`m=JUO~(63nk{!fmOE40xB8I0)LV&iQe-*KqaDec zl`pm}G)NbboW+yeSr8;ZVhoSwbPEWMo9&AkU1g0c^ z1M%~fd!qdHXoBN9RVDjwDFy556d*Z#PjLi~slM6;_{_`BFF6Xr-LW`<>Naf$f&N}F z2F|&he9@BS8?R~7VEw1WFv#1=y)P^E4!!TkAE?_i1XytphU^}c7VjHU&+ziJx;kJO z@@?9E)j~s2Xh#FQ`cBogK9EZ!PWnK;FR=`9(L3o zY^+ghO=@615`CfN_}E`kQ_^Of{x<%sbV0ACt2fdn7os+uw zx#0e|&l&Z?&tGava$1e{J0m5LZ+|ZfQC`ePFc5{`o>G3vMT);amWF$Hh*jlwEV2Vc z*@reUFu(h}8}BYhLZh0wt_?jMMs*WIpoah)Ot1HY^~D2#1;;Z8TOACvn9w@3kPbX7 z?QPoyXn++xnbvh{x}Sujeh%&E15)l}-U{Yn`=sIIP7;AgnzXv%(r+Tlr%dq$XI4$( zOuR2F^SyJfS#B+?)F|GEoU|}s01W7Uy#KAB=OL+-v$iB>l2)VVSDc=kq(3Qtwrfs{ z=PaozTmF75e{acOa-aOn$0T0#Qk0kcBzFZ#?po4$LgFi4rssT%u3OOaJ>-d>*E*m) z=eNcE?96o1M=8#(H07DaHK#ZO6ZCh!=!<@r4!vprY))}JHIxsFcx9n73i|5}NH6^i z!tb&DBdRM7yc|S1PUC!JE&H4+7rC-K7H2Veex3iT``Jt zavDi-xqZ?8ySgReIn+Wu{dVs(D&nL{x1@S7HECQqWw9vmSyZW?msf)R3aS3QFu-%F z?{XR57vTruld!XTkpj_viLkBeon>Cfcfdujfre_o7~tNP7Nw&A>&S+|#|@Pl1m6mQ z%f?O)9?MyLeR8hXb*2Y_UtG}Ros|NhDNhHt)WQ2*-4?&8tfTl~4HPZi?|H}HeBjX} zUWt?Cqlh|-EZ(r!a0ac;>`{aZEy_6$)DRAz2Kd!B`I!HN!#oW4CHy^h?j6SwTeNWW zO405GhBCUF2UuyL8&9iguY=w39Ju|zoPQsZoJl%_`u<+zcWH&6bE{2qR$f^sEV~jvOhxr) zUVy-BnRgY6h=9Rl@TGg4Gug z`2HayuO2zZrgI#^cE!Y_YEW^WSlh#J^sbY#`65)5}m2rkPH%gD$aqT@BvQoNT$^7cb2klB}!Lc~6Agu*69`3kTOL!DI z;8n4!MCl-=&XEa@o&7+?U*FsBOlbd-G??0TNvU<1yZ%|fJMBwtl14N?br8Vwv zwga$#pn~?J4|S-XM>V`|+7so}8p!@z2jwrc0?Au3v;)WKk)P!y%u4Yns@rvq;s=s} zhx=}xv=?Sn+5FXU(M}NV5&c-`KZhgUMFHk*g~M+fcnND($3OrMg}^>m9q0PgAUL*| zN9+cq;aFN<{^H_NO$Qe~@IuuZ2oC_Y5epOw3i^T;=nX3O30%Uzi`oTpbBr0+K3 zK9e*b=~nVfYB%HcrX+Xwow0qD_P27h52|TgLDECbH{yQhY=?5@GZ#OXItK6Gdl1@j zUD{B5fhp~GRZ*_;T9SQNA^ozhvfG;h$=;zj`pv~Dzbr=@zAW=mUK}>-28i;44yVZ9 znFMch4|BCB5f8n%xQvLUI@?iTT(ITa@^!<3SH~GqClrPlw|PGDSO_qu3-|m(0g*@+ z0JH7REG@N=*;h%=T=rp8vtLX4!Z=ariS1)N{Jkj4uVhu$HdJDAX~rEl>H1_Y7BysA z$QPvlL|wutV#hgOIMa<6iqit(G-+Wd0|o5!QJaPJCzf*RU~)8BzF!pE{^s|6N=AsUyk>CG!O6Jcr1QCVKDYPwZ(H@{`weC zP^1dV+mdc1r|I&E@jMk#?{qCeIL|o=@B9sHc)}AqQ^1q$f{CS0Q{BQOV4z};y;{ao z{i|5;=iAfLD2j&*2lSVQz` z-{iy%ySDn}G@!8EwABzMmu9beuoM>v>(8^Bh26wEW-Yr} zX6=YgGR%i;{||P=yLiigr@dbO{gOXq`{(CFeLpMChqTw*=;5z3-bak^9WsT+%{lO$ z0n|N*#)Vz_@i5nJLT-pUx#g<;e7l%*Wnh#=yKhW?Vn$F zR4?(CVL9p6MCN?k{Qs%-J$sWp|JEfvThfcvwkLgCQV-I*<3kiNJ(B*c{{7&hxQ{K- zZ2kK?arFy#;(VVVyZap%Jg*j{c%-g#@cMnywZOI;0AJ56RD=8TFS4t7#MvZcz zeQNnz!Xdf|yAw+>D4;VS;1}K|rQwJHFLQh~77jCQw%Uy-V4z@Tr+bCNtrqoJsgh7& zAPRby9)*Brb5E1ls30(TZW~dve*omJ_aXmj2OV&9-3`NaY566~%toG?F-rq|jRxnP zX36l3CI%SYo2-G|*FQfW8qL}buiBn0;tPiM7hV^b_4LRERQ9Wos=rC#kOG)?tvAtKypO)lI(xD`0 zrU5|>f#h-S;(cjRt+Jj#{C`_kWxa`Mm?`v=TIsv=yfJC4{6 zSdTQy1!;S_UYuy0Gs=&?J|lDol)^xZ7+^|n1?&5~qq;vgp%gc(PCo-!AABFh`^H0PrD;m+L&DjuO#@P`sr-s1IC=L`91P`9icmML(AIe%Jr53HPgU$?P$T*J<%8Nn| zES`|+IJ9oIQG~3S9hEMN*#1Smsjq@~%?12HUC>J%BqFlV6f!7C+ zoLyW;^0F5D9cllQxI_Uot{nAGScP&U{i7)pvAxU+@6)U!#S@reT&~&};{*aaVSN4W zcIbcE(E!K0RZ-9`T3W$)+C0ht5up2?;skDyAMz}SG|hm0U?muS7dY|q<)33gz!aT} zf=Gs%xwv*XtZvzPhxLmPU_rbr10@8C*Qu|t5-=ddXSpxDddCA0n17|)`GA!+TJT)n z>PhMamj6k9Qdy_sy|sJUGsXPCo~PkV_k$W@qxA>j_B@+Od5fU@BuS-QwIsQc^enZf zQarb$e^lP=t|C83nk{d?FK<_<&U%8hw{O9DGhK%J+^p%y=L|1Q8b4|qdZ#-N4>FUC-ke#Y=rqI`6^HzX9twCE$z{dIjJH3 z9_y+YSKqdbBEofA0YK??T*vVofyfuiBn)OA8>U~@ue$O6=yZA8u}fQjP$wOj zb_SEF;YyHk{*y+cPO|QogE?JW1+B?Q$evl5VH=pa^`f zlKxS-BRR7@a82AtQdKrCfsnssL(G%5-Isv#+>!Kt4|nA6adt>gP<~0vX*kYg4a#9S z<&|}@$9@}%(`ipQ(9hOje{n~=j~{(*^b-91)G0`tji7unbBYsatcWl;yYjT3{epR& zzqTiNGb8z{hWEEEi+ZL_KFVMD1T5KcIF4)F!T zVWvm~3nBc>5Qsg+3jhW|FAmEU!S)%34!cBO2fXUajVtN$*KqaSP5s)X>zOC^dd+$_ zp7}K=OlOd{Ebp1>0?umCw98OWU(hir>sw3@1E#RfNoe#@zfzp~}>8@W`?A-Noi)V4j**VJ!H=W?|N$r1T0-3j+Hp5)7i z&UI#FPr0@n$1kV2*!Cos^+_KsK>DkR13s4&igUe9c+yED@H{xOGwP!@O-NoVQ2ulo zyuVvpl6#WZ->Rj>xiSzBhUAf7K|I<|AmG~%-T`(PykPj^9POX`q3ggEP0#P&1z(s| zq`}Qtl9wn5`*L#l?)WgM(>hneoZKN$-?`Y6TH0V39CMx_4Fm#@zByOXAFPWXtrN0> z$=5g@`Ns4~vb$(MX)pWVW_8FojOI-rsNHBj3oT&#=`kwsiZ;^sWM7km|JQK!p7b@`F9*HtZ)(=m^yl%FWe%Leti@NEosZV^VC8(s_Rq^%YA&3Qlq@A@nVpX|8RkQ_ ze||oa^Wl6L=zl9i|Ezy2v-6>3m=D?h`T0mGg!7@wQu86#;Q;z^VQ;>DS{YkmG%oM^ zj90#G$I&>W58w7@Y$edR8Xq%W`Twvj{huaer%1>^C9d}O^!4MUpG!KF^zv63sMO2b zM=4?!C7lwBZH4SD|E2D|(jY~Qmo!^X&ww+wFA?W+c>>a%WH-)DZ7;Gnnt9`OQ>yFj zKzS_nr(nM;;Ura9f%@|HZ0EhtXE%`_a}dla6d*L5X(RWtP0 zRH%sY_KnM+UF%^R{GQT{QJ?p%p$Odh?&T@3rl2Ci@D0C$CH|WjFgHF8w6$CK# z0XW}}NB-@BT@E9*>(0f&Iem^O77Zp_PIS-RDICrnu^aaAUI+|nRejgw1tD;*TKA{3 zYXrei+l4b89`T3!Z#EvV`Gw`fSTvYc@{JZm)M0L45bE-57I$x1{C=^y={VS;T16@5 zml1$2EYFO?kg2qApP_KuF4O`?bIA2U3)Ui`Mh8c`a3m!bk61^)S&iX-z)u%#PWl12 zY9DXx5BY7M7cBF^A8zt;6X6-o(^4*9ZNTRv>3@^YP?AGQr5vWxxmi+olFM$*O$o8*J)n1vC} z=7(i?{Wkf5Hjo`mLvnd(Io`*N&eM;UkzbYMv+`(+Pl&ZZd)CtCSa0}BMT$o-qIiA< zzbCRW*`;eL5uVRjMrZ|!V<<#&s-igLd+-}*Z|p975vGblMhBsVIH>M{m0$AJ?RPc~ zM0#e!Xjpl(@`x@~SpG~XPaqAUepZKK!KQnO*8-s!s{%f4xjzVa!JAr$+075KEbicX zAMM_Y>2V9v&MfPOvtli`W`5SB#p^38-@nN5*cW21)HS;>Ha$OP+#98~>8SMl7y-Z8 z%MW^W_Oz|G*AKA1hCgr=MD{0x$=NNx{o}IsUu~+XLpJJRX~XZ8g!WWbdxK4Ewc9(0-NQ9zS1y1-2L2;pd%K;<$Hq zNR4LTI95;1oQE&0>`u6zEzv%vUs>^U=K&F=XkP}*8#qaIppMnT``8*Qw{?7AVdW@? zL;I|vyt(@@f}!j_YkCQav0cpi+(D}UlmL8ldopZ0aBEY4muRY^6b|k=^{2mjm!2Qz z<@$O=jp{5vfTz_Y`2#aPSD9B@2Pr-^lB>Va0`B*|uzoKutf+x-%kf)gdoX+@4ySaR z`I8d&uP$T3?AK%u>sg;?4zT>v-J`p2-mit{T|dVMn(N?9{O#hu9M%EW>F|S4<}GTn z+Yg$)+{!AJ`$L3nK5j{aoy?kc$2tJSw^R=N)}CFCCuw@LAh6%-@@N1Y4*PT?Y_N@=1l8JLf^Q`U!+nxE*T0*?i4u zWi&9OHY=1E4)LMG!pg>kfH(6_MRW`SUOchKoj~w=F^gsD{pWr!^JS$!i)X3#TEz77 zhIPy3J*#La(_***q!2x%S}P~UoaN;`y5nH&W`Fj zc3ykTY! zvlB-a4K~%PXem;^b;(VC?EBT{ihWU`~pPvs^HJlHtEHxj}uW^Wk zx#=g#x4%!5Xf(G@W+^%QOW9ZY%k)X_wX==NWzA0&D^?#}F@41b?Q_}3ZS#88~d%G9s#gz1YUz+FH zr1zWCc|f85BsvGIBz=Dy*>&9MJQqmz-Jx34^N&%SxB=<)##8ZJ(AkFU0AvRYm3sdI z)a(0^ecDaZ(Eu>Bx$m7aLI)>zt?at4 zg$}CQ?nlYQzyzKlwkQ$)GzTsbj4?|5LGSY8hb$V;4zD9y>ysWmezrQ1Ifkm zp$rc0h+CT;59hPJJA;g8^{rWStYHzR_jeQ(g+Q1w2L{5vVNQ)6Rt!S5Rm%Sal7AbL zb4jbvIsZ-uD&^c`GxCq5-}9?W`@aS8jigfUYY=})`iR<*#BY-RLAg(`y(4~~q^fMW zN4XE$htKT-$^A-NXPeOZJ{QgV7kXaP=y}{C`QAWu2Fbr#hxZ%eh;o*X?0M$XNbbhq z^%|5PQjD&r&cl9-rAYNj-aQth+_$9uusPWF7?1JnUwWav&!-LH{Zywsk&+7H^%O%t zK;34zK76RBTzSt5F}$j*2(zBfknDc?%0Y$#+s5h%l=r`HUia%YX#e2}L!vwi90hjn z!9-}wyaD&d?S#$KdBueoD7RpEgTAH_;Hwv(%dq%12zGqoG|oE|V(0M!DM4U+mOJ7D zfJePX?9f5UH&IE$oET0QFI$x03q5%mlwQ6de$NbFIB+%EZb7myuq0&WIMl*kiz|B~ zPBVN-S}(40Sl3rN=vHU&qcwKSkH^z24gKLk;Spo940CyLyWED>Wpi4m> zX%GlJ8X*5GChvTCN=%n`q2pLm9J@^hdXd(F_9X9;7N`CGaR!>U4!Gj;^C0@))&VK^ z3h|w!Es4+C5dTSkp2S}zv(VY@itVD6{&pz?IwmrL3ZwMKb&UWDsM_^ zpNrQElRa?zV!WUILL8Ti^odP#v7goJtg!r1&$tu4sXk|0!sn}jc`vO>VI0mzL(2bb zf}gWBrFIp@weR$}0618;Co5I2gAsPA zHn;6{z_X?+?e>M;JZ-G6FI4KWG1P5=FYwXkR~h~@tM++(0L!1Id^H)D)8X~`gns$_ zsBTQUzkei)B65uL2V<7K)3$X09P6;4y7jIA@Zw<)O#-1KZyE)`rpFia?t7HJUa0s> z=D8OnKY^qzsXZYB4WM&`8LbzR7N_&VrwmlCBWK-k{w2-qI?{spSJE!TpQDLyC9O~U zn?C8clK#PU#J0$9;&)1#y^hd2VxY$7C+oC&cwzgMFVZR$uP6PymFaot(erCg>xetG zN6>mQh|U*&bUmEf8^ZCvztTFB?2qjlg!6fE3AU{$ejsox`V-rY#&x9iLZn0I;kZMy zk($VLiPoRe3$Wje&LMATyzL~^Kl1b?zhXP|Kdi2$2)T+iQW%(K)S)*;@q7FllRi`* z*PmTgmE22Qm`9;F;r8YxzeOtHU#CD}R@JGN(Qy!AfiZibs`x-WtUp+=&w)L$u$RNC z_(s6tb!(OQ&}}f9c>`xQ*#?MnzZFC|!ZQI-p>DIM&V?C{!t2F_U*6RL;^8x#cpe@Y zm5x8So5R(2&_dxWC0R)}7MH|YqXy_8$E8QtDm>G{^9!Ta&a-FnemhL6HLl>FR_|2e za3A77>tq1oxCFx5r^d7Vk^_OGru2(t_`wPXEZq9-`9rQFDYTAA%JvBrf0<{X&U8*O zrFBG7xsDVkJw(zOw5~kfg7=d&v+GMsT2CakqIJcQ_*c^5w5~MUjo%~bA6#D&uHpJ3 zY4-X;>&q-!Ujk@darD4>j-&O(G6dV(X?<8t&us{;BaWn(Y^U~PYR8d&BHOmKo(v7c z`<112#!BK<(E2hU175|HaVY2KMcF9aq>Uzg%@EwFe#pJjTkPtNPLd~OOZ&s4LTj(I=$ebP#7!>AVJV< zF-NWl0G~0D52r`7xI@a<8Duvhd|_FSct`zK!&l|ef+7wyEWViHlqG(CWi+sW z&ZJ$rN4NZ7C#znbXEuvdEIfvlt_uK;LQ{2P02Esk=-s45Ae0k;`+?9_0Q(2Q*NFM- zPsp&oc+k2MOzViG_o-cl^c6`L(D|pzW}H_^<@(Z;))Pt9w64hd3X*1aeQ87MiKO#r zT}h<%L{hnq%%F3-r2i$)AT91O{2u;2l4jdO&>kY|$EVP`;^2<+)KP;}vlXdaAG(o# z(v9>HKUzN`sr{bTA04eDN9g)$T1Sj1F1j_XFCK);p>`tufb5)=XQ8}59*^ru?vbd! z9G#2fiq6Ki?+m2HNKZ+n^=8Q&{N3hcSDiG4^o)L}XP9-sxTsSmn8$FV4BAD^3^0B= zhZ%m}mT+)NS63$9`xyBAbP2-s$V=!rB*D z@D&TRO2^r_e<}FGkcliViZ2&=8BUu>WAK9wMdmf$+S(sj(U3tg=l#>pT_Zab@*Dmx z0QxNqjTx~x5T*#A=k#+>0{_~5TQeOmA^lJ}RQQtZpZ7}{v=BBMRo3!L{d0agEA}k- zD{IY%^oQTaAEz=XL2r-^{^S2tx@=SG$K&47i)gry~9y zNi%z%=t}1YN!{rj@qz3TlFIYMkOcf5Ni%z%kmm>c`}lcDEuK*OCDL*41hwsvqm7Ec z#^n{%wud3Ljj(N3p86|Me+}xdE2x!eJwcNV8)Ex*W4hj&t^@UVr2g)LCL8vq{(gd{ zMh>FuLj|?3HC)hS!;yj-eH5Mh5gBdKYM97oh!({I%+1i3lkpT1M;^_qVq;LecqJvI?hhP@%{Rvzx6~1 zq@l7dSy`%+R2cIy78Ry^=mz*ZY^y7-lUOEGmy&e;$cyJ1x96}d;Rwsnx(ULadFcck zX}|%;($5qBuf6Y%k11{YpJ1_s)kTjUWo<|hyVf;Y)K!BZh~A>Nu-MT%(OG4)kxV9i zk|@!GAliyv7OSnbdhhIepL2b$$s^ChhWGt_9?vgN{_u65xz9Z_=iKu>SHGmL=AIVK z)`t1TP?N+=->zsAMOF4Cwm*F=*}w8_+&9VuYN1?uB6Sd{U~N@kHi*{c)dH@Q`>)j0 z4W(X+9t zp_Tttb%9y1$3nyJURrRaSqR0h(-J(B?VrIaD*WpODeMz{V29|jkgtoN1>;*OQ23-G z_(M?m!yo)1SPDF{1^ZDjwS3Y6{2?g3Q53uZG|k00-46a9f*;B!J^sXX0H$?5i3{NC zVZ!~pi08;T@Qeq3@BP3NPhp4H1^ebu@J0oUi-JGg!5pnp=*oO^4t9%0!`Lp7ekk+ARM;b)p+3U|^!Gi6I!yR=A zzcN00ttDU;F+1A_>3b`8tXYcG4Io zj3Vz<%71^A#nyjN@~?h-cGjgyxQ>{9 z%F~8_x8jxq$^)zF4ytawJ@%SOz3I%?J(Ke_R;%u#-IrR(T@UyQA+2larKk{!FQXOE z4y7+|^_{l8Oc?29812Hzp#6z_k_bMj2|K3Xbc~OuK;aW-Jf{SQgGUna+!9p0VW+*{ zt7O!MhLrUWCGJb&JqAMl5Co5yb~0}W8mfXn!cyP|@<@E$8%FJ82^!KmkBnWTUk}4z z@JJ%;+e`fUxCQ==w%~_}uwzC7r-3huf@g|>H_YIT26#Om_`(Xi=5g@TBk)HY@;H^9 zQr#Yd!GwCjHNBZHzK0w#a5%qy8hK-nz$=sSdL_Kx5O&a#;F(_SVb^TH=Np}k)ADJN zUqw+jKL;N-*fYQEbVj^bInE!=?f|(Z3+n4XB+d)HLT-A&mLw{xJ6IEFw5O%&p>MXq zZW%`h$7mH3V@P+a?uxScO|&W~eZCGOl^ot6VUdL%H`WR;29ea?Ws6@yo)FTyPR?!|N?(t8b-_$wq>9D#9$hpAe<=UCLHjGb_-*vTI6egmf4G7#1ebv~ z?6KbjMegv&z7zaA{t&rin#3&v7eMal5554JE@Iq%Cw~vY59AS((@p(+!8Fe!g9G_` z`r*EwTf$=}+|OTvc{~mL&|oc(PeU%53B9Ktc&0XZ;{oK2ZQzsa7!L+7O$MJ_UBWyO zJC)bXF^1)j(Svy$--YE6d*Ec$ZM7eY*9WmYVunAqCE`PJ4&&DiBYFK}m~ZO9`J;a{ zgkPjG%ON|8A&xYsk@Ll~ARhHA@JV?i_q%aG{no5Tz3R%Ihcqp?olEBOEj_=T_$&A{l^c#Qkl8DDjznP>p z%22#gdDNH1iv`;HDc_npU{;kv$^AA~sqQwvmrd@op|>CCXQr!PWL*9I7fKH4v0#3z zRblx(&sVgtP}XuvW{F#pnrj}1ygFk?#ufQ7)J9#~z$&bgTth@B% zf`$*|6JxGh`uT!so=<#X4~YYxNdKyHxX(4gCpWtgvb8gjPad#-VEq9g^TS)ovOjs6Ms1PIOG)F^Ed` z?%MpDSE@fW{H|laSXy{g2|hl&dG9#-o~U=DD*GICwk*A#U~}q^>S?u zC-JXkPG`{n99?w{ga`lP-yE3dTkZYx958Qs_4m(n#Jpv0@1OTCKFTL`D7F3*@ekl@ z8)W?P{(>YgzH~Bf}&4&`15f=C+HCw1NgXLYW0b$ z@qC_O66+JD#GU$aL7U!?u#4vlHa*1S>Bkt+Z^;(N3BPOnU#v$6rg=TWZ074_551uw z>{PFyM~YpcGW3WGF+6`L#%9b&rFV+{LhqC`AYlUS>8T&Hhw zqsnG4@$2Z^jL$PeZ-76xQbnV?)<*4FK^)J{=MR32x|Y9F;xerv*-_Hd+_UVVbM9Iu zd_0bCBHbGj9ZPyTdJik5%dTajTWNvUo4PM zyVkvk7+80Ici3oOxk6>RW%0Jdtqlgnc*W0hf-jUxF9bppN#+e*Ag;#SenI*oD9^ zi!fGwz5`X4s8c1IB97aI`zz)t0Y7UVqsl5#C89rveIgThssiFS3c*kEg^?3spU{-C zn&un`1%cELB zE9q_EPK8ko%{}`?2o2JS#7xJ2dfhN*0hQ1C?6<}veN|G z&g-DEMta{)?$5m|+p|~Bb&~l*88@cRn;t@ywKj33LT&kJT9)tPx?v>zv2Uc-X!0HW zkp!L)lzfnth|3UEJfX{6{>f<2{^GzJk>Cr#xdBG)Wy8odZ)0ltBl}K1Pf&QGI(R~G zBgWtC=J|pjz#rQ8ivRPD{{1vmLsj=Rjq=Ae$OT5&v*Q=ydSIXA3E=Tc*dfP*N1kEq z4D5#EvoRLFFk}5(kUREZe*}O>;(la)*glhauY@NYSjCyJ}6<-SR{WDfX6o=}0yp7}Q%xrVwL$7td(_eD=K?xkkkw2YE?THU#vR_cc8EdKG^ zq%3QxsmGvfK8vliP<6Rz^PQENtEx8}N`_M=z42p+t&j2T%_S6bsC?CcPE!{B6iBLh zY!SEKfs`wcR(c~jzH7NwzWZ!)oejO*u4^#4>J9W%y^Rf8B1>|@CFqTP(F*%QFqQma3g4;kQ$d@a@&bGzI2OE-eh=?!!4Kt+>vx$y6dJU@w8|gZ z@q4cax%43P+?p7d!2PZRIl>qENkL!$_+bIY7cdt2;thD?N9>D{0Y>d|8|q-csCYrW zk8{K<{{27sfzciLWox>#U-OM8ULVAKv&)BhA-p#quY~@XGY0bO$9(ws>cRYej^6CI z7}N~)d0dTp8hmg`_&4&R?nYL`Ut~nyniH?ntB6q*>MpJ~UwW8xrk~ z6guwfqAo2$sNzq*j~+Zq)z$5#b=*<)w6|qfrMd!ba){s4(8Zr`Pwu0x$Y6P$Ka?KR zNDBZ}`ip)bm|nzH^O*fG*jAq>sov>(j~a(iXk9HbEQGY`UaA2t{Je@cytGdoMp|=# zk{iRR)`5o??}WdDKMF%H5%j}21N2=%;g3?-KZ2D4_BtU?cEGQVRS4{xGP&_!0N?_Y+LB{Lv5khcD!f%-|2H+f)Y6k3{GnKH!^? z&{JCBcrL6n5@S#7ht=2*KG-kv{HYS4Uw=a!`sE%$ANgRZqcMW{VQdecX8?bUgx<1c z03Y{8-5(`y^wR#WAx^y!@?S>|;(77jT*oJ@IoJ2OP?e9HN+FIspOJaM9{H)CF|REz zYt-Y(+ZTl%l@oF4&uFl!`d6Ukb?O=F?a`|2aZ*LT%7^4mBt5*V*$!F~W3O8OY^D*1 zuC5p`a4l(7wH*#yNsGqZG|5VO+xPC_6nivKRqP8Py)5{0vn^g+`gpX|+HNgS_4|L$ zdZXhtrO!-Wzp}}YAX3iPLdSBOY5J+hSx4?lj_1VYQ&<1seEPHIEg(~g2a1n&g! zbuWtjA!v_%k{Qn%!Bp~xjW6yqj~E0igFj}1F9h@L<0ju;mPZ#2VxsfYa|@jLQdid)6M|I2gu zciRYc0UwRv&(WeMAFtubV?#UWCE$@H)E6y`IzNg>4BB6d{)|cJKd@&2pXb$;<(h#l zSRUC_jpMl5mFCwE=0}`YHpG)>K)iVwV@WN%yk;@uFg?s7C*sea)7nJMaPTLdEm2kGL-%gv;Ug{lyTwN$;jxopMzR%*4vZR2l0TdCusGmZ;qhm$UkEf1kx zTHE^=GcBGt@TKoK+J@4d%)y_1^&*r8-yFF5hF=(| zW`63Eh0_|X+?Z z4}1T-cYU5u<#=j;XR?p*wf_(N=7wo}Xaw3zkZ>+z5AI{x!M#q06^^mi-q z7~)SI3;jV*;w*+kUKA92#P|TC{`my0Yk2$!@}b~$=o8l=M+$zZKCvW$*B872eZmWR zgWz?HMwPT zWeMQxTM_%l9r~%G#Etv&>q)CwpU5?j_5Q!0PpCXqUB0f-i(hwd&-tXK+p^rR>TVgd zKmXp0h8{exSARb5Pt+5ggu2D<&@&9(cwWAitXD)Xzsz1KRDl0_6f@l zSNmeobfvXncJCxBrPmAgDfRna=ew44RNZ*&(zv*5iQ` z@l~#~Ln{T*^(}i;W#S-mKAlZf1T|B)heHNf98}%$25w*a9tkGRn%RD;g?v@D@fkHk zDAV#DS1aC9^~aGHB6)9Hyb?xrI~hitj~h<=RF;1BA0kLE3pyt9U&m+ogEs^fZ~UeG zZB9o0{tU-`uY!FbI1>9L3HwJ-c;jd=-&a8+{PrpHRweID{HWq56}0T(bt-``1h-?{ zW*@&U_<=lfG4lice8DtpmzaqAQwqw7*%X`?(x>rw@Da`#W~vx_YJB@VrU=_;`z6%oqC*zp<+? zANTIV^K$oQzESZTIvgBC?8GA z2&$;Yir2$!d8GL#7FMPQGqL+Z^&!#wDF>2nHcbkoEK@SOO)ei~>x1yT(Yoe)yC?I9 z)}`RgFTwQ158s`%OtX-l9{Y{*hv<_sk3;06NALDBK`F#XGlt&hUM}TRTN9J4kI$H4iR{idEf5YXud|c|q zI>Y{15c{D7cx5Wa_X2qS1?-P-jIToexP|@YgYinpIS0XC>V5S0Z8nL2_sda-=N9A* zV@KwP?jGO;@P@K~>O6F;E5E;~52JH$KK`&Lk6U2AV|VtKG-(AnL;9kWMLudlj^BvM zfqqHp*xot3q_MM}SXLZ$%5x*X^(E|yH_^}hBynGa{WMN5ps|y@M>a2-uisXh`mI)3 zV;wy!kYAO`h^9e$Kvy)S*Aqn}XoBj75}PTMMm_XUDeu9yy1#{*YMElnpL#aBpf^UaP{Q&qhwLVXkluc!YA88y zqg7j;D*J)M>GHv-Yla=+;Kc|!Iqc}mJ*VEyCkE|r5$>ztWsF_%{17}FVAQ4?La}cI zg*O5re+WukaeN4$Czx73Sp+#m&;|0xP4I`{1n@}V1N=UMAIK-(B{|MZFwOEw%V@ri znejVU_3QQL>DoDb+yU~5J^ZqU06y*n9$A2KI^cfnm$i^Xq;7)n%4zV21LA|+z$?!4 z`F#zOnJ*@we%DF(Lz;R)&S=d%(HgvA?7_#+c{5h*&F4QCp6S8!)A!-y8@sdL^-e49 zlhOx$)vh@i^|oqBg%H1y)2O;1C`X7L>THzYeh4FrBMv<`;yPZT?)WWIN{DK`e2RR> z1gQ$(2k2hDhN`v7F8Zma)`58|>UgcAsO4j056y|D_Hg@B{h8#vkhtKPs5k;?Td0;p^B6&y|Gv zEGMj=&A-2F(0kkv7uF5?#V3HzJDB1*bQ1d}9oF*&f4G8QJi#9c;HMkU1|t2=+r&%r15i0io3gO8u-&hszA?lQ6m+dX1iqYwKx zET=p!XFQ>*epSv~1aYOgSPvRwXXHL-c1FD%`Jv*71IuI7+O3S-_8NKbx2d=4-d5qs zDbnN1{y0E-8>hv)$YXe^rGNQtqb>7CtCAP%>Cd|Pk8d~_O}hJZL^Q2%yiv;#538?-V!qgDL5 z&Q~FH>Byd({qKiRskOg*%?wlVqcO`)kM&YKlBrIN>GN>1KQy9X-;i+9+rZV2peDsO z`ETnPNxE63-{0`b9o&CGckGLR6zGNLOfBd&g2E?KmtHUiJaRmY@1J05`Q#P&MDP&g zk{#d=K`-#gs}%T=d=lU2k^XrF4QY)}egmJp!tdX-fS>36XYup?I`ke3;?|`e(=F_y z96-^Vw_zX6#Cp>I!YFzT_(S4WUBN3pc;8V|@c!TPK8i=*l&PJWAIgDGjJ+VINZdMJ z-`4}@p&!C7*j3g*UP<4R<4c=+Kt8Dny`?>@qy?)f%48L@idh3f0+OYd|W?p%we z-^-q;Rb*Z?U8s4Yc=cTo6#uzq7zrir`xH$58=mX_XWd|0o=XeVQ1vtP#JOaCJKO(WZ*wS(?zG*zx=$GOD48Sb zW_}e%iaHuARP?-k`((T8%&`62Dol(@bv#Es1OKH1%J}2`_4iiT>1%yj@eAUw?}IOD z9L}E-U;Ce6Fy<{A`Tlvnm}eRF{&^2`lsm!Tr!9?uGMuf~nOfHo+bt7zn$>dgu*;qBp!e#P1{c zH}y3nKVIw-{+C#f5KMD=#JpI(erC}J7P38{;7mSlg}q=4>=MfFX3+ln1Q<0i?D1z@ z1iQmptWywp4tjzY?3(wmj?|mALO;2;fZsRYRMr=2jAdNjpVv2aVS7UNme3~q+a@ zdh)Mm(z?Ir4UDF>HNQ})r4hD#G0kh*XSnJ+Frn7`vW_R1P@{?d1ivb6OX_jg8N~|OaNa9 zUITCBLwu27YWc$#_6I@X2`hL)u+1-g-(RG_59ANS*vGt%U|O?7v|r2LU-%%;BK}=v zn!(4D7Bg;v{lFVMQ40G-{2>|Ow@m<_7+`ns249FBb29is`gwX|{n-ooeU4A%_jf@b z=q0GX)7qO~k8i=R*A+g1KlWhsmB_koJ?%GpUFuXRhJM>sdhe>E(Ud2$!KE>IBI(Ah9)lNc38m8wf+NQ+ zQM|Fry+ipy!PG%7ARI)I8vTq{eOxQ;ZqxI(AUfRWZlNyUn<-t-)~fYEavjf0^|i<@ z3k7ItY{Nrn#xl)-5lZQ`vLAOsscNY^(T{#s`K%>03r09;){kt@ls)sn#yss7N6_8y z^CN%W7)d|P)f!p<4SzVnPARwmV*}(3LE(?1@Pi2Mf!uKsd?Ao8<#O z5T2-r{gMTI5eI*V0pn$mHypqheZV6vu)YuYq8{?!jNq-PY52X5XFICou^J#})M&w< z%ha9G(23`Tw}(7|d1t&?9?1wC*^TXy-fiIzA%6YZm#91L$bFpeNS&Sxh~IckTJ)F0 zi(+gVFhZ0k!s$0 z?fU)XIvRC~Yn8?#$f~uIy&XYk=jEB-HE$H{cGD_|S>NwpeS-Th*ay6EE(HpI9E7|f zcmZ;U^w$y;o~UZ&`z!c&@uSb7rwA(EursK?xWhKoh9(#03&X4w_<_9Py_tDKFwOCX z*LtINJq?n_UULcm4!x)Gaq+ufhrF>x;wi8{PJ$XAUbcYO z?}NV1`@mC)Kk)lkL?42a$Ulv2&V2GG{2ik^^ZEXs%opjq^4J-E)q|bEJ6&1M=-Q6; zo=HZI`|SD^%NIWuM;)HLh$sIHa>Z-JZxx2WH9z7u?ASl@&3)AC`jzf?*E%2_MqlRL z^vz(cBT76q(}`jOEz#N}-;ALj>-+JO>sCcd-*A8U0M$x5i zTEpiQd~tmqU)P3yZ=>Y3h#h1P_Jg4KOC=suPa{X^=oi4xQ^5En? zNu5`9SD{}j?p#&#SzFygmz=al1Ic|F%V}*1mWR^I(XHoa=^jSjH~#ppe)(|9p%(xP zC##kQzbb-ASzS&wjU+9~^}C0Ws28GmWsi#r*|4VCpX$yR*pLn=s7{jm9buB ztXBo=Rl#~yuwE6cR|V@;!FpA&UKOlY1?yG8dM;Sc1?#zBJr}I!g7sXmo(tASdeyLAHLO<+>s7;g)v#VQtXB=|Rl|DKuwFH+R~_qB$9mPVUUjTj9qU!cdeyOB zb*xt%>s7~kZdlI^>$zb)H>~G|_1v(Y8`g8fdTv(#(|HLzX{ ztXBi;)xdf+c)iNbydM8{m8dF&pzrek{`21}@V~nPsf{~mjJJ{PC*zNQ9>xE{d)oh# z@+<%Eu4XFF7JsUjpHV-f)cm(GrGBx@%&!ZI9%#l%M2Lzz&H_l@ghHd zze(`ZctgJ^fd20zyKz5(aZOq85kJW2F6`&%(T3m0AN|v6AP?nM2iQHkaGd7#c8r;8 zqd$Q&pJyz={acUagZ}Ut`ai$mdY{<}AYa}Aeil3UWAD?(IT=+7$vHasgI4)Bi6T#T zQJHCb&_^Jih78?ax269EVl!t9rLWs|!#;qx}jloU}dFQGY+YGo+k=B zE?W~!%J`sCJ%UIjxvEox==D###*VgA^+J)C<~*UJa=8Z6<~Z#GSSVI2ANZ?<9L9a7 zS|_N!XLbB0ts-y+}ilm%-H>nj< z*e?>nFBul_^%6|UFTy7&(PZKC1joQ`u^RpoLE)D_z%PQ*=d)rAuOs+({33kP5&l-7 zDI@cXX&m@O@Dz^?u1EQO1V7LoDt3$bNah5$CHcVcnKVL z!tqKvdOLrmV^S^`9$%@3*K6T*!f|&TZ>D2XuGTo-R>ym>9q@W*9jBM?s$){F?m9ZW z?4hIgm0mhp4)oFSUTi-d-B$YOXqr4w$LZyV=$Mpin2v^B8~FO20MAGt)*#%!ukjq} zvy6ZDS-~?EVZU$!1}@=sR>OXAeFcvlaXb+AQOhDe{w&2lH6c&VbuPbdng_eXOde;R zz?cX-M$PU#zSso6^CmXC1?R7$e#Zz8<{=+1{`|SoFJVu6J|2#|l@p%4UVLrDeLJI{ zK}p1m=STn6jL18ALAxG1dClAToH+i&0sh-`=u7_qd1mL~mrWwQe0ky?^jFz|zKe16 za!HL&`Eth+r=`YFUEL8IMO#K|m4U;k)XIEYu2)d~o6En9%=|{>o8;WNW==_!XWw{X ztEAgOwD-~E(nVXCX?{g5(A-SD^+J8Y)Wtpe+|ho?b;{FwYaLirp992oH1DKU{0O53 z8&*o*Ngz zLxS((p*r9l!6b~mA>Ro4K+f@x;QJ*gd^9tfj|)ov)Jg1T!PN4R6XHAsComrw?7=&N z{W0F30za0IIz89FM?5gi@X^Fpz)zl#pS&>7fc>p_OV86U3VCQc^6++`K55HY%rmo4 zSEwoaB%JFGd8jG-ad$W7d9AxJzXW?SfAwt*IjAF#H@0X0%y{^HT)<-s(bv9Z9sVAs z%Elf0RpE_TllQrSt11tBx)^D4^TQ?VnGHqdPaX0_(99c1vrh5JS6%HVG>} zB+8cmI^BI#j&E99Na_^__0kgaRsWU8C*#|GZl>BlB`yr?ZzjDBOqO8yje=?HiZ{0! z%(PH5r+V?lJwm8m-OELr6bq#-nmg!PD4qF5OT-DY)mMHva`Lf>`Bgq#`>AdxcdIy( zeY;f0nMf*MB>b;o%c7{Fp1^Ch^^3JYk@a0}XL#SA74wYI@1GZkc`f4Kx0JssgGTrw_}T^;f4tqm`xxH|zfa0;klO3fG{yCh@yB0}d#Ct%#Q*ntsB4?v z&U)@Y`Fb=-u|H(|@%v-&CB7ckPilWi{zf8*A!!4DN__o)f+Wn#ocR8Er!miA-}~oX z#5||{@1OTiALJ7|kXqi6eT1)VknzXwE5~bm)t8wBlQ$R?iaisFwN*G?i=}g&%=Fr3O%I(^pxVbe?C^8 zw*#0*>}$|Rl92!Y(?p(UfZic?lRnT#{9w10el8M!YJxt~5O$W~b9o&n*i%YDA1Q`D zPL0uzrTRqHGu{kgJtL$C-oF{gpMTMqkGJp4c`uvV@#}uA_`GXkKk?-Fj|uHqpHcDf zx_&gZ1KUxu*JXRj*Hw6(nWd0#U%<$9|1zUa;VaV5>-CJ5)Scz?$M^ZTE|70}A2({y?>T?lJyCW%u$UZJOIMW^(WmFMoZ$YY^Df81C7*K^7B zf3hx%*DNFyyd>gr{W{6@I;M=;IQ?!ya$Qr;{Z&+_2Nf5xFHu!Gjij?~-M*~#nU%Dv zW~CRTu%}!FFB!o{g1!s+`X^vN2)>J-GJ}T%12J}o{X|fBXk`?yFIXRXwe&?66n?6; zmgfup9Y38o%si#gp#7PS=xF+yc}P(FM5fm`4m7-@?{&(&JhA5*6PcF;(+n@205A0g z4-Ewm#R2={z9!;%)D8DN2E63Df`7+@rt0^{5GQt5*g<@-|C|=`>r$_I5ys*-Dhxi_ zGN0EQ2HvSKmw9FDEIz*v;?%~^Wc$dXiHu(8bF2Ec=yr>B%@CK`5d7oCc9aoq`S^&I zu;+R*-}Gz8ar2IC@&4`E|1`WVpKk^qWhreOr8T1LR1p1KKIeR$2Cr#|-l+B|J)DH^FFl*QHF3Nv2;`i8^ zUqsSGKh0pQ>?L}jS`@w6u%Jkx2`Tu@Jpg8jb_cqpxn?={z1fgZam2 z5`TWTVLUF}i{+kzEfA02kY6`?G5WP+|MvM7EJyv^4)Ro6#vyHZeX}RmP3%^W{rDd9VcYhw9zdQriGg|S8iWht2({|iN($jSoAEmr{ z0mOZzRS)_x^DatUwt0!)7n>+yQNw_xE!Wa@RXuI#`B+j-bL;gTWF_4leLkGB+P6P{ zqeckn?Rf1&$a!_`x_@N|rdEwsnab=mQ~sXS5=)f~ru-H9IcM7yOd@v;zk6xHm1fCy z-Wp3UhsN9urE97R+Wu`}q`S?hh12M$4$JQ~i9mg}NQx^u6i|)cs?*X|=Q2gRou^$ATuumjkRPbH=6%W1=90=ZW#Qqc%{yL2P zD=7RWa+Kf$@K*Rb-Y0^8$6pEHFQARTMlesA1kZrCj3=-kfgi|Uh8nNuRIPQk~76p0h4saXft-{zRp};l$Jn^V+ zpU{%Wrsga+$JnfuT^ptH zk?xmM5|f3Jw6vp_N-xywMTU^%3#=~H$>n*kF!UD=r?m@rw0v+QoGN9g;ndX}LFF}Z z&oh!NUtjGNXBS0Qx0TON6t>dM-7QqBjAZ^YXn&H2^a}hXn1Fo{zlg8Dpzzdi*aZY7 zZ%N{P0c$?&{!PTbj9bi~;~b9rqkc?@IlP`1Rm+8Q-uar8~%^**Kr@DAB@APhs`uNZr7uM(FdtH%#T^4m9 ziy)8M9`*1IkiQ8C@XhflO!OGOJ!@a?o+^%c%A@cg0H5mb2nm-~JE zLTIqcZv3`u2%YXRXUueWRsTgN9y6_4Wj}X8yI^u|Fz>ftAE`Q#_xfC!zRp5(vT6lo zRei|YUE23;6-wR5?9BCrQy7i&DpY>Y^Dt`ur?(1SRQXM+O60An5%lz%JfFR&7m0e7 zk#x6Bh1}yqqUc;X&B1P^&0T*DIG8b-4BDUA*G;&eKJ)qh3VL8a490!8A@+k8o@0V3 z`D!xwNw6MxDhr;cg2Gpaux|y0uY%U`=MwA%p868M1Hr%HE5n*2ypAC3=ceu8C7_|q zNuD<^1%9BMW!myuU&n@NMNcfag|FWX@RbYhQv~>`A^0lP%JZ%Q8{_%r3c1QwulOhC zlW(BEmHUCu8;<9x5&PB=`#1^vRqX3car_zTPL-LjhwQAa-U| zCrt>p)seZ~U}mCqX0R<@XthRPe^?NA6oNeRP;#uTWuAvoFuj~#qKb-B9Mw3}&A-BJ z^>q6V&lK%4C6Y3S7gm+FqHK1piI?{(P1Z_FEf)%&R(*_8>pyL5;6L%7%J}2`r%4ak zf9n59`Q`t8*JJS0|K#h@G{ydq@yG9v9IyC##D7xzgX1vRA8^HTde zlQ<3G9T|W8zDhz)`KJ}Pk=pCgJjL%!#vgw@OrPIX4@?1RRK0&d!L66ez!`Uuc;h{vXCr}X)PA8J=I zO@71Y3#Ju4B@uc`JoJ=oxIg2er`(R^`%vA=*dJIBdaVKetC7IY{;VgwhF{!%JfG)< zecA*2HgPT=w?I!3|CP$0GH8Du&{OWrWxXdq>?XaSuf#(?*^a)0^I=aJID*fQ?*n^E zE4Jf~tIzgeRsThg*E!i5$D44xPv_Qb_x+(2>nq_cQ8(8EeZSlp(^X?X{R=0fu2&2x zhB%+>#zK0*Ts!1hJz_KP#GFRGeD0i#MlGws(fWk4t}~C`TkbLqch@?49VgaD4$#Vg zk{e%5OQ4rt19rcO*-T&Rj;?jo>*>|5Kc8Gn&#JvXGcQjx?W_49Q~JXZq;EwDAiYR>G_4z zrl@`!r8I|o7#+%}H9l5-d#v-Cb$T97UoSYIGEP+;D%By^I2@At_7#gn_e!3JtYzMIWzQA!BUI({+Pf|g2GSX;3vU<;3wTL_E{9i z<0w?SMA>hEldvxZmx7r;HS%o4{8g3 zI_}TwDgTn*e{|S59(!ZIx+4D2aSk7Ei2XVi{mYcy#Gw66m*?(0=8Nl-8J%F~b%VX+ z&xuBDy1@$n)Xfq6zODPReRf(a=AFa!P%pC`ujA#x{I#qRpPvJ~_`jnQ)_6| zFWJ;fWqPW5y}DsX@m2X&U44pMNXuFu^P7bZ-q8}GR6XABD~`K9QRPqFJYVc?n=mS; zrGYtx(@&~9(Y7byq=yG>N$#KA(@u*-iXttlX!O-&ev&@W2JKJeCwU$w(seTsNb24nF^2IxOL@w+4LaaBX+xh!tHey<9~ak|D=+^Bp}6Bh5vfxao};dg&P?&}}EKJJ^{ z*hg=`oY}~In4i%V&CTO?h4yP{;1^Diz2&RSI%M_2Mh3)=u@Y@+ic+Fn^RU&EY*LzL3 z9w*l$uG;;`@R3!lR6ZnkQa#02iASeZvJ6VWSMl@semLU!BPe+yfsnHWZ9E3vG6JoT zznoX`dV=rbv3%ew!9_BL+?5X~JT@8oUQl=}62A+<5#X=dh)WTC2anl!>kjw~XfU1R z@6nfe%OH3P zJojqB|2`AA2l7~B)GrGKCIm1a8Bq79jEToh$MAjg2fsNXo^i(vKCbMwx;*r34)fxy znam@~-lcAv!Qciz{9@GO@Sn(h`W<*J58h{siO+j7ieKN@pXH}M==+qV9^x76GXIrr z10HO~JlDdV<+$MHESI&x@f_bVzlB22?AM6j=b9VhA1iVl_WQ*VN0bwNXftsAGp7g0 zL&*lYEHmo0Kc|HbnuGl+Wm4Vj^Ugj=4fMk2zmQ@6_?KS~-h;a5TM>u7o~q~xaO-J8 zT}!LjOfj^z$fdOx>qn8#tz4=Sdl*HPzB+03p-|LCS9+u_2dX|V_ua3}Ckt(;zdT!& zeHO}kcWBPvW~hGEM-DhG=@?2{7u)hx!>F~=Ry*ZPu8%BwrG2rStxL>^q;KrLo1CqA zl&x;rtC0IDv&Krw7^WIGM57*BGT-TTss_kAX_1WjeXR)oTblweE->n^doF$(g$J90 z_XOi1ztvmK^9A3N!GD6ngL(0L5ELGq0p1hz0RNT7?@RD+c+k-IIDc-S zsR-+#rb^6rHe3n*%l0eJ7yLjTyf}(^P%zE#;CbAiqrg4D=eS>8fQdN19OG6S_e8hR1QxQzwJ3s*Ps8%7#3ZxEVYiF`ehfL*APP zK2&jA^BcZ|-de&~T@Mt`X}lM_r%Cl0 z8Q>?l&;EgQ_K08l9C^dfX}OlZGUh5NBfF{$f0F7dsc^{5U&zu)CoeY2Eb( ndarray +# Table shape: (N, 2) with columns [Energy (MeV), μen/ρ (cm^2/g)] +_MUEN_TABLES = { + ("nist126", "air"): _NIST126_AIR, +} + + +def mass_energy_absorption_coefficient( + material: str, data_source: str = "nist126" +) -> Tabulated1D: + """Return the mass energy-absorption coefficient as a function of energy. + + The mass energy-absorption coefficient, :math:`\mu_\text{en}/\rho`, is + defined as the fraction of incident photon energy absorbed in a material per + unit mass less the energy carried away by scattered photons. It is obtained + from `NIST Standard Reference Database 126 + `_: X-Ray Mass Attenuation Coefficients. + + Parameters + ---------- + material : {'air'} + Material compound for which to load coefficients. + data_source : {'nist126'} + Source library. + + Returns + ------- + Tabulated1D + Mass energy-absorption coefficient [cm^2/g] as a function of photon + energy [eV], using log-log interpolation. + + """ + cv.check_value("material", material, {"air"}) + cv.check_value("data_source", data_source, {"nist126"}) + + key = (data_source, material) + if key not in _MUEN_TABLES: + available = sorted({m for (ds, m) in _MUEN_TABLES.keys() if ds == data_source}) + raise ValueError( + f"No mass energy-absorption data for '{material}' in data source " + f"'{data_source}'. Available materials: {available}" + ) + + data = _MUEN_TABLES[key] + energy = data[:, 0].copy() * EV_PER_MEV # MeV -> eV + mu_en_coeffs = data[:, 1].copy() + return Tabulated1D(energy, mu_en_coeffs, + breakpoints=[len(energy)], interpolation=[5]) + + +# Used in mass_attenuation_coefficient function as a cache. +# Maps atomic number Z (int) -> Tabulated1D of (mu/rho) [cm^2/g] vs E [eV] +_MASS_ATTENUATION: dict[int, object] = {} + + +def mass_attenuation_coefficient(element): + """Return the photon mass attenuation coefficient as a function of energy. + + The mass energy-absorption coefficient, :math:`\mu_\text{en}/\rho`, is + defined as the fraction of incident photon energy absorbed in a material per + unit mass. Values for each element are obtained from `NIST Standard + Reference Database 8 `_: XCOM Photon Cross + Sections Database. + + Parameters + ---------- + element : str or int + Element symbol (e.g., 'Fe') or atomic number (e.g., 26). + + Returns + ------- + Tabulated1D + Mass attenuation coefficient [cm^2/g] as a function of photon energy + [eV], using log-log interpolation. + + """ + if not _MASS_ATTENUATION: + data_file = Path(__file__).with_name('mass_attenuation.h5') + with h5py.File(data_file, 'r') as f: + for key, dataset in f.items(): + energies, mu_rho = dataset[()] # shape (2, N) + _MASS_ATTENUATION[int(key)] = Tabulated1D( + energies, mu_rho, + breakpoints=[len(energies)], + interpolation=[5] # log-log + ) + + # Resolve element argument to atomic number + if isinstance(element, str): + if element not in ATOMIC_NUMBER: + raise ValueError(f"'{element}' is not a recognized element symbol") + Z = ATOMIC_NUMBER[element] + else: + Z = int(element) + + if Z not in _MASS_ATTENUATION: + raise ValueError(f"No mass attenuation data available for Z={Z}") + + return _MASS_ATTENUATION[Z] diff --git a/openmc/material.py b/openmc/material.py index 735a05743..239bc01f7 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -2,6 +2,7 @@ from __future__ import annotations from collections import defaultdict, namedtuple, Counter from collections.abc import Iterable from copy import deepcopy +from functools import reduce from numbers import Real from pathlib import Path import re @@ -22,8 +23,10 @@ from .mixin import IDManagerMixin from .utility_funcs import input_path from . import waste from openmc.checkvalue import PathLike -from openmc.stats import Univariate, Discrete, Mixture -from openmc.data.data import _get_element_symbol +from openmc.stats import Univariate, Discrete, Mixture, Tabular +from openmc.data.data import _get_element_symbol, JOULE_PER_EV +from openmc.data.function import Tabulated1D +from openmc.data import mass_energy_absorption_coefficient, dose_coefficients # Units for density supported by OpenMC @@ -409,6 +412,190 @@ class Material(IDManagerMixin): return combined + def get_photon_contact_dose_rate( + self, + dose_quantity: str = "absorbed-air", + build_up: float = 2.0, + by_nuclide: bool = False + ) -> float | dict[str, float]: + """Compute the photon contact dose rate (CDR) produced by radioactive decay + of the material. + + The contact dose rate is calculated from decay photon energy spectra for + each nuclide in the material, combined with photon mass attenuation data + for the material and the appropriate response function for the dose quantity. + A slab-geometry approximation and a photon build-up factor are used. + + Absorbed-air dose: + The approach follows the FISPACT-II manual (UKAEA-CCFE-RE(21)02 - May 2021). + Appendix C.7.1. + This method integrates over the photon energy: + + (B/2) * (mu_en_air(E) / mu_material(E)) * E * S(E) + + Effective dose: + The approach uses ICRP-116 effective dose coefficients to convert the photon + fluence due to decay photons to effective dose. + This method integrates over the photon energy: + + (B/2) * (h_e(E) / mu_material(E)) * S(E) + + where: + - mu_en_air(E) is the air mass energy-absorption coefficient, + - mu_material(E) is the photon mass attenuation coefficient of the material, + - S(E) is the photon emission spectrum per atom, + - h_e(E) is the ICRP-116 effective dose coefficient, + - B is the build-up factor, + - E is the photon energy. + + Parameters + ---------- + dose_quantity : {'absorbed-air', 'effective'}, optional + Specifies the dose quantity to be calculated. + The only supported options are 'absorbed-air' which implements the methodology + from FISPACT-II, and 'effective' which uses ICRP-116 effective dose coefficients. + build_up : float, optional. The default value is 2.0 as suggested in the FISPACT-II + manual. + by_nuclide : bool, optional + Specifies if the cdr should be returned for the material as a + whole or per nuclide. Default is False. + + Limitations + ---------- + This method does not implement correction from Bremsstrahlung particles which can be + relevant at close distances. + In addition, it computes the gamma contact dose rate only for the unstable nuclides + for which the radiation source specification is present in the chain file. + + Returns + ------- + cdr : float or dict[str, float] + Contact Dose Rate due to decay photons. + 'absorbed-air': returns the absorbed dose in air [Gy/hr]. + 'effective': returns the effective dose [Sv/hr]. + """ + + cv.check_type("by_nuclide", by_nuclide, bool) + cv.check_type("dose_quantity", dose_quantity, str) + cv.check_value("dose_quantity", dose_quantity, {'absorbed-air', 'effective'}) + cv.check_type("build_up", build_up, Real) + cv.check_greater_than("build_up", build_up, 0.0) + + nuc_densities = self.get_nuclide_atom_densities() + if not nuc_densities: + raise ValueError("Material has no nuclides; cannot compute mass attenuation") + + # Collect partial mass densities ρ_i [g/cm³] and elemental mass + # attenuation coefficients µ_i/ρ_i [cm²/g] per nuclide + nuc_attenuation = [] + for nuc, atom_density_bcm in nuc_densities.items(): + Z = openmc.data.zam(nuc)[0] + mu_over_rho = openmc.data.mass_attenuation_coefficient(Z) + rho_i = ( + atom_density_bcm * 1.0e24 + * openmc.data.atomic_mass(nuc) / openmc.data.AVOGADRO + ) + nuc_attenuation.append((rho_i, mu_over_rho)) + + # Build union energy grid across all nuclides + mu_e_vals = reduce(np.union1d, [t.x for _, t in nuc_attenuation]) + + # Build the material linear attenuation coefficient µ_material(E) [cm⁻¹] + # as the sum of ρ_i * (µ_i/ρ_i)(E) over all nuclides + mu_material_vals = np.zeros(len(mu_e_vals)) + for rho_i, mu_over_rho in nuc_attenuation: + mu_material_vals += rho_i * mu_over_rho(mu_e_vals) + mu_material = Tabulated1D( + mu_e_vals, mu_material_vals, breakpoints=[len(mu_e_vals)], interpolation=[5]) + + # CDR computation + cdr = {} + + geometry_factor_slab = 0.5 + + # ancillary conversion factors for clarity + seconds_per_hour = 3600.0 + grams_per_kg = 1000.0 + sv_per_psv = 1e-12 + + if dose_quantity == 'absorbed-air': + # mu_en/rho for air [cm²/g] as a function of energy [eV] + response_f = mass_energy_absorption_coefficient("air", data_source="nist126") + + # Factor to convert [eV cm²/(b g s)] to [Gy/h] + multiplier = (build_up * geometry_factor_slab * seconds_per_hour + * grams_per_kg * 1e24 * JOULE_PER_EV) + + elif dose_quantity == 'effective': + # effective dose as a function of photon fluence [pSv cm²] + response_f_x, response_f_y = dose_coefficients( + "photon", geometry='AP', data_source='icrp116') + response_f = Tabulated1D(response_f_x, response_f_y, breakpoints=[ + len(response_f_x)], interpolation=[5]) + + # Convert [pSv cm²/(b-s)] to [Sv/h] + multiplier = (build_up * geometry_factor_slab * seconds_per_hour + * sv_per_psv * 1e24) + + for nuc, nuc_atoms_per_bcm in self.get_nuclide_atom_densities().items(): + photon_source_per_atom = openmc.data.decay_photon_energy(nuc) + + # nuclides with no contribution + if photon_source_per_atom is None or nuc_atoms_per_bcm <= 0.0: + cdr[nuc] = 0.0 + continue + + if not isinstance(photon_source_per_atom, (Discrete, Tabular)): + raise ValueError( + f"Unknown decay photon energy data type for nuclide {nuc}" + f"value returned: {type(photon_source_per_atom)}" + ) + + e_vals = photon_source_per_atom.x + p_vals = photon_source_per_atom.p + + # Construct list of energies from (photon source, response function, + # mu_en_air) for clipping to common energy range + e_lists = [e_vals, response_f.x, mu_e_vals] + + # clip distributions for values outside the tabulated values + left_bound = max(a.min() for a in e_lists) + right_bound = min(a.max() for a in e_lists) + + mask = (e_vals >= left_bound) & (e_vals <= right_bound) + e_vals = e_vals[mask] + p_vals = p_vals[mask] + + if isinstance(photon_source_per_atom, Tabular): + # limit the computation to the tabulated mu_en_air range + e_union = reduce(np.union1d, e_lists) + e_union = e_union[(e_union >= left_bound) & (e_union <= right_bound)] + if len(e_union) < 2: + raise ValueError("Not enough overlapping energy points to compute CDR") + + # Histogram interpolation: each new point inherits the value of + # the nearest original point to its left + p_vals = p_vals[np.searchsorted(e_vals, e_union, side='right') - 1] + e_vals = e_union + + mu_vals = mu_material(e_vals) + if dose_quantity == 'absorbed-air': + # Compute (µ_en_air(E) / µ_material(E)) * E * S(E) + integrand = (response_f(e_vals) / mu_vals) * p_vals * e_vals + elif dose_quantity == 'effective': + # Compute (h_e(E) / µ_material(E)) * S(E) + integrand = (response_f(e_vals) / mu_vals) * p_vals + + if isinstance(photon_source_per_atom, Discrete): + cdr_nuc = np.sum(integrand) + elif isinstance(photon_source_per_atom, Tabular): + cdr_nuc = np.trapezoid(integrand, e_vals) + + # Compute air-absorbed dose [Gy/h] or effective dose [Sv/h] + cdr[nuc] = float(cdr_nuc * nuc_atoms_per_bcm * multiplier) + + return cdr if by_nuclide else sum(cdr.values()) + @classmethod def from_hdf5(cls, group: h5py.Group) -> Material: """Create material from HDF5 group diff --git a/pyproject.toml b/pyproject.toml index b2b5ff5e4..bf4011344 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -69,7 +69,7 @@ include = ['openmc*'] exclude = ['tests*'] [tool.setuptools.package-data] -"openmc.data.effective_dose" = ["**/*.txt"] +"openmc.data.dose" = ["**/*.txt"] "openmc.data" = ["*.txt", "*.DAT", "*.json", "*.h5"] "openmc.lib" = ["libopenmc.dylib", "libopenmc.so"] diff --git a/tests/unit_tests/test_data_mass_attenuation.py b/tests/unit_tests/test_data_mass_attenuation.py new file mode 100644 index 000000000..0fe12a7db --- /dev/null +++ b/tests/unit_tests/test_data_mass_attenuation.py @@ -0,0 +1,53 @@ +from pytest import approx, raises + +from openmc.data import mass_energy_absorption_coefficient, mass_attenuation_coefficient +from openmc.data.function import Tabulated1D + + +def test_mass_attenuation_type(): + mu = mass_attenuation_coefficient(26) # Fe + assert isinstance(mu, Tabulated1D) + + +def test_mass_attenuation_spot_values(): + # Spot checks for Fe (Z=26) against NIST data: first/last tabulated points + # and a mid-range value at 1 MeV + mu = mass_attenuation_coefficient(26) + assert mu(1e3) == approx(9085.0) + assert mu(1e6) == approx(0.05995) + assert mu(2e7) == approx(0.03224) + + +def test_mass_attenuation_caching(): + # Repeated calls with the same Z should return the identical object + mu1 = mass_attenuation_coefficient(26) + mu2 = mass_attenuation_coefficient('Fe') + assert mu1 is mu2 + + +def test_mass_attenuation_invalid_z(): + with raises(ValueError, match="Z=0"): + mass_attenuation_coefficient(0) + with raises(ValueError, match="Z=200"): + mass_attenuation_coefficient(200) + + +def test_mass_energy_absorption_type(): + # Spot checks on values from NIST tables + mu_en = mass_energy_absorption_coefficient("air") + assert isinstance(mu_en, Tabulated1D) + + +def test_mass_energy_absorption_spot_values(): + mu_en = mass_energy_absorption_coefficient("air") + assert mu_en(1e3) == approx(3.599e3) + assert mu_en(10.e3) == approx(4.742) + assert mu_en(2e7) == approx(1.311e-2) + + +def test_mass_energy_absorption_invalid(): + # Invalid material/data_source should raise an exception + with raises(ValueError): + mass_energy_absorption_coefficient("pasta") + with raises(ValueError): + mass_energy_absorption_coefficient("air", data_source="nist000") diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 0b1b5fce2..58cd4d563 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -826,3 +826,58 @@ def test_material_from_constructor(): assert mat2.density == 1e-7 assert mat2.density_units == "g/cm3" assert mat2.nuclides == [] + + +def test_get_photon_contact_dose_rate(): + # Set chain file for testing + openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml' + + # A purely stable material (Fe) should give zero dose + m_stable = openmc.Material() + m_stable.add_element('Fe', 1.0) + m_stable.set_density('g/cm3', 7.87) + assert m_stable.get_photon_contact_dose_rate('absorbed-air') == 0.0 + assert m_stable.get_photon_contact_dose_rate('effective') == 0.0 + + # I135 has a Discrete photon source (lines) + m_i135 = openmc.Material() + m_i135.add_nuclide('I135', 1.0) + m_i135.set_density('atom/b-cm', 1.0) + + cdr_abs = m_i135.get_photon_contact_dose_rate('absorbed-air') + cdr_eff = m_i135.get_photon_contact_dose_rate('effective') + assert cdr_abs == pytest.approx(6.091547e10, rel=1e-4) # [Gy/h] + assert cdr_eff == pytest.approx(6.102167e10, rel=1e-4) # [Sv/h] + + # Xe135 has a Tabular photon source (continuous distribution) + m_xe135 = openmc.Material() + m_xe135.add_nuclide('Xe135', 1.0) + m_xe135.set_density('atom/b-cm', 1.0) + + cdr_xe_abs = m_xe135.get_photon_contact_dose_rate('absorbed-air') + cdr_xe_eff = m_xe135.get_photon_contact_dose_rate('effective') + assert cdr_xe_abs == pytest.approx(7.886077e8, rel=1e-4) # [Gy/h] + assert cdr_xe_eff == pytest.approx(9.488298e8, rel=1e-4) # [Sv/h] + + # by_nuclide=True should return a dict whose values sum to the total + cdr_by_nuc = m_i135.get_photon_contact_dose_rate('absorbed-air', by_nuclide=True) + assert isinstance(cdr_by_nuc, dict) + assert 'I135' in cdr_by_nuc + assert sum(cdr_by_nuc.values()) == pytest.approx(cdr_abs) + + # For a mixed material the sum over nuclides must equal the total + m_mix = openmc.Material() + m_mix.add_nuclide('I135', 0.5) + m_mix.add_nuclide('Xe135', 0.5) + m_mix.set_density('atom/b-cm', 1.0) + cdr_mix_total = m_mix.get_photon_contact_dose_rate('absorbed-air') + cdr_mix_nuc = m_mix.get_photon_contact_dose_rate('absorbed-air', by_nuclide=True) + assert sum(cdr_mix_nuc.values()) == pytest.approx(cdr_mix_total) + + # Input validation + with pytest.raises(ValueError): + m_i135.get_photon_contact_dose_rate('invalid-quantity') + with pytest.raises(TypeError): + m_i135.get_photon_contact_dose_rate('absorbed-air', build_up='two') + with pytest.raises(ValueError): + m_i135.get_photon_contact_dose_rate('absorbed-air', build_up=-1.0) diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 9d07eda0d..aa8bcae5f 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -612,6 +612,7 @@ def test_mesh_get_homogenized_materials(): @pytest.fixture def sphere_model(): + openmc.reset_auto_ids() # Model with three materials separated by planes x=0 and z=0 mats = [] for i in range(3):