From bbdf6a3420155d29d6fe00c99ae9a57f2cbff844 Mon Sep 17 00:00:00 2001 From: Roman Krivonos Date: Sat, 13 Apr 2024 16:23:16 +0300 Subject: [PATCH] maps --- data/modelrxte_ait_3to20keV_flux_1deg.fits | Bin 0 -> 308160 bytes data/modelrxte_ait_3to20keV_flux_2deg.fits | Bin 0 -> 80640 bytes ridge/ridge/config.py | 4 + ridge/ridge/utils.py | 2 +- scripts/01_bgdmodel.py | 5 +- scripts/01_crabmodel.py | 307 +++++++++++++++++++++ scripts/01_crabmodel.sh | 11 + scripts/02_grxe_resid.py | 164 +++++++++++ scripts/02_grxe_resid.sh | 35 +++ scripts/03_grxe_map.py | 155 +++++++++++ scripts/03_grxe_map.sh | 11 + scripts/README.md | 39 +-- 12 files changed, 701 insertions(+), 32 deletions(-) create mode 100644 data/modelrxte_ait_3to20keV_flux_1deg.fits create mode 100644 data/modelrxte_ait_3to20keV_flux_2deg.fits create mode 100755 scripts/01_crabmodel.py create mode 100755 scripts/01_crabmodel.sh create mode 100755 scripts/02_grxe_resid.py create mode 100755 scripts/02_grxe_resid.sh create mode 100755 scripts/03_grxe_map.py create mode 100644 scripts/03_grxe_map.sh diff --git a/data/modelrxte_ait_3to20keV_flux_1deg.fits b/data/modelrxte_ait_3to20keV_flux_1deg.fits new file mode 100644 index 0000000000000000000000000000000000000000..5614c109e3b29c01ce210064f76e733173e9c971 GIT binary patch literal 308160 zcmeFZ2UL|?maYwgA_feYBS{P(l61`piXtYIpkTs`2r5YgR0IXYKoUVk1p|nHV2(#c zNd{EJoKO+W0dr0S{}-!A{axKv-Bo?NZ`VC%IR=fu-rqOZGvDWBLwVZK+IG0jKn;yP z8h_=_NR1vEZZkbyHKxsS_0w>fa&i>9CJ~IPcJ^$iTzdo?D9Bb`}zx~zu z{_=XiF5?`4N&kaCzpiKWSFdM6hM3RO)qA@CjGFoW<$5#yH1PPp_2<|1jQ{HOj1B+R z_3-$=^=IJNkpt}p{MGAeXw>ZE7y6pf{9OGtyqxFPti|7XrniP`fWND^-^^Lwe{p^* zdt2Lqb|Z1SUr+2Z$kR1o=2WcI+RJ&mtHwxQXKz2BS-$=nj`RKeUA?+$)U4RqUjwI} z=IS=n+jW`-uHD(v&);{J_bjjZ|Nr%{e2rh9_L(u?Z>Ec1XMK&ivpnZ`xoVh*Uiun7 zHJ3N($>!MeRF_%PT+K89#K=7g;UyK?}B@t^PG`ma9!H~;$UdcXF8>sgNc_t$g!2e0Qa+S2CFc`-CF{PUlIdcXGT z;Oge;>+0>|`Y#Ed>+Ct_KbX|2=6Zk5%YVo9pxdA4Gcw?RJztF<|8~B=N$yWyuyVll z{yd-2U!M2Ba=ifqZASj9`Tjhf;lJ=3;5yyc)zz=Z-|+w6j%!tOy}vbIjc@ebzlx*HjdH0)(&YGP(0I%B|Jh3ew#>g*3! zYf8pdql**%bhowbK48Ge88ghhyv+RkG)9f=it7xqb{uK%;G|(y^90B^b-sp%gTXN0 znRB7zUzq-%j`r&TRzn6_4Rai2i}`FUeXQ)sn{6z;Ep3E1*Lc@MW5|F}8ZI+jT|E5e zcxlY>shMvY%>8#q`Ip<_yq36J&3VblgiM67iP66jxW}JbbdL8-|G)35zsm3LJBoai zBp?Y$0{^Q76yzUCKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW z0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf4 z5|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDy zBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1T zKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i z1SA1TKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW0ZBj-kOU+F zNk9^i1SA1TKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW0ZBj- zkOU+FNk9^i1SA1TKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW z0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf4 z5|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDy zBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1T zKoXDyBmqf45|9KW0ZBj-kOU+FNk9^i1SA1TKoXDyBmqf45|9KW0ZBj-kOck*35Y=3 zrNZAR;(zcK^3sxkB=83$P}3_l{qhGq`aj)~MQosri2l(*M7*9UmS)uvVXqsD(B9uf zkbQG8cZ~Hv-N14nNk9_#rzU{j7=->%M~{3(pEUZ_D>eO6(=+Isf9m0X$HR%uHyp&q zs?j2G=N_@P^B}Rh-C_~@AW6iu=pdq8kBab1j$-ME60x|`Y!P~Gj#${kSIoCK{yRn@ zk0uF70{;L3^g#}K!3+K1hn|>)zPN+la6o??L65XYpKL;}g#GH5#GSi;^$mLGAGn2l z=zjwNad1Ya*xTch*s1$MWFBrQ(v?C)N=ZABbVW&QcvK-0o8A^{CxnUBc~8X3U=0zg zrYDxa(G*dc-9$uYx>$CoKrFc$E5a(ih(*TF{u@S;myiS`fq!BG@ca>c{|N7+;r|Kr zKnnVxGkO89+ts?DC)%Sg9-=ow(I2{is+R#=$ByhObhf)IeO<}P5=B8AN6-N znkam}MHJk+ERN6IByx7&6j>Yhi9-g}V*i8vV$XqPVrS|Av2D)@v9sP{cM)6*2DKBD!deh-%kbM8tg&;c7X* zD{%5KlE5E;0K7c{f8WF7X7D)DlO^btz37)J^vquL&5d8ZQ`0}_p+5l5|J5Ukhr4%*vLWW;_6#R+ zJ-19;e$-7|(DN0AYe$KK6jgCz{ePk#vd|N^&=*6|8*9)X6VW4m&?kQA z6&3W0BYI{w`lc&-=PvqZ?5`d|AN_&X@IQKVQJwffd@uG9pMqk<+wpzGOM`gvtV?h4 zFu6dKy{Hr=WA2IKCqi6bc}rX==_W2_j}t{lYl%X0cX7JVQ=DoLAdY8`5P9>%M9u&w zacn|kk=3U_95$#Uew*Md4#YW%eMODM-i8~*?y(od&RrwL4y`9*+a^ts*`Slyy5^%u zH+2@N=l6;&^R|c-&9fr;(FBo{enD)SJxy#JsxLO^%@XUI?-7ahJN}QtkQbK(e#Zph zPfK{T6+TUeSK9Ea#h8Xm8K&#Lg+ z27Zr+=L+yW58l^@|BmQ^h3Eqt^nxGyVFr4_O%!ESqc_H&KXlO}H_#`h=#^6R%MbL- zN%V~wdS^BIr#E_NIr?Z;O)ve9SMk4cG-{+$m+F5ROsb1jsm{LDr0lbt6mLa}pHV@g zBC?VAx+7G4d>1R;jj1bMm7WqWS{xTouNjHQ-FJxx?c+sR;uUdMD?;2pVk?TruM;;K zoEBHlcM_NHJ`@)ZDvI-0W{Yzli$!6(cyY$&vM5-(RGfTjDoz+Y75VGUMPB{ZA}8#M z$ZirMj;<>gStg3&@QrrjP>7xQjSh+fpC*X?xj|yz%5`GzY$dVBX_(kOOhxP(lpuE6 zT8JGJ6U6p;*TlAso+9(y0g+iNRb(J@TVn%6`Xl)o;eRDn@}~bc2_WZNkasWSehBiP z2oEyh!&G=t2tQW96I1xILhK$|1%E8XUZ3yqX$8FUgu+2#?Rg=jHJFGW?z?E*|&+-%H?q2l$_l9$1M!7>r(6jD9$e zo>+mt$V6}Gqd#)cBl+kPW%SB&^vgo@Oat`IeDqEP`X?GawCYzM)$|hj>Hj9n-}}9& z{iL(h=2$Cgnev{Ro0?Knry^?H-jW)A@}c@=yGXTUG^vzdpgO97q|8NB>(X0NY`a7J zT>Mm24mS|ra}JBIz3z!m*^|WwyX)eejjMS5v7>nT#afgfyeFQW3>Hsr4iS&uoEHz( zCW-r#j)}79x#I4_fuf}MZgG2~zqqB;NZefPBCdbS7uS~d5?5QLiOV_t#ic2`#l>cG z#QDotMA7>F;@r$?Q8;L)INQTboYCklPPaD@1v<&%ln4+fM;3?^3(~~#J>x|FD>spE za#Q3jwHLYfw~O3?r$o-du_8xbM`Z7OBeHuf5Xa8OiDOe=h-0_E^^l7RIA z|4)rZhTkWXf%a(Xb}W&)beuyvmAy$T{xfOxUqS8X52v=HR#58!qp4-~8ft#_6{$a+ zO-%#xsPTRaYUEf>4S)2d2I<$R{+NzbPc4q>o)00Fl|*%7iE2-ZA?3*-q%`Fx)taVG ziZeHpLX3gxqZYB)qUdIICJr3U$FSxDPMftG*f(NvsS!Isupkceu~#c z8^kN$E8?ZLf+&BeES_&y6wema7f;8&5KpY?i^rrc9&xOAIH;3&FiBV357QO*_AL}; zuX>0wV@FZCc#*h!Z>zXF$VZeMJ|Id=x`{hSqQsrSo5k(>TH^Mii{iGvi@0?!S=`#t zU)-A9K-}t(wo_)nFXE$+cWr?_T+Ux(RFvwy5j|gDjld$(9?0+cF0j4autevokPy9A#eB{M|C~qZ@0KLsS4kTdeV4mtNl-UlG}8pyvlJZJ?UTEdH#@FNtSD2pfktKiLD@hqrH zJm2~8Tum|J>rc%DMzoYLccho zXO^RH?9n^R(LbNiL+0osfAmrU`e{^6PyHV)M?P7B2CH5ntDirpU;p0JXYOn=Q(8pA zXEvEyULvE1S=3`)DRrN+ne=*lQPU6Y(bZqiT>p>W4Hv2*yvl>$e=bhBP zku$YDv7g!;KToX}WKgU5ov5Y7Xll_roSK_0p=Lan)F&oT(}{1XNw5YrUSdLx((_2| zY&U9HWkL66m_H>6~EmTLVLLbdvxCB>`y$YB8~s^pTwj$%?6-I5gQ3>QBGibeIi{i52h zKvXL}7FFv{iYkkiqVnYiQJJbNDyKCNl?LvjQt^(cC~F}qPTUq1dmf02v{#}cB}-In ztt%?_trr!ij6_B0N>Nd@KvZfU#_JXpQ5j||DzmCY>$=r>|7Y zIFxGnBvY+Dy{Xn4b5iQQgp@q1NNHybDZOe-%6cnE+5G}3Z#zNC&mE|?t{c^Mc}}%A z=gDX*D1bv2csu}?Dc~~bb%KXq_i3Z9?AbpiOAfU_xhuLJi0@W;=uM_OPX zb+MOw*w1g+QyTVl7JD$V(b>GZXm{WJ2l6*)uh_FO{l&LDrtqvBcQajd9x`7SEAd_Ybk zk=NSD?LOr9Dsp@PdESg%7a-q*k@H5#`$gpb5b~c5547P!KD?L*Kl0(pQ~06-Z#uxA zDe%ac>a6VmuT0^WH$1xm-|(8PjxGG_0uL9#$4T%q8h+}-Q*HQ~1#i2<-!}00IecCT zuYKTmFg#xb--pBdj_^MPJ#YqnunWB~0sWAHo@ju+_<-IBMt@|ZM~cxW?&y`H=ofeN zOcMI01if<@{c{05Gz)zcieBgYk#VfR{U|Ivxso0U_$WvxtGiRPZqKRNxKwHu)0>)|h{t`S@H1&O zHJ^Q%n(w_&%|A4v7Crk=i@;&j;y8YHrn-z;I?kn*NoT3$%L&wqn^UXsYt*WEF17Cb zfLaH9rq<_HQJZ$_sf{;&u08HgZJX7jwr;Da?UA?CuD&w0n{bfY?Rr4%s%KGqyZ+RE zlLfVZIfyz?KkBfg33a%bL>*gRppKK5QpfaN)bT}k($IfR8eTU@W9KE(C_hb_I;o^N z(UCM)pCHXMcS-Ym9BFmAKw6{1NGl|nv@%+dR?#iedVPnq8#W^Cu3n_wznrw40!iDe zHEA!*C++Czq`kHtX>Zy}+FJx^r{?_W{ibw$UKL+=!{c)C_*Qsc0G|I0*Kx=78sfUe zxc&{?rz!5Y0rw5Y{R=V2a?EoJbLC(@JItAcdGj&%In0mQ^`GH9opCO7$>d<1^EA%8 z7Ux#Q`3W4rVpA`0F$16J;B*_jc2K)|mEiXj93#O~5nMgM_Zc{^1n&rNe+mA@*n=nb zQHZ_N$9^2JryT6d6MI{L{i$J({jg6H?6o%bdlP#uz`j3X?;+TKBywPcJR~3&rO1aC za8<}-kqyU?r9_$mA#I7)p|ue9SzAeX&9M&Ttp@|X=J?p zA{jS{Bcp}K$?%;X8M^PK9xty@5C7BDL-{BfY^Wdu^V-z?fe&?$ok!gV+@$V`8%V!k z73s(AB>kyJN#8t~^xJkPeMJk>d!0^t_XX+QI!AiN#-vv=fb<^UAia-xoJtYtYm||` z1>QT)ApM9Hq<^qE=|3t)u6I#)d)y=HJ9R&&Lk0~B$-vm>R+#j`nRi1 z{kzVl{wDRPKV78$W?6WD0UxX4^FH{xG9Gshk3WOwDdYJOxXuJzFB;d?!u1VspKRQ( z6!%S`eiK?_j&Yc04(3{f`Br1j6PQ;QbKk=J?{SW3oaY42wGZc;i*p*_ykBtcH#mO{ zhv(p-2QKZw=N>r4gVzLb8xMX5z;PJ~hwtFp9ejs?vln>hfx8R%FTfs#4tq?%KDDveL)h^TDaZj8N~V*jm>gNMk&VdNqk`Dlrpyh2`F zkekiO&pG6%F7mVpxf+jrTcv2IWAa=e8(c^=a6>~4FXPw~N19%q$|8BxVSNIqSFVDiyeem=Re7!+^+LXZG8Spq5KCgz?O7L49p4-Cr zVesAs{wtsddZQ1Tq8Dt?4?WQnZO|9S=#BO0kCEsRMxRVUuXv$f%F#1U=$lgXjsp6} z0X^h~KDv%x8iRgXf}TqJ)mJsWh5q{Ew)a1Na*95bL6LnfP#DXi&EgSFeOHU7w#=s~=SR@w=WA$El><#`kW3R>Ine}{-)Q{wku)x4Avs-JPEHD{ zG`4G38Z+CRMz2{+qpqK#QBB&?$Vm^$F;S5m9t|Xi?uThaU=JF8NS^!>NY0prx(e3Pa#={ z-Xv?MlVoi+hpe@(lXcS+vQ|Ax)++DGx}iE*xB5uddQ-@HKsT~>(IV@J)?~f^9$7yc zOha31(a@1rG&Ftz4ZWR0!?c&tu-R2K>{JQawEaXj^V*Tk&8}oiZ^$;eAKBIYOm=hb zliibBWIuKv*S}VS zRgYXIyd{@|Iy6lsmZptSr)gUfY1#`Fa_v5ZTs;Sn>(+8|E!816wMXREyE(Z{+eK~> zN#vGMLT)FTaZY&rEj(`np1&Q}X^QLJ z!*w-q{U^9jKiscB?pq7@UyV5$Vx9=h)dBN8$DB%-HyU#%WB!Ra$5WhV9L`mO^DW0Y z<8a;vICm$UUjZCWg2y^=NdcewG`^q`ypq7J4frL3<3RAV2G@(=s{qdRz`HrP+kyW* z8a43~_HhJz8AFa|DzK-^*cWoyI}7{kggwr}K96Is*4VEh_MC!!7h>;Mu>U8>ff4dh zhFr8pJ{lq?kCB%m08dhwA!-Yz3|s>mNc9?%AP9FAN*LOzv| z(+9}w4&-(v@>_r$yCKg9k?S+aH*)+m4|!Ka?x!RF_uxS;d|-Hy4?oVplV|W{G`u+g zf40M;5%39T&^rUa9>Fs$__iJ1p@$Pb!^4~Ku@SsH0zU)bsUdvrKx2Pa!{0^l*bhEi z!|R>!+a8{;h40qzJ_h~|K@S|Jsc(Ov7fR3%&CnAs(HBqA8#mA&*65MB=#w4jmA2>? zOY{t*Z?2$sV$nYp=%LN%qj2<+2l`0^J@x%pU)A(h^r5Z)c;y}S`I{K{~Y30(QFb!Iy_nAVYqA9eS3xz~jQSkmQwD5W%EyS)CjC)8y z`|nU-UJnYaQ$_Rf`-*_w-ZbyQY?{|=5zSpPh2~s%MsxHE$=_!S`JFW&KYb_i4QoKN zAIH$_A!+2Z*OYuTt7uk&GR;!YBkzQd3 z3cGui78^FE#bG^Z@wG0rMC&RoncbF_qa4>4_BnZ8Jr5aiWM(ffNz)i6SzNQ3SFP@ur$0)uvHow-*#S zAc7*t_oPT~{5@m~MaGPwNbDps8LuPKdQ)Wj$v?eM-ipsF;p?OExD-6T9M3bs^OxZ| z&v8BE>rExDpO5<-!TsLjzH@Q^QJ7;R=5fMY8!(?g=KP3x^)R;q=2ycx?%_O#aIUR5 zUjfdEbLTw8xtrtsMc|+b9tPl|3O+^P6a!wV^iVHB!c4bG9^T?+2{ z;J*WVSdM)tVJ`!)pT!i|`V;nb4tpc)&lh{#fqgc|UY}sUw%Bth_PrN-Z-o6{K@K`2 z4+D{l@yN$U3S7GvHG-P-$9CP7*DYtdnmT(ImH=hQrxOO6!)zIt(=@eE6--pDicAgl55fG`jxag zP=!{PD^vW~8x&tSp4OPnqcuD0(c1Q(X>D9z-y%!q@;xD11hVX4ffk zhbFCWvy0ZFN7jFBKpR|AXv4({+DO&3F(ryN)eWOf=%r0p22zso4N6)ynUda~rR1TX zC^_{UCI9F_o9*h+=430{{Q5Yhus5ZInNSK|!)~d!kG2f8q%F%5XiHuiZTUQsQadM5 z>S$+5U9_1}w^>u_1$#<;y^Yf94Wl&e?vzH>lxDq_(#ABQv?&KE%`KSHW_bSTea**X z8sqD0@VLQvygr_XYrNiy>)7IYgK*tVxc)@kCm#13j{9!L{a0g--k4`6=K6y93^3;m z%zF%T_rd(iI7e5U=OoVc2Io7DbB@J%FW}sLaQ+kE5C|S&;PM%KK7vyMcqxIK1NfzZ zV@vQX0@qUT^#yfLy$d?Lo zmWRB3MD9)?e?yVOOp40z{OAr(%;3u)cyk2)%z#Iy;nQw-H57ii!L#S^%@^K{fPcpD zFd05xg_kPuQxBeQhOg%Ewio&>v;!k=^K%!|0VM=$D!3nak*#cIcf5^v@>rP*e2LSM*X3^wUN3 zlmhyyD|+h)`s>559z&n~{q6V1dnBDI2%+O?8I+rzL)lNm>8NH99f|lvhpMK~!6obI z0Dd>I58t=|aS4!JjCj{|MSz-JiA|G^4Gly=bdVPui+GnbN;FQhII&N>6J` z=@YwCx*=0qMQ2LOSEsbltCVK{h|*f46Ye#n)Lr8#byi+%t+0m0yQrc3=MpL|frKIXJ#W{Ni4nB-%`@5;T1@-Z4QYS+Wja8ybf9Pm z9h`WB4!&zbzb(2+hniN?q10(~*!VsjK5j=x>>}vMotcz1t%$N-N77OM5IXujoQ^FB zrejt8C_897Wq(bi9G`N^c@{*uEymKbn?+cI@P2eo#JFv`FwHC2RLsm&Ygnu z3vieN9t*)`4)}zCQ$Bc20=I?WX9bQoz;i6PegxklaIOvBh2U-u{$bd|N$f)hdpU>w zlweOm*w+>8O#%Bez#c=e&z9J0TkLl??fCi$`)-WA>tg>-$UzbEu$eMH{h-VR705{e z@)C&LSRg-!$dNztq=#J1MZN-%GY8~N7rASJ{0%@3FC&l0V#5{4XD#H^33)w$+!`Ri zFDd0?74kd>xt@)Dt03pv$onh2bwU0c!Go^w;V!%w4L|C@liTn`pSIY4fj>6zs4aZj z0k8JMuNCmj6uymscS`WDFFd>rAJ4)|RruK+o*KeeZFt)h{x*ll!{GC1c-J>z<HxvTFM=j zNx3}_Qf@mp%2jwxIZsR}=Uf5h9CD_dtzRi;T|DK)R!~k<4CO={QqIcll(T6xM?PPd4l(>qSmnE`=x=6*Px4LMGQ znm?)VOdC4q{hW$gTT)R$1)cYQO&2=7q6@b!(#4pWbjfNZU8-nHmk->fD_)s&wTmNN zeYA(JCCAbA@#pA9D+9W5J(h09jb?r?oO}FXVEK*C3x?JkGtXXIrzFB z9v6hiAH?(S;raD&oxZqU9Io2{*Dt_*F5`YyxUUKBe++Ycr00I$FxNB8S4PjYt1$00 z%-sv~ufsXE<2>;=S1`^OgLA&Xd2?}YW1Qa?95TT}0bI6%&lzxv1uuPYD+0d_;CL21 z7t`&aui!fhoD0EwF}Uvo|4{6qBld9=ds&42Bx6rj*jFOV`gK|SQ*Epo9J`Iv{CxFatIkedMH=OJ>`3VG^>TrEbv>e7kl@C*M2 z@p$EC-zc=u3C49_= zmp$NTb9mYfz9zw2P58SQ9^>Qbc6h^Uh}3v^{u_LU+jsiHe{;<2j6T5gt8~#1!_X6D z=!g}VqSySls_eb>NoD!AytV2)x@21B$2GJw` zG4!xwKHV>#LiaW_qkFOAscgzsDz%tMrQP~bX?-6msd!5zxBRGNS2mR_RHBl>u2hmc zmP!JSQOWt8bQgV4^4mHpDceMMJ>Sq>e}B5`dztQ@jG?J!O5Q$uAF8?>z%B4T$O8ye6F?Y5G!e2V5PO?tgQN)m4k9w`Bh)8J>HRPpS{O* zOwVwgq_?c1(wkMhR#|;9Oa)Uj(+~DCLZrJ2DH|%qn8@g$5!>BxNxMMRnJe$Q0 zA1ZRguimUyw~EzTWU`uuC#&fWX05%lR=A%Q?z;&0pT!LnsxVJ)%+&|;HN%`YFmF8OcEkKragHRM#|Y;##rg7Z z&KEfEWt=-2=kEy)&%q-fTyBAn4LAw#Is|S{!0$ddo&wKZ;Ccjn6~H+aymP?a1N>vL zhx^!v0`}4w`&o`Xxnf^0v9}S}-zDrZ5&L|Cy*|NyFJjNvvG1<*{&6+-A5QPOeMTPU zAs79Sk6`4)8F^`q++0I`@{ywk$P>QLqR*f0`Hq||N8X%~yL9AlF>=@idCWyFy^+sp z$mt2>bw6@D7WvIVj@u&79?10-`e-tg5J-oApr7Vvl% zeEtBh!{9gGd<)>aJG_4k|KrgEE6@ka(F>{QhyCb@;pmHrRDz#%N_PE*9`QtGh@*=7tlMS(LasRLle+P+2|!l^ph2OY6tpC8NH>5{(6WWoBXTKz7!bzv6}op zKP}hGAH#Le&Ez`eFS$r?k5l;jk#$sopR+4QlvBmSYg9RD z6;+YU5;vu`~8+-b@RO1`YnwLdGkC$Pe1J61S9ofY0?vZ7it zW(j4*Nsg=-RL+Vi`mC6DnH6tuXT^7exK@3CuBCm6YxOnaT4RrKt=TzTD`r2}+TMX{ zoq5Bx%GFt^p#v+KJZ7bd<5(#wl9i5mu+rPhtlaq%D^F|7%4wFY{Cou0Hk{107mVZD z=lgPb}ifuQPGI+yHE^_SJ$ChiKi?dHR6 zPug?4@k6-X2QzLTp~xND&F2n#Lb+qV_uTRFcGj2_!5S~7v!;J1)~rZotk=)635qG+2$hubZSa+8hcUGUko#!0q&e!d^i@{6ovNV{x6u;rF zy1Lxer!jZ^Z5?<0vW>f$2=3-#%H5LUx!VP6?)KG!_1aBgz256tZ%hg6d3R>LkU6Xu zUBG&)JFwoGFxFf34j;R-UhsR?^ID3>>9d}BIqS7O#ofM~=5FV|;d;|?-Hy1v0q&EI z`>n=(8{qz0n4=u?1YoYGn9mw>9>ct&F}Dfk_rW>3;yi&kR~XJW8Rs;{d0XS$0_Wcj z4r{@q9k@&apAc|L0k3D^HW&P6ar^ME;Asf17T`MooZZ0tG`L5A|0V1p8~gCXURq#3 zXRxQO*w;<$Edcx5i9N2yK5e+^v2WP#Z`gAx_I)3FkH-E_AO|VP!yx4167n&Y8){S{ zFSn4Jy2wuka&#GaN=B~qk*|2<>>ctp4!L`U{N*Bt`2J-Sfn4rEKIb8)ACOm53$er|@RQ{d}hcpC?Q&EYY=FYeTZ*URAdDR@2{ zzHfy0FX8`5^guZJpb))q5B<;xJ+T6P5sBU?LVuW`N21Ut3(zZ@&@WTbGn(j|y6Bzh z=%0P)AzSoOU-VKk`pF8jEJa^wptnlUUxDbcv%mVRrq|GK|9K1kF&@EfUk~Os?`^o% zgKBP3K8l;y8OzNExwCrAaBfoO!cEM(apNQ>ZmjZ(8_m1NYH##d&HW6k9advCGi6p= zSfABYkXx(atQIhz)iO4*+Ldrt`|^<+wOYlE4m{yTg}1n|tu8k{k-<&$-f@$t0B&+` zH8-9Aotqv#$?7CnJ)v{E!#VD z%ZabKWu^|dEI+}mblz|)S5I!0{Fz(b4CU61mvigE+qm^YZ*G0aom;={#clLja+?`- zxlQIUZu8;5mvZ3Fvr+cVp_U8{%O&O4CXowecix*6O)(vRD}yUHCV9OVudT5(6q zm)!9fdVw~v#?eu%*{?fmo~y@NW3I4P`E=F}-px8KBU$IbZSLfVzWg$kbvL3X2M*`X z@7r^id*^TC$$+WY*VT&iciY4KF)&e{dThb@Vji^`3Bp4OJ=*11KBS8Ca&wkw$+um&mr6|9rt~J`-fqU z6wEUobDhV0O)%$B%$tb0JMfUNRXE2uoW}|0(!}}h;GElV-YlHk2<+H`!1p>hp8)Tz;Jy?56|jdQ?86LunT!2wW7Fs=>}x&t zwhjAp#2$08Pj&2d9`@S=dse}|7h>;;*nco`FdTU>K`w?M9|w?=c;rP5xv@lk?2x08 z$dey(RfK#6B4?YCHz(xoBx@E{B8MH3$HB;DEb^&_oR%W5^^n`0$Zr91?1ntAN3Q!K z-=4_%4dguux$lYmKZXZ~;DZIXnO*@u>cf+6@MRLbnF@b4!J|<4m-dvPjx|GO+|0nqrYaN$1q>{Nc372`fXTE&;2nP{Xae>8|ZqnUinY%dU_0Ztu=r< z*X_l+6YjE3XI>VY4<*d2) z9c#5Y!&(cSS?hK(Ya5PYosPJU+h*=09FSiZ)=l5bodq%E7YieQtF8*K6-l}+sqv+3>`+_U}`?&;~xJumIzUPjNjSNt~a_2U~0mtHI` z;U6CN-OhBto*7x^jGb)O)q>6TY-V%w`fOfwnR`2L=H8WA+-Ku^w&>rFEnaQlz8ju# zzu~HEiQhk1o}R+}BgV3o-3}hmnt8ypQXY78E)R+;P&yB4s0%03>xu9`>033iAwWXGk)cx2di9u<0tM~8gpF(E&AZ0Kfo zT9n1(7OV64B@1}M(tMt{Y&%a1SK-OaQh3VJ{_Gs~8&3^B$u9FvdD`r2>^j4Z-6kC2 z>2|R^qt8K}+4UW}tJh=qZ$<2Jr-D7RT-a;%1@?BY#k0&Wuut7@>~qDQXU8Y8ufsw1 zt3QeTa=qB!V-L@1`km(-Y|nG;-tyc>gL$4$70>%x&H?kDbHKZsJb!8m&%ZE}1C3U3 zVC-}be5S)e`2Q3T6q?3C$D44_*EC+B{goFCy2lHqXz~J|`@CTOFkTQ?pBMPK@d8%| zUSOZi3yfoVK|?&}{!c)SLe9CmM7 z2~H2dYb&^gf?p;$jsedat~%g*0G!jo+Zf#Ac+#>j*h3=ru?KtczM z`+JT(=3t)**y}dzw*&Uv3HyG9z2ke?E_dW$KHC(1K`!bcANP@yWaMQQax)M4xrZEG zK%QKYtMkZLXXI=>@-`T`vq%1Ja6kKR$m31qG79-LL{9NNzj7OL8;ktrv$Hq@c(@%tX2MH?pBv$6CVV{!Z@0qV0(fi%pBKUFz3}@oJl_G|ufqET_@9p+ zc#S^ThF+L}ei(+HNJL-UKyTQiKh~p1#-mS8pjTd^Uml}p@OSMh^o}ZPD*r?e;b+w+ z=g>><(NBr!sgvj{5A;?n`s)sQtS9=+3cWV^SHGd>{^KV6{eP9M;wN(dqwlz%%LBIP ztirvOY}wqaGqc|=rZ1hCLNuAw_OM7RVPUz4dp$AXUa?->tKU=>YKh#dz@Npk{!He% zEPezqWi(?Rrp;!D(%IZ98o6u4eH8S#&%7kISf0UsW3{=aMY5q0K7{Qymau*BCmwE4%){TD zMdo_%8G^Xvo3>|2_` z{`Kv7PVf0V*XJ$I%M9j#ho3l5qdf;rKE(@C4)Ve$hd4xUHHUHphn80IqQ$j2tV;xk zU3BBc{vo`i%|Tw0|BRQqec@#dj`OnpojH7r3Wrzp<%pzR9BFC5kqRMIhT3mk8!*zXcDjbuEVSS zqImTSUyh%c&GGoX%o?iXHJgv{TIH#{*5w|r&1=T%n!n(6GluZGL+yFp&nQk9Xu$~) zW}HyimlNwY;zX+iPMp7&6EjCuX_@XB#R`=%hs)f8>?KIvt;_pXzy&@j_X+0-C zp3aGvbvW_xVoqE?kP`z;IB|p#Cw4O7ginGKa))rjVjoVh%ESEz;=V0#{{fie3Fc{p zxo%^=37B&~=KX@XtuX&@IEM$$vl!=UkMoVeIUR9cBb>V#&fftXe89s1T!O$S5S;A6 zs~xzhfS(aK=7Z-Aa9s<&I^diN-aR>dbQSow#~u!2A1|?&yV%cm?8zVd8jrm# z*yDTbGZ=eai2ZKHo|CZe-PpSd_Fsk^%s?J)BNr~n2O}rZ$V()0D! zbL49`ayA-yQ|0LuZ;(Gd+sWg~#vV^9*=B1b&Z#=dIxT19;y6{(nLbd`2I% zMK2_yAJ(ELbkG+n=#3@lk7V>nd-RDKdc_z0l7*h>gT9%C-dThGnTQ^mi9YIqUfPL% zx`m#?&t3T0j=pH4zjmU>I-$>O&}((jZ}-u2Z)^JQ_uuM&`-D6WuW!fdJ?2rz(UGS2 z+2Kbmc8FBv;X$=|xOOqyUmMEyk^XFt-|yJh>%(@}@iY6@EVc_-$o3hoY`^pa4^I+2 zVp}LXY*S~)Bn=)Jd6q|cwcyd#b$LwnD;`^Lg~#Fd9OKU^@(*QO^27~nOgQF{*T6v06q&*OEXA1`?Sk{905=HQcSIV8CwhxTvEp?gR1 zA}u#wv?`p#>NMxD`OzHqZW1q^5y*=l*X1QsJMof7M!a--A}@XG#>)bibGT|R4&UI) z5k|E+;+!!@x?SR^IyX6LJ8|^zCmj9bG%wF6&)v=&hIQkO8jE=2 zmuTK}L7S7(v^Y6vDsOhk;goKwyhTxiw-n#x)ZNOQwm6v6$F||EdOLY*bzRQ5q{Nv? zKD^CSledd<-l4dWcNEs+ohwv$*Kk+f-RLv#J|D?@mVV~Fz3cPdH=B6hrV!pgWC`zo z8^#A#+3-QbBYf~&68|>&I{)_0m=6W5=EI6t_;A=6K3qMLkNA(`Bae4*mOXLSK?OeA z%#M%xh4axHSNK?;&U`Fm31>Gf<806CoPG2=XMYUioYsdqr|U_)nQ>10A)KR9#M$^+ zJo`!lXCJK3*=z9T-<7k6>vMMJYCiTOhmW1f<6|q*`Ix;0A5&e(M{}n0QI7_Ev}q`3 z9k_%08)J@i%wvzaY%yOx<}Ad#(U`j%=0A^ftigHq;#^iZUjWV-g!9h9xmV--3gGY( zJc7aH1^77e&X`KxarP&;#erWbI5r1Q7tXl&mA6(^f^#x>Hw5?DoVx2P_E3m@sB=o! zO6+G3_B0CnYKFZv$Nqw`M}O?IJ@&c@`)!3i@5jEEWAA6N|8V5s1@h1cxoC-ebVN=D zAuk!o&1B>!3OTYup3;%4ntT-^XW7VG7IJqB`O`xV3z5f@$Ym#vF#L|3jzwO(A-5sO z?>OW*26;xNr?%n6k1CP#S;)I1avy^H*MkQ^@L@c>7z#gB;7L9B5(aP7;LlBXGzC5# zfmd_jSABT40={X(J8SqC4-eJYcgH(;>B~O+8J-@6uiEf-EBsA|$1ULVTX_8(em{ff z#_;_LyjO((uh0WM&<9~WzVIvhApkvb7Jcyuy`g~qXoMcY_nM^E=#>WOmq7H45&8x{ zA7xBN{{*0i4x*2~qn8-{HGTKD+wk}Q z752Zhk$t;o@~rKR*bA@AJS*$7`{@97U+>87{SA|zrz2Q}-uJGzmJ&x~d!12$g@S5Z~ymnwWUi;dL z*R56J1iZ#dxV@4SgX;47Hrc#BTbDOXbmol;WxR3g2Hs>_&YQlj;H30coII*KZ&n@0 zoAd8;3g+0-dnl(WwB*z?m7Er@#Oc#|@>W{R8I4OgpUB%yv`iQfe^hC*l5#rqG_u~Auz2d@i2XPVKMY-HqQ4bm33+9N&Jo7MD73MRSy)G{6d_&HD$lDybZ$kc6 zxPzBC8~F!!>5ltYi!*(`<6f_DH+$T#7Iz$qd&Z00{ExWrGTeCv?p=+$H^=?kp$97T zVGVlG8vW>wp1ed~Mxi&^=#MUXv;ln@fL<*`zs%6H(dgR{^v(kPyM!LjM<3gX{aZhw zpRdqUP4tz~+s)|jRrJ^!eYQZawbAdD==nYLeIt5*75yIs4vfJ=rr5Eh4t#`ylQ|-G z8(xdN?g4)Ax@7Nk@H7})S%WY9Jo>u?ygd<7s-NKRE;#H89wWfzKJba3M|Ok3Yddf| z8~io~$1&jfF}VH)zAuCG2jCs=S@C{()s7l?z#TroYqQ@|;Rh3V!V|uj18*FKKZe62 z_u-QQc;zOyJT(ixvW2&9!e0a6v4!xN7Q7Y< zzv;qr$?)B-f4|lL?+r!x;I(3H%@ncrVk@yG__PR9RfsV0RD>Ee5TUK2M5ukV2vypM zkUmRA=ya6`4Sgp<_b(Hn=e~%r+xTqZkCj-fHbSiJG(@bkzbeA5u8Z{sMq>SsUSdOJ zzSuYwpCjC_5Svn@h-%bNY#A6Xwnp_4(SH|-ZLh6GEPl>!|7;_6e6tk0{Q8L9AELw_ zn~`EKKHrHy<|7i+CW%BmquQT3MkIY}Ba*wii4=UVCN<7Lq+Mtu(yPJ4`&S}EJ5*#g z8!HZ59}!t?K8b8o4ROTqfH>MCBJ zbkqCd^tUb|uj+uvFO3%kc}qoM)>u)L+E$!Nau&r0Y{cK`Z^hZ;YeY$jsW^A{pg8}@ zQd~6HBrbKkDK3w$7NzTx#FgCNqU?RSxN7%PT$?dQR2-}j*FW|UHsTl5t=li|!&Q|od&R@}C{fj*T0EXQMLa2eBA#|{E1vEBB3|g2i}FS#~>f^srs1s;`>E>J$YVyZ~a8nt^Ff@)Ls?8e42>g=`*Na z!+umRvn|yhZbA*}H=qUwa;TwKQ)*cAf*OT4Begbxq;~oWsr%`YdbJv9%#9?CUm>I! zT0mNwW~8-t5otHbAnnC_Nc&Yg(iwqY#jB~YjVU!=dy^VJ-buRbLb~xAsY#m&)a2$e z(mQp5^pA#9)2wr3uyr9BdS;Sg(HAnTokoV&CXnIE+hnMfO9l(BQ`7rt)YKz^^oxg* zzTXnk`w~KW>Gw%5;48`s&x;?${@85)c^^gYkCFdA+(C#B+HY|e4cw;$cbbWN z5$-k{_w&XbHE_==+*Kd+=^*gB4cra~zrDb*F?b#buGfKY1)SqO&Sy&zyYxHw zH-!gI!3WOp!fN=z9iHe0UzEZdgW!*o@Q4DRG=*2n;FnT(#t*(30`FMBKTqJHkMPkF zc!_zbTlRoCFR-2Zs)AK$@y|F;qU-TzeVxYtK)pR6Wg z>x~ex27|@6ONJtP>md*6O8hJ6+_TwZX#{j zIC1cBj>za?Ee`YFBI`hjI1)Nq9P=3_avJ7{lT9XzQ+P&tddYQ>zgbHZQg2aoElU*d zye7_i_Y~($&BS?pNAsfV1aYZ!izuBjTU>dsA<8Fx5?8B)sPOJ8t{)g8ZmLZew}x7a z+iC5^oe!PFJqi%_mxPN4r=E$*FJ|Hq=Bt|7MLbTi7Eh|3#8c}a@oZA1c%B|8UVK_2 zYTUMp+HFI{1YfXFSQk4e!dmoLK8)uL3i;ZX|4E$ z_andWKc{*w{#3tNDb=q%OAXGpr-t!&snPPkNo`Cisk_`EjW(A_Q#+hA-+v^nyEUYJ zb^z&QHK4|ERiqnnftoCdA-x%GNPk3sYU=fh3>coYe?iU6KTz`~SE%_PA2NFDPAw|e zk@4l{)bdnUGC3g0bn{g*TUKKL05j(UTq*5Jwpd_56Y-hBmcbIe}mKg;I*g7-}DRo zo(IPn;CTeNJ_Npf!1*BXo(Jx8!M_PSkO?0Qg%=2Z!21dWRNC{6C$=)|ba(a+RY4TVcs52JH_%k`l z)xSl`xRoMJsSs(WUWjygL>$sI7l$@Z7n%JWMOMAtB0GJDIND!R9D9)`a)MWjliDpr zZdiMfS6^KejOZzf_J)b#YtiB?zNd0t(?DD-z-JJL*NCzn$>M4fi;5qy;--8lZofXV4bWePk+g^N1*(AOc4-wx= zGDO|Ym*QvHHu3xDHL5p&sD7U)YGC_~8X7xO!(SVy5nQZx+Jn?%my^c4`lQ*f6=}6I zBkg+6Nc)~Y=^UR2c(Qq#I|WKfIOs#QIx*g0~`qfBk+hPv2 zb8k%T7ig1Ba&PKzFNW-z-Xwd^d*l!jOdaz#P^Y>FaGV=;Sr|@TkF6x9PoJr~(-G=1 zKZJT7?MUL|QsOQhNzT?K<=|aXy%<6+R<7haVmP^N9!0(WR;S*d4XKZ{B#-_?eHWCF z=Z@{oBNHcK4lm5p2XlpCK7l#6 zVcx5l+ZOZR#~wD=XCU_4js1MF=Rf;uWAFXg|2%RiL>{-0%Mj!fj+}CkmmPAmM}7^d zSMgWm*$=sTB40=3+y{AIBsuFB^0&YpX5b#?xXWDJCjob=#=U0YZi{fg&bZ@c+;b`J z`WNnd5qCa|dxzofU8(*2@6@jM5AtjU?T9Khfu<==DhS+XFqnjlMTQ z?+x&N=My-X10EP$+yfsg!AV!}@*3QPfuBy`=rnj*39hb#uN~mb4ZH<|yJq0;0yu09 z9v#7DEcjdnPBp-5eQ>J*eh-0Tyl;8A3tZm^-|NA7C-A--+%E{rSM4(ywU)E*$U6Zz&CdA&Sdz<7#_mUi&NF`QW^Xd3{U-l zuQtJ3c<$2V5j?gKKDz<0CBbjo|KquTe0MM`?*BI8zxSVtBG)KUkn&#Sx5*ayaSlAbuZfOZ9T^QG=1b)Ub9isioolxW7MXo_IsraHP(g zrlbpg^whtQ{*n)5P;WCC&ab3q9|Oq9vz%Jw{ic?B-N96syY?jd8hl&*J^GEZt#_U9LL^mIHO{TM<=lH4h)Nj@F& z`HNDnHlqFYM$=x$4BB~gHN}*FrmabbDN1~yO~Gjt*)NDTu2!cFUn^+C<_@%>M?;Eu zlt&S3cThwJN$dZv#2mvgPd4V#!hFjyXIISYhq=`-|03*h8~Y5!UN+coDfSG*zDu$9 zW?B{X13A1!9tp^0GxD)VPAh4N`FrH{4Edcwj?am~Ak_!Dt~pyzn4kv|{3e}Vq1gM&BVVKcaB4nDSlld0gP8r(bp zKUv_&20Xn7SH<9KJUF`s-a^6MbMRLS4#$BUHjXx;v&mJirggJ;ulShgtJ)K20z6bx@q=R^QtG9R+p(EaQUn1V$T_rwO zJBn{}Q^Zfd=2Wk-6V-oml^PztPHOY5Nu%R-(wcgnbefowZigz;>xTEiT_qWA974@! zo}?D>Nz}3kpHGwzAq$HUWI1>^S7S8s*|nqrb+}*dq^V`~tidm7i&{ z_H>%^AeE-&7}JcHESecSngXY_q}jvo(_Hv-e*fRJ(62oOjkrfkroExy;BK@$dKIl4 zc7;|=PoverITR9gi$YUsD6FWA*4&>=Yrn0ba3ft>?^Hn%gEMKvfooS2Z3)Kv=nVu&6Rd5vNHJJ8{zCe2e zeJOtF9@-Z+lM*5VC~>1d?O%VBl0p*bz@izHJozD|^v|Z$&e4?Cqz9!{_owtD8|mQE zgLFtvq(h%)Q^tX4%JkhvneSZa@VaS)=OmPU;yWGpbE4CmoapSV{Z#I>f$k4}NH6+M zrjI8d)34(_xWU$D+~~apt6eE#HR~JPP%V}Ij4Plo3FY+aPyjvd5J9(U-RR1uNmR1w zAr;!kQ|=5cI$^S!j=j4?N4K@3BRNrY#59w#caElPi!{pGt4>)hYB9%n%+nKdZNz*L zm~$%T9f7$MFn$MLto;$qjksBDX%s z&ki}_XHos;$Tbo9#vJ1iu-8cPVaE9>A2fI-0v)fH~)is z;^)r&A-Hc5?tCBjzJR-@;{LnQgF^ISGJ0Wxeq2LO3elIw=*<)KXBv8Rh61PlMz4CI zU$Hdp_#gCb6M8oQ{kxAIjzb?y(93n`XGioj7JcoH-rAG@rdsqEKVLkv(d*ghw?xlv z$@}~(^gfmPE_($IMuP`Wa8V6DZh@0Z@R9*;9KjF%0t4jSJ@XB?N&sIb;H(>XGbQ`m zb>QzPI9yNdmf^M2s!!na3^?r%UVXu>H~8%Wj(34)Jj2=O1HQY0^J4IR2i!LU|Ap{C z6ntO?FW|G!Czs%fj_`#xyb%X~w1Y>2;FCG2Vwz(4-*P(g3q4AYo74iSa{AEzB>c&_4vns|GyRgz5hgfoSG%xpFJwx z(QWaj#9zF9e^R_j94X#b_Z07K2a0$4%f$P?ev1#KPU4f(Wbql#I=_4W6m|K=;%ARp z;!hbqV?5+c4FX$GqmJ!J4PUm<*y~SP12>XRRa??sJCF4KOr@r~=TWnIZe%p{E*bB4 zA(Ll;WPaIzT77v&R;MJjo$-Tgw2zQ2K38=xt0u>$w$$bIP3m?`mwHYsBc45&lwb46 z<>qbb^?E4viK!vaMJK4Avpf0X=Z0Us0vdeMhlWNt(g@d8H0qH$jorAHCRkh2q{5># z6(49!KXsjE_VJ{^v?Q9NUr6)tT5ZAQXbQ4oS`sjyf=>*f6%AZyWugPEI-5ePpQ}@- z)^Q56zesBa{-U)(=V)ER8w$VjoYw#DK^q)@(Z;D(6d7Mjn<~>N%5ng0nW95mGs`IY zR}#heZl~C!ZM6MQH0>O{mUb14qPVv8Xiwx2+WVt7?fWZ(5+2T>{X=F_Qh5+1`^=`4 z(hij7yNJ@S+tZ4I!Q7jq5ha`$djddQHVns=k8Q`%5{*F@w_)##R57>?IqXFqPRPXp`7}UI*2t?Ka#KToamcYF@-#!P=aBD45kmw?WUlqwgcodn@#RDL6O+9#X)? zM)1)EoaBR-p5W#messW51dV$50bB)wFLQ9l;4K5(tpR^d;4lR|-UXL^!RG~VIt9FX zfZH?Rw>~(&2%as#HJ%CWybjJ&!TWn~-xd5nfCseUgDLQWfFCO035G9b!W-k^j{tbY z5z-nq<8r{nw_{^!amYW_FdDn+H$}pK5BkcnO)tUm}Y; zcs}^?AX)3Ypmyj?2fN*5*Dabl+TEtkMhD5MqAT@C{F_+xC*{*>a-H^zdSB{CeJ@*6 zKTTWmz2!oFpEPKwZ5@p`YD8mf&1r&8Ax*wIfTq``QozzqG-p_IS}^AdEpA#!!7uG- zkZ?b@3{Et6@Dt^uvRa+AW3e4&VWO=x4kL$pba36C;q%c{NR#nZvH$&_K2Mwu0# zC~Kt)9qHVJjy~6=;|czBVnR8cvN=k*pR6hGtTz?xZbn5xp;SC#8=ZBNbgokXU9eK8 zOD*3}X^Vzb*3yEmT3S$r-Br3FwdvM?t#oI4H@X*goE{|IriT}XP}SR7s&45+&)gm8 z#k?L=lh}-2U7ta3)o;>!@rgbL45!cWvGlc~C)G77rk_qZ^m{@S*X#3`>yLQE4Q77h zhQar^(fS^&7JY)%cdN6;-gB(EH;%P-*|7GOMywNZnH$eGV%;H2xJma7tf%*o^&V@o z{=r?`be1Cb5-Y zN--~ZdzHiN)^p?}TaFphl6Qwa;P}lOIl^gQJC8X^KW6jsz2E0CiYs5{bpm&#@M%rwWGdce{zg!Gu(GE-Qah)_X*rx4fj{4QlsbS!wtI7`VIQwkDjzaUzVdc0aUQ-7kcD` zKIx)Yap>16^vs!#KKn{X9KWM~cIe?l^f3m#+=_l~pwzy<(bvZ4?N;=6GI|_A`+Pp4 z*D>gKAbM_wzU!g)HRyj8I9LZBe8GhPA9mox0le%5H}%2KOK@}nJmIyE5jb}N?;pT@5cuy6 z51fJzY~h7O_@NA**Z^Nd!yAp@k1BZN0DPhiuav+q3*ecD@Xa)MrxgBq0uSMtc-Cup z3D3iW@O`qdO8810-YS8=zQJR7P8XODueFBXdc$)G|MA^_ya)fmga2*Df9LPWxTqPm zsF^`6T04=kRTZ_EeVvR~>QjsA5!AB2B^fVmL@kqskx9i-GE?^=bH$4+@H(=!OA%SE zT1{=;hETikoz%fSoa}zKAP1lQ`S~HpJ_WhJ9)AZ9Q7P^$V@~TuPyo zN^8#Hy|Fub_}h;n$84g_pY>^L&J&8UtD)`Zb7lF9v5XIv=fQinQlyqYbr3@!Z zs|}z-GpZ=_<5bF?{*;bB-AFkD4%5l}sdU=zE#D-1kn5W#aDx?}x#2T6RvWO0)$>=dreiPG zN;<$g7VEk3&Qfk-(vS6aS8`LUVQi4Lk(+hv$IZ`-Gyfa&W4du2Q9l3q$2W)e!72Eb%!*(}zv%~l*?pSl39cR7d&hIvH*O?CN z^w^ra54^-ZPAygCSPbG_Il`W?Hb&0x2z@7ycx1^14QWB0&a z+^1(G_xUh}`}UsA{Z-3&#K4st;OxL5%V+VH!J9d8(jz{4w~R}I$MMZCOZe%fC;Xv< zG5;74F6#{(F6-Us!@nK$`HPn&*C;#q?jt+CU{jxu`^RwN$tApH>t_ymYs<4s9C@78 z0ru}5!~NlX@AWn8d21H;ec74&sRPl9QoWvPQ}Qp9J#GWe)Ex|F7mvHTrH6AT;yzz zybmGw{>Xm{?$8JKAlxMa_X)7ZBP=vO*=W{kd#MDLDLmC+aUFbI7# zp}Wogpr4`WDc-mBY=GV#LVr7;$M)#+T=Y5|{l0~sZ$;np(fd8<|6On}2t34rizVQr z0-X4RmojiO3jElDqw(OWDY%*nzV?8#67ZG=?jC@@Ti~z}cnkrT9^f+?oZ@{2#e>^o z@LNrbZ~X?(J;AjB`1S%n_haK8ckpMwYL!w2)=g}(5^GM!+Ku z@X1hkWj6dW7oOPz-!z4HLgAln@X$c`XgR!;2R~K8Qy%b@2fT%64@v&;SOfUX0$v*i zzv;nq74Y37c&{t`H~k+E{yPo--&>L0?qOtmIe=`nUy>cZ+iaU%M7AyPta7zK*}o4a z`^B2%@D<;=Te_7xH^OIkTRW3e`ylEO=|-Z%Pm&RNq;gzBuA8<{@0wxM*ENc~W9n1? z7vsrKOeg>F^E6_(294g|fF_KyrK$Q2=&z%dG<$$L&3}A?mSlaW6`lUj>OfrzBWGIs zYa^{c@tih(tfZ(`DHN?7q*z}!+UcJ`afR1tp9wy1p7@QD_l~FZ4xcF_{vu_oeWhcm z3+QC`Vmf_z0u>HlM8&o1>D;~(bV+rTuAI%GYg5ALX3NQRC%XYXaOz5r&g`RVE~n?G zpHi*sc6yWDm_8Ud(dV%D^zCj0{Z#FvKij@@{kbOG&}0HPD&ENIGXhxCP@A>#0$6AA z2i7&5!%fbFvHpB>yq0{=hEFDP^LRLq1oN8E4ZK$439k*_#_M9$dE?$D9J^4953DWa zQ~fmfTI))F71B&L>U~2R9BUyh6Ro5}MvfE?eWcrhaM|}+1L^&tp7icILH3DgER~;b zvWwnSX}!3eG<5Zo4K_#c(;h{9L8BjM>#KNQ!*t#jUCNOc`*K8gHx3WA;&mr4^E&rZ zUR$k!IW}XS2be1i^DV%fFEFn^<{pLl3$TY5_PLF{T(Mss_Eh7gt=?hphS>iwau|X< zoRLcc@;QK|`hnyE6@5{)22J*j#J1oOJ?&2ezdd@e9DOiBFU-)7QRvAF^yMge za}WKoM307}Pk!jtS@i1#H;=DF-v*#}mFS-~dN>b#9D`myLO(~cronggbrgCVjQ;we z$J*%gUG#c2`W=X#-$37u=}p#q^nW=xxC9<@z(qLtFrzzJwcy1U+{^|)3E(ImJSBmv zAn-K>oDBqTN5I`=@b>~7=7Pr;;BpZ7tOci@;I#zY;(3Q}4{)3Sp7((38{oSVoQHz< z?%-Yx{J(<-YT<)4n(bQ$KU|}!dSBrSOB$W{jYbUn1&UTMj%|2;X&p_g?(RfAHXc zoAKZIJ92qqOfD9)$a%UTRk0?yEH@sb_HfX@^0eereA`%uKC4YX--Z`v}rCvAI|N88Va(eA-f6u-F%?boTL&`@6+G6PSAz@*Xhd0<#g?LTe`KWDcxW4gsOaJ(lfo8RBJqh z-f0e|&s`0u&NUvdm&_Ue**R`hK8!Uc=CRhxbKE$%g!PP{aMJ?|xLNymY;<@ww^YWl zY5rZd=#`AuZmrqcgSaiYw^7-$ZPsyi=$Oi#lJ9VrR{hy&XA$=>%4HFo%hLQ7tCCIF z)#)7fI`@wIOt{CMYDd^Rqlo*DE9Zfxemto191l5Yz{6KBU* zY5f}WU*kMEaQRQ3lVZm6?``3rmTtUs+!tP+9mlKGKk(|4o*X*6Bd@WW$7^ezb9lx; zj#x0CH@fKXCf$p?xpFjb$(qj5Yu0kiXiwhW&53vFp5=H)Va!sa%wnpxQY+Iapmj<{rJeMI6me#oR9C1S}nwf)I?`-XDCfMs0r=?33OE#ul+BiU$8f^2@#Pj>i_ES>$n zOaH|mDWGNOH=+;|_e0v2R++$e$H5dLp8zvGKWkGR-*9~W9%bAEL_&L8c_`3;9)j)9ox3+5V! z`Cehp$Cy`Q?f}f+4STG_KAzZXBleTn(-HgHWAEG8KOZ^xB9BVsasc_9Lry-(%K^E0 zAio*NaS-x6j9e{|uM={fhrIofyC?EL&KsS-;vVC07X|mZjysLRy-aYo?zrDW+%X9E zyo$S8;J$%8C;1!heHV9c$J4xP(1W(O^pEHbKHD;NLXWc1r~T+v z4Ej}oo?S=Z@E$C=9{QJu9_~XQ6VS^P^s@(gdIo*ng5JiXzs=BN5A?Y$daa9oZ$;0~ zq3`|C`!VSMLvSz|JX`}8IpE_2ILQPrUf{+S{4@kdYw2^BI&h_cFBfoz&lr7&fx9xg z75N(+8iGfBPPxAre5QcYQt;{nZg+uS4RE{;JdXp{rr`S#IBx{r*Ma+B@P8B@z_aef z5?%;}ADrO{5BTCAZ%l_j9?(L^&+v%}yiy3iw1Q{;f^RC|ofYuU7I^3beDpWG6b?Vl zhNmvVSA*cKbMV(Acx(uK_K(+G;kUi;oB@2N4)67Z{{sK<;J?%G|GgCrtH$@K9ZP5^ zp2rNSTSmieI?%A}Sv2C)0UE7xq){I((b$36H2y?3O|ofCQ{o5H^wzs6(AJe^pBh37 z&hDq6De<&y(H>gW+mAv@2h*CNzO?RfE=A1NrcEZBC@S?1ZS(C$v5(VfS6DRdZTp)N z{?4Z4G3P1mn;&I{^`@iRVRRzggYp_2rJ_YORB~nzU2^zN$?tP!-E;z!lO5vR4ihPH6Pi^ zb~(2#oX#B1MCGsGnDLmxn z7aoz(oyW{F;0e7>@)UzHJpI8_4#?`lb0RMAf+;_FvAYc~bNGu_Zp`G>>ql^CggLL- ze1_M>dGPwoblz}&El0kXz)||sd8_k6-nPh@JmZJGEtXmR)eIv;Vh|MQz11DXh}`C3aMEzR%)4S zkXnO1q*lNSsl8lEy>owL%eQxA2ai!wseeNHANVBa6t0l#rXH3%y4;az;p+0VsgR{z ztK|KfVX~&xL-}#*X8BF8Tz)ys@~u^&tUCEkmd){$Co4M2gq`@=aLP>voQRQLdBxJP z$t-D}Q7R3#ACQg9%BA+8H&QDkLuyS(kXkciFvl~@n7wIg`5u~@1w|l1oH32C4RqgkNUXF8QdoxcY26>mE&%I<9>(uu>CvSb0O}!3HNQt zNgqGq-lK4Ld)$8jdXR!XG)6DH(T@)3$yoHo554h5f99Y^cpnwl1-;5ZzckUaKIq#n z^v)gqYlt4+Mjw;VOFi`SZ}e1!zV75PfgjOdWAykA`s{^XuS35B*(3Qk`mT=N8?fy6 z4IJ2mhe&YI418E}hXtR(%TRC=0Dc~Wqp9Gj8@SpFzP!O%5o?8i0(bktpAR@Z2p%tk z%V6-i3Y@xu*K>4l;h{MA$QNE>_$d~i>JML;!&@2f zR|Y(G9zJt`*An5k0C=uFd>0Jwt^bez;KBbkg*~m&izVU%j8ba-Q&jcH7xMEQv?l z`o`l%uHY$J^ZBm>zj)5Ha$aPr#=);<@T#FVIMkyNuXWwR>$@%Bja@hMX8Sna+NmSQ zbX?3ky8qza9;bNkP=ym_4dtYDV>o5M38!CF=ZsIgINR2kk4xcr(IUtju^Z##YBd%s)q!&}w-^@1)?%LAZXuIrp&Q11r9VC61y_N$O4v_=fOp*gbEakvW8FHZ4L^Q!ZM&QLf9qCAaUYk?D&k$-?S1dHa&H{IIf>qFysdY4&@QVs%qp z=`<%u5sNbw=YSH$MHgSJDV?Y|>OWLkH|VP9xqHfQ?}p0j3-#oYJvMUtcyAd}(oX)h zx=fBpJSY839c2HTB>NA&A^ZRAA$_0Cz#Q6`XUhM~6=H@tM`KNAM{X#9L zB>q7Ud!vsEdWqK%9*O8_HTpUcy={s9u0xN5(dP;1^$HH{^ASD&8-4fWMJ7+t|7>vZ z2s{h}7n8t8B{>j*^KLdE;7(Azf@8aOSh49~ne?0i_H2i;WMd7`YDO_y{ttykdudcTLXUYw?g z&xRE7I*B&sl+vbkgD7h3BHF4fqHR`D6sy;WcIKJT?wrXK9~ednL2c;(4WTqV4>?rY z8sCX-K*tu1qmx4(Q2w2Tbf#n-oyYfPN-NW;Vs;4K(b!9sNm^9xd4X!@oTc||zSB2< zL4U%#al`S%8ri#8XYhX3OKHG{KPR!V-*qiFczCLUd6 z(*Z+dvnejJMffdga&ne%~7s)Pfe#q{JJS43$m#XPI zq}vpG*(b^JX_9aNiy*DP&x1VLm5;SC70c6 zCs!?vl%e{|1 zu1butP`3AIt89JPO$oohLJ1lWp$sovt8lvyip7IiMSZxCPi3V%-|Ct?im!vGZt;-& zr~EAwTSd#n&o5+B<~y0Z^ET$N#XJiz*DB1n8FQY&ympvd7xPcV9`~?MZ|r4>{XDVf zWf|}J1A7Nz|NF=x0C_|smzK!KA32ReUgMG5apczpIZi;H^W?I-ACT`D|`Wt%M5`A?+Z%?AXDd=%q^tl^)twO(b(eq;TeLQ-fjQ;Nc2XnxK6S!CaK8Ayn zGVpRA-0TEDIpF9Scq#%{cs{iGEI9Mz{?}fEyA9y46F4jdk73|)Iry9nPQ$@#Q*iqj z{GJEL!&xu+E4bbTzPEt$F5tZ_xVHuWHSoX*_#lre_y2+)qTz|f@WnnlKj1g~u>u|m zgimh3E5j)Z&v6fxn!-26@XkQ^#}Xb|2_KDzmyF=2dhk>vd^I25nh$>sgvTzyXNmAy zQ}`_xo~wlK-otx!@ZaKpJow*c{P+EiwnvxI_VGvYOr<&P?CwuHTUXNVPG$JjoA&fy zN_)l^QT&z*lyIqzl62!Kd6*`prZ1y|nm%-RkppEvi=muJYv|+)KgwT}N5y*6>D*ys zx;!F+uIk}C0%gT?e?v>EnmC$XbSa^?IuZ2w&3pP;DY^dZVXQXy9Bb9wWZe}**r3ro zHrnjVrq*}ZGCr5v_CCe7Ww+UJ!XEDS+mGd9M|P{rV$aa!Jiu7+kepB+W4Dker*`51 ze2;r!pJ}{2aSn$(cHnj9GkBvxd)}g5%CTDdyt~P6-WTwlldY%o!O9VwU5nRAzu$3z zNerL8*q$%TJzT+H#kGK|22SA@$p!q*HHE)E3+F$jezM`Cw^HNBHrcq_SgF5goorV8 zL$-7rEzNT;ORGUKvb}b`v?~vhjt6RGw>{uG<*Rf)zgW8K^pIYY1ElXAUFkn}tQ=|n zPL6wjQBHX}Th6RIE$4MzAs0ujmMeAKWXJ>yx#mhoxqg7V-1xGyjM|Yaw@pZrJKUzr zIH|~eekC$#^FQ2|1C2O~Lm+uEh$*-*{ine)%()eAdqL*2w82Gy> z&2DQcEjmRjEeEbvOgq0&%!+Fi^Y6D5i?dovE2FK7#iP@Tg>jf-meXA^$=6hjd+t-3 zzw%X@FK|#=7>`v<-TzXo0-7oIvo!D*m#k1c@7z*GFWs&znEFRq_v4+ibIMI6Gr&(d z-%3Myoc~DqzQk5lzw=S$>&C}Q<>s48k^XXJ|5JTsgW6UlP-mUu{kB}`e9b^BF6OboTzxR#Cd^rgd806QL(E@_J?>+l zo7gJ=`z=$N-Fk<8ZLs$d>^}`TL?Msn$YnV4$v{rWkk=06HdNNe{z8szkmp+Dx&-+a z%W|8`$U8`$efk#pGwu+Cdo-5W$KS||t{-tH9o%aG?zR{AYk@oV#ywZzuG?{6J>0ns z?yZNr*TenWpa-7lgC}~i7X47rlO5vrAtI*R~=&J{MyBz(!iyqse&lAz>Q1p8@dfo?pk4NuMqW`PGfi8Gh z4leG1kImr39lZF1n`huBf@3xRfTuigWemRX4E(VrcpDDx273OlTG63p zOozKT)8Q}+I^xrdj$1_0iSK(UxAGq4U%}_=7v|78Cu6$kUQT5t26XM+H@Y1jL6sfz z=t-3Yz1;Vm-p+NU&wiWf*UqV2|H@fbZ#ajwUFx{WIvqB+lfW&yd$8%=X56~RVs3lN zf$e7caA*7J-2KN;R=)akuV;VQyWtV`>(YTotlz@pABFLZaR+$rmug;`WWcNQ4)NMf zF}xwS9!Hhe5oD&IIbgDXwP@N?4# z{7&F=+ZuIQ|G^Wfq2ER}_FpTTJ~fdo)-S^6#Eqr(@knW_lOi1t^^iSwwUDahw$lCj zEa~H?lKz(GIeF|L8F0P7To`mgu2{HVhE5$S!_7)$q~TS$)wq}3-v6PDOR<-U z?fS`-OVeaV$~t-Ua2I*1I$IX?9V*W~SR$_+SSD{Iot5_=6v^rVFXXHG+4AFuz4B+f zPD%qWH$^>do1$aaRMES6TWOZGP-(e+fnss;x6-BoDIFHBS30%GQJfwxQ|Mlx;@Y4{ z@fa4M_}q?B2CuPIMvmH}Oc4yD3y zt5TDxqf$%kpfa~`SGjn8R?V4VrwUi=u9}_arBa=KuQJisR@IBtSMGhXQVMnFDCuTr zmEFzfD_fK#WkXqjvi`$jCE`^VCGy1@Wvg8iWqZAm%I;fnn8OzHbjDmBnD2p--oM`e z%c`S;XY-!lLzkQt>|5StLWIi#2shh zp62pr+n2cSDBO82?wy9aufhE{pa&l4g9m!i4*dv0Pt?(u(ddl}`qK_Q(m z-y5O#^U(jM;NU8MfB6(#+y)wPZM%|9+pgQt zc7ubo-45T8=|7NmrY6wNhglTYe=EgpccOUHiIkA?l#*NxDY?22rER%~zw^YFvU?4r zV@A*DRD1lHBBN$>rdJ0#w_+<@zTAbbb$Cg)cI>1FX7*J5n5ZVQ8GTrhNp-G8jc+p&3kk{@&0G^IsLIQABkMgxfA+v z@vD}6IdL4{h*-l9os9YAzE=FnTZ{jk87$Qbo6E-9_oTt{`LgB2d(tw0xU})!A|3VI zWX~NvrK?#(=@sx$4%Rp)N58ourv<9yya@PRV~GrlDUlHmhs!M?&*YBDmNI_LeVL-K zEipzqpYn^+@ZoVq>)}B~zx6i7D5XR( z4?3i@iN2}W*L$jTE!?i~tusn*ul|aUc7Za)MPC_vX_7MijE^$k=8>}eeYO(T(pXvl zwWSi(EK7-<7p3fJyhhpoGF~~@u%mKh!hR)J(^olD%gQCgK}yAnhRS_wBjstM0_Ba- zY~|~St*UzC+Nsp64OLoiPpX=n{H`(#KBqFaS)ww_Oi;B}zpQF&8mY3mK1OBVvahP+ zbt{$Qi-D>xT_32r?5I<9QG23t9ABYw$k0`F`21AWcE<))Yt>bi`QUacli_(PtGp+@Rdw$^L)E@dg34^i303oPD^yL}tp0!8o#|VS{rdLPs7XnQOhrWl z6~b|R&nr{O6iTKt4-qAlp@?QBl?F1EL`X=33}v1(CP`9MG9`1y;@R)@d+}^*dw$P9 zaKBjFTJH7Y{vP|We?G3^yw1~3Ywbiktv`3|bY>UYwI7sh*Qu?Wo#B`#b|#A!*_nOQ z#(DbSTn#wiN1U?-=Y52855xI2Fh_69V~4qxV7?5@>5qBGWA5FU{|<79u~X{)%1&<9 zZ{%Z#oVp;drLw2pZ^%m5pG1yovP*B9kgGcKohLgY{}Fl5NA6C@-wk^Rlx_O)2xe_Nj@zhGV}ovFBshcRKdI4*TDZ9xO&5LeYyY=*Kbi zBmjL0mYzm8pg+UVBX#u20=>eW!D=o@hqUU@x0BM2`7P+*RrFAxkD2J@2=sFYdfFC! z-H+bxL4OOSp`B{b=Z@&LG5Wn4Js*O;UqSC@qW{;xfjoGy7GFod1|M?bS-)oR;stIT z!H+yR8UUVFfh#4kzy2dQdja0U!QC40_Yxdl0goHOr9JrU4o=sDR~9aRn!&FjI1U5P zi^26`@SO|JAA)yXaBlnNxPkO*B@$k!Y zc*X+0*$nUes zj+-9N@f+~@eO5c(`YfHd<9D>Xhm7IWct=jZrp*}^?>W=2j*s5ZfZNN?fAs%H>-Mx!Qa=zj%0z8@3$g&$H6_xAF!d_t#7)dDIHk3jD5fN~LIb z{*&l1<*(>caar{63>Ie3cMGeyeqw;hA|c(zHIET^?cY;FxTv=lZu`x|-0&4*$(S-w|ekA;n-@b|E!YN9W-3BNA2YpIhu#uZE5X4^<6K4T<{ z-Fv0JYQ@r^U7FJ134^2&qc2KM(IX_+CyS)%akHiQ)2t+~WfhX&Id>`0K}iZ#)R#8K z8%vu9_m|?e+elk~HA%aSrbuZEmPr|<2c#_Tr&5k>vvjieC+X~zOzHB84pQ+DL#aev zReG#aAw3%sD7`GWCcXDsDSew~EB)D~BvXh;m8sOL$~5K`$=WIpk?GyOEHgMM$+}&F zOFqqznGe1yv&t4STd_#SFNwQx zd|Vu4erW|V|D9^G^`=p>V0~X%$n|5gFuftNu=kO&@F97!uwN%+VOniuA-?spAhQ}- zfd3Xd0$pk*H4yi+D{g9=AkSw z_o~eEx`k}!@L1Wjuo9V@I?HDAI@$cgBV}Hll4ZWzKgj|D^l+ZvIM+{{uM^H$i1X&+ z+#xvs0L+nwd1hm-nKJ(!m6$UO^J-)6j+oyUIY`K3Aab!sK1*dI*1tzySCCsG^6Q5j zm62xxaUZvFCHx_igMw8~YDH4@RI5R_Mh&DPHR>dZL5AC`h3SKcqnWo9Iyw^l6_o ze`+)OWsROO`ql}(+m8OVl?Ls6jXoYkFN4reOY~G1ea%L1Z=%21=&>C7oP=JtM!%EM z^Nr}cBYIzr{-=NgJ@60>E_#9wb#Q`fi0Z$IJ-3^{&w6mQ8ayR}t6AXd3OIWV-jcxG zIq)|F9O?-dwQpjK>2L7)9-N*6uhYP-3HXfz$1B0}DRA8czF&d!2jG1^iX*03w*W~URw;mWx;c=;X4<25BCFv|NikHeE5Iu`0x8)>{V09 zUJ0|GG9P`h4@m1ioi^njdw(#g$pxxi-6j-<&(eA1jve&zT-VZj`G~+*T=Cznv&F?N14v zv6F=Uv<{-P*BD{6v9B=A(h-(^a-y$&m>87bP7F4?D;y5)7h|o~iizWOguC2!F}JUa zSiHtxtgzP+0lV}>aAAdr=;JLmJ#P?k+uMkwNk_zv2z`;dQCIA18YeObtry4YvPJ&Z z5OHSh6>(v4gt)dDuW?TO6{Y(3M8$!x;`wl0@$!S2XuN$_e0{Q4{QjLLwK|X?DY=zM zs-5a2jcpx5y_^wcrE}WEx z=IoG0Tqu&pd{`+tgYn77tEA~8^rSg|KS>KKUrFAz7bU+ou2R4}XDMj83$9a%lQ!C% zl{V$PlVV4nkmA3QwDq*Fv?I|^N(pU|_HNoQWn_GlG9Rf*$2zn~xl36pC>nLA19@cbzh_} zlYU4|v|IYNW`OkVSGM#m>$=nwuvq#sU0eF-GGBVBwg}a@vf%+>x6;@@tPAwUFlnDYNVu^1Y3m*CFpm$bB60ABR2YVIT3> z%T4TOH}>R%ef7oObg{n<*kdI2nT@>`O4A*lW6xQVvt2Fr?tuNrq6gmSLkxP6fPNUD zC+E-?4fG}s{c%B$9HkELzoJ)<(62V=8Lo%zguDB=okag)(8Cqz<4W{$Kl<59JRkNR zeZ7a?9zlPfp~p_>vnP6e2>o7!o^M0nJ<$6>A{DQnc0?$EhaDm=xfy&s0VgxSOR)&p z`9rK2{v90622TvGg2lvfE#S-pylnw@Mc_{!9Qp}M-;dx@5qz3}Q#zIvrUd{DFH(^^Jf}v^d?d z1%5dU&jiCajquJ?_{RbsS`HtT!Ap(slRP{X4PQ-$w|c=}o8hrH@YzRrO%;AS4bLUP zckSW58u)MTe>?~u{`)Nd-xKkqj;ic>Y8|`c`dC+EW1f6Yjwf&X%F`Tec$!G(>GIEb zX6av^lRb~;Z}sMd5$kzzNE0ukmFzpli~Y4e^E%4sV3m9h{f74mueWj3v)#P8@(9OO z4do=`o}6s$$h!i^@t)gPIQ`vlKD25KA1V36IsN8w!NsXu*!weIxHO)xO-tZgU*GbB zLlyjZ>{foJkj(W(mi%t@NB%m%i+}x$5v^vl6G{d0LUr3kq1oswbcSyc`YPW2g&x^d<4Jm$_)NeI$NOUo8Hnw36i4O_Y=dSxBwrH%MxgXC=)mVN9sBLlVpC!R_f!vLbCbgCD~^BO0;OXBpVPY4Q*#4+2eYX zkzX36u}!NbXE|4ClFnFZnwf<(V`!5!XIiARV95)~)9-}jy|h?b?%G;f*(Fz6bNRls zPPSUw5Oh!q9;zmV?2LoAI!j?iQ>E~gxl;JsV^W0fL@B&cOA23}C55G@Nul9oQmD;d zDWXTcw0Y-4Y1_B=QbzJ9>16Xc>F&b$((A}nnZn2_nXdL^nQ2jkOz0KKCIo8A7VKy% zThnH;EX2uG7N*i!7Ua`OwoEfYHg)YjnJi|ath>T;S!?6&Qf=K$>6}Wsl%Ad`ZLy7# z!sAv*L4GDu;JY2t`tm9%C=cf}cg8ttaGp$@D+lM(`JZ#XJ&beXoJHd>hc@O3k%EW( z#(WW&b0X%=!ral)@`=xp!(QaE8M&N6J~NTiE#%b-xphZ==E%_rd0s)T4$@Ga7UVn` zd5=c!PRQR8dl-j(w8dWTVn12f(<$t0nq;8(8~Y2u9;>iVPwcf1_M3)1&%?g8v3C-G zci%-12B8lJ#K)e$(GPX>WIOuej^4PUKYP(52lQzKdbJAuN=45U#M#>4=v^WD_eA9Q zeL)|WqL=pQ=VkO%9(~=3-uj@wtI%U-^jQnNK7)SOqUWp7_uJ^bHu^sa91wo<^$lF) zfDbir;tXEW!Oe8=(+(Uh1W)DQ${l<~g0mFx_8Q#H0e^46VJLX?0GB6((uwcjv^97Q z0Jm=7_Zc|u1fDm7>j?0D3!GPh_Xu!r5B^`l15WV4G7B1-f)CJhQK36 z@QHv|4!|$T@Ju{>(-Gb|4*wj1hu*?R3*aRm_~|1&wFbUg1aHlRzkb7GPvNsY@LC`E z?K4ltz0zHcyTg0m|KmS+@W0LY@B9{#l$nC&PY_b__EL%Lt*!QTpb z*ual$A5_UBO78L~iw-;{MwOkue_)rqRy?WqDW3MBmZz_q!gEHBWeIQO-(`PP-o6d#BX?(GjIqpq3gKsT}=aNa*TrS<@s;)h`4)+9kYtg_T({uTI+oSxa z{Rttzae`1Du~MiO)CkS2v7)VOXQAh3Av)qUMwiJygi)*{OioV~W>4ZpAEO_l@21gW zKz*eUGZMw%o94oP>=!Ze>UrTbajS6tq9Gdh60yViutPyglFek;a&1YtVp~g zR!=-F*6C)7pcO+zNX;-24n8+t-y}BKe-fLoFB5SSI*Rz(%_0fkYbAF*DRvZh7rTSI zh&@A!MY_fhvHwM)IC!I{IDE=V96ixh96vur$_51xDq6;^mr$( zC&Y_eS~tbrU9ZK10qLUb`B+i0eW<9K=_hLXMu_@$wc?HZVDYYbm-r+fEShv4h@Y0R zqGjA{Nv_aWl7BEjQhXaAsVJD^_ej4bb?ZP$bC`;x76_d2P|_PtWqH{Oz=wxMK%Ych?BO>iB^da3)K{ZjWM6Q%C{1yXnQ$x`>X zzLJ^zRmoO2MsjSgDa{J?kphPi{9G#SZ?i!<)22ZxiO7`RO7^l=dCfA-h74JUmuWJi z@N${y79&|t`C3_bPfb}z#d|XKvp1yg)xV@iir=LCdUI)8bx&#SQ3YvI<_i2+HbJt~ zPLqsv=ShZ*Zj#}Ij*?-~dCABN=TdjWc~1S$x%TYEIsI|oE;#oBoL?7n%)&egnCqYU zT4T-&n70?^mdE_(ki$IWks)ag{e*lLBBy%fh5JFgosawsk>d>HX^ULvA>V_@c?$Ag zj@;WL|I^sR4D92dy%=CWCfL(D?CUc2wjcZR#vX5CpEI!6iQTcOVy!uj(D^xGdj-;cgmp!YiH|5|Wx1w8Zz7j@tx z5S%!G7jtm)3j8oQ>I|MV!Id5OdI-)Ifj4__mkR#2fx{H=SPm}RgHJthihIFLG6c85 z;CDPY?gE|@!S#6X{Sus8f_EQqp9B8i!2`IC^6(>g!3=(IfG4KF7Z2c#eE8!KJaQF2 z*#WOif?wvtGfnVKUwFp>{+R#|Rl-NP?2_9IKfQ&g6yYmfc&ixxdJK;>!e^)8wX5*k zeRwVozVn3lWdHFWJoxW4{C`iyz1)Uy&%1lM=lEA_T2;usV&k}1zn*ONs=SIHT9Kg2(?lJlG zBbRMB&rdd8=V$N6b3NDb+XuJ!qbt7GzjU5|<1mmURhcHH*rf=!4K-rM9eXjy_=K1r zHCinEHBKzS?|OX{wu%+I_ls4-3x)sB4kECqK?GgT5+SGBi?HMJA|iW)h&*vsY`UBx zHb34XV!tMdc>QpZIMhsRT|QML`LnqUhRdabEeLxG2_$%ky-_)qr2(de9dfNm2osJ$KL!Wn{*H_T*O7#2!`u+~RZ$SSK2`%%V-~q3d@b`$ay9GE2 z1TWk8x7R1|a{wF-2Tyt6>NfadaFz?+vcO#x-{Mc;&=Wka0+)Tk=NfP-2VRGO+g;$d zFF5`Jp6z+3{a^514$h6iy9Ky!1O8va16SaKXn0`}{O}2$xDQ|4g*VLLkE!s84tz2R zULp8J1D-K}Zw|A4(O>vy8a&h*KDq@jZHJ#0!&7PP(#Q*MqDZ|C4pKzp&n-Tddo>k=xHzW_?_DuOH~j z9hzpd!HOB&>1!}|S*OU|R1LWMkPPlo_L+OC|KVPTU$aHcMz(6-o%@cO&$es7@}OgW zENa4du--Wyp8JYNy!7HRYvg!Bt5kLk^WbT%mhlY#VxDJmiWglt!`>!&yaM-JUSo2W z*Js@45Ii3{+-fODI{)S9C7~Sq@f*i~|HE4sjOXpw*7EM437ppJ6z{vZjSr5+GnBfz z@Ubdq^Vzrp=L8!M9d(^&;>%Gl07H)L?j%q~K1&vqfd>4s2TeF*n;?ikB5}=ll>8!0kZUJ@yPE{WvkCb89_ zMI`k-ArdbJiNrPEagIee5AGr4l85ts$2mQ4-fK9wEzYlrIf^mQC(NaW`4U8w;}6Wc z4RhyW{`JT~7kR8jE*Z#Y1#+5&ylj!%W#s25mN&jap4rHCEAm~3oNbV|kC=V`1@iBW zJyc*Hp4f{k_M?bBm1AF)*qa0Pw*Y&5gneGeUZ-Ner?6*5?AsN4KZ*ShL=Q^QhcNU) zpda1QlU?Y`ar9<8`ZEwc8h}0pp;vb3R}y+=iM~w{?LK`IZGV184-cb{?n0yIC!xk4 zMeDJ+U)mfa^fnp&J&7J~{fHQCKRs!yZgTKAtFcLf-2A9LZrw=%_2d`hiZ3FoA z2gl#Q^8#?40lqJQ^P}MXAMP{2|511#3qCjpFU*G@+~J8N9<%y8yfF{{sD?)h;FCaj z#T$NUglA^LH=W^~5csDK9{LF%HNs1K;U|4~sug^-9^P_;zh=N=&hXg}cuf(0i-6~f z;k(c9-jDzI4<7t?8vehhVx_*{SusAA6>x2ilH3SZ42)z&B}-Q7n#9T?n3bxUSw%&O zTlf3Uste>;ZKo5fKUC*7tDmse$yuy3#ev(U+~D>il{-ArXM>*Cx$~v>Y?wTOdzcSm zlcHU0hTq>>Do__l|Y2A5owJ!T?kz>DVJd>oYK5wva;t>4;4lD1^8#m&4dS)fOxjBvF&e?JT ze!sibaUXBjc*Q&4_vVzsdQMvt%lk$VAHXvR50!P{tgMcF%zGW@beqF@<@329Xf2=a zy_^dxFY&o`b$qdF4PP$U$JfUD@r~DO`4;X8b+=6j-%mZrr8euh>~b8FuZJ}akCt5!a5^AeMg+`~|Li6%+p|#9Hv^5e!=h0fB8+Tc>ANyK#zd;_ev-XeZ z=aeICJD(5(Ub~Ay*{=n!*eWEeSRuo;b%T>`iJ`(!3@d*o>^&QW!^h6T(Q|+pS^8Cs zwo(&g{LRJK93|mYHeHOXR~6%*wG|UC6$x@o=WDfc9JE2UuN+^RqC0RBz1dTF6j+d zkXnzwBEEbq6p!#e>`dS+vG<;_h`BB&0zVXp#g;wA%%I_7s)vf0v}vet6)j?7aG`K{ zbp_{mf%Bx|Tsb(OiWv8z1?L?v#vcEU^LG-X`~1c{FEE!0=2OR(AKB z9qea1_JrS+YNlgv53#>&>~SIX*$aE^fc?(Io;zXR%dz*R*uN!ua0Y!aK`;E!k3;Cm zI`m~5dZUc~SffWt=u;JXH4*(%LC@x*Z%5I)AbxIDiyq!aAEVLB4D@pydO8<<-HhHY zLVpX;%=(EI7=zXv#&03Mvcg(vt>0VnF=%5t0bWah-^}5;jqshrf4ujP z|KP#@HsinZPt@40H#L;Jptn0V&^xn*)S%l$Z`<9ah8xDz7}AB_Y5w;CGB3bPAXasEkG4nD}OcR90KF0%WB<|)>eInV(ZUQ^^tK~^;ukh3-&OH5B4bT4bm*>AZz@F{bvbXV6Uf%vZuN+jy z{^A-3&OXdRx&1jLBbUSW&)^6=uQ@V0oj3X2;>~mL+GE&g-eR1_@woSAV(BB^nkD08 z|Fyis5#KMj!8NC4=Qwry15TU%j?;~P@V@#a&d3_e2bXW-OrgbD$`|>_qeMQI`H_!@ z^yS=H3YC;gQ9Y4~TZuG+^n^oZ+r+H?J%G=7bHL%l62{+e9D-!(Vz&(Np*`j=?+KNH zWYPM?Y@yZ$_hGlS5gJn*g{I$H(Pp!w&`Rwlv@=(TwnsB@-*$7+E_IsFjoK#kmJAo| zhv7cFS{sG_{Xe2ZM6T#)JydkO<0TAc8i-EIt_#DsWMQGRNQlE9gv*2_!po;hgyVjr zNl7JQzp9SNUD8XOANyF`{CZfFoXZoBdaMwSE?J3E?_1*bD0y))CPd`@=_dBCoGg+w z0!1jkTVE`qgmcHe0xyq+)hHWb+FgdvBeF$TRV&f?KpWBN)pF5k_8nnR`xEDAhx07J zxo{79+&fACeiP0ch;!qaEK3M;L}Q+Dn5!G+OTnBTn719~7MR}$IWY1tMlQ#Y&qU-D zg}m_GMHdg`SA`rSk>?fU`V#quBj@qR`#f@=i~NhQhq2g41ol#e{hY#{)?r_v*jq2` z?+o@BiG4oeCw&{Z{Ka4Fc|7*5iM=nv{$tRCJ?O&)^rAibu@61jfWF*7ZziKZlhC77 z^hpoBdW3!jp=WnFcV;7c*BSk*L=V@ak4Nx?h+pVuFnT%&eGNiyU!cDg=<#{<8TTt} z&olmM=>!Phx( zhUX0#4+nP(!CxXcYy^+@z~vC|nFmfw!Rt_P8vuS|!EsOUyc}EyfNyzl-pJNoP2m0- zo5kXJ?JD~4ffl?l8GiT-PvDsxrC#uc3H-4G9ytx4 zYFSSuPp;F0Z%?S?x-#9*yhtUlCer=4t?2>xr4siXdT?$UJv7y)vP4TNS1PB+>s{$- zt0=14;7iXG4XAc)52}Cll3uOqMsF*B(R&_3pW^H3Yuy9-S!O^jdM+&YaU08T>cEQ6 zC04QA&#KETSYxs`w+Y|E+Q(*aJFm{%-oS@DoH@ju#wl{wsSDWfej@icIh0MKquK1t z1-4u)v6c5g?(62kwo$Wq(AX3fV`|y1KmJ}AJ)Z3y#<63s6drA!%1%0KczoR^cDb>L zCv8pQsY@%_-TEod_&I`SACc#IbEdF|LI^K9V8BZT9cS<2r@YLqhF4UV@=7;dUj3*Q z2e@?SbvN(uhT$$8^lS--bkX8a_stxZvV_B*+ww-!o*e13fTIdebF}&c-aIP__vskI zaqar>7T+-(U!Kp2LpE{J(Q~}5n+qpz9nL$nVtHq*KJV79;FRP7-qSsq({g5T`tbL> z?@15N@U7s3+OPRgP7glpI+L@0dGgVHTAV$uBOm`gj&l#M;C%OCT%Z}qCofd;>6O}C zXm*#2s;}_5#GQP>@f%;e}kTn9q;bmT<*eAFdqVz)!RGa8>gQem1-vKgVmKnyfUgy?u%6p3mYJ z^-XxzhAzLnx`nxwg@%fW zP{40Jel5!5&v!2I`{m~RPB)g{UUlF%FWvY}FGqguIU46UhV#tCx$xVZ`ZqY|M4a~o z&Yghs`(qA8%<~X)#bLfDm~$rPJ&(EfVg5|y;DbD3k;`u6V~(7z@NJb($gMx}Q|Id% zACTu{53%QK*f)L)H2ewnKN>vb;rSx!m(agV^l%dTsE=NnqMy^z(?Il94!xa${@S9)Iq0)8uP%GZE2sTL&rhT8 z*U|e+=s%u^dq5sMECv@tcy`u1aAFN!u7R6t;3o(ixq_$F;A#Q*+6>O3!J7=+*?~WA za99c+W5J~}_{4j8Z*%Z!4Q}(nZxA>h0G=bkwL5om`UTGSfOijY-wymwf(N4Dg9Y${ zG5nAMPuRj2I`Bq3{4o+9(S}cMz$;JTmqGB%XZWTAyz>nH=?o9)!beNtrBL{3JUsOP zzUl>U>B3)=;jv@znKito48J+Ua~I${C3vp_{=4>%2mjlQ|K2}Q(VNXwn4wE&R&1fu z&hzOE9ih_$I?)-cHdJWYmd-ejpu(B;RJ801o!z*P&h0a!i)(w+rJa-L^7&45ElSYM zk=0b}*O~5Etf6~P?$d+NR#ayEoGQBF8Nx^NshV7=R$&FzpL$1c+ zKS)1D-Js^g0G6BGljS|#SkYrGD^E@1*899zt@0IXDx6`h&c58XXE?W8R>yh=ZgU5X z&TP>6C3l{>iMs}EWy9P^@@zZn z4i6mlk_q4A2<=TQz0+g6qBtJ1*^Gxx#P>f2we0ZpH9PJ}aq;ZF>@Lrka+qh?AL2O!kMlem4PIb#KR)`FP%-glFa+tl;2#t{n2?Glv<+ zb2y$?9T6c<)KHhe(! z0w2UZOAe`M@}XwjE9q@pK79WMXJy~vBN12m=%lH9tW$qJ_E?>>qj&Lf+e3W(5uO*G zFog35IPn<|E51_uoJ-q(;TPd@{G(%pP`Vi}wCwO&`&56?b?#cxLuZ6A2^cI){0@m8 zN-Vlg`yuonjTBlHGlh!Q32yGOjo4@f?oqOPJ#<=27NDzgsY0 zE6nMOdGj&13FbFN4(*YL3vwBRe2kFO0_5d}+$JNxOyr2q2)895*9*vZ6(=?SLf&!6 zeKp6sea0RH_7Q`<7-2t_*i$X`wF-M{kNthX9<4YW&p!?`4#9rwuxAJCyBK@#f&KfT z2ba)?q3A_V^dkj5xz9eE-=H@t=#L|MRD?eHp;w!Ep7m$+EE#=^L+}3S-wE_^4f^PV zUY^-Z$Hj)#AP*kafQydc;~_XH z;(lvBftyj_rx6@o22ZcRReSK&0i0=rx1-?B3;acbLwr`=8PBEG&<39uz-a||T?=kK z!S8W!Tn(P*gX`_!+ZUYUHQ;PRaK8lnpMeKd*`K0HwiUo3<-#=swk;SneJ zEuY&}<+b{BfNQ4@;yh-2}>Nm`KM=<0!jT7UlHaL3wT&bRys&olIOzr!(@X zFjt~;5697k`<>`=$!xlI-kfe^%%j^LUUaYfB`PU#r!xCgDlh#;l^x@#dc`ER{fT;{8L+2^jpA6Ek9Xh*B(~w zmC0(C99Uzt7i(TV#9G7iS^IlhevyO4a=ZMenlN0zdFMi1`reJ>jfOXE&4$GCI( zPwv{qg1gOiXT$xc+32+j_po$k<7IAaa>9vCTUxl6Bi?&%J;vrQ6WEgfaPPzn?(@-y zt(~v4%}Maly*Kxd58wexi+SLBLmu=yia9WqCB@+^i}U8eUH|Zq+)5rcx`uHd5q5~E zV@LA@9$8Y#qu2TKSo6v3R7pHO_5(XR3}cs;U${27nkTQ!=BcbxZQHharf^3rendD)S#ykgcL_S2fmE6?B_ z*K-`%U%rI{_6+BBwg))y>QUY>u8xBid2q0aKL=OieN3VphYtJAq3;bjEU7Pt_y57+ z_lYAGoa2ZlSKhe(CP%iyy#?22aOBsE9JQ2q)3IER*}01oicay)j=sFF?N&Z&?ae29 z4d4seYJB5ccfO}x&kvptIzGtk;H<$k9OX?l@blzt^{ymouMc{k*>So?K zSC_ZiRdT%VT8@3xnPZF}@n(H(j>hZ3O<(Y=__!wCWZi+I?x=Ird^4P*5a-#0b2Z?6 zmvBxeocAQo{m=P7VGaw-(*bjR;gAG8ySsWW=3R!l@%ne%edMqcc_bhgE9CPUIrT+e zLC9@5@_U0EpCQj!~#Cgi;qxgSLS_pyfw*hc{Nl7{{0U{9T}FJyXeJJ^y4;qat(cXgWgOdI24*Y@Fo4NhrlOyoTJoqIIp4kK6B*Hrf;Ga$K(0KSL240GTpDx2w z{ot#!@YX>1%MTt4gU^n_YZEEM`42pI5x#p5@45ZszyH^Y|K5K>Tbnyj(z!R37{7@U zKMbXWr;8|YdL zN*+QL7U@)(7EV<^uh4U@qT1vM^x|o+tsFtxyy#`M_`IKElhO)Qw~I{pOgW0N%Xy4{uuF z!5bYkI5goAZy4E)SL3>B(7=Pbc_@k~_Jn>fD@=3vaT3Ue*Re7=~o z6!SjE+zpuD5ji*@k3q;~8S)v-!v}vuUUkTAGVsQ2LM4!H6KU(>+Z1Mp@C?rwpe=u<3#X04_q6A@44U{pDEaH2lw{iza2af1s}AB7q-C<#qh)e_(Boh_yT{_!XqZ| z$xL|V9sJS*o|ypO6u>)i@J~8CB!iERz)J#tnhH-Pz*p|@Ry_Qr1&^JE&#uF3KJeQ@ zcy0)M7xf?S!GHf=yZ`T#QCN}*g`FNmp($z<5?e|kw>wdA*&+(*7E2*s`zUzpbqfA4 zkwV6NppdKr3hSy#VH@zDx_pWlX-pB>VHBxmOi|MdDeC%5if*@$Vl0MI?ED~#+y0T_ z!w*uzk*Sol{wi%7dXbX9CDD%56|{>3Y4@$`lv=)x((de`^z2gFAEZDB23OOei8_?| zDU!0DRMOFG7s@UiPB{x#Q=WYtozNIX1x?;`Dro?nnec^*jONhU6WVm%?GRm58%38k z9j7aLs&sYV8Mg-R-%T?uE{!`;YzULGL$I8fr%m%cjvI zQ;Ev`Bk1wHarDG=J5>fp(bH#pse0H*dX||>HF_7RHsT%Cwd7NMfH}Qv9!jsn%;;^~ z`_z!Z)M)2H@2e#G7`&T4TlmtKCn?kv`<=dzoJK#Dx6`kiQPjMpnf^?_#B%K?aw|FJ zRt-N{zG5XS+=ydEJV!<;%ZQazr?E=HQEt8Q7pt!C%xZp(tnTH=8uJxc)BPK_nTWq* z#yn;1;j6eU?_wRBj@-`l3hUzgKfN|(++MDN+kYC&`p+J4hugEbW5F>tNR8x9p#|I- z*G6}7=*?Y?l(_3pC+>FVHXEjnVxz^;+`aEK?%w3jJ@WF|*sGLH4EM0f-C8zXy@Y!j zD09z(Q0_H0pL;#^VY3CB*zDCaHlK&UnIDwT>5A@8aFqtW2%K#PZ>_-HMDSM+4ts&e z+2B%u&vJ0O1iTIew z{IwSz+YFyAgx3tVJF(@BQya{CEFTT6*U#c}2db#k24;@$tp9X#IEc*tMBF zYA?`2t0lBx<0x9tpiCaFGs)xfdh%2;AFT8JBdNh-KtRrYSX_IgN zd$iK{2d!4irhqT^v~%(o+I7^Kc7H6URAVjL zGk+MRt@Wh6dAlk7Mr+zvr$!la-gH3w7#(b%M2B?OQ)Z_MI^3?AvNVU#5#JLzMr4b#!zngMygi-FLpOklg4dtJ^OD8U^rh=Qp=;Xr;I`zVsPXAD(Lha2| z)cY8n9Un#KR<)<|dxz155*50n7E{ zSvG`f>+9%+`+2G_NuXDb2K2f}kKS72{qnvs)TsM|-iQ6B4_{i-r+!-c*U*V6}ff&LsmUs!RiVV zStC7`+pK~IU1GTHuXo(eIgxey7O`&NLe_Qo#O-#^W}RJ2x$UAr*6tY2ZOW8b^N1&F zrtO=`?>#$njeO9ym%BtV+S>0)0Rt;FsssmHG^@j>>o&JzpPgsk8e>px| zXu)|DaV}+?&k*PI!FhFYZY!Ms1?HH9c@APOdCa#7b9!RlYnb~c=HHDR?jesX;KBxcG=h_2@Ztk*!og1(IGRoSYMUwj`Y-Sm z2hL`Iw*lZz9{imHhYsNJHMraZK8JwQ)!;Q7+-?QG>%g%Oc+LXX1ir1nc|Le&aNhy^ zKZge_;R87eRrmux7{C+#;EO)+h9dm&79KeOpL~Eql_-i3Nb{0O{1g|y0Z#UpM0pB(L<*iw%wc zeu120cMZwe_<`HU2aQ(vsTcCsPhzbBbb7V z>L_HtHid4eqOc{A6h7(;McB-wjcQ+MV?`&5OnFC93;WV0yk6P#<{?Gzccaa-R#J>^ zAjMoANwKT-Db5Ukr&K7>mW@{^evm}*Z}I*wWgsQG_)(JTPD;9bhPDQM!RrlkN>)5U z$+v#e_GCHQvETyj95j=5wJE1vZwzSn^+A-9sY|J`S+r+$9i`3BqP?ySls+1LEnxfbDdUyy=kLXUcR`pu^6=lr?h%9r3oJqwDqQ*p|nXoiU1z|MNcPK?>zI zSWv#=WjfJqAQkX5IytGCPAz{=r(^!onan>_c(s6vo`0pYzeDJ}o&{a#eTpvHf22!p zyXo@MV7e0IL07kE(Y2Hzbp7CYx{*DOZsuiEasENNmGgscXU(TOX-Da9{5HB5SWEX8 zOrjD;T#sd3L=T#oO7DE8hhVg9R&RRLbvQkGbdt&=r_kfRpXl+UsZ_CICOuJpOivQ` zP-W*sRCyBbg*NP^s+coW{rVj}U(i4`uYXhB(697jY8$F=bAal{6jFWnPxNAcU#e@g zp*rSKX%Spl?*&X9iV08$?yH)2V8ZFFmc?M^9I) z(Nnc4RJk>PDmz5r-@hH_xR3Kp!nwxde6Mj%C!F^v&P_Of56qE_OduhaeK44GP*jEGgHVXSQ#U699&raBD4EEa(dnWArI`*!D{dYqT2A~hiD09RY z^rHqnnSs7UqBp_lk0*LG1AX#GuPV{6UFcaU`nDUrQ$qjD&_gx!@dkQ19{tRr?b|+~ zuMOz!LG-s1dR&h_A3?7-qu&eBb940F3B9*M|4&lf;~(I`7F?)-4*^b=gBMS5GY?2j@q?`!;Z21OD~kfeG+|GrTYcesG2-Cc+m+@WxE|BOM-@4xgmMD{J7F5_l#N zzKMo+3gI71ct{66GJ%(z;HN?Glq-Do5Z*chf9-_F>ftjbc>omZ+js~;|CEHIE$hKku^`9?C{b#o(+uX5ad)%9B z^M{e`*(WsM=Rz8&5=Mh;`w~rAN$fL+MBq=7!upXcMLk*C!oTHN2C zme^R6m$e0Xo2AfF6MgdOl1|GyoTTOLm1#x02=di=PJZnV(aNsBX_bX5t(Kx`jjJyC zFV~}h#6Pt5tOKok_mI|i|3(|oio6j*Q86{NY0`X(?%GJvZ|BhFllv(qBAH^RbOmo6D6V-SUJJdY z_ya+d5U5OvuI`j%iSJ7kYiMiLP1;t_NXc7z(soa9JYU%rBB{V`?{6WzLz^`|KT%~;iX0g%-!ig{RTRiGLjBC zM^L7`7G)mGqr+21(Ba=Hl!fhP0hjEk#zvLDlrQEQglshn- za$c;ZoM>yxu_(m9zdz1Vfb-12xkll9J8(`Poc9CHZHMy*V~*#TClqrfVLln=%*4D( zm^%>j%aFrL-*yA`#Soafq^~HXxuxB^yTOE6Ei~Vmx50cRbSM=gJ z`Y{DPDMMdYqBkz+PXu~ojy_F6uX>>I zelCL}Ra&Ui1g?_6mp3?j2i`)#-Fon+1`boeqXoDO0iQQ%TKF&US`Tg`!0#+@Y!9CI zf$O*6I|iH!8f*9++(&}{^YB0|eDHLmmEj2aj07Ctu(dclgBz zp6LzWY=C#>!atATAv^fU6J8nxKPgiGc`fi2{=ZL7@K-QARs)|^!)syi+j4lW8NSo| z_gnq{KOvc&)TEvl;>hIr7&7s7ryi$9P>_!)G74KoMh@lF?WrDhJ3gEYGdGZ7 zP+Kybu0@8jATsLSOh&C6$f(wVy5~ewkJYEi*!CP5zgtZvy9>zFaWM7V`HXtmwk0#q z^<*C1NEWFEWVtt;dT&>vJ|UT8b$aVi)T7f7(0Cu&TFyT`LG8Ap(j4Dq;sH zh>T~b7+_#y5O!b-VqgIx24W*2U|=gs+SrPK1sIqpHWt_-DhB6q?aybQz0bb(KJQs; zeBk0`GUtDcasQtCnR5cm7y8mpw2*BYmq`2mCQ^)bmlzu*<)C4*?bwg9-Kb#MzNe#f z;LnUY7#YY8Z|r5qW5KdhXtwO!?XYyLwnaK#J}#Y>>y~E_)bOmp#tTl09eU%3du-%3d$yW$)e5vd`phvajtt*|+Sz?3WuLUE+Jl z{vpfdfXPkdzyY1*ppF-$tL=8_*6f0GZ&FWsG>edfZOh3a9h~LRe%<7-v3KS0z`Jrp zY+X6>%u+duKZhApMPGV4+>v8_s>^XZw#f1P`PIZ~*W{#uHRa^+w{ps@n$o-GWI1(Y zqV$Q2m(xC4$rn0a<-+moa2)%=bkzt{q!QG|KOT3V7IQE|EZx|;LuD4 zIZcs4FK5byD^|)yEyCoYoH25-*L}G}H%~5!lQOvT9T|LYtz0_2vkcM7lp)cja#`b( za@nDma(TzWa{2l9az(GDaz*xUxw4m?TzOtcu5uVGSMAr4tLf(I)fMIHFJ^Mh;2Ltx z(PuKW>S!4{#b1V=t|8az3%S;fKY!cuR<134D#K)58Md6SA8(J6;k6&haE})Beaz;5hH+mHxIYu_vl{nXp8Nj5{cG?X z19+bDJeL*E_m$^d!1Mk&cYmJ$5jh+rk5A-sihMv$A>?(1+**;}HFCU8o*v{HLcWX0 zc@lY#CU;ZvH9LP0vzG$cPo3CPChV(Q>@8pRmtc?0 zXPYj`O|m*`!@Xhe@9l&w3OAJJ(JaX4V5Nq&Pb!@+ojQTO=`$m$6#WsTM8vgU}z zvKHT8)XF^|&E{Q@=6!r+?aFs$?X$kJ&e-L$&M$LWH))u(=zCvU6jYV(U2 zXPK4Y*dTSIiWXYW9 zv^nW0ZD-$>t=pu_*55uzyR&?+vvQ-fAM#KNt3rtiTO}SolJaD*Y#Zw%+b#5$?I$jg z4gdC%eH_3kAlBCObOWB`)BRk;hcsZ~*TMqh=E?tXErQ5UB z(mn69^tfUq2cOv|ha|6;LpMaqVT-58;S--qN+VORsiJW+Rqny+?UQWudE+>yUAtyiel~aaV z$|>1J(hG;B*On*JTQgXCdk&M{r#{H3RrtFd_wsUTn7y2O@|B#L?<=RiA0?-L8!4xL zpUU_4_2twTIdbaF;d1JUCvxhhr*i7-VmY-(nw(m-OnTp4DZL}zq_^`}=~bL2z1A<5 zUiN3XUIf?W9-ex0ABVUfL+)!E_ZPu^>T#a8NvMcGDlaKr%>i9kNL`B&Q>ySHq6~)=IqyXnciNZzflK$NqT89tmTgc(PX_*f0IrGuzlV z4cI&7*+1s&p;hdoHtZ!6_LC)hYAE|EmA%!1{k4xhrolej#$G$kejCD`d-|vE{{F50 z@1Iastlvj!P0*2=r+!M!hNootO{KDY^+J>lnTy}PqwyG1t7V1sGO5)*O==~Kk=hOmq;^bztZ3X`Rt$=e zIvR(h&ae$q=fX)@$zhVs>c}cy>t&UjS<rmr8#>WeqY8Wvk* zjk8x}&9R2EmTnhW>!6D?8(Bk|8_bvHr$@`$(|5}{jjd#zhkU;jvqV}9GM4p>#>jft z21v`bXJq{W%VdLU{K(rR8*bhv8;zYN8#j9{8-Khnn;dH@n+CBRFoYeKg*8# z-(<%H{CImqc5?HQosR96ovV+RoqaaS&ga~uqv1{I=;16KH|9&ndpgp|=#_MGnkAjQ z!=zK_C+U>H=e0A&(ka_ZI^E#&+|8cS>8g)($~YyR4mFidu{)&GLVkb5B$>9}w1d>Z5^3fwFck*gOZn@;=PL3bRGmBiSkZ&P5 z2a~rhx$h?bAkHC!^H|Th{FzTI=VZ=#<#29%jXKtobBy3TLpfJJ&exA~-p_e=;oJvt z{uP*m+04Tw=Ath15yPAeWL|uj8(rqdi8-=op2{;^>u47&&GqG-VLXJI`pt7eLP4n!{}!MJzYs(?dYv9 z{hdpXhtlU`^!h9PUP{lC>HBJW--`Z!XAhK<<=&RE7k0BB%Cjf_^u;yy#&GtB4tvCm zed5Ai>B)Xs%br=vzLD&mf$X2<>>+dZQ6_t7Df{V9Pr0(MMzFVLvAxc z+;s&HuD-{D@r7eB%L(39B9>2uTbCLM9_XmDQ@n08M z-cgq8G*Xs3y+xMq*j1K4*H&tDA15`gBuP#8uTt|_b6H`UwXE>_xzq}^m)iBTr1l{p zD?0C!6(0_iIP;}mJ|psTDr=bJP% z9wQACJ!O@y_hglaJEhU=ZPM5vSsL?i=u9wOnq)ndRfpe~Ro_WjE#RIs<=>;2Mh=(N z%}ZtV9Rp>J#$9ENJ^N(M#xrHj9aCj3^B`F(e5EwgV&BX%lx9!*Nb{ch(mcgR*5==8 z)SlW=*1k|I>(t&O>-adyIwwNTn&>(%cq>zT0w4LeDTUnSDwpk};xz7#UZ*%TjoBQv^ zbLj9q6L~JacQl{JbB^bEr|{f+czy?R=t&-*$Yl`uM3B>U@@hqHx#X8ejuXgpE4lvR z8%@p|$ond}4<>&P&fzWR(T;N&&G`gzPDeSfW}Mq+&aVyUXux@{mpZ=RIA3SZc^Bv1 zLTX!-a{kQ1?|aOHC39iHe7s>!wlOb5n46Z&Pg~|Ff_WOmTs39B?lNZu%v&IHcb)mu zXATE5kI$IPT;{VUb2^=Q9l_kT#utN+%<&QCc`bA8%zSGy=Utfh?#z9C=D!y`XiFdR z=*2MlafY6_(3gSqrX~IHrblb((@lCck$!pMjma^V^1ZpubQy8V%cAT z?6GC+vvz;=+Mj-7&;2tR{r`MQ+!^G7n?1VV`p8+hrX7l$K}E=}^A(r2`QYNJnz#_O z8RrAvB6DdmGGewP{ooUvySo|bIsEz%ATLfYj=NXzSwv?s@KF1!}b zE!%?hY8`N{%m(Q`UP!;%4jC3n$QT!gj4e-*X>|mdy9OciB-f}u6X$no;QT{BTyWZs z3lmClVSfu;G||Gv&MvqZ+zJ;@uEoW|F1S?X1THz|;*z&1E^XGrr3=4t>1_-yS1HD2 zxdWF+_rT@V6L9&sAubm+L{^>g$a2*{R#*jOUH3tD^>xVhv$k9th&cx5i$sUd? zZEoYrt_8Si)(%&rb#Sf5C0yG&8P}~YIgYX)w+ zSHc~m-&%Oo{5u}&PsZc-o$=(R z3!a`Bk7qjvp2@5ix>dw$do2_;T7aVJPI#kZjyKN<9PlGI2K0TPz zqqrLNmAP%o{7z+#GnwaJ%=J^|dpC2wg?S&z+y^uNBayYX7+Hfp(F-^FQH`FIqc0xx zrY-%6q({%`(+Ya^0heOg1>QRJErZ@A(?3BE7t_c3^wJ3zjla`V3;Ozw-bT}3ADrLu zl|E0R*Aeu)Cq1Y0Tk`1rX8OOAJz&c|7|UL0%6_nCPsFeOI)9vq05*uQ2K_64j)if?@+&$2_(f*nYVa7W_G z03_@!MZ$*Gh|d^{_>_T&kJy6v8T}F8+XL~Ik%$jWM!fSv#4n$L_+Fh6KgR~~QFjoZ zP=xr@pGc_G1qt0#kzjuU3B$W1Va^jI9Mwa@$!PWSNOX=sV!z%<9B6<<*X>9g_8y628z9jq5Q+0{aG&*%xcNB}kC-Ac zS0eG#CL~oqgCxgzB+YO}Qru=FJ*|dh%a%xuxt$P-!^y26+ z2ORyKjbl42;drlWIR51aPHZ2GlP=G3s=PZ+9gM{32{UoV{0h$82*BBueBSJ0jkI5U zUOnoB^qEbOVLcfcuk-27K%5`76Bo?*al0`tuGPS$9y+*OzA`SSzC@PyL}b_VMD~qj zu_hj8}59Vin}9D;_gLDs;aWjCkE>UVl2*Il%R%bKP@XpU&kafNdk#`&0WPIjDEF6UO4^XtkvdUKv1BY*35}7^Rl10(PMsgGe^6bCoAS^6!SHjIUC5l4QB45 zn7`~A#+C_QLS zA9U!&UHV}{PtMbq>GY-s{dr1{qUck3dew=3h0wFJ^sPU=qgR+l4_)ZvUV5oTKXd76 zCVf3dZ`;#f4SIZ+K6jwk2k3V)J#SCni|PF{`hSf*(2RZ1mc1~A{ov1@*u%b9!`>Lc z{^-RX*}y*O%U&7Dez9iHlyIMO**gQ-KU>*Dt=LBv?4_0Lrxole!M@tf-qK}%y<(58 zXP>#V*ZlukP5vJ{E$1{Ed#eRv_eyWkD+G3Ml6gFyf z!v>8eSg&V>m^+^lT~vnXgEok|Q3X*^?-Av?2~kxvvF`jgtn+umI;&_zKHPxF;(SE@ zNhq*BFFe4(rZ5=dmToklQSZ1lMrc^ zg-Daeh}4Kc#2Y(ATz`OwGoul)wKgI`?;~QiE+R(nLPYOJh-mMIh!$fIQKumys!BxY z^ERDph^XX+h|0l;sJ;Oa4fPQrMj@i_GDJ+~=am|W*vH%Ny+CB8--vYJdb8~id5~)s zufV#tow06-4c6s$KvaX?i1KUCGk-+1O&>&u+akJT5n@KvMa;!Qthdj^`nYgxFl&ts zVIkP4I~^N?bFoQt1vV}AN37Oj#4cZg&4$CVnas9W)WDYbRBW|R!q)WL*fzif+n#j7 z_UR7TQT{o0MAX2}1_!b8NEmiG^W$D9c6&|0?q6RK7g`g0>Ly@M$};Tjum*dxqY*zi z2=UJwBEiR#9`ZSEff-)D>+vt?;>&^M*d4UhozhcKN`ewKKh(f9_Q7HbF<_8DsqmOInS+}YXax1$2mXc zy!&(Rb2$Im%)wOVA&|K^&wPwyPW+jdM$FAp=I0o5^niIXX0B#2UyYcv&dl3$=B^d< z*NHi_WgZ(cmu;BOK<3npd2Pbno?(7JF~|3q=TJnf_{w|_V$SC?@2<>!dFFo&J@BIs z^hYO?e%z!dE$K@wdgDQVdefs@^l2l#Do4Mr(=!wL=1%W;-#5+a;R^a#hhCnepS$R3 zC;CcHz4p-G8T9x(eYWJsG5+u6^xPDYCx6lVF#69v*l~h=u$Ozk&3;(No(N}ObYO23 zBg*hMdnAT^Qp7XsvtO*(Ge6lkAJ{uJ*gp-~Lr(0YNcNIoKe=Gb$lvU%7woO(?5|#b z_1K?2W3Tk$N}zC*BQWdt)nOIkI=;@car*v|=z&2D3HZ9gnZwZbBo0xW9#4U5bg zga5rP76y01!Zo+C(A5bGYgEC)3f&O&;sS!wh9PKOEP}>v2b7z7RK3O z;jhtHY-M6|axqBmAT zbfy=g-}FF?VH9FoIwOY99WkCq5Hr^kG0WW$6SW*MTS^hLYZqeT_`m;r{`2*p+m{dG z{YLZthxoXGe0(=vCz99e&g)L%^|iT9Gp<*_bvtqWe%yy8_j8f^a^?Q+bDvYW-}>CQ zA@|SM2z_55;`?u&E1u`m;5jXLUUQz?nCCYmhk@jAi(F2TPh)ao)<=9Iw^ih4Mvf=Q zb2+*0B;O|FY)jsc$h|N5Z{!>bIFD|e%LC5mBIjhudFgO&(VX8K&apn{>CU7r^LBu_ z+ra#pGKWK%$BoRTKJ)p6Idx`Uncai2%Woaxm)`c+8JCeXKRdgn?1<{HF?s+O08|`o03ZYlgzBa~LKszK=yTF*LNf^E0Hby-U!lC(V2_RP_*xqtQTg!TdtUd)yWt*s9`4<$ z!>#Nq+|ujAZRJY14K0S7T~oN3?1yXVEx5jZ2iM0B;reJMT%RYw^>YZ^bn@WV&<}1s zoZ;qM18%#{!tG&IxSJQieT+8TlbqrHdjULL{o!$-G6q-G#bAF&41T!`Lp)M2ijiFM5bE4HGfuaxy%p zAI8{P5g40U9pk1~#rQgdF+MvL6a0T;Vsn1vhhWm`NtoPWCMLh1jVW6@!^>p{yvupR zJ7p@S4xbGlU03)VzK>}m^DtefE~Y03V#a`rm{EKMGb26WYcm+WSxYc$Oe$uTT*T~< zRhUyX4|6sLV{U`Vn7ex<=J7S;yx14;GqHx>;u83MFoVDQ0{9;4YCHmS{So+bBmztNA+W4F7HD++>vLd9 zU*0}~-=D|(?cx1D@^PaO=v2n*@VWe+Hm|#e*AL-3T3nAE|go#&e9~dHi{jNYyQ$2dfcxk#%;OE_Qj_^iVNNSDuY$Rq#Qa`jj%PE^vCOq6^L>RmXU0o2 znfq`eTIcPkP&h{sv)WU>SX0K(D9L?^*Obo4)U$_g3_OJbS>9eUQLjxWRsi zVoyxO1cRUKjlS%U_Uw^m?2}vUm6hz56`1<2jD0hly>pxWGnPFxjD1v%z4VIx^rxqE z*jHNYtu*%6o4Z$--dnkCD?~fg8k~Fun%7Z z``BHuPcDP~`SGxSwhu!4Fht8f5d0e-v351Y#R5RP3E+c z&OxUStI*ju1)WO`;TW_8P8DCkX{`ae)W}7b?Rx0id@Z^jZbyz!(e2`JbnhRG?)Plq zJf9{)(Rb9d6Gbv(PJG273KCiQe<9(YtsX`gngtpMoFgJ9q>7W+kCt zyDI1xw;e7r%3zz#g=-;k0`bR%U|BvM`pxZnQh-r-hFREi;D}4-{^#cP_uVUaA zUktLojX|EzFeq{*2Aw^NL9b@Ol|TD)ZB`GiT`IwKa5uP4%7yE!IJgF!gX;p`{^#?b zuP0sS_q+0b6?y-o$^IR2qzMVX07|+Wyd%xuQqsgH@c^n{@3FI@GoF0pb#hW;@&>Z))M0VmtG(32oPWXI>sNH*=Vu7tE0>^Yogz%4WVA zFlV)xw_D8J0_LwHbJ&4-WOkozXFlDT(~rz+G;{0A{KhlKLz(9j%r!H8GKV?ez`So^ z?uRh{Yv@5!`mlmt9Ht-7=*c+xl1Oiw(w{_nWJ;fQqaANPlU2BdayUnvOg}fN2ah(3fU{0*)IdxGZ)!6DeRp{_K!AuXchZtDSK%p`^lR< z70tdn#NOJ<{(8e6GyZ2b_<#P?@W}d(L0_(*KYu>j_gYuSY zxQ33OG|=I;4cdQig?24$Ag5jh#`!_;XP<3OEkg_beyfG!3p5Y5L$mW*Xl7S{rW*yC{BS}Ox8G>IuQnR%)JCJ(o6+#` zYczEJj0XEW(V)Ry)L%Us^~?W)<$}(rSF#rMW^aYXcQ05h?hA`5^HFzO3F@}3fjYOY zqt3ht)UiB<+SiYw_R>eF?U;<(6@J1zCm80NhQWM#IL!O#!rXcj_;1s~T=N*rzUsp4 zLphken+LN`17P-JEX?`7!`yru%cITn9n`j1g4%B7Q9INM zwX>e1j=m1+^qYY?Q5vZ8C>nKJ7^1GfGwSBlgGIf$u$bEd7T3*DuSFx&TiyiqUTMJ6 zIUJVpM^L|VAJm^|kNSCi(V){TG}yBY4UNX5q2D_+EclH^{e02rxHTF#G(h9fw`g3N zj3(o((d24hG({nrZYw}D-CAfi-3!g`oIvyTr0j!rSfc2eRv}}@wmUHaTGTj)h^c~S^Kr&iIbwsP%k6}|a05*NgV6)U6Hpgec zrl=!qYox%|@gQs`Ily+6BWx26!}h`&*giT3+v1*Rt!07M#tvw0wi&JKbwca<{QA%5 zKVKWS=lARIem!{qvwYktKAzWF#p_OL%j@#mH4?Z^JFYjC>ki=hrQF9H?&m)D)t39) z!+o~re)YL;YwrIz&ymIR%;mWr@O&mbr{H@~w%LDQ`Ku(_IWkhcG z$nP;Zb|KFOUN#y^Tb9lyi4CP$rb3T(eCnwIUIp^lc`6Y0UA*lEA zE9Y99^R2}>SK+)ZICn43zZ!GUo_Ww@E>@z>x)SDOJnHEEU~cl5pL(e6`h|I#!(3%C zU&hSYH0G@m%w5WuKi+P>j(J?bTzW8{-I&ve`+=O{v3-jK*|MU>% z+=Y4fW9~mN|2ycxP5O{TFOuoUa(Z%}z8KJ(BKk9v9y!sc(ex^cel4YEg1!~gJ9qkb zg&q#3kA?K|G5uUePpi?_YV`IC{cT5&ThixzG|T)&zgy7rNAx|5-Y=s6quB#7>;ul# zIF9{rhCSiTz8K2hXuW5Rt1=U@6PVrj#@=qQOjNz zHSJPRV{S02r#^vcv4pAhEmWJBhpLw=qN+v*m<%(3G58I${kizAn-VoC&ST_Rtzq z0482b9bFtZJ+RFkMOb_%K-IEX4`B{1r32BYXoFnSyZ<0cPb>=Oy&Lx*7e ztr1M*Ett&P3zO8fF!?+aRjoRr>e#2Ky1^M$b4H=+?;5Da--TE6utK$klTmGt0jg!| zpxWDhFx9&Q)5hyy+G!z7-6CK*=^RY`bYL1X8m3|AVY=QPrm;yd-O}veKCf@{Z?}hx z;{7V}{%85PFg|_(ufxad-Q;z9^ZHs`N0;kO;JQt@{s``)D)%##`*Py`9&n#+x!G!Sh!ohxg=>NG@Z^$A+A&$ZICKnUh}^a?~cz z`Q(}bgQefdng8O)lW=mkA^+x_LoDa9gL4_e`Q&p>&YagV&aF4+XU;h~bDo8qYhBJ) zi*rupya#gb#hiaCb8wt_sL5P>WIjeRCqc~1bmqo|`Dw=-xiC+O%#}Cub%;5OW!`!+ zcXye;Nz7qw=CLYsIgR=3&z#<4Uagp0elL1IbKHP=E?};8nD0Z(`6uSRh`Eno{_D_# z5%ggjy$GfsS@fg^eYr$$jOou5dbFHAy{A_h^ve(W+kVnFcY3#k{&k>-FX>}Gy?jYO z*U{6C^tCd*^`*bs^jOg6^7Q%!{k}@iZRz`cdcT1F?`02EU>|H`FT}DR98sr28T$h4 z4F~o|PxiRr9Jz_g*~&DeKU@|)0X`c#2(tqK5EKdn#X>!XHQLMUupf-TkNlY zZhQaVog7w)X=vW0D;j@mi3WvNV0mN#>OP-}+7DjA?BzJrDnAi5ng_siTnki9t_Tw? zy7X)(stmA%p==0))6JltG6s6vv!T0Q8e&%Ze+ch?E!OnQsoC86Ti`qAQNhb!WT<|Og`{66t*Woz+upphuC&=zs#GCfxMbXSqQF!R0C@gO!Ub{aK zua3W{ZNQ<3z!*RidDs zmMEwfE(#34ivq2mqCj(uD5&^R6qrmE1r2wKf({;{V7QwoSTtS~Bn=b=4_=GsCZELf ze(l8b(6!?E&C}vVopABOi|d{76fdj&6fdXE6fdttiC3-c#H)4Y;?=L(;`Kyr@%q*- zQP?$G6drXJMJ) zT^GMDJ{P~TI{(|}v}?Rwo8Mo-`+4&IF?^g2AJ6}55rfmbUKX!AgV#^vIyqc#FV~&S z^&4{^-?*Qz+*b$g?=|;Xo%?;xeb40npYj}oc%CAjD~IQ+&2wt-yxV#120Z_Ma+pCL zYsp2BPbfK^C$Ea+7D#>`)zybNGsRTqFwG7Bim#%;_NJ^`t2H{DJwkVvggP z=Q!qiH1l1KIls)j=P>spng4P0AcsEGrWZ%)M-V-6r!SZ3O<(#mj2@-Xrz-R+kA5ws zXD#X5CQ(xNH~rJ3hYjiD4thC^esbpP2hrEX^tK!Qb)d&}>2n8qy%Cxxe$(?1`ff|_ z@6!Kc?1A3w15frs6#GGgJ#mqJF_OK}nEheN9;w4V*~VVk%zn{g&opJ<__24!vwu9< zLoe7zUhJhn_EQdf%9?%E^1=m45y(I3u8U=M8kr zU4f1q1MP*aP~lEGG}|XZ!?Y*Lh20Rpiw20FmE*+sH#@}FYmws1iD2LtOgM zMO=27EiUh$BC?FXi!9$#k#(oL$d;)hds73E{dEHlNGMQ_EGr>(@* zPTnfn<;szRCwVj?m#5@&f}C!W*FbU`L4GyKF@ZeWk?VW%eM`0?1q%l9%%uzD)G@ZGMW4;3Uxe4=jhqs?({`5q;?^ovW2y@Bred2YRl`*eP znOiWwp3L!c=J_3SeVq9YX3nQF?@`RXGxI-_9t6+_b9&K)ehi@}jp@sAdh?0?JQi24 z|DaDj=v5B=Do4*c(zjRi?jim2q=(Px<7;|ZoqnF9r*ZT(k>2Xl-vjhGgg$4}>l^g@ z1wEfb-;dM#t@Pi8J@A-)P{>{w%zg-DPpoEN#IiTGu|ImVN3OF^y0cf}*)M7A8Ef{< zDE3Y^`=^jSWY0eGV=t{`KN+&8vi|gy%JzSa04nw3&pcMWM+M(TDF0@wDEpZ&es=yO zzGird;z8@gM~#`{ZBBJjcyN<=HIA=W=B9|JyOPD@1>3~Kk-6gj#TgTN0aK(Lb zAi+WG&ovVJH5-e44tXMFVRez5V=0m?XNjbsl_K%wB$4ReKqOq9A`)C|MSLDa{N&*x zUUQt-yR(tl+k2we^W%ru!~c+D&$P~Bk6oyUD;p)^t}PL9ahFA0u!e}6@JqzG@b(VN zM4ZiO5!Xyh#5J}NaZPKAxRwzjuI+3Q*JHhi8+J>?`SSj2zKggd-hOqjh@&TaY9Icu zO-DWY|0e>R_n*1zasIK)K?w8kR>aLLVLmP~C##v4$;=IJZ)(6CHRNaJ#pVF>6~LS= zX5K!ExFx)Q9J6+9F7r5uxm?VAwqZ^;Gq3NN+hNRaE#`PI^W2}gwqm{;Gv}k2_ZrOo zaOQsjJ(x@%s?!U7`f-Jxn9~<+dSgp}uF<1S^rs<;7cN@sH;faT1zCjRTq1|PZxW#a_NSN*zFc3cG>(EJB_x8?InA~wx`Bo z>%9Q6Acu5eY04<=eme_Ra3;YY9*qB28*b>@gho& z66>~C6YGp`i^!G3M5Nvd5fRZ^nbpjuv6A6GWJPyjXkb2md6; zO00D=7Hb;>{wvQ*?d|VPfcf)f4t1Ev`po4==5s4^I-GfRXKqh1zx+75lzDbxuG=x+ zLCpC8<~^Ia*JJ+O>A^Po(2rg;rXPm%Oi{`tTjYQaA0&t5A0 z(@%fziT}@@PCR&>FY=n0ikllpifcd2M0VC@acP#FINv@=q+9$JXX{ytQ|%jwX}O7aKjejVCCmxUo2lkY)uh~Ni{|MjGJQ5`XOTXr3Ye{(_XP-iM7~v^`Y46mL@ho zb`Y`N>%_(ftHp-scVfMLw20346;WFTiFM<;h{%=;MR-|d5q2|8tUWPCgvLG>YeEl+ z)vMQvRV&+wm65~5iml7V^5bq|+1>Ucq%2y5Sho~Q$GVE(9lB!4YbllpEwOmj9kHm; zQY`BGQY<_HvCukS1nph&KTA`srU?Ac6JX8*nD?8^eJkew5k06yAFhZc_P^+d4L#AI zFMa4uWBL7EUwXBOe!0>!Yx=f~-X+n$hxBkLeLO)gE7DItdfI`$_NBMp^f!zi zzoO3}^g5J&gPxzH?|JmT3H{&69>`-KxUm3^Q+-+M)oxyV7JS(g*1_PC2<$qmJkY462> zi=#zKj{gKY=r;!y<*;%`C`tm)?&8NSK-?(Ma-CzC8nj^6h2ygg!kCJfA49i(hkHl&Yt=%t{ao9JnG`szY&N6_Dz^!Pn} zoZ<`~`jYruUxoza@L%Ec@U!d%=tS@Pj>3gMG1^y)l;kQIkE=n|<<{y^_s- zdB~nIV&9Bp@4RFGH2ABB{`3)h>A#P}Kl{oeShUYkEOfje z7I^Lx0qzIHylw5pZ2rv0x9&?ZJ@lLK(f=yEq92OMF89QQ5(hCZJzaQi+$Tma7$ipe zoD{?T#{aV?p-!&|{J$eW|Ax}T1N5;cz1&SdkJHoE^z|COy+wcD(c|~@`4PQdMZfLn zc|3jJL+>-_e+TwJ7W?2hd*K`V!Jj=*i+xd*y|ID)agROnhJBLEUMa_Z$!5<~WZ$%B z?`-<3e-d6f{^_ItcNG2^D~nBA+llqMZX!zOfCxVwEJA;Hh}8o;#EK~`M98-VB6!mg zvDoc{SUAO11V&8}e)I2&**_kOnT>x4pSwSV*XGM&(!3pF{G@o{IrWtowcK3{|4>H^ zZTVJs>^~sfI{J!%_u7j7i|+n2h^kX60)HO@^!7UaZB388>2m_TenG#h({mU4{*>PD zp#QP#feiM+7xqF4`{6!&;wJlI0()a1`{NyZA^acAEw`Kx#S^bdRJ z??dfh`y^szL|d^eVWbFNQ9&$vG)4rKYa{$?l@qhaM2qQndkF7&-NYoPTVkBOix}On zff)WaKn&?&Alz%U5QEB%5-vv;i9Y;WnV#Q<31=_$9R067ZnZl_;J-|Oo?oEvZuGu3 z{h!Ys7|K3q%wDL(ewfRisK>tO!`{ece>ky6y0K5z3ICeE*)MkNnbGVUP4>=+KmGGx z-jKS$e~N&ZxA>iyxw?_?zBf}${J2$k9(55T?XAU-NAHB&!P899HR19uBOAPX=Ec(rCAbOeGiEgX=2*)s6(cxtGfAQ_B z?I{9^z`qy)_C{CsM=SP7FZM}m_DX&Bi@{$#!@l_!+wMPdJEC*Cv2d6q{*i}QhgJj> zf&VrF>=k?Vi_w4E0_wtwfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2Or zBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs> zC<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%I zfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`e zihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhl zpa>`eihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2OrBA^H; z0*Zhlpa>`eihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2Or zBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs> zC<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%I zfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`eihv@Z2q*%IfFhs>C<2OrBA^H;0*Zhlpa>`e Yihv@Z2q*%IfFhs>C<2PWKY+mh0*&kniU0rr literal 0 HcmV?d00001 diff --git a/data/modelrxte_ait_3to20keV_flux_2deg.fits b/data/modelrxte_ait_3to20keV_flux_2deg.fits new file mode 100644 index 0000000000000000000000000000000000000000..196810257792f6f3a493a085261b007a0d1e4229 GIT binary patch literal 80640 zcmeFa2{_jK+V`D8lp&Q+3e5vW#Q!{h|Nn*1N|F-IwTd*O&>#(HMzgV5lS(v5^CX(* zG?z50G-%W)@3(c|$A0d8KWpv1-}iXd-ur&utK(SPa=EVa{GOlBSF5$6-P6;%dsh!j zOU3ece)w6ovm7>ZOpxW!(4cV3z|fFkp<%(65uuiyJpFoFhDQtu89F3v=wH3Cqo-eY z&py99&&Iwz=7o$84hRaf3>{_}Fftp7^`p z&G*;yej77n5|aLnAHU7B``vkVWQ%-bf*J22F|x3JD(>8uC~7yZLx~ zdvx=|y5C;dveTHLNh1TG)H8U<@E}XSupuGgV?)CtEPF9vU=k zWJu6ZOU!QJ8XgfA8WI{D`M>5tzU6NZj~x*iJ~A-8g|+2`&@tnKgDmZx9IY+K{v2=L zfdo7U1cnX`V$1fnwsx+~UE8*8i)Y|}I#6dzH^08!J%01i;+L*I-E3Sv{aXANHtj6? z*bJG3FcFbsgZ|6ofB3h*&HK#<=DGU)@ALZpzntgW%eCuY`eNJG_Afs?+Wh9%H)vQ; zSWrk{(0?{~!jLiJ|3*@`pY#6Gm;V#<@aguKeAKqXf6ZryzWn!me@N~>KH%nyd4I{* z;ji!g|1!_rqpROv^0oi#c|Sit+=GUP1qFq-`-9K_{~p)v=e)n<`#&(x*T>KGFZupA zAOC&*-+kZR)ek?L{`~xH(cae1(Z*J?v9t5Dbz)oM4)!f@@NYk%28IO;0H2Hjk(QReZN0)qPQa(* zUwx+k-Dtl(!L750n^#YNZ{+J3>fY6TjBD49|M}O?&%1V(o!$K{14je}jt(CmY&m}H z&wNAinft$;@}KK)uPes=yw}dr!O70aUUG2yU!TG4{^MuS@gXB4{;{dQ{p9zLouYo0 znt+;sn!vx5Ko#{5H32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mN zH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mN zH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mNH32mN zH32mNH32mNH32mNH32mNH32mNH32n&|EmZ{ll@fxR~@82=1(Pnz9@g{f`7trX{E+w zDY?yPY2lM&QryUU(yU%pr05&VrO03O{|W2WEB?(1pa(kWL#i~({|EX}gr4M~FQb0; z=HEQ}AI6i87d@5^<-C%1JD5n@_O6pQ?|LAux7aVOrXXqAnif*>%umviGZ&?W#+Rh{ zOUI?T$<6+7nEC_%;|aj8B|IO6?{n~e7XCk=2Tjq30Q6!n`f(3ENk?Cd(VLe)`}05T zXz5GQTB&SEsPwdLq;&7CEZul9Rl58jN;-SAw{+s-IVtbhZs|~$2Kl`>m}N}Fo$kT&pJY3($7Y4!6mY30OnX?c4MDa~l-|1=K&f65C#o2BJ;U*PLC zygh}#Ht_fyKDWT@2>9&^&vEd59o`+_{|I`p1ARD!Ui3siCZZ>k(HBN3aMwsc`;$>C*Al(gpMD(%E$;(rKT%(#e(?(s3IDDSx1o zl(*-SbX4vx9e&HwyZ*7(LpFKF#=dV*nS95Byt2H(45l2>L=cv)JMbzNrS~82-KqhX3$=G=X)d@IAMl*cKaLXq$s2V|f z1GA~dg=wVQa|%^|2_DO5lUBfa(sWLuYK|SKY9}U*P(f98_)?Xcm!u!x4oN@szDgCI zEPc;*kiOGq>D%)_>DvL7^ewBo^zD!d-nW##x9cN)pZ`Rvs9G*nB;-p!TFsGu+?Xp> zCK*r_TYIXqyo@yTN=PHtjWoX5Q`HG2RP{qVRSWM;)!t1b&2dbcuTPQIsNtmbXyAXn z{h$5!*xw5|1|v@%a@ix_MC9~A-WkY!5Bb~R4uN}$a90WLYb{kKe#gBNaQArJza1PN zgGV8_cz};SIB9`b9dK)d_j`b2GI%xy*OuV>44gB;`y{w~f`5B>u!9d9c$p7B$Kh!j ze4U23@$feh9y8!G7+&ANZ)bSUh3{l|?+^cL(1VreLp*wM3jLUkp0q$;%+Q-@=+98} z$ndv5p;!NGzW?YxnyOQW!Y?_|n91KMXm4K{9Qu;_Pn$=*iZ@bEcYE@=B$Ib-JL)pf z3;*|sF642?o!nB_lIyDJhPsZ*q#uAV=N983ldsj&9Y8Fu2jwWQ2yo79Yj*xZAEwb)-gRFI;sm;?_)aG71wRy>8ZK#lSCr`3o zWJcC6#*s~Lyyo&PYTIuvwJrNh?H2tYTREI;-?k&W?3&bm+(WV#3&=tD1v%U)qz>Dg zlH;s{n4C)|lS@!D za=B7QuJ#_}IxB-*Pj4jGAJ3^{v*y%MwxN!0MR;qC_nYDIr}+4KJZ~MI?~2#G!0RnB zXFTS`W3DmgU&o%6*!K(e9>D%r$Z-RCrXyDw@;M=ADDrMcZd2qxi93pMPb}`bPIg(9 zxU(_t#m`5J-s1ij;1C5K=fI@~_zVRnJ@Asj%@+J@z;POQP6F3(@I3?0=fOJ%+-HM- zO?cP~A1C3(6n-M$DH^^S!dnjft$@dl@Hr4($HDJ=c&-WG6X1Or{43}Ip${JDg*E!& zgr2lTU)G{Gv(TUC=uvO<>FaO3LcjjeF8(v0N6V)7rDS18OS*2Sg}q#_on;Fd*wTWVzPotQ2rzpCu2Tj{=KvTO{rKkf{X|kaUP1=2tCaS*E`0J4r z5$Q$Y9ZYE4t3fn&Up)$4^NB*z%PII+1_f73qcI~^(&)E(G(X_Kbqh0)sS zR{mnl1UC~doHLfgyR)6V({l;ixB zb_MmI+pY`8v%gKlT;nzc@n$9bKtl)=Mhb9YqC&yQrXOJl^Ku{YX4sijS-Cyc2l- zbG*)qjz#^zoTix91#=1WPhw9I_Kn5f801n<`He;hri zk3JkiF9OgHPxPcS`qBfv$wq&)(W5u$Q+M>rAN`v9vuFQEAOBg;q4&0y^y-AKBt;Y4?GTGO4b<#cCSINhy1k?uTROLr&ZuhaLb$aW$X z^)I2KMHN(ZEQgAoAESHK3h7?$5p=Iv8r^HLlkV|My64@O?uGBAds{Zqz0V0${Offp ze&|H^qvPoTzoUn>zSAR(e0r>BPfywzQpw0_^mOk|dTzdlN|$b?m*(~9)vi=3>%WEG z8Wz&Ki}UG2T30F`9Yvqrr_sFR>NUxYiJ;hz4uD;@WZ!JQ*;FX8UVxc>+^B!kCH zaM1>zByd^}UaP5O#2b3j_8T}_g6DQ{eGb0G;2a6w&A`0?{G;F@1wJz1#S?zQ;b|>= z#lc%Y_&W}dyWn#pyuN{7bisBEd|!k2RQOLt4?@s~QRu~M^rH!SVu!v2qBqyjpK|nQ zF8cHiy;_Za4M)!`{#l*=dwaP4`&(S^#Y?X1X~aekOu1IN6&p@Izy`s&Y@pO&gBk&B z5VegBu3l%u24C6me1EQWqB9#U=)kq7UgbKn4cFE8;JWwI*mzbauII6a>)jf_CIM5~ zjJkvI+(57A7-1am$;qHM7GUI;Py%} z+n>e%S#`hu?D+f#OOviMYb{{0-~uakh*gV2*hNR5UB~a_jyEgW-8O|iqO7>np6>k1 zqr2R>GJ(73w&yODx$N~&hr4cbV((#r+)eu(cU$y`eX31mp8+N8v-}a$Ou&FoWU1b17w4A0+>*A2ky?J!3J^Rh6v4d&0lomf^b@yq|~vw&+1e^uZs!_<(*~Mo%))mqF-_Hu|#@J=%*tS)f-Y=$9dS_Tp#X z{$3~l?2ofs3bD(uLYB3?nI`116zjuIo0_uYvm|zOyUUKJ^;vSB%F>oM%=@>nxY&V} zs8V)bn#Qgpma*HE9qf^f-YiI9&xs!F)vGFdH>$04M|pI; z=^QLhb7*)s9(N#%BkD}z36Tyw=|wM&8qtxbzV5-%%43e16Uwny@J|&B!Lvrk^X%Pc zcuwUTp6BtB?d_f&f_*#<}v>eY1yLI42qilI`R5xA{J&l)6+Qo?hWt^mvI9Yci zCl?Inl;Gu@s`-dh=cjR+{sc~&ZpmqdvAnEW4lk?Ayv(o_r+wAov}+4FZFM%M^-AH? zFUvS}?s>d!C0?J5IsGs%9dp}b{xj^!!oCUE+X?$yB8N5dj6trd$X5e7(~Sw=>pl4pc^$orIv%CGb=5xfA zCmfnJn#WYe@~D-Mc~swrJj!l4k81XjM;BOdu-_C8xjT% zGLBhsl4r!n@~q@1JZH2i&+{70@vTyMzEyu-aFcn_qz$}8<1{awX~~H%Lpa&9JEtsj z;nbqXyv(8rFCTK3S7f&0mEYU(>c01R&GB2j*0nLOJNJ|~jC9}(<5#@tN&;_LHH^1T z%I56B(Y)R7G4J%d%Da9g-W^uPdlwwz{X4?=;H@ouM7t*Eb!g1_L7{v+X)vG6@585G zZQ^sa)A)jIAYbwzzT)G?SG#8Lbup7~npNG&>Oy3U46+g{@5<3{qc10npZWgmXJxg$SqGLTDV zKISJU4fyf1(|G+J%qhma2+Tc&`Q5Q+1@`s8-qzUv2011oPa1M{K)%+<>4&`YkUJIm zJL8V-xMvXVGU6jzUvcL$+&ct!@5KH0z@a60q7YL4bQdU`yRYU z!G9il@CbeIMK8vqA1%?7U(uIe(VGeAj~ROO27O9KukN8=574t{^ljA7-uN$H&Uu@xsa9mI*-O*uK_IHy|v;IuR2d3nf0Uioo7r+2XB^)3fF!?K*W zL~rJ7PUan3hH!3=hP*Ukk&K{=&$|LDaEx5yo|#M7{S% zh3VxmQGd%#VLmNHH0rWWSn3xEs}qr;N!SI^%&47czWR-5*>I$2HRp?HeNSIlH)$V{ZxeH${$P$P1%Hj;F&jQ5c#VPINATcm6+-8AZD$o zFJ@nI6La3Q5Ocrw6Z2jriMX?C#eC8biGx0g)rU8U%z6vO&WZtI*Ebijqh4pR>Be)h z+A&8YUpE(vRTIU6?khyXw@V^{---D686xgf9_B2>yq=hAiTN`2)W*Ix*lUUXNyuS~ zJhsSn2l*Bw=V;^|hTH+jZ;m_K^jZg!>)9Arm}`!Nm-GqQR*y zcx8dxW$?QQjCcW}=De+_sT2Os_5B>{c}JPm}eEO^^440X%l z@eO=t!E0~$y#&uM;oAq^ZQ(x(b4PVd;KlCFRJ;_90649I1=udO>s1y36hhA+# zzx2_wp6J^m^e*>j|Nf|(f9B)-e&iW`qiAv22MsPOS;lW_f97`=FZ0Ls3;65NO0K+e zT2u|)CA3QOgl^A)LVw<1VYF(hsMk+h)PHtBSbXU%ngu6{)+<*F+q!i{hrLaN2u~F* zw>yhY%i={>HxuEjSs;2Z94-bt4iy0{_KV?N&BPcl4>8WTRE)p;LQD?(B&Jt?C1Oq` zikW)5#cY2gF?XY}h}Re==8qaH7Cx*gmW+%MNgC}$%4U19YzT{$El-Knl}*IjM?qr4 z)fr;bWiOF=Ge~4V>m_znMu}W=U9nfz6$b|P6^Ey0io9jhM8S@x;^g6j;>=M;aeikv zacOaqDC{>@T-R+RZl)E8J313ZQIDRYcx?|HJSwXWPcW|C~2_DR;VsOBk%&Z2j9SNFgON)r#84&1>X#Co(WJp??xIaav9Q(bD;&263dJd2xTW_HUbUYK-@c7Rzo@k$ z(BDmrellN-*UlBwQdWtX@=7u9mmINR@CmVWzq3d^D~VNAdWvpAAqO_nulofOq zA0BNHUziVzIn(kn(dNXcL|qmf*#AZ5mjaTvHfJnUV1WFj;#9Sfb5bMA$Qz8Pj*j} zj%z};Jp{zBf&ou9(3Tt6kcAzj|Dsx!!|Wujjj(X(IC zHy8A->d*fDQ8)k0$Aycdu&3y<0zMCaE>Mc0xKqG!EfqHp&eVn|Z1 z7&mOljsIX4Gvh<~~~{7EaC;N$#h(5x0{L! zgBpq(iyw#wzG31;gGBMsFbLbatM^0oox^3F&TIk#7PdC!xU^8R`qC}WytD-IT!igxQ;@{I=7<$SjS`S8={^5N<( znDYeliZHho<~P8eCfKKqy^pa!5II7S#|60-BHu3Le1yC{$lVwD*WwOG+_MyS72v*b z+<6K2y5sJ-xIYXW?7+hZTnxa+8k{1*D-7IJ;MZ1mGJ6Z2L&0^G+|uQT+{E)0coVo^ z0)Jx<~{qPhHUsK_2ApBi`M?Ls7hSvq~+X$Y!!*>?Ew}*dM^dJj;Sc+aW zML$}hCo=ld5WUGlf4ZPYgVCol^y(n`)gC>2fWDnZ?;8E=-=Epxr3JxniF0d9k2+iAZ$%D3*S*D|Q4$h<#1giTvRS z;{3)+am%2#cpUmcynWwNR1W+hYu_Fs8-z8K>)o-GEm{Z4Eo~OacBik&Vo8qdG4qbx zy=A4`@AeXTn4X6`ZuA{_O4$y1X6KV~oKqipQQLGmS#B(^=u%ByJMNggY2#)&`;)hv zJEB59P;pE?cBqwnI_s`{>DoT|CVu`azLOz8&G(VZ-aVE-kLj;eNz+uc8XGIRmsv5` za6qZEGEXtxeOGDt%0+1+E-NjU_fy(dovYXneXmG6qZQ?PFQsGdQ{|U-S&CQv!Adue z#!C0eO_d(AA1FP#>Eo?~GIY*7W$s}sC38w+<>rM%mG;X3RpXn+%C^o| zmEkYfm8+5OmDHill+kl)E542Lly3byD?ZQ5F{c*h+{4`Zm|q8bys<9@doN(Wj2uD8 za~-+HAfJMqxyaiTx$}^JBkoAYJywdARfSR|^%d@Hi+lIt?&i3EIyejmk0fx}2R=IB zv=+RI!Oaf*7J%bsd4=Z-aMc6f1aS5M?_6+S1pWu%p*4I=f|qCTa{`{$z*hmhO@lu_ zc&r1Twc)il{Jwx^6Zn1t?^^Kx4n2rLAJWi^(dfqxv1VTd`f?Awd4v8uK#xwNPdm{o zUG%FLdiDiVqJP`|sGEP*<03VwwpivLAeJ?D5G&ek5UT>L#oBY@L`I)@u~p}; z*!k36?0Z%z@@usa=lTy2*Dhs<2Sa;{GQBn8d+``qD?C`%kJgv#EjE=c)6?YEm(Ix@ zDu0zbenYp-V&p;A#d7e+A@Z~h*7Dp&gXASO@5(C-PstfgcgZ_n8Ouk^t>ts7LGo?4 zZ28&d`|_t}50q*TB}LEilv4NjET!S0NTr3|V#Rh-iNcxYiiiFj#ph0((%*80GVED5 zCERhb60KLF%)VHr%zt61Bn@4wtW*t9Hq6|wWIN7K_Vhij9DUbFIrD9fa&1DTa=+hO z<;C)3>v}TdEp|JX29}xN6|k1FE?dZB*Mr9;qJYwR6_6@1VN3)LON^BvTb)psljIYN4ug zc8aRz(UmHLuWl-1_gGc^qgt4kg1OsN4ZnTFo+{Xvj=j#^)a_&do zGsx|)9DVy5cckK;KDcW!?sLSQ)0O$9uW|Q9+?Oh_Z+hDP*_@r2w+C!|}Ttj3$4i{P8gG8=L6LF|s zUvYB!Lvf|Gn?w?p9 z2e)`AM<07H$9w%Mr}*m2>r3X!ISU)h`7aCPYY86m^Yo8$Wqv0`uV$cPHgvqwva7Qq zrF~Pngd9)?mK7-B3HOv)DMOScGnXoB=B!e-Jv^oyP4ZH%6wFXcQlBW_Rgy~g{s2|I z&9hZalQLBHc~L5NFLzZhS5wu<{1K{Y@9wMOKJ8bfM6haoX%AIS;6T-}7hhF{(@&_L zXw6W4>hJHY)nuHr!G2$7)0StQt*RS2x4gR4x$V`C&h70>oZD*^I@@GDb#9i{-??GW zN@vrpJ)F(HcXGCxGSAue(p%^54K6q@wVLQ+vuBR;lCrnXUO{c0TXeeUY&Yp}lsa&$r-4diN%e3Owg2YJnr z+Zg$);f}euXEyFKRrPYYggdL@URT^5f&2S_gM;#2c>^vR!N&`n>VsD*xJ80rI5>KM zXEwMl1K(tDz6;(>!F?e3`@lnE__zTtt>I@YJQ={(L^-8L1^oHK<1_fY0I%2JHxr&q z;M)q`C&GUgdN2%q@Io(q(2pqeWG?z5(3`vHPfzry5&F~*y$VIYs-tH<=-UbO?$_V? zhaUc3KmV-9#F&RoMCh^xBDBjY5jNRZjDM6SrX2PZF(rq^thn(auDGXIk(NqB~{kIjD{rPB@~)It^3gGUHAipHi|if+STrPhP1is_qT#j0n3(z>?0V(&Obk@L-zPWhdc9YL=|F ze(-Qrj(hM4Ue3YKDR?S^uSxLs4E_inb>XuKyuOFuY4E%dzO&)|0sL=754MQ{&(G*Z zee{FSlcVTMSM(+Z{n>#Y`Jhi>=+y!A%N;#SM&EqUyA|kP{hvMjvpW7qbA^om?uTE* z3FTwDaQ>7g++)rQPiu4G{Y6{&w|y!GOkF94nXVF{$6kxb8EztW_aZTO=XtSkS1*yW zwTf81#!hUCydrk^3>OE^N#cZKeQ~*WT~TD|Bud*_i%)eeebasb*(K*P5ZivVfPVA9q%N?q+Tth{=y$h!{v(< z%R0vR?I=^F@!Ew-<0%W2#;;tJ4ko*lU+2imBB`~q-{^``>ikfp|Efe~H8)<>LhHQB zta_aC=~0xj>)BEzYDcjm4jCv-O&k?V-&2ZJ$WzQ&gn0>=YlQi_*b{?&FBKEx3ha+Y z4qv5a`ya@ahZmUgE9*`DD@)+_@e1*1_GDa@L+wa3}HO1&(gu=>V?A;H${J;=Y2n1Go%|Gw<0SXT5BPrmFB zuZkkY$IUsSA|P0tN4o?3;nTk=XBx94h2#id+%M7l51_k+&Oin#8%nS!$sui)$r-ihEo5d4$jL4gkqaiP^~_&EVji{Z-(-de%mC3wt%&mQoa1HY#5Tms+m z@a_Zu#$sSpIr=aLy?B6rWS}Qa(U(Z{Wky*n> zWM`#{odGArZtL-4zkV}uSUXl6Gcyt=JZgwDaoxm)XSSkn;6-uc`9pCx?v{9{iM5Bwqudo)8QocK-h7MPNbi%};6)v|eo2yScBumMWX#pU{QcNdfPK}m z_vikD$Ps`%(a3cV`JN-^VdU+D++&g7OO)ll!#!1T*H+wj40j&Fy(e+^F5F)a92~?6 z_qXDh=@;-B1x^OwH3;0cfnO{*lE{eo0IoN|*BG47g0}?jM&KU{50Bxa6TIAlpSNP_ z(3kLa7v7q~pMXaf_$-50PxuXjXEXSo4)0#@-vT`_6AhRDKrcF=A8zQ02l{dhy)i<6 zE}%z==+l1msv-LI3_a6B-`1daEr0g!kGlD1KF$YU*zo=c7v3+zc)#sa-rv)n4~Bi^ z!yCHru@xzNB7Oy*nXrm47De*47CZRP_8t6CdnG>~bb;TD3*yh-OGTA{HbQg6DxtGs zo6y&`5Js)_gz@HyqQ1o*(P&13X#Dk-Xwi6-XfrB6w41(2*yC#noV@desGA^EM|X*i zz26Ft2i=9|@G#N!Ww_`*ahLF`Wh;8`jT8OHmx*6JpNT-n(PFr)FGlyTB0`tsi|~hO zVv^^tV#@JdB6`br5&OM`nC0&)=H&JkarNJbgqeHAg16RU@t}2L=~-Wq+@Y&TT{c-P z`@Bx9^bQiMlNFJEX}VbVB|vP{uN0eV_7|I92aC*I9Yj{&^CJ6hoY>B`Vn=eJ*gqP- zm)5XQd~Uo|HgwOG8;ullvzQff6TNM+Nl{0+y76xDLKDAV@hMas)+iIZ$LNZj!&+kJ z>@Q-+I48_$iFvV@n~(WzuqOcfe6aUD_WL7881l43uHML(h@7L5S3zz!)SVz)v(h^@K0{zU9}m@HZVEGvRX#ylTL2J$TN5?__vy3jedv zgFfg(K6^zjsD zE!@L9UbNxehkEjX&YFBQy9yt-v*9!THhghX8@}4%3*WMB$wfc#b&V@zeyTH%Uj((| zvbDwhKI0pIT2+I;4d2R@4K;*DP=8TPS}wHo8wu^FhC=5|OHpHQ3!%63wJyQUVo@)6kAZVd?kRjr z-wVGv<)RmUf4?_A2l}^}Er#Vfib(gxV%Fmnk?>1zG5^?35jQ$m%t-AnCeA7rqZ$;7 zA>I*UK*3(o@4#ZwcjCmxN9%&!>=8t<7)&ybj95!xZg!IiYW&VKXB;WxOd}9^XKr8$e&iehnL6ja|WIoz}EzLYXyH} z;L#L5OX0O4{93@X6(7y|$OnG;0{>mmgU0AXHhPhVeng=sH_?}==*?U7Cm%i9i9Wfb zSMAZS%jj7o`lf^4)&JSQKkDY6{Wy<_Kg}b2_wcY)i9Ec05(mBR!oyBH=8;i%c}&9? z4%rdU8E`^QQ(m|tpO?7ja*~Z6r<^^) z%lggZ6(3@G^@s+Xz9)&-J+@IN0tn2$chp%*vMkAdh(YxJcK zdNULKi9(NVpijlx7qO#JXYO#i3jI48Ff$_BM8L^X87HUD@M;1A8{N z?{XFQJJo{+Mm^xccCR>K>0%yQx}Jxd`SM7|PCVKvpM!0EIkf3S9%ufQ!_D6E_~xy6 zqT6+joHUiA3eNDf)@dB=ufZ|t_&$=awmfqXzW?Gv4A1GD$8*me<+#2(IsWxgo*!SG z7uuxpqH9xlNz7(Wv@_tOhYL6*v4YcnY0Aq!<@1US3B1auHLrdzbNZrzyslXzUccRl zH?}*<8Eg0R=IRT1%YZ`8OmD_p3)^zm+l`$4tuJpYYtGxxY4VQL#=LWmE$1AL;9cb_ zIrnG`=RSJExeNMnZZ$97HLVWsGTO{JJ1RKG%a(UuisPNiK;E&qDsMmg19L1e&lhu7 zV}5h&F~`0$*gFgRwo9iz-JdYMS#~XaBB*FfgJ7k0X$>CH3xi$f%8%Dt^@9Oz#rdx(KG-) zRPeGLe(d4N1-_QRTQvL?!Q*B4>=8;vyCO(Jn<&Ch_U6C!PmG|Uee^2L5pPadKxx`&QH{`D0#&fqHd$~uIp6vT0i~T-ZaIexI+~-;f_uJi` z2gIB3Am8~sxJD5V+3C-LzjWiFr&{qa``$c!PG25z>?DtT+>J+-?&i_g9C^$pM-Co+ zm_tlAa7fx#4te~NL#tfnkZ()zJr<`q1phu1;@gaa@4w*S?uU6y`gI}iC3ZrJOI{g;qqI`ZfuS5xGhi=1W1TMM~YBEK8%$iO{AahE>sYk@o4 z;@(BLdk5|>0f))pf$xKxvl)C=fzt-?IuCA*z%QQd^2))}3S85`cNsXFgSR)hr-T0u zcsL0kjp5}i{M>>kP58=#w^aBmg-2)j41m{F@H-uzSHX98c)tt(-Oz(<^q~iOQHFke zLQn8LsgeVF^9KF7h#p--pR&=bcIa0f^z7qreM9g5>~8N2(c6Qg z>D>`$dVlvZeXKVK-#7i3KHYszUq&paZ~Ddb{p@J^vHCt&aW7#Fx7%E`s}5Hiu$eU@ zUb9xBCTs7j%hew_v#!xFX=<_MmB(zkAPUdZWs7sK@w!!beLChW!8}XM zt%3RG*wX>~1ok$;{^!WC1bJE@*B<2SiJbau6jhGg_#S7I<+vjh_u%{7C%WQ374A&L zy$-lrm$jCD0f&d+v68ELRe;Y$aJmCtXTZ%L{7!@8Z15}u*9+k53C@?myA8Pa27fnr z*asi8;l+qvUHbt~p73=H-d0g*$q#tk4xel3S+$Swdm5f6!*>k44}yPZ^k5?Tuo1l| zL_cnzCr{CrDdba3Zs+B0D&?fS&DV@W&O zHt{ZH%@{{p7Y?P&RSvWzXAEtbFo(7@FQF}6-_jPJ3$$hV4BB$AfHJ$Tq09g;$~5ZbyJfA+^gGurw&o3c8erL3%r zl-)9vviF{)ZG#rmb~7W|{&XVkICY+O9@L|pqs?g7RZYq*UqpMFXVcz5746G3r~MU; z>0loVI<)@*9ckK*jxIS($FxRL{^a9S@HU=Kj2K2Ii@ww8?r-SKkq30HRVbZLAi7X; zk1k17>GCvhy0SNf3ZIzKwQ5Fm-Ka6$Fla$HzSN|fm(I|w#0a|W){bu9b)-9k?$ez! zE9q{dQo7rJAl;q5i|)o>#M=#c-xZIa!^h|0c{lKUJG`zHULS%v`IuJ?b1g9c=N=vG zGsNCN?C(gIoIWAX738u;z5wLBio6ez`!(|S!5vd^&u-k+ANTFYotC)w3hu6j`)h-P zK6tbRmsQ{sL%Rw;f>%A-d7uLP4uj(e@T>-|Gr;#CI1d5uJa9h@{^{^g6F#cLiy!>t z!qYMMa)7rN@b>{8P2jTyyncY+Ch*)5zEj}+ApA$62Pe^o7wE-(^aFpkOU`ig#TLC8 zjQ#|pN4L-?YxIiHuR!$dU;69%|NaR@soNXU(y|y@6w;RFJ9tstHw&6uvYqB61=Fkp z63x;(N;7WPpcxkZD0crzidAi(nDTuTvt%p9H26%>rTG+{mPygwR5bl=1WoriNYfTr z($trWX=>+FG-V6EgwZsTqE=>7lt`e-?_bd5b9yxScrs0XFrA_-K2ucmT$)mqXv(_X zG&OJzO>;4$=}zY;y6JR^-ls(|eKIJvY9htvexVs*)o5mido-&`4$Zo}iDqX!q&cxO zY3|@Gn&%oxaZNQTzWQ*AFWo~4=WEjZtVmiAUyl|J{X&bFXtC})T3pzXmLz%8QXiQT ztGH5P#!^ahpFm019VmG~OG>_$ODP?@QA+xFN_n@BQac==)KS)yy6`-uZVac??5h9q zcHucZ-Wng@hv)6a^TY5uFT8#S<}|^)c9=UB^XFraz&=Op?S%b|9Hq!(hg@%v&j2~c zB5zgX&PM)$6xZY%?iq`_*5SU%xbr6NU4^?h;r^>MW84?;Ci;`Ex0xV-w1Gi0^Vc5eFjZF`vo3u!pCHINrRu^@MHsDd*Lkx{*2%;3qFJ3 z)d_yr!}CG-eh2Rx;eQu;um*iNhFBIMFWRsQ(tzcUaA`8J3EhjE{`Mc2chKUSB*U1wW7{p&#BYH z7v!<>BDuSMCAT^g$+aSpTrTz_=kw*H%6&-6#$_b0X+>iFZsLO}M8)`rp80W-g1?Z{ z?M!kSf!`sA3TOO6|=k>mNxJ`P?Sy!8W3q2Z^0d z5Pw-p;?xV0*QbyYKZsOQbjW%90CJhokX)x4QOCKT@aNbqA@?omMm=dV^~`f8 zKdA@#t8z41C1ygmYRmSA2l z%#Fu_)1n zueiHA?w<<|Mc^?LT=sxZ7C5Z|ue0D541OoU(F#2C$#MNha!mY2j;kucTMOLroFUop zunIn2!pmOxSqV>Z@U;cr8o}Q}c)S6hFX6R2{Kms`2z=p(MEq<(W6ZCX$E?A;Ag-7)gAx6{p8lGH7Tv|XEaWVp!O6-Z7pt5YtOUP zJg0z~lvI+HR7@5o_sP8C8`XakNM`uc?Mx1yqMB0tokoLIIqfRXts_Zin-l4ZmZW=r z7}W?epqd&A)l4@eJ-2nF_sp2|=g%es`+8(><0%O?Q1I$sA<-O<;n?$v%|Y#%_zbH|c#VQs2c`!Ut~WeL?AdyVQX7*F-qbftP5 zUGUZm??>SALVSD?p7#*XH^l3L@%mperx)g3#@zmxA4j!+RAS#y?ESg_CUOizo(ITf zhJ3}ysgJxm*Ea=+~a_|v~b@5+}Q~C+TiX`+jm)j0_XAItpfKD@LvuOCh)=VVh2Cx;i(IJorAY@_?r%ob>K4( zUTxsl6rQ!lOO-d%gR2 zc#P_YbSL9BOQ}{*Evk7Zige^}q4wPm1ZYN zPv>8k9*3Qm9wZe=_a0r4?mE7cZf$QS-RQklx>i$2SMK$fE?wCwT`27$oo~@bIvamQ zI%BdzI(7PsbaGu&>BQOw((&U7(s8|-QbG7_DgT>`l)q-Gls}|M%2#rw{PrWIeAOr^ zf3UrjpWIc-FZPfMC|4?2FGq#dZw@D|i){st4z95}ywOu+@93-7ywn;kE zuf25E{J3=Xft7S_qk(jOR4?g*^;zk{OE>9a#;SiMbL>w-4nm%z$mNB6Rgu#edG{fA zbL3x+JNn?Bb+{`}I#Ku?ccx3n*ZsiVLAbv)IBWzDCvZ`LPcv}p4_^37$pF8L;1~*? z?ZNe&bi8H-I8Os_eQ-Ac|CjKv06vbxOC$JE;Astf`M_H{_*(;y_2Kg=yynC2et6D@ z@AL4U0ROengGBVf7`wjdo^uaDkdeNnY z^vL9bbSG}Abj`JjbnzQHG=8*n{FIKAmzN_QD()%mZ)PIxE_^5DTriWiU#*a`-maB0 z2TYSTU2%~%4%C#^o0>@J&lX6li!G&QszFsG#%r23Vo!>}FW+SD< z+E&uimgUlt0mG%mncC9gI=!StOYcdG+Bp3qdDYMPeFETT4USjA(@R>sN$`a4P%i%IZ8rQFN*f1!hffJ!C&2G=crJ(U+VK7o{*BOs zZsdq=RA7U9DPds{e=I1k4ky%PfPoDzL$2Mx0klQ-6v&C zJ0qrYrGVv%xslGLQ(KzmVWN)m;G9~oJP zN$e^q+b|>?ocsYOL{2UoB!waeza1pCn>5A2L-_-gGL!ddPV@bKKCkC}|6XgX>{PMc zTt(-@8k`%|z*evN13BuOZP4%P1@(NtQRKs(x<1#dBY3Q~<&SD^IY_*la-wGe@3X14 zrMUOn@otsgZ^JuYd(Ro~TI+o?ymQ2RN4e^X7t zq$-xP^>aT{B}ZQr2=Q&t#F*$!PvH*o-NPq$k?uZzxRY`B^3~l8zv;nV;pGaf2Vdba zJmQ@Q(91g+nrk9}00IagfB*srAb literal 0 HcmV?d00001 diff --git a/ridge/ridge/config.py b/ridge/ridge/config.py index f0a3bb6..b641462 100644 --- a/ridge/ridge/config.py +++ b/ridge/ridge/config.py @@ -5,6 +5,7 @@ """ datadir="../data/" proddir="../products/" +mapsdir="../products/maps/" """ Определение областей "Галактика" и "Внегалактика" """ bmax=30.0 @@ -39,3 +40,6 @@ SEM означает standatd error on mean (~RMS/sqrt(n)) """ sem_cut=98 +""" Координаты Краба """ +crab_ra = 83.6287 +crab_dec = 22.014 diff --git a/ridge/ridge/utils.py b/ridge/ridge/utils.py index 2b9d7f9..3de9c5a 100644 --- a/ridge/ridge/utils.py +++ b/ridge/ridge/utils.py @@ -63,7 +63,7 @@ def plot_best_fit(X, y, model): y_pred = model.predict(X) resid = y - y_pred err = np.sqrt(np.sum((resid)**2))/len(resid) - print("ax+b: a={:.2e}, b={:.2e}, std={:.2f}".format(a,b,err)) + print("ax+b: a={:.2e}, b={:.2e}, std={:.2e}".format(a,b,err)) """ plt.scatter(X, y) xaxis = arange(X.min(), X.max(), 0.01) diff --git a/scripts/01_bgdmodel.py b/scripts/01_bgdmodel.py index ad9cbce..8f9f8e2 100755 --- a/scripts/01_bgdmodel.py +++ b/scripts/01_bgdmodel.py @@ -7,7 +7,7 @@ import numpy as np import pandas as pd from astropy.io import fits import matplotlib.pyplot as plt -import math, sys +import math, sys, os import pickle from sklearn.linear_model import LinearRegression @@ -42,6 +42,9 @@ nrev=0 bgdmodel={} ignored_scw=[] +if not os.path.exists(proddir): + os.makedirs(proddir) + for rev in range(revmin,revmax): # if not (rev==341): diff --git a/scripts/01_crabmodel.py b/scripts/01_crabmodel.py new file mode 100755 index 0000000..b6fa9c5 --- /dev/null +++ b/scripts/01_crabmodel.py @@ -0,0 +1,307 @@ +#!/usr/bin/env python + +__author__ = "Roman Krivonos" +__copyright__ = "Space Research Institute (IKI)" + +import numpy as np +import pandas as pd +from astropy.io import fits +from astropy import wcs +import matplotlib.pyplot as plt +import math, sys, os +import pickle +from numpy.polynomial import Polynomial + +from sklearn.linear_model import LinearRegression +from sklearn.linear_model import HuberRegressor +from sklearn.linear_model import RANSACRegressor +from sklearn.linear_model import TheilSenRegressor +from sklearn.model_selection import cross_val_score +from sklearn.model_selection import RepeatedKFold + + +from astropy.coordinates import SkyCoord # High-level coordinates +from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames +from astropy.coordinates import Angle, Latitude, Longitude # Angles +import astropy.units as u + +#from statsmodels.robust.scale import huber +from astropy.stats import sigma_clip +from numpy import absolute +from numpy import arange + +from ridge.utils import * +from ridge.config import * + +if not os.path.exists(proddir): + os.makedirs(proddir) + +enkey = sys.argv[1] + +sigma=3 + +# some static revs/scws to be removed +ignore_orbits=[352,834,912,1019,1021,1028,2275,2405,2493] +ignore_scws=['066600420020','066600420030','132800350010','090200390010','269500190010'] + +fn="detcnts.{}.fits".format(enkey) + +with open(proddir+fn.replace(".fits",".ignored_scw.pkl"), 'rb') as fp: + ignored_scw = pickle.load(fp) + +d = fits.getdata(datadir+fn) +df=pd.DataFrame(np.array(d).byteswap().newbyteorder()) +#print(df.columns) + +crab_crd = SkyCoord(crab_ra, crab_dec, frame=FK5(), unit="deg") + +plotme=False + +npix = 50 +sw = 30.0 # deg +pix = sw/npix # pixel size in degrees + + +mean_map = [[0.0 for i in range(npix)] for j in range(npix)] +count_map = [[0 for i in range(npix)] for j in range(npix)] + + +crabmodel={} +rota_arr=[] +a0=[] +b0=[] +a_full0=[] +b_full0=[] +b_est0=[] +err0=[] +rev0=[] + +totx=[] +toty=[] + +for i,rec in df.iterrows(): + obsid = rec['OBSID'].decode("utf-8") + if(obsid in ignore_scws): + print("Remove {}".format(obsid)) + df.drop(index=i, inplace=True) + + if (obsid in ignored_scw): + print("Skip ScW",obsid) + continue + +# accumulate full data set +for rev in range(revmin,revmax): + if(rev in ignore_orbits): + continue + if(enkey == 'E02'): + if(rev > 800): + continue + if(enkey == 'E03'): + if(rev > 800): + continue + #if(rev > 1000): + # continue + df0 = df.query('SRC > 0.0 & REV == {} & PHASE > {} & PHASE < {} & CRAB_SEP < {}'.format(rev,phmin,phmax,crab_sep_max)) + nobs=len(df0) + if not (nobs> crab_nmax): + continue + for n in df0['CRAB_SEP'].values: + totx.append(n) + for n in df0['SRC'].values: + toty.append(n) + #if(np.min(df0['SRC'])<1e-3): + # print(rev) + # sys.exit() + + +x = np.array(totx) +y = np.array(toty) +x = x.reshape((-1, 1)) +model = LinearRegression() +#model = TheilSenRegressor() +results = evaluate_model(x, y, model) +a_full,b_full,err_full = plot_best_fit(x, y, model) +if(plotme): + plot_ab(x, y, a_full, b_full, err_full, title="REGRESSION") + +# go over orbits +poly_x=[] +poly_y=[] + +ntotal=0 +nrev=0 +for rev in range(revmin,revmax): + if(rev in ignore_orbits): + continue + df0 = df.query('SRC > 0.0 & REV == {} & PHASE > {} & PHASE < {} & CRAB_SEP < {}'.format(rev,phmin,phmax,crab_sep_max)) + nobs=len(df0) + if not (nobs> crab_nmax): + continue + cen_ra = np.array(df0['RA'].values) + cen_dec = np.array(df0['DEC'].values) + print("*** Orbit ",rev) + x = np.array(df0['CRAB_SEP'].values) + y = np.array(df0['SRC'].values) + rota_deg = np.array(df0['ROTA'].values) + rota = np.array(df0['ROTA'].values) * np.pi / 180. # in radians + detx = np.cos(rota)*x + dety = np.sin(rota)*x + + for i in range(rota_deg.shape[0]): + rota_arr.append(rota_deg[i]) + + w = wcs.WCS(naxis=2) + w.wcs.crpix = [npix/2, npix/2] + w.wcs.cdelt = np.array([-pix, pix]) + w.wcs.crota = [rota_deg[i], rota_deg[i]] + w.wcs.crval = [cen_ra[i], cen_dec[i]] + w.wcs.ctype = ["RA---TAN", "DEC--TAN"] + #w.wcs.set_pv([(2, 1, 45.0)]) + + #header = w.to_header() + #hdu = fits.PrimaryHDU(header=header) + + + pixcrd = w.wcs_world2pix([(crab_ra,crab_dec)], 1) + ix=int(np.round(pixcrd[0][0])) + iy=int(np.round(pixcrd[0][1])) + #print(ix,iy) + mean_map[ix][iy] += y[i] + count_map[ix][iy] += 1 + + x = x.reshape((-1, 1)) + model = LinearRegression() + #model = TheilSenRegressor() + results = evaluate_model(x, y, model) + a,b,err = plot_best_fit(x, y, model) + + b_est = np.mean(y - a_full*x) + + if(enkey in ['E10','E11','E12',] or a > 0.0): + filtered_data = sigma_clip(y, sigma=sigma, maxiters=10, return_bounds=True) + filtered_y = filtered_data[0] + filtered_min = filtered_data[1] + filtered_max = filtered_data[2] + b = np.mean(filtered_y) + a = 0.0 + err = np.sqrt(np.sum((b-filtered_y)**2))/len(filtered_y) + + b_est = b + a_full = a + + a_full0.append(a_full) + b_full0.append(b_full) + b_est0.append(b_est) + + #if(b>2.0e-4): + # print(rev) + # sys.exit() + + a0.append(a) + b0.append(b) + err0.append(err) + rev0.append(rev) + + + poly_x.append(rev) + poly_y.append(b_est) + + if(enkey in ['E04','E05','E06','E07','E08','E09','E10','E11','E12']): + crabmodel[rev]={'a':a_full, 'b':b_est, 'err':err} + if(plotme): + plot_ab(x, y, a_full, b_est, err, title="REGRESSION rev {}".format(rev)) + + if(enkey in ['E02', 'E03']): + crabmodel[rev]={'a':a, 'b':b, 'err':err} + if(plotme): + plot_ab(x, y, a, b, err, title="REGRESSION rev {}".format(rev)) + + print("ax+b: a={:.2e}, b={:.2e}, std={:.2e}\n".format(a,b,err)) + ntotal=ntotal+nobs + nrev=nrev+1 + +print("Orbits: {}, obs: {}".format(nrev,ntotal)) + +for i in range(npix): + for j in range(npix): + if(count_map[i][j]>0): + mean_map[i][j]=mean_map[i][j]/count_map[i][j] + +w = wcs.WCS(naxis=2) +w.wcs.crpix = [npix/2, npix/2] +w.wcs.cdelt = np.array([-pix, pix]) +w.wcs.crval = [0.0, 0.0] +w.wcs.ctype = ["RA---TAN", "DEC--TAN"] + +header = w.to_header() +hdu = fits.PrimaryHDU(header=header) + +hdu.data=mean_map +hdu.writeto(proddir+fn.replace(".fits",".crab_mean_map.fits"), overwrite=True) +hdu.data=count_map +hdu.writeto(proddir+fn.replace(".fits",".crab_count_map.fits"), overwrite=True) + +npoly=4 +if(enkey in ['E11','E12',]): + npoly=0 + +z = np.polyfit(poly_x, poly_y, npoly) + +p = np.poly1d(z) +poly_z=[] +for t in poly_x: + poly_z.append(p(t)) +plt.scatter(poly_x, poly_y) +plt.plot(poly_x, poly_z, color='r') +plt.title("Crab detector count rate evolution") +plt.ylabel("Crab count rate cts/s/pix") +plt.xlabel("INTEGRAL orbit") +#plt.show() +plt.savefig(proddir+fn.replace(".fits",".crab_rate.png")) + +indices = sorted( + range(len(a0)), + key=lambda index: a0[index] +) + +coldefs = fits.ColDefs([ + #fits.Column(name='OBSID', format='11A', array=[obs_id[index] for index in indices]), + #fits.Column(name='RA', format='D', unit='deg', array=[ra[index] for index in indices]), + #fits.Column(name='DEC', format='D', unit='deg', array=[dec[index] for index in indices]), + #fits.Column(name='LON', format='D', unit='deg', array=[lon0[index] for index in indices]), + #fits.Column(name='LAT', format='D', unit='deg', array=[lat0[index] for index in indices]), + fits.Column(name='REV', format='J', unit='', array=[rev0[index] for index in indices]), + #fits.Column(name='MJD', format='D', unit='', array=[mjd0[index] for index in indices]), + #fits.Column(name='PHASE', format='D', unit='', array=[phase0[index] for index in indices]), + #fits.Column(name='CLEAN', format='D', unit='cts/s', array=[clean0[index] for index in indices]), + #fits.Column(name='MODEL', format='D', unit='cts/s', array=[model0[index] for index in indices]), + #fits.Column(name='RESID', format='D', unit='cts/s', array=[resid0[index] for index in indices]), + fits.Column(name='A', format='D', unit='', array=[a0[index] for index in indices]), + fits.Column(name='B', format='D', unit='', array=[b0[index] for index in indices]), + fits.Column(name='ERR', format='D', unit='', array=[err0[index] for index in indices]), + fits.Column(name='A_FULL', format='D', unit='', array=[a_full0[index] for index in indices]), + fits.Column(name='B_FULL', format='D', unit='', array=[b_full0[index] for index in indices]), + fits.Column(name='B_EST', format='D', unit='', array=[b_est0[index] for index in indices]), + fits.Column(name='B_POLY', format='D', unit='', array=[poly_z[index] for index in indices]), +]) + +fout = fn.replace(".fits",".crabmodel.fits") +hdu = fits.BinTableHDU.from_columns(coldefs, name='GRXE') +hdu.header['MISSION'] = ('INTEGRAL', '') +hdu.header['TELESCOP'] = ('IBIS', '') +hdu.header['INSTITUT'] = ('IKI', 'Affiliation') +hdu.header['AUTHOR'] = ('Roman Krivonos', 'Responsible person') +hdu.header['EMAIL'] = ('krivonos@cosmos.ru', 'E-mail') +#hdu.add_checksum() +print(hdu.columns) +hdu.writeto(proddir+fout, overwrite=True) + + +with open(proddir+fn.replace(".fits",".crabmodel.pkl"), 'wb') as fp: + pickle.dump([crabmodel, z], fp, protocol=pickle.HIGHEST_PROTOCOL) + + +with open(proddir+fn.replace(".fits",".rota.dat"), 'w') as fp: + for i in range(len(rota_arr)): + fp.write("{} {}\n".format(i,rota_arr[i])) + diff --git a/scripts/01_crabmodel.sh b/scripts/01_crabmodel.sh new file mode 100755 index 0000000..7fb25d2 --- /dev/null +++ b/scripts/01_crabmodel.sh @@ -0,0 +1,11 @@ +./01_crabmodel.py E02 +./01_crabmodel.py E03 +./01_crabmodel.py E04 +./01_crabmodel.py E05 +./01_crabmodel.py E06 +./01_crabmodel.py E07 +./01_crabmodel.py E08 +./01_crabmodel.py E09 +./01_crabmodel.py E10 +./01_crabmodel.py E11 +./01_crabmodel.py E12 diff --git a/scripts/02_grxe_resid.py b/scripts/02_grxe_resid.py new file mode 100755 index 0000000..3678652 --- /dev/null +++ b/scripts/02_grxe_resid.py @@ -0,0 +1,164 @@ +#!/usr/bin/env python + +__author__ = "Roman Krivonos" +__copyright__ = "Space Research Institute (IKI)" + +import numpy as np +import pandas as pd +from astropy.io import fits +from astropy.table import Table, Column +import matplotlib.pyplot as plt +import math, sys +import pickle + +from sklearn.linear_model import LinearRegression +from sklearn.linear_model import HuberRegressor +from sklearn.linear_model import RANSACRegressor +from sklearn.linear_model import TheilSenRegressor +from sklearn.model_selection import cross_val_score +from sklearn.model_selection import RepeatedKFold + +#from statsmodels.robust.scale import huber +from astropy.stats import sigma_clip +from numpy import absolute +from numpy import arange + +from ridge.utils import * +from ridge.config import * + +enkey = sys.argv[1] +outkey = sys.argv[2] + +fn="detcnts.{}.fits".format(enkey) + +with open(proddir+fn.replace(".fits",".pkl"), 'rb') as fp: + bgdmodel = pickle.load(fp) + +with open(proddir+fn.replace(".fits",".ignored_scw.pkl"), 'rb') as fp: + ignored_scw = pickle.load(fp) + +with open(proddir+fn.replace(".fits",".crabmodel.pkl"), 'rb') as fp: + crabmodel, z = pickle.load(fp) + p = np.poly1d(z) + #print(crabmodel) + +crab_rev_max = np.max(list(crabmodel.keys())) +print("Crab is defined untill orbit {}".format(crab_rev_max)) + +with fits.open(datadir+fn) as hdul: + data=hdul[1].data + +#print(data.columns) + +rev = data.field('rev') +mjd = data.field('mjd') +clean = data.field('clean') +phase = data.field('phase') + + +rev0=[] +phase0=[] +clean0=[] +model0=[] +resid0=[] # residuals in cts/s +grxe0=[] # mCrab +crab0=[] # Crab count rate +mjd0=[] +a0=[] +b0=[] +err0=[] +lon0=[] +lat0=[] + +d = fits.getdata(datadir+fn) +df = pd.DataFrame(np.array(d).byteswap().newbyteorder()) + +# BKG +if(outkey == 'BKG'): + df = df.query('CLEAN > 0.0 & ( abs(LAT) > {} | abs(LON) > {}) & PHASE > {} & PHASE < {}'.format(bmax,lmax,phmin,phmax)) + +# GAL +if(outkey=='GAL'): + df = df.query('CLEAN > 0.0 & abs(LAT) < {} & abs(LON) < {} & PHASE > {} & PHASE < {}'.format(bmax,lmax,phmin,phmax)) + +# ALL +if(outkey=='ALL'): + df = df.query('CLEAN > 0.0 & PHASE > {} & PHASE < {}'.format(phmin,phmax)) + +for i, row in df.iterrows(): + orbit=row['REV'] + obsid=row['OBSID'].decode("UTF-8") + + if not (orbit > revmin and orbit < revmax): + print("Skip orbit",orbit,row['OBSID']) + continue + + if not (orbit < crab_rev_max): + print("Skip orbit",orbit,row['OBSID']) + continue + + if (obsid in ignored_scw): + print("Skip ScW",obsid) + continue + + a = bgdmodel[orbit]['a'] + b = bgdmodel[orbit]['b'] + err = bgdmodel[orbit]['err'] + m = a*row['PHASE']+b + + clean0.append(clean[i]) + mjd0.append(mjd[i]) + model0.append(m) + resid0.append(clean[i]-m) + grxe0.append(1000*(clean[i]-m)/p(orbit)) + crab0.append(p(orbit)) + + a0.append(a) + b0.append(b) + err0.append(err) + phase0.append(row['PHASE']) + rev0.append(orbit) + lon0.append(row['LON']) + lat0.append(row['LAT']) + + +indices = sorted( + range(len(mjd0)), + key=lambda index: mjd0[index] +) + +coldefs = fits.ColDefs([ + #fits.Column(name='OBSID', format='11A', array=[obs_id[index] for index in indices]), + #fits.Column(name='RA', format='D', unit='deg', array=[ra[index] for index in indices]), + #fits.Column(name='DEC', format='D', unit='deg', array=[dec[index] for index in indices]), + fits.Column(name='LON', format='D', unit='deg', array=[lon0[index] for index in indices]), + fits.Column(name='LAT', format='D', unit='deg', array=[lat0[index] for index in indices]), + fits.Column(name='REV', format='J', unit='', array=[rev0[index] for index in indices]), + fits.Column(name='MJD', format='D', unit='', array=[mjd0[index] for index in indices]), + fits.Column(name='PHASE', format='D', unit='', array=[phase0[index] for index in indices]), + fits.Column(name='CLEAN', format='D', unit='cts/s', array=[clean0[index] for index in indices]), + fits.Column(name='MODEL', format='D', unit='cts/s', array=[model0[index] for index in indices]), + fits.Column(name='RESID', format='D', unit='cts/s', array=[resid0[index] for index in indices]), + fits.Column(name='GRXE', format='D', unit='mCrab', array=[grxe0[index] for index in indices]), + fits.Column(name='CRAB', format='D', unit='cts/s', array=[crab0[index] for index in indices]), + fits.Column(name='A', format='D', unit='', array=[a0[index] for index in indices]), + fits.Column(name='B', format='D', unit='', array=[b0[index] for index in indices]), + fits.Column(name='ERR', format='D', unit='', array=[err0[index] for index in indices]), +]) + +fout = fn.replace(".fits",".{}.resid.fits".format(outkey)) +hdu = fits.BinTableHDU.from_columns(coldefs, name='GRXE') +hdu.header['MISSION'] = ('INTEGRAL', '') +hdu.header['TELESCOP'] = (outkey, '') +hdu.header['INSTITUT'] = ('IKI', 'Affiliation') +hdu.header['AUTHOR'] = ('Roman Krivonos', 'Responsible person') +hdu.header['EMAIL'] = ('krivonos@cosmos.ru', 'E-mail') +#hdu.add_checksum() +print(hdu.columns) +hdu.writeto(proddir+fout, overwrite=True) + +with fits.open(proddir+fout, mode='update') as hdus: + hdus[1].add_checksum() + + + diff --git a/scripts/02_grxe_resid.sh b/scripts/02_grxe_resid.sh new file mode 100755 index 0000000..189d3a6 --- /dev/null +++ b/scripts/02_grxe_resid.sh @@ -0,0 +1,35 @@ +./02_grxe_resid.py E02 GAL +./02_grxe_resid.py E03 GAL +./02_grxe_resid.py E04 GAL +./02_grxe_resid.py E05 GAL +./02_grxe_resid.py E06 GAL +./02_grxe_resid.py E07 GAL +./02_grxe_resid.py E08 GAL +./02_grxe_resid.py E09 GAL +./02_grxe_resid.py E10 GAL +./02_grxe_resid.py E11 GAL +./02_grxe_resid.py E12 GAL + +./02_grxe_resid.py E02 BKG +./02_grxe_resid.py E03 BKG +./02_grxe_resid.py E04 BKG +./02_grxe_resid.py E05 BKG +./02_grxe_resid.py E06 BKG +./02_grxe_resid.py E07 BKG +./02_grxe_resid.py E08 BKG +./02_grxe_resid.py E09 BKG +./02_grxe_resid.py E10 BKG +./02_grxe_resid.py E11 BKG +./02_grxe_resid.py E12 BKG + +./02_grxe_resid.py E02 ALL +./02_grxe_resid.py E03 ALL +./02_grxe_resid.py E04 ALL +./02_grxe_resid.py E05 ALL +./02_grxe_resid.py E06 ALL +./02_grxe_resid.py E07 ALL +./02_grxe_resid.py E08 ALL +./02_grxe_resid.py E09 ALL +./02_grxe_resid.py E10 ALL +./02_grxe_resid.py E11 ALL +./02_grxe_resid.py E12 ALL diff --git a/scripts/03_grxe_map.py b/scripts/03_grxe_map.py new file mode 100755 index 0000000..837ec89 --- /dev/null +++ b/scripts/03_grxe_map.py @@ -0,0 +1,155 @@ +#!/usr/bin/env python + +__author__ = "Roman Krivonos" +__copyright__ = "Space Research Institute (IKI)" + +import numpy as np +import numpy.ma as ma +import pandas as pd +from astropy.wcs import WCS +from astropy import wcs +from astropy.io import fits +from astropy.table import Table, Column +import matplotlib.pyplot as plt +import math, sys, os +import pickle + +from astropy.coordinates import SkyCoord # High-level coordinates +from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames +from astropy.coordinates import Angle, Latitude, Longitude # Angles +import astropy.units as u + +from sklearn.linear_model import LinearRegression +from sklearn.linear_model import HuberRegressor +from sklearn.linear_model import RANSACRegressor +from sklearn.linear_model import TheilSenRegressor +from sklearn.model_selection import cross_val_score +from sklearn.model_selection import RepeatedKFold + +from astropy.stats import sigma_clip +from astropy.stats import sigma_clipped_stats + +from scipy.stats import norm +from scipy.stats import describe +from scipy.stats import sem + + +from numpy import absolute +from numpy import arange + +from ridge.utils import * +from ridge.config import * + +enkey = sys.argv[1] +key = sys.argv[2] +#key="ALL" +fn='detcnts.{}.{}.resid.fits'.format(enkey,key) + +d = fits.getdata(proddir+fn) +df=pd.DataFrame(np.array(d).byteswap().newbyteorder()) +#print(df.columns) + +#df = df.query('abs(LAT) < {} & abs(LON) < {}'.format(5,5)) + +n_bins = 80 +minmax=[-300,800] +sigma=3 +maxiters=10 + +modelrxte="modelrxte_ait_3to20keV_flux_2deg.fits" +hdulist = fits.open(datadir+modelrxte) +w = wcs.WCS(hdulist[0].header) +smap =hdulist[0].data + +sx=int(hdulist[0].header['NAXIS1']) +sy=int(hdulist[0].header['NAXIS2']) + +# fill AITOF map indexes +ds9x=[] +ds9y=[] +for i,row in df.iterrows(): + lon=row['LON'] + lat=row['LAT'] + world = SkyCoord(lon,lat, frame=Galactic, unit="deg") + ra=world.fk5.ra.deg + dec=world.fk5.dec.deg + pixcrd = w.wcs_world2pix([(lon,lat)], 1) + y=int(pixcrd[0][0]) + x=int(pixcrd[0][1]) + ds9x.append(x) + ds9y.append(y) + #print(x,y,smap[y-1,x-1]) +df['DS9Y']=ds9x +df['DS9X']=ds9y + +mean_map = np.array([[0.0 for i in range(sx)] for j in range(sy)]) +sign_map = np.array([[0.0 for i in range(sx)] for j in range(sy)]) +sem_map = np.array([[0.0 for i in range(sx)] for j in range(sy)]) +cnt_map = np.array([[0 for i in range(sx)] for j in range(sy)]) + +for i in range(sx): + for j in range(sy): + world = w.wcs_pix2world([(i+1,j+1)], 1) + lon = world[0][0] + lat = world[0][1] + if(np.isnan(lon) or np.isnan(lat)): + continue + ds9i=i+1 + ds9j=j+1 + df0 = df.query('DS9X == {} & DS9Y == {}'.format(ds9i,ds9j)) + if (len(df0) <= nscw_min): + continue + + # check coordinates + #print("***",i+1,j+1,lon,lat,smap[j][i]) + #for i0,row0 in df0.iterrows(): + # print(row0['LON'],row0['LAT'],row0['GRXE']) + + grxe = np.array(df0['GRXE']) + + sg_mean, sg_med, sg_std = sigma_clipped_stats(grxe, sigma=sigma, maxiters=maxiters) + filtered_data = sigma_clip(grxe, sigma=sigma, maxiters=maxiters, return_bounds=True) + filtered_grxe = filtered_data[0] + filtered_min = filtered_data[1] + filtered_max = filtered_data[2] + + # final error on flux measurement ~RMS/sqrt(n) + sg_sem = sem(grxe[filtered_grxe.mask==False]) + + #print("Sigma clipping: mean {:.2f} med {:.2f} std {:.2f} sem {:.2f}".format(sg_mean, sg_med, sg_std, sg_sem)) + + #plt.hist(grxe, bins=n_bins, range=minmax) + #plt.hist(grxe[filtered_grxe.mask], bins=n_bins, range=minmax) + #plt.show() + + mean_map[j][i] = sg_mean + sem_map[j][i] = sg_sem + sign_map[j][i] = sg_mean/sg_sem + cnt_map[j][i] = df0.shape[0] + + +# Calculate the percentiles across the x and y dimension +perc = np.percentile(sem_map, sem_cut, axis=(0, 1), keepdims=False) + +print("{}: {}% cut of SEM map: {:.2f} mCrab".format(enkey,sem_cut,perc)) + +idx=np.where(sem_map > perc) +mean_map[idx]=0.0 +sem_map[idx]=0.0 +cnt_map[idx]=0 +sign_map[idx]=0.0 + +if not os.path.exists(mapsdir): + os.makedirs(mapsdir) + +hdulist[0].data=mean_map +hdulist.writeto(mapsdir+fn.replace(".fits",".mean.fits"),overwrite=True) + +hdulist[0].data=sem_map +hdulist.writeto(mapsdir+fn.replace(".fits",".error.fits"),overwrite=True) + +hdulist[0].data=cnt_map +hdulist.writeto(mapsdir+fn.replace(".fits",".cnt.fits"),overwrite=True) + +hdulist[0].data=sign_map +hdulist.writeto(mapsdir+fn.replace(".fits",".sign.fits"),overwrite=True) diff --git a/scripts/03_grxe_map.sh b/scripts/03_grxe_map.sh new file mode 100644 index 0000000..b3b7863 --- /dev/null +++ b/scripts/03_grxe_map.sh @@ -0,0 +1,11 @@ +./03_grxe_map.py E02 ALL +./03_grxe_map.py E03 ALL +./03_grxe_map.py E04 ALL +./03_grxe_map.py E05 ALL +./03_grxe_map.py E06 ALL +./03_grxe_map.py E07 ALL +./03_grxe_map.py E08 ALL +./03_grxe_map.py E09 ALL +./03_grxe_map.py E10 ALL +./03_grxe_map.py E11 ALL +./03_grxe_map.py E12 ALL diff --git a/scripts/README.md b/scripts/README.md index 0968533..0bd67c5 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -3,41 +3,20 @@ source /venv/bin/activate.csh ``` +Номер скрипта означает последовательность выполнения. Результаты работы скриптов помещаются в ```ridge/products/```. -### 01_init_events.py +### 01_bgdmodel.py / 01_bgdmodel.sh -Создает начальные списки событий и помещает их в ```uds/data/processed``` -Оригинальные файлы со списками событий задаются в файлах ```uds/data/evtlists/*.txt``` +Создает модель фона -### 02_merge_events.py +### 01_crabmodel.py / 01_crabmodel.sh -Создает объедененный список событий и помещает его в ```uds/products```. Этот список событий нужен, в основном для извлечения спектров с помощью ```srctool```. +Создает модель скорости счета Краба -Попутно этот скрипт унифицирует оригинальные списки событий для последующей обработки. А именно, корректируются слова OBS_MODE=POINING/SURVEY в зависимости от типа наблюдения и производится центрирование на одни и те же координаты с помощью команды ```radec2xy```. +### 02_grxe_resid.py / 01_grxe_resid.sh -Для запуска адаптивного сглаживания ```do_adapt = True``` требуется запустить окружение ```ciao```, так как нужна команда ```dmimgadapt``` +Вычисляет остатки после вычитания модели из данных в единицах мКраб -``` -conda activate ciao-4.15 -source /eSASS4EDR/bin/esass-init.csh -``` +### 02_grxe_map.py / 01_grxe_map.sh -### 03_init_obs.py - -1) Подготавливает списки событий в разных энергетических диапазонах. -2) Запускает ```erbox``` в три этапа, чтобы получить рабочий список источников для ```ermldet```. -3) Запускает ```ermldet``` -4) Делает кросс-корреляцию с каталогом Gaia-unWISE ```do_cross_match=True```, которая создает три файла: ```.cross``` -- все пересечения, и ```.ref``` / ```.src``` -- входные каталоги для последующей команды ```wcs_match``` -5) Делает матрицу преобразования координат и корректирует списки событий. Для запуска команд```wcs_match/wcs_update``` требуется запустить окружение ```ciao``` (см. выше) - -### 04_mosaics.py - -Создает сборные изображения (мозайки) в разных энергетических диапазонах. - -### 05_scrtool.py - -Запускает scrtool для самого широкого канала 0.2-10 кэВ, чтобы спектры имели самое полное покрытие по энергиям. Список источников берется из 0.3-2.3 кэВ. - -Вычисляет ECF для всех диапазонов. - -Делает принудительную фотометрию в выбранных каналах (параметр```forced=True```). Внимание! ermldet из eSASS4EDR не делает ассимитричные ошибки на потоки. Мы запускаем более последнюю версию ermldet (v1.56/2.18 esass_200412 Jul 2 12:04:46 2022). Для этого используется параметр ```local_run=True```, который высвечивает какую команду надо запустить на другой машине и ждет ввода. \ No newline at end of file +Создает карту остатков в единицах мКраб. \ No newline at end of file