From 0484bf5c157d6286f7c2f1b25ef6100031e0ee90 Mon Sep 17 00:00:00 2001 From: yan Date: Mon, 15 Apr 2024 15:46:32 +0800 Subject: [PATCH] debug --- .../__pycache__/csst_mci_sim.cpython-311.pyc | Bin 262679 -> 261046 bytes csst_mci_sim/csst_mci_sim.py | 961 +++++++++--------- csst_mci_sim/mci_data/mci_all_9K.config | 8 +- 3 files changed, 469 insertions(+), 500 deletions(-) diff --git a/csst_mci_sim/__pycache__/csst_mci_sim.cpython-311.pyc b/csst_mci_sim/__pycache__/csst_mci_sim.cpython-311.pyc index 12d60db044ecf6362e20dd28e52394913b24d0ed..6c1399c8b032d263a2aa6d974a1c89c0e29d82f2 100644 GIT binary patch delta 27690 zcmeHw33wF6_U}~p%w*q4LUyu47RUw!WK9SOYr+;}Ck*LH$R2to2w_4%#Vg9{00(dZ zRPZXe-H4;Ofy;Hp8z72|m(}ZniUt)ITzKbHCo|~`BJcg*|9kIy-H~EIL9#8l271%Uuwh-Uk|8~vnanb;k#Mx(^4YVXIQ7y_H3ydQp$98v7**&JwCKT_bEnW zZOI+>A!8@4p(6@)?HI79rEEb486Yt5TBPF!F!ke^Z$?H^*R29o3X{jYN`} ze1&47NRkCvKnIywTStZKJWze!`t7KPy-7;gtASa;@y>3lKZCF@Ihlauv?-%h@a z`O{GK`3Ba&dt=XeS4Irr;zr5G}xC2Qdy z0*kU{8jY;L3dRL&s`Y2%P+cP$D_RDZeI{aiiLOj!&DQ+NXx$wcyrQL|awrS(LKQC+ zuMsL_vR1y;`dn4A?k+TJY5A-wLC|ePb+^^8_QNoXg6>DPF^bTFyat2hChH%y;ktWJ zeYPd4Zm2GhqJ5Of6!Iuu?6eGdjkU_$U--9R-DNJ+e?<{ATfa9C(sTb_8^IeNIXnp4~e}b9zkIjYHRF5QK1}1xx3Xm$&#ag8pFR5TDDm3&<3m^ zj9abIARA3)EP0Ey@2bdw&E&q4Orl|-NzNtTd+AhHT4OLTlB+3z$|UEJMs+Fgwl=NG z35z9_@@97{tv1vwlni$9mNKKf)B546RQ(3BuD80~nx%Ue)tHtMw}P8~5``YHe!OOE z%nOR=B13JhvAWT(FmGJU0Ae3OzyMB@}kCv}`GfW%7%3s=LLo#8|83S1cn_8JgZ9lYnOe zp1e-J@mNZ`wydJsDDR*GzhKQ=TcZ03O~ox+*5-QXj}Yp8>%X@QX2-0-TgOKKivn^8 z+e+F{(vFhWQCQ2?tbklRw|mr^%1TV728{ZI^_8vB>;vnWts}gS6Sj@S3)UnU-F%v` z+6uF=hPP4q6q&3Tt6f`VTxOD=C%cVk0mMmR&-SrSzBcC87Eo*dz< z8~Fm+z9gGFRZDe+$xL-zZ8pl2qp>b(<*x~Gkpxwd{0)h3NxW|L+?Esc6{>cjqC#$` zu#SqFziqrO92@x~Est+oqSyCD^{UvCymNy-mfv{ZGH-1B?&Ft`5OO@l@TK+Sp1%47 zI#0BE?d`>STT}K{2Bczu-5sl~&`@QR<*%4^&)x#vXbkPy(!TeA-qY^BxWORF27JCe z_IR~7uV$@ekgMcptZ`46CfO*TSpFJgk%0-k;Mw1M5%Qv8oe-fTw8!!IW?{}Rfu zTU|XxYLe%m!z;b{1y!P?e#7orZZt2ftE)B|YUOe&@gfoxBq~W%k*Fq7L!y>M9f^7p zi%HNUQKtSxZX(e@VhM>`NGv6>j6@@e*ySWGBzP-w??~K)eI+l~3hT~8liA}f?T2>i zb5Br7KcZ5LGQU!G3iy};+!gp43OG&TEQxa@&Xf3r#DA=)5lKS!v6uvnyEUuXLra@EnE9MvCcgB-T6Plt|O@k`eSj*0tv|*#hg!=fm?I zsoQmQqneZwxs!tj3M3xJ>JnpFne2;a$(~BZm6FBXRM#Mv8Y^n$W7eooYS_Y-yFc;J zHy^_k8hcgT2@?`b(oYR-kCEQi(OB^`W$woQ*m-TLeA2?PH+ zzjlvlRm)0>>MBe|rJws*uonLjzwiph^(%?rNc@4q?p;*BsLothTT)ituvFHF7|uut zBt#Ng5=u!a<8E~7PJ+`2Ptv?dc%y*W1uVw;%O5Fha7)6~Izh1XgG+C*82F`*dD#vz z*4G1Kyx0ns1n0b11RV5Yz2N~b774$3v4>d-?Dl4}A>D_)5uCag`S3pC#QfLP%hGFC zFPI_O{-F%0vC zq$mu|_%k1t2^alYvXC_l_PetPPbs<4OQVtcH+n*gKZ}L=0n9th5vI{tG?krM&`1Me zdk;+X$p97|m`#2v>mWE6z{Xv}Fqj)e^Z-);8wNf-SX3WZ{~?%Bim`8Jf*YA!2*U%} z8a5Or1hRDaF_8J9Kk0sb5$V~9l!o&_kA5_QGy7qIH}z$~urY|;!$!i>wJe0?D2H6- zFiJV(@dK>+Lhxq!@JKLg@hWJ_traYOO}Ww-F!p4Vbe{9}HI0VczAT=Ng;PE8kPR4L zFdjle*cd#N;F*wA22a(pP*@+ryfN(F5H`z;*JC4QQtA&^Ls%bGKxQb5jTlGPD(dYt zrWngOX?zv7G1S8?5s?I^ z!n>HFB3RX*n0}35eIWb|(Hj~gnKulHWWMmuVbKegMKT{)K&N2~1b!O9CP7^!3xWL@ z$|fuNqFyYtqYZx>$!gdXNcUq=;2+JZd~4sa_*zC={CEbF;rl4oKVoW`5Y|~@P1BSv zfN7B1i-oabnA3|TUOUQr;gTNjzAMpeJ{D(I3`<{drXS-Ja5Xsol3 zw@Y$qD#Q3Zl36;N0GGn?E-OlAec1w7p3I`yIA8&IYZ&-_6|m45i>1s`vA+~fB(r2T z5m*Wf@tW9~0c>=Jf3i?*QZZsS69inm9aG2x;z;wxw%cODm^# z`tLV(dty{uw3ce8oyu*gE^**79!GgBBF}yGpVWxZtHbQ9j4ly(G)7;ktfSM3uv9xs z7}aZ|cCF_|&$W-6AL-00)xE+N{JbE2l9TXfi2! z5$R?yWUzivmWB`K$1ngkllBp5uac(EM6ErI>4LEHy}7ew(lV+ksSz?W@orcSyS-TG zJ$LAs1!^)`EcOCTuql%@1g^NreaT=XW z9ygN*;xP+nN@4pz)(cJ!WJU%sGMi0hkHbCL>|)lp729_}lJb zGMXpJ73DZ+qKRBNgv3`GxMK**W9Q)hoh)>mouZ|OnBAwbPO32AG)oE1ZOB1bi*`Sp z^)fFqEvl$BDxrZ?x1l6x7F13n?{pINkUErQ_qOXK@Xlh#fU1W-#I8s$rRY3~m|5;z!*}mitiH-O3l$naqVt<#0+m z0)<^?USu>Ia7F{0N~ytY7p1Z+68+_BO1Fl@2ey9?W22aUHhIl~gb}Ph*hqRAiESBPzW4r|mR#f!;j?+BOOw-&H_gm!m$au1sW+A%Clxab#Yl=pL09k*FlGiv*uZh=wWiSs<7ySbz38 zJW|2ZQobe!-XUKk?Hd#q;?dzk5_x?~VeKR?iEy=oWe0yiCzw6#wh8qcX=yONk|l-n zo6wB{s8dsNL9AbHhP9O}Y7D=-MQq-7u^K1WW{`>B_Kv9c0}^@W1B4Qg2d zgMfOhL0@HX+gQ(fGe6i}&vIdYJ>veTXW^V)vH)e=^`jQCRWC9RTc5>jjYz>}7U@b7bzhcPtTUfNB zz3I>{-@@iGf0(uuiMv0nRt~nMm|cI+E#rZM_`z1OjD0S+`8H+M3X(6(u;7TV9F|+y zB*C{4HoYM9g!2~mpm^SFfKBg+!Qj6FquAfX(j4~hnsA8k3yW7EJ$V>*u3$@HdtV#| z5$6EEFf2(*Gi%^sIJCvVO0S@`dyqy9#B%c=67&Qdss*e>LhZYf{RZO~V|hJSv17Ub zk__X1NP-wo5`JukyHWr}-g~98$;kcJq$%ZhdgV5XA2sF*O9mjOO zd^F_?a65}vAaCE!<4RwRmEjA9)m=clS9ezEzgDy1XuSr=)_ezR7SML}P8=EVL=hYU z^=mNfOPv3JHES@_e(>lTHi`M!ep$m_Be-E5PECBcmJROIcLXE$D{ri&(d!ru91$_4 z)A-&xC$QIgwvi`fducu1am?3tX#;y(;ECDZxQ7iAkV^i(8K+-jSKvil0Cnq;QSXDp zZkXcWE$j~XdJ8s7`xb0a;agF2--_B;(y~cgLYhR{!=&9q+A-1&p|%xJ`G$hhNQ(s2 z=8`s^w7W=KO4@Uz?IrCq(%Mi%%DQO_5?o;$>a`HOjrBpcI(!?8a;xX`aWNS1gvp)s z5$@Z@67V)Wyp1(vm^-MWbAFdz(w1+UaIw$5w@u9{-#+?@OQ!>K4;9G0=zfldvh6Io znR*jGsVMV~Cyq?|51}6=QB@^AxJEV6cVW<`jVjeYq=k3X(ZA~ zWRS=dN21IkHI^{(C@ONo7|rY zU^9CC1@znw_P5AEl-PDk>?YJqRHDw8fI}etJ~oUQ;qbj|>R$u(Wk?I;Fj#ya8~E2~ z>?*F_cO%gjE~027ahgUROxYYt86HMrIEfJ?Mv}-OkxPQ!Y+az28O8EAZXv&960_jP z``Oe{^xj1tOs4T9CXgUfDNjVrE-o{gH3Zl{}mh=frfjd>RJrWC;oU z4M|r%iHOC7rAM=DB4LK+o!I;pirMU-m>th2LO7n?*CdmCGTT9FN!3(OS|bTg ze>vwnNCpRAmXA}8XJeRMZ(NEm=uFKr&AP}7Df~;)rjsCAs=OKMNvBIG&SfN=%#|L; z$~(GdGV~zaO49PU2fr)?@u#rmB)+1sD6+1iQ_eMRGa$U8D0! zh0OhIV8HziyQXemcHgq2{>g3r$*um$uy{Wk5VS|PYvFE>hdqw^r?mN}wDNz@vY+kf zp=}#H_O0-?!Nsl9=eEu-YhAS1HuZ5<=NW$1-FvP32A{ipT6+$DZp_i1qmR1hx4GxH zy64*t+nAl{d51d7w)=VZu(sLp5Pno7)?KRHw71jMjHR0|jzm0)WS=gb41ZgwojkQ@ zI)4vm*W;^nsSXz&P($T!+~_TwU0+Lsj%sB{?|9dUWB(Va^@=#m7UV7P&zIRv?Dc0~ zVc&^Uf+^uXlqL;1@RhIg)qerSO$pd_*aj;kyIVcJAT!}Wp&DCGa;P!}8rg|brbvn) z?=|+caLUW}<7;fWIF4tqTfsb)_r4|Mxq!r7)Ce7=kng5bp8K_=@ztcZ-$in(v4y|KtnLZZ zC@IIqk@-RhX4jX><7#mvHIaz4Om$YCF_51OAK6wk&!()v-g z@}~+tg5_WmlnuKd&LvE*FRm+)4a%|xUIEU$@-+Z_$pL?H&L1aV(1qgtF+G3qb#|?} zf0o^+w63AnTvA(C+h~+w~N7Hk`Jg1si{pFxLGJ*@hH1=*rM#xbt@D7tJ6`CZ<5 zwuqgC{pVRPJA1Et+tGV<9eq$$`Ip!s7j>(}lU2#0(NM z;ov7Y@}6ib{tqkBHS_C$m81-3Wvc8eUgePvh1sNBf+Bo# z&#Qe4MaV02Gir9v4rn2c!=1I;tvc6X3MR~}vH}pYzGPW!lg;=gnqvgVGr8!E;5T&|2>tX6yCVRHjfCVdZi{|_o%Bino6oGaH$NB&vBzE zg_uZ%aKzk^ae(hw18cGQe9uxnwM4Pt!WS$g{>OQJRM`Fry66Xb72fz9Xxoo6uxSvq(J)sDrzMx4!54xQ%_>;y5 zGx&_@7^&>RD!#Ci>X~Mf=HxuIUcpzBLZU4yQt)T~e5DU>1e`6# z!|*5}a+>qnjiKl`H{=<0R=P8z{Qk3%Cv6v?nZN~p#+f?>UcZ{7XXnUASGSMX3yFT>4LkySC86;nuSI=rmiezM&7g-!&0Nk z?uJ_K%=~FM{dn&{A%&&EKL=v7NQJKl3i!$pdS(lK^^>u{jlH2cTfm`qwh)XjejBre zf%vljsca#gWx!cFWJ2yBAx&H=XOS{=3{7||okg;@!q3l9#eGJC-;~s$W@al~{Sm53 z-e>cg{G33%>-_>=8zd|iPWnLhU}1o^ur;R0Ry$aj<34sWHZ}Py`SH;dte{yQbg2P@86e~nr?HMGv&BXok0!{*s$O%Tqyhx3NQubuOG_Z-Dy zgK5P=DsOTZ_^^cT6Krzpi-pjFKA5Kq)OTL+K~aSCgjbF;G8!D&Mbk2qU5f*uj>TQP zoltI-2jNNhx>$%v=A_{jD&!#&*J_lsU=`az#&jWxCBv-g!esV0+hfy(yF}rgBev{W zLQ6F4Ss?WE3G1?`V}|WNONBI{*-T@ij-T}XUX055Hd^|jOt5ryy38jJujlb&H;1u=jl209SgE{H)KG1!5)$ptZlBZfK< z!(0%09g=s*;? zASQ9dWCvo33t}oqOmiTLT@ceb;${b8h6`dQN6c~{@Y|zqH^Ur`nCn2yb3x4Kh!O{4 zfeXUG5eprNQWu285k?21%mq=-5sMs%3Kv8rM^rfw)h>t{j;M7Y>Rb@@9I@DekX;Za zjxakA4K9c!9C3>SvD5{zj3XKyh~+K_3r93@L=bjKE4m{fGXuN7^y5M>^{gu$QLX9@ zly2ow-R3~t?t)m&5qCHe`&)MeVnwJ5nCLHtu6@Qh;0tUb{E7Bj=0x>xX%S~KSw;^K~lfv=ZME0h$mbSPjbXl4#d+gh<_;0q868D z(K9a412B1^5Dqb43OB)K*0LOQ)VteH35NerDEgJlHq-N6ynGEBsD}Ru4OD}! zNdwi8Ytle9#7P6y;Z7PDbI64Ts-bEcc<_)=7)XRq4eF&ZQZ*-3L$1pS)lkj}VfA4l zGLYz?+NjXMQ)DDMs5Wvs2m@b5BQZg>QDK7X$Vg03ZRAW4J|!cOK(&#RKp6fS8j1g@ zjhz3%7BUj`QyV$;gNtM&)~7ad)(82oqmc-o+Qgp9<< z)JD$8An{!^5*bq)IT?c{G7=9{8#xby56~Du6ij7wQZPt)PlygA>ZJlG)C*R>hfx#j zQX4tzf>US=Ai||GItdpHJc8cDwNwCwYr#4)63tQ@In9Dk$w&-KZR89KhQE(SB3Ei7 zCs(kAjKr(dM$W6?A{mKNsg0acLH@tcNNh@Nc>faZZ)Cp2H`gH?*I39-=YqCf#fJRF77|I}QhUxWbJ%{(ub+`5}Co4E%n{Ib0tP<<8_IV8hTw~q=u@Ak=mvv zMrvDEVl?9FHHeWK{wKsp4Z0>VQbSb4XfRe;Pi*RHh?5wp!yUxP5)Vt8g?`;;$kbpp zHJZ;vg&L_r*Q7>j$aSfa8cNhCpR*!$cp^MSvt+A;j#)CbHHobGoD!)Gsft^s*fBe% zwx+=X5#5LnsSTNm+c|P0LZr410DmpI5gAe&1}Sb?&xm0(u}OqT4HyE4$dkyB+Av)4 z%++Cn#D>(?9H`f!8!;iZAy09;L~cZd)YbyPHRbu72B{6>6t~@abR!<5woZgraw86; zHcV37rnsRSks!5oDy(!vHzGi4!*s=sxuY8~AhmTSWVoXnu^+Wzj^g$-xe@(QTj#@B zawF=aHW(B)g9o}1_fcCV*yw?7#Cz0+a>Xsc6WxgTsI8Td!F?8i1DZmCdG~6)HA0%YU>i%<%MoUchrVuid(EVrbc{6ZMDE;Z!w>99hIRA*LmMt zRJe``fLzB1(<7Rr4qvULx5`IUc#g{I;5j~U!bj}Id5#Jo-L1G!jb;3HAg-eZtcL}@ zm@x4iwPBOuDa0$D#C6nw&EW5cp2Tz122ectY&sDgwRHy^A~zyAYQz1CTdqIWEU_H5 zbr;n8qZ=_Cwc#Pf?fl(JK8fY10ec}P06mG}s11)Qo_yM#D302?A6m(c=#AR&q~bQE z2j-LbjoNx$exrtVi4dS4PP_J)3^g`T#Z38JT%O=~)m8Q)T* zu2wS1)9JOPUU?fWe@!nugYbQ!-8*mEwBiz4Vw{&*BL4%<^bv2xRr?E*#4`2`yqzR| z6Z;LmyS(C6x3Cg7CyA-Wy-V!+Nu$Q+=M^i-R7{0ijlq!CW>-LBB6$j$i z=y(W&yZeffY`CHyCVd2)?kf($kI(TC2I>98Ubqs!sGs;8ethkgCgRuRwgGAQ(TW#; ztBje(chc}sx;Q}iQM8>*7b}IZFEG&Vo;Nl>Z^E?cg?6{mxy8At$cZpFOH9*Wq{8{a z#w>9$u2>TKi^aImeMWzAct~mme$ykHW2B0*vIf**@arUVg9#q)kK0Vtn&F%NVkW*; ziWwmO;z93>?JPwug1CWVe;wbP0p<-9GfChl+5^QMY%h$-7Ki$z;%8mv3S7GrQ!g9g zNVXUYPiKpdMYmB#J|Mx@et$?>8odZUL0T$){w5~j>wY}+)@9*(t;z8IATiqgBede# z6u3M{jM1MYY&t{?7W;W+kaH%9EZelf;sc&+6nr;IjI$->i7iYI1hvDt(PB2MvW4c0 zLj*id8zXuN%@P(z`B06oc-2$A@)z6A^=N!01z$D9*9`FGy!@56o#Oq3UR3kN0(?oX zmQML<*R$M<`UYg9O8B!sULrv1re5B%hoFG{fM8lCS1Hii1iKfmCo3tftEn;KdH`ji0!7A`(D4-@eBCl%M@MVn zl6x!b0d;i?YhmtqF%etnt>eYM5g$`Bd|@dq3&mw#v=f>iWi}b!9xo2^`ZobikvI)r z6U1lTKKI4NQ%g#0Z%q(0{Iv61V@u$-+2VP)I7ti_+V!@xb41)C=F_-7f7k9#eK@ss zz?|0Dxp(T;dT*S4)IH{uIvmzDiVN6L_@+^e(3ZBwO5n3xT&b+G*tA^q^w75!E@+K4 z*mmD8{^(}=y;j`f;hH93no02BGV$O49qDcea7oc!Nihe+lKu&OcA4Ay40_J3txwL; z_}sSm+(Yq4daw~8!Uv6Z_=A1IwM1~ zR`Qgzv?`m{Tr0oPntAJ++$h4%5Xn2NQ`Ja56g4dsfhUf!P+y)Lmrd^GXzQ++pkRx; z`b9^R2PO&a;jE{JhvX;utrjFV`gAQEb_Oy35%WTtJS}>Qhs9HJ>$=+n|8UAp*{DO- zO4=suBx>2(Pz&l+czfaf6krx%(yuHsjaxQYJTcd_1c(24KI1#j@S}Wn;uEX52|odg zrnI|zs`Jw1g&0iCqh|dQCq0W&NZ*2qZDqER{8MGLv=k9 zT(4*M%8;~h^eM6M|9AyNb*lg~J#adH7BBezA7TwIAThLvx8P=3co#7(eA^=4#7~0g zTR4BX`iwX<#9QiN@y6P~4?Flx6vZWO3v>XV-BDDf=#37~h?AoJ1WOxNH2K2Q2k`0c z|BM){m!egE0p0y#;P(S!wCE?r0((}B@bRzB$D7{=|NC|Xb+yI)nXTF1;-k9eanurU z-gU8zzaCPW zGiltFDd3Qp0l&Y1m6y>aJhMx9R+sSpN_fQ)(T5FyRnLq4*g$ylc`;KA#BV||y-2wF zyx7NWkkXq$$q~Giy$^~#;O2v3obD#-Unq*0P-zJ1mOxc@^!tax;g_+!47&!fYyLWc z@Wnx~cSvVcU2{r^^jE}WKXnQt90d!*EN8&s7qG4Zr7Q@0Q5+^HW$1?iX>jI6F_}Z2 ze2FVr@bODXHI!;i#}%UpNX12<+}s0WK=W{2t|Q|;J4&djq^$(_PjX~C^kuPku##;} zCG9fM&3;E$JAE1Lv}2(LveC|Bz!R_F&PLKJ;@^62)bqZp6RwfQD1AQm*T=E#2RLu1 zW@aelb2pe~euZ3 z#X<+t=dc)pWSd|;Twq18>981K^}NTAb_3J>DQ*%BeO19JOn_o@?~Dyj9TpQQj%(R6 z(GPaLDtf8`u9Z5~QDA3zr(Gkg_!?o;uMu`L4^v;H85CyGEM=8`&W1LVDpcujU&D>A zR=g&L>biv;bO}ogp4HJ_1UZ(sm)V`sxE4)mFSGuCd+BJr{};Cx^+n?Cg>UBsU%Vp* zg*jX19A`(*@5kKZVh|+1EBdi{Q1YhO{|2G|x^d{px5U)u`Hn`{gEzWZRhv@U*?d*l z|JHsf2A|f>j>%aY|L^om&Q5LBJJ@lxzax6z;HCU;yXxW1lKPVJR@v=Q^Iz^MYxzjC z^KmG3-7F5058{Gu8hPT-^5|26aQ;0p2ll=v20;HKVjjdD!L~oSRqWwMU%>O99BKiNX20<0{n7#>?q_1F@ap%ek}x^dPVdQc8%gY>4CqG5a8qJ+dDB{4%(o?IYJ#c z2Yu*@gWj!}?{Z%BXIru8w+6r=MSS77=p)2i;OxeZIfUk;*pJmvd65uc&lY@B zXSs&`Ww}xj;jEejxO7|$=sJKYg%VAH|n+& zc77nni_0tw7O!VPddUy}Mnxl>`9OSSf?ocOCKrDv@dt^kBs9p1WrhN`Qr40#kPu1y zIUL^pP>e|B-;Wl{G{>TBAWXCU7qaN|N1P^7_A(hoG37z%2{Ed9D24Km0{PAweCrdw zy+$(m@I5wCNE=DGRMK)#Yt;8ktm&61DZ3Zr=kEC7GzOQWg1?{~<)g3|@`)wk-1?xL z(vCx6cbDsMf1qbhh@+L=_jBQo6XGi2H!ZCCNG#MVJGMdFN8;}$W<6OlmMU8XUVCE|^y2%fG+@r;B_!e~Ccd>x3ThR1;%pGF-HMn<)STal(Kpdl z_Yh(iGh-Igr7@yt494l_ygE}&Md?(-GE+=hMJ;ZFu?VKv#VM(cQHJ{Z>SZwttBppf zqP9E+H`}JyB^@C70L5Y>@fo zzaYV9>iH+4C$S`Uy@cOF;{-d_XUD`mpNZ!mp&g_!nnh8j-@hb{ZzD1W*@ZG2&u8)Z zJUne?W!Eli#_i~ooiJ!}RF0qu;G1^v{o9@2e$rNWaybdUbBFUoJid3yak@w%eF1xh za?dw=&mpbbm&kmo|86Q4|6(1VEzKR9-SLmK`4;lL)*A>zRj=$zzl5|p3VVn&zD@nT zq?sseJ89HO$e)o$U+XFx*7K_8ACmq;hDtI#Nm?ywZ<03t2+3n4`1iq&k+zu3Uy@c% zHFuIuZ%2{blfPErw;XqwlXxam`5tKY4$kZ@Xz$vBLy({u&DrC+IO57GpUYz-WN|L*iW8!)An|GU%u zte4-~q1zUIA@3^Dh32oj9t=~|1DJ2OGzvvLErV`SN!y?DIPwD zZQqM2%|@yyzAgTz)UMal=~fc_R(XQ7^Cb99^A%})6PZwI)zm!YFG!e3@RrH9QCL99 zIc{^AH?u7i#+wvxVZ6aTrp1SV(sX#6;!Vfda(JWXZI3q^-q?7{xt%WYcFHf`&nTO; z0?HeSQrP^1_*0L0p1P2BO&3`J_g@xcd(3yI@Gi2*)_Pf-!V0Db>N49kF4Ek)z)hFK z+BJWEP;$Z3)>?NL2HBSXB!;rAk1;dJdS$yut0RlNHhJ$urNHld35I`p0>v?S<7Y8N zw+6%4!Vf=-8{Cu)?7d*!FJh=qg}j=|&AW^$_{T4zq(4LV!a<0*=qJ1X1`2`0|SQN;7W} z?G)gAH|aq!@ScwT(2??YjCk+HciE&rG9(WpF_OeY5~U>g2Ag+~b{7f0zYy;<@1fJj zNYG!h;eXCV{)x13p$PNkV=AqvUzTR7H<&AMdy+^hV_y=ZN$_h<)zz&rTyY;(Qztc4 z8%M~y(TdS&OwG3Czlk}bt@WySf#pb}8Sb`#n?7Fmh}%JO`w!i5Q{fk;?e#wZQ+RQ3 delta 29325 zcmeHw33wF6_U}|zW-{3)`<_fzn1qCceMul8YylDiL{=dV=}8D#>6swJ2?0e>Kmi9W zP_MWzhzdj<+0&t3qOajwT9rwgM$1u@cZWWgqpOJEx|D)Zco2pfn4JlCZu)D7Kf}Q^Gn5{ zYc=zyh{Hfm6^DbICXN7ECXNJIE{*~@T^tRvLM#9|Lo5V2Qyc?wmN*t=8nGDU zY;hdOIer>(jyN8?&lM-Y`8;tV$ob+VkPF1gAg{G7?5H7PDL@v9bHyoetb*fII9?}C z1Atj91G!i%2U!j7h$u9Ojh6Wqq7mnb(*ai_&KGCEu~yJb(ugzHYQ$NO2^!5)@L> zUm{jQ>^f1Hr7@>m5F^DX zG5U&_7%^6i6XV5%U5US&=Oi&%Gy*2&cQL7As+cCGiy3{UrIDiS#L;~r%F*>%upXKr za~4*}o>5R;py_~0?$?&>Z-ifk)Ld<`S}j&%b#r5*xk)rO)HPX*W~sKd(b8nI=5kb< zkW<|jn3Gf8+*DInTT|Cyart|~ZAs}@* zd_;@Y*5UG5Vrgh``8HZ?=H+Gyyj$D+a&oMeYK-IZbqCt=di!-N)yZ)G#s#xVN+*<; zOjrn>mbLlkN?#vA3cjwYMY)@-!aHMTZDFEcL48(_>IU>r2S zIC!Cp)-~C{&5AZ%PR?0EC_OGZWFPbXfDf_FHz!B5EN-oJY2|~z)z&Oos%%oLr7d(kT56$c zjg8GBv_O-IN*}=`;L@3SRZ1>6fvngtD>kTlmk-DeTU+4F)+S?XQ=JW!0I0!qYU);4 zM3wyse72w||RSvRbOl9a62;B)I&zX)^W%9k)nU91x5~X&8D5Xpy^SF`%WdLE|7A zjddG{)tVpmn&7zG7HwR#dg4hR|6@LpM|>hTi`&X~CmfF-ePHfcs?kQB(-7@+@&lJ= zWvSC>0|4A1qo+PmcRaqVYfhB{H*>g|ARMF{E~BSC(biQ|dps#kL1&=v48&%rp_88& z)-|H;_n{p;7qZ){VS7PKpP007+VS{YB?EbwfjrDWo>vC;4LBZeQjpn*%tmB3Am^){ zc8zNMoxGm(HJsH@ZQ5~t@G*Vt5q)g8KK`gazN_co$&konA;u#i#_o{RqamqX`c%g& zweg*w46Y_x=>#^}grRlX^AOnA**Ww>LQXi=4=>Ydp$!jremwj=qK$=Mcjq%B2a!3x z5R_oB(nx+Nuq#OZ7z#iUh$0BYDbN?Ko<6Qv3Px1~2&)h1a4?fQj<2n*(K2YMC?$aH zm}B_p&nNZa+eD-gbbxV9f)tJ6D0J4l59%A2ii6Odlmw;|4s*dHff$0`D?kF#j?ltb zZ5k*(a||r3FvUYi+KvthcW8xkjXPI5bOmuTL3hgWL}5{8MlW36&9zV|T8^*{7;i5|# ztI>SXP`pnkaX@71G8gb{+!ki2TJv@jYmN0!B;-2sOEopUYeg7#4ezTsF? z`(f-B8J&aGG7iy#bU&!1I~^@c;O}N&D4s|00tzKf$%&`?P`rp@KZ-x2 zcnO3{Xlir`v z6E!ckO2rs=H=b&%;U>34YCr{+Nh-k*R<87d<4jvoEY=1SsydW7SZwI*`X@vpVwoUc&(ZxlyBxM+3F3KT1)ZamdqZ(eR` zlCwo?&{PAaf1pW#GXYLM!nkQr7?%$c7wH}>SC`}Ao2#^eFfOd@^u49nPxl$3K6jLF z8bMAu)@_=c`~@nC5c?&DMq}uo7;=~C+f9Q)vA4r5&stMutu}+}e>?IwCzG!nGd7P2 z_!nYNqv&$nvw3{S_lWhWvsoH>-;};XlLNfEeCjPLtSfS{qcu(C(m(cbb+Vog04Ryc`ppF}C7O5f`>Y_;^P+&Q~jsd!ncs|O}z9Wf@cHFmPX=o8BTzaUsGIPB}k~D%N`hkhsd7vHLx%h$S zbz12Ys24}d;|+nlhJ8e{R4={l*#5Y6{JU5@D0-u%TpqFX(lKnBk5R;+xElr5wzL&P zSjGIl1((1}7(S+PEO;V|Z0o$`iGK?MoY|f-lh-;wd}?|isqf6(mq856FoV7_+p&Uc z$AT9#c>~HWaeQ$3*II0go0}Uf<|b(+maq-Q4Jhm=+EJ`R(Sc$$iZv*1L~#>}wJ2^z zu?fXG6t|+d4aIsC8&GUSk#alAJ5camqjx86g)xzrYqKL|e+4;Daw^ji=> z9(l{0LB&~A=w($oGyV?+zt99QVIlNEUZ~ z{BAvki)?aNaoTS5w+97xfjM>gHZ+?>+X|c9a$6m{x>|^}GqJlx50UQpXf%m;9R6r4 zY3;0Xts?mbj2r@@+!QBWk0;wvJb;4xT#g}ph?IW8KJlv~_~U#cIf_1xo1mob(&BY# zl}nU}5t2Z-d~KG7Doag`lm=%ef4SaT(W19Dw@THPx+dwAW5>shWL;;_Cw{sKryzy4 zu(FcUs=B6@R@->9%{;JT*2IBQ0z|5mcwSR5E2$`OOfvb&x%>$L($6SeE=NO4 zSWk^6zjQu(Dq07(VuTg$Q{Tm(h$~B&1HSr7I zrGw`Ki_{|DLw0h!yufZE5A7TlH1QZkn zUjGtIF~Bf|rzu!(S`4LPNQa>`4EbV6kHQZ{f7C%Njg$$(6yhm@D+F$cwk1}0x^ry! zIcn(DOXTvKRll;Tyt&S5kw-Ou>iF>IG%*kp3PKT#!T`b*Sl+Uv+1A`tRnySALJC8D zIEn}qktm{2$aR55i^bD86!9n$P$Z&A0%47G>}(y0?Hqhzv>y|6Bv}yqGxK2@&OX^D__Gy0#Kv;9B!TI*B$H&ZWm*!?rfA6{BwL2H zrITd#o=AdNOC*V9&4;O$p-!;Ju#r0Qx6lEP0f&AHc{KmEC$ZnOB$Yk-45VT7C2=H| zjrJwkaH#bq>t$opQo+EcMw1YB$(Q(R^Td2{5DU_i7uZ(^i9ZX_ktn@5_%SH%SMO6e4LqBh) zV+tTYkj#LdY;7h|Ew2wj94l9L`xnHIUl?#DA8>8Z(3!~XT_^XG<&jE2xNDJ5`T)BT(*j&v8ACT5Ap=z zRQ7TxsUjwp7)D0=ny_zxyMc5=9!+Ji2AWI1nBAw;2#3{J_U+eL2 zOdzRb4=YI^Wn`gqO9B*2YZ4djs$zdmB%A1>UDvU)BvPm~+ZR`ru#*Q#1i#o>ze6+< zhZSQ;wSDo@z@935U_5*9CTL9mKp{NAwxLH24Jcj)TVp4wj=;NvS?N@_|W? zEqfyq7u#;@(Lh{_cigmoH8=>Ur=@q*s?Xc2ZChA&GF*!CK)OV}XoKITeylo~gpw(Y z{0A<^sgPzg>k^57a-HWI+O!5?wi)Kbs$Y@f|#-@Gn{@o&ovy*|V{h=-vCQ^{Dg_U5Jz1wX z=z-Ov=~WoJkmNRdiqp$6Tm{g*(}E=T;1_Ri=vCH)sLehbJ2wYxT(7!7tjv;5hEf|` zu;*H!ckNAwdqq+|xYJK=fq|g3A5>CN2HZ>AP}c+@vVyJ2AUR|u9^%`Nk(bOS_Phfa8JC_LwtdUit`hn_3IlKhJH_73< zAsoBbeY2e1Q~To0{HAgXd#xW_AM4mM0d72bZxAhvjTifo6c}4>V=4VfYxw$d$$~U+ zCklMYB<*5P_9wOE6bm(xRXIreOun`EY{%b|uS)-cByd_0W|oH(^*$S1Bwb`LnMh0r zLZo3xws>ebhAyIRgzW3Rckvf_KsH|iIGF#Ut`6k6i%pk-9al>G*uyP_{&U>OD*G?_cC?KU{3*e(bTap)aXz!nEpbip-2XL3JWajMR zY5DHsApXdUCl#f0xe|_pV8)r?XbW0wf!8%q4Gl(k*y9gPl~2Na8}j?RK|uK@Lh(+?eRKQRHWQD4Z+L3W z4e-V)w}L-36ibqMWo{$PiAH@_f()^nmmdCIhX!&808crhJkX`!rVYGay{PwM&|ENMwT8#6tp`^_xMyK%$+Nj?Yv_7aSs<2x>=BR#5mwwDy9ZMNS4o4l=2g=6H zwcYxpqxz&SebSFQ_RT1AXY|P)TUXS$Zhgs7eMy(T1R|~(O>A-SGR@<6P`7^YQT^a9 z{opUdBiXa}kj5B&FY~Z&{qUpu;a&RSXLZgc1!Mti!IC)dq`m(OLg6_b=qk)B8 zfrZDT;+$_4l7513rH^01IZdVWr?KRWPjB$h?!aM31BZ154*N1B-1*gb@)gO4tOSM1 z-nI?w4$L|lnAH`S^`#+Vi*Bp#){-wHv-%1-G;)WhQcmG;n^5Ab_%&mKb2yRf2?qRJcT6q>Tu%he2SwxEu4kU zBx{I+IcAcfqxj^^JnVLSW{@_*0C13BAK1y@zebvb7sai3idVELQl6yS;S^@U!wQ+| ztf(Y!>Lc-0oQt;Fs^P7J4&HfVM|H3P3n80)&MuiJy>KBshW&?aSx9oTE}~1$3oc>k zHxQ@OAi>joYDwb*I)Xy*VW$_8VUhnt6y!yoaXQUomarj(E&yb1fqAIzAevvI~pI zR|z+RzAX}FAl(Zgnx)*H(eNfGxsH?xj=&3k+Cy1Ab4V;{oDk@E;;umF?Tg7M;<)Kj zjL=RT#g}4q?Sd@`VCQ~?=ed7Wlc8*z2xZ9=VM-}fo>GE$A;19T@N*)Gqr@&aFN)+Y z8XsH$xwXftCSXEPXDehW7I+p}T@6&P%a_DD=3NS6LLCWnY`zrW7zif=@TB}w44Yj? zk|3?tIx>`mI-jZ|uMtSDu^#AV7`sRXmd_& zB&P&lL;Ij6K{T+DEhG&lVHUQK>BQi;t|OWqY9SBNmuzPC={|UFUA2siaT|^FZ&h#=u^_p2W+lSoYoZqzH;IcmXFm1Jt)2RbuFV4Bd#K zcQLdZL+3Gc0z7x0B?KAs8OYFKbAH&*{t4 z*f>B*qtK9whS3<}1C5$9ySV^b2cZ~@q5wr9iZLj@MF(RsRD_}!#W)lt6nsnLF_?mw zG!S(~u#kgKr?=7iL4|D-S>7gSpt4P*pRTlQQk8Uu9lwLjfhBxB`gIAV=wUm0mb zAs1a~u}asW&WxfOg^0p}q6Wom_S2mtHA`y3lV%hxD3+nXyPRZ2VMEc1VmT8XWLWm= zkVaeBKzR15;%}PZVFKPe!QIWyk(dL{#OF)jv%fekGA)LBzM8ieYC8G6CJawSh%6r3I7?gQDV(!LMhGHAXpc{qMOWRR*4+>=P(!Ch^ z6(L_>-22#;yGU-wQ>fjA;%OAmpm-L=ZWMb^>_yRu;yD!0qj&*?OW)9JmDd``)E^By zPrlf&-yPm{+WnUq(ugpem@=+`)Cttu}Z=De=`6PhyM=ye2YH{NZkh(T4L& z3B^bhqfm@Su@jS{P4$CZLj53V0aORPHI=HZRqW-vNpQInBk$vw`ucj8V9SU9gQOQx zy&uJ&QM`oWWfZS~08jZYEfi{|*`GKc_p2Hw6pRihbo?t2+XKvQ9k=sat?hy1*>!a+3ZRAbiY%o5kv`3zdj%V$6wm&fQ ziSQ@HUCrImV~$3TalXHuJg+Ti>yB*u+YPJEYBWW{H~|hN!g%2vp7Dd7Oc3BbJ{~|B zH31P5g))9p#t!Wumj9lIx(CP%UHVy#HYlAPcz_HI*d}!wQjZ!^y9}vJ_aGS({D9DH zNIPms>*D`d`Ge#ht?z-Dw?;b8KS-MWoj2_wE)u{;!+NNmYKlw4oc3qP2A@I7{qaeP zq_b7G$q~O()iXsh!DM*q>2{y8mBn+N~-HX3##`n)IT@&*826?(ceuWNdn3s47 z`y9n7fXCUbon$2W*7ABC}S(c!rG8p z*`E#ot<7?Ne1JIg>HNjI(hPiIhl#j!)lwn zoBVC@0Yvd^yo@3Kq{43ixJpKf#8fdGt`L~xxTdADd7|WA#>^|glUKg7WhJ`c>p}R) ze-T}L!G^v={z}F=ecmMpsqSsiRqJ>+`t*HxcGOp5!w_i*(;g#B$O+bT45s>c4)~ph z&+H^%SI}8o8YD}L+}fwuGsj3&+_g|)E-e@=%V8NcEDtlKxSN=rI0lq%0Heo=DUlDc ze3%~s$-4wAj8@h`F1i?&YF)!hkCQY*4OV?EiX|xOnEg1}PL?~PJ|b0G1FvZ)ocwqQ z?Qru6VjPFJ=WQyaRRXfkIAJ|)IX?(?!b_UaEc?9iw1nDMSt`x&wO==6{u$lJPqlC5Li(2?XE6@Q+#zK1!GADwT1Pj>qG za%Z1tSldN1Adb7jxmZ{i*<4atT2(7&{8J4(aFIj+?fKh9k`ci>Al7*gOp@+r!I$8* z=K(h45Mruvm`KQ>b_i=%lK=*PN_O#j#_A(WUNuBdIDzpSn#R>M)nsmV7 zwRhTpAzY`@YmxK|P^qU^W6kn$ieG+6Fv>iXDr7kSVGs<&z!#hG zX6Gzy1p6&iNSNhmt-0uKzmic;HG4A32f`A7^sIKeyjP=(F{~*}XdtVde-9JR=nbBg z-I!_VFmuHV$>lsVaz3zCj-crS6ub?1E_qgzw4c<;7hc~fPD&8QkOR(J5`-Ug$(N@H zZn0iOSddp|mfmIC(*)ySe)0cqvV9f)6ZY>kVHi2-OiLH^xy9;xY(JHBP`HQ7O{t( z!+6-Eq;J4wx#Ylze6f^P?!pJLqo7Zhp9N|HzX!2Md*QsNYgDjThYFdxchJZ8?CYVz zG$A$EIdYf~><916Cyy1T4B>6fv+RRkj7TqHqWp2dQ~x}j5=>0)3ha%sLJ?W#^eqxj z>f?Bq;jQSoO}&j?-eROw$mX5zG$;P<&R0+>#7rCsc{+_l<>??S9PkHTIUZ;D_~D|~ zl~$JzJZiX?8uA{4xs_Mno?!Qs3h^d>ukXS_cBAO4d-4U_~YSW zBCMyv{qH)zC=)u9ov$_uCZR+8-X7a)?Z^)Wkec|TO7|1$Y-2&Vf8x*TB?jUzHG(R} zuUCz6-W>xg3V-4S7W|cv9r~ZkJY8f5z8B)s#Y>l^_?sHHnTRl80+Lni6~e=bmsLP( zVg^;afoMKnQ0?zQbzV?k_U+d~Rus6=d)fRHo4=PWfSsrp;@O{m5W>O(y{tiOmn3A9 zV0kjn9KCBvdcK`QU}0zSSM6qgHa?w$q1+}iE2u$UsKE+qh=LmGg&L-yhAXHMUZ{}@3RbVVJIrV=RDpsjR8V8QP-7KT zk%B7rLXA^UB?@Z17ixlnny8>Ad7&mNs8R(r#S1l6K}}OoWnQRq1vOnkRd}IhD5#kV zYL*wOQbAp#pk{la<|wGS3JR9jU3QJkS5ONS)U{ryg$insf~xXDU8kVT3Tm+zs#-yb z3d-Vzs^KVMGhx@w5#n&5S>lDMQ!s{iWK6vmra?)qQAw@I3)QT6Dtb`%w9E@5DW0r~ zCz}_lRq=FP5as#4-V3ur@w8I$)aHe{K|$FSRJ#{ym6G1chvb~E_QI@DFnM3dn47#X zYn9Y)R#Lmg3$;!`-O5oB&@Q)Kj$%c(3-I2dRR|#)6ki)JM~Sy9zV1*^o4in)71S05 zb*C4~p`f-ZsJpyScPl8SptgCTwkxQ66x6+5sQVPupA^*nUZ@=k>H!7ypcm>P1@*9k zdc+I$sDgS-K|St;dO|@xsi1axp`KDuyA;&ZUZ`gj)UzCgtqu1-^$5a>9v2eu{<+ud z{@HmsR(y_C_|tf{>3(4Z`}-}zc0I&<9#L$wAD!iVZ=G;Y5D8${-7SpNz794XER;2c zBuO?L1OpGcHHG9k*>DI9Jm}UGvgD6o8FZCHNO!)^griZ;cb*ZlIUP;vgO2vr{V_V) zJH{WQqZp$vI@;UQmFQ?M92E9)I;yr_k&ddNSEQqAsG5$dZE8BIwp~s~)mCphsttjD zk;Lwogppp9R2{cZN~(tRNlDcZ4<%LWJ(N_fS5Z>z$H*;{*ch7ttK*c5iYq;R4Pd=U zX_AMOs(oLLl&aA(DP`wh5R$`@jH<0N8D%q_U_}zDw#p=w?L;e*PqkGhpG?>XRwSKj zt4uoCe6%9jR9j`T$@ZcZNv7H=lT7CSB3O}Js;x4)WQ);?q*85_NhRBlRwR>Zt4t@Zr9 z9ICA{Ib^9XffY%i+A5Pm)`nIjgKDcx2H6p`A_-JmWfI7;UIr_YKebgRf9yuIBI#3G zWzxr7XhpK8w#sCW4SEHvNb=NHndGtcXhm|Tw#wv=eTG(~b!w|j>)5DQ!HR@VZIuZf zyA!QQcL1zN-PBf@y0LrEilj|#l}Q`>4y{Pm)K;0Su_>>C z6-k=fDw8z!5L%I*sjV_SW9QL|#7u3Ki5Z*uI#`jAsjV_0V>`hbibPCh^$;H=!YrBdM%c;Ycd1iX)wV)5DRHWRCRrw`7jgbyVR<5FCYb5&^6y zl*ICm!BkFCFGY{(In<5R%3MjM9OZoAxDXKS{NRk6B&B(7Tw>3hU2okVBT2nu{4tV* zG5R7&y**utB=y4eL6X$gE0QEN^ok@&4ONpQwM|Wu)V9k>lG^G`l3>nm_YZx{*{S38 zNs!c#J_(W<;vqVuq_(=r(Gj#FIZ|8QL@Sabwbf0I zzCtUKBehi~M|X{*GB;8iWo~rWJzzv?q&CXb=&tX;h|EZBbTgwV1gyx6)K)h$dI+t^ zjMP>)Gdhn}WJYSMn;Fd%z>3UBZFMuFooGd7q_(=5kwC$U%t&o@Go$%vMP{V7x|z{l zuySUkvMS6dhz0o2B+iUffF5SF*aw^=Gg4dK%xFJakr}D2Ze|pw1uHTmwbjjxmZBAz zk=p8JMz5n4nUUJ+W=3&3up%>3TiwhE_J@#}k=p8JMu*Xg%t&o@Gow^rup%>3Tiwj4 z4XwzG)K)h$I)YYYMrx~@8D;6AA&?oVtXE-1Dy)hbozv6p$a$W3ZuFzD4*<**Mbe&0 zqoQbvlScdiGeu~SCym@1et?-GG{{8At-*slX%rGe(-}+}#W{bArEBRhTw?-@gm6E{ z<|ff9)#0mDo&2s`x`?l4P59C{uw>S($|~Fnr6j*fI>cT~q&JcjR**z%$Zzc7B>HX2 zZ}2|#d_eQ!rLgCsu?cqcap|Uyn=+xKQeL*K(K`1g)8#~mZ>igI*_L!_VsldHE3lox zv{V|)2By*k_|g}`L&GqtGk3nOJ+;^(`B$k{Nw(#lH{{| z6D>5P*TEM}lFcaA)zq{?h}~?WlLcE7d*4L!AvuyofAz!H*DjJJEoL;E4%VB5^FDdG zc|+OgY?|*2%d^(PcgCzCo9-b_`1+8Jju`;!)@*h7QDF;g;8ksIGA^;y)h@BImI3s& zq7WKC%*21FLU6p<-bP1zd%Xk1rL1TF04_9ZBtN8KnFLm)njQy$4eOItD5;s zeZF`YxAc?Zpqlwg3YV@`UPZyzDBOUaabbeX2e(y{mx%GD6}X_n6@Z&x;$Et-d#ZbN z1v=x)_4uj9OTzhy3b@%Bk-`Xl|H$bT^Vk?ihaVpocV)y|@;<#ixm z&!YZ*y6!S_SIT1Ng=^`CP*o}nQt8ijH`9;*$5Z;r*ZJ~|bh}?4>8-wf5O;l2xkd;c2HvCF{p?!1-tLD>GCZF7cS?HqRxdsx zZHKCAA{%4uRHW?FhiRlE{=*RXD28h{+5K3^BQ){yG^Rg71B$^lWWVPbnD0-cPO9A> zlGjw>5`X@)LvgFZl`}tf?;|h>Q?~3;n#le>UGRq$8-JQgnUOK42#K&pb|G^T4^Ag!` z1a|-I{%q`%w2>sT)t|#A6t6!?7ZV@HgBKzM$ld6jbR_T1*v;A0ICy7q@hLhwF+dEr zLkGpK!FzNPm&v$6OM4(Yx0B{*ld#Ir533IpBH5*<=&0!5iN?O(Ne63<-Vq`%kC1}t zusyqIJPn4$Di{a~X_^4nV<7wwa;tjVQh(3Z5n>NiwQ3q43xZU-bgm0l=`y%3M5W6_ zo!C#M>(6n$N1F`KMP=8pmnwyT$V`t4E3@5TH!MwMUDJg~E#LWwP2Ww^lC$iAVm7bh z0emB@Ak~?H1wBKrp`l_9+tCRZ#*NR=bl+V0y(MnVl@<~z4)mlW=CNgaX(=f4dnpI? zC^UdEDeG$M){s(C!eKpz@jcZOY>Uos0q>@$SixI*kNKm zdwDM{U?I;!XXuSB=D8OjlTa-KVNam;-Itb%oVJRbI~(A=hILR49*P zf&h<%fx)cbb3mq#JV!Gl<$F_M671X>62bOY3PBlVp6=z9(QqvN_C&|9x?fYnkufcsi;HA~4({6NHbb$I_9h&l*JADa5?n5_+>r`zTC(cEsU9*Dp>XP%YKKU?w;YAw24!=lav?^Vi z`ZR!z+fU;n=exU?Af@o`wV=n%?{%;F2vCaq|GSsFDth`j`@gY!soGQ!=du|uyLlH| z`Wkd5yr9WKw&i6SOBS)GUZ#UNDeKD^|H>7N8(yK=9aTz?3*$X5Rn?*7F|10Z`ycIO zMz}~lgP5l$|GzbcC0=n1``_69aW$g4arszw+0E#`JzzGaG4nwR3sjzg{X2V}3grz` za}UxY_WD7%?c9i={)cD;(3-+SG@1Q$kPc#ZfPxhqqW;=6yXJw>Y|kM&fp3AqFLY@L z_DuF*G~DK5{{n2OkZn0g^V!_L(1`poJ;S1ONY#HAdTxbct=#*iPVOh=O&S@r#I1S= z(i_Y7E404}x4a_UX)cv*7KBjN^%ojQWZQ|$Z3%BFwxYK#i?`}68b#`4LI<17`mui; zq<$gANBRmktzXPNm10#C^o8ExiENtHpnt9b~^CTPnJIONhn>2>+QCoG8hLC3V>tUKq zT3GBmKwMkca5zlh8$!Snr$Wad-APKNR9rO|V-dly@yHVP>^n4)E~^*Vt8dZKu$u&D z&#)1tzcvD@UhTtjRGNb@Pb+^K_2d~mPV*-RcA57tf?O4t1mLoLZ&x-Z890WUd zgwF7l-oo>@*{ly}Gug}zen9KV8*JcFx}>9A;Ks=d%0(s>l0Vovzf zsZO+2!FzySo2~PCSv~{1A5nY~mru>|30c1N#5|Y*mS<}DM{H@BRou|X<+HK_cIuXQ zPq+mOnun$h7~zg1G2fw|SKT`FiuEUN=l>jr_}>4!Fmy9ocVdWt0qKu5#y<(Qqu~G= zR%7TmhEAaP90mWJ_0Jf(39T0}l)4h@=4(8^8${3ER(NH=*TH-Pdba388gJm+iSrrd zp;!w8*li!u>~v3G@y3Lys6$m^JPuNyV1{R;H|`4n70(e%H59n2W-L6~!AUPNE1D z27tsrjRHFxR@b$x%(b?_AIO8PDDttO{Li4^-(iqD!F^2s6HWAg4k`#DYph$H1An4L qG}1`Y&T77ub#_20:# - msec=int(float('0.'+strTime[20:])*1000000) # microsecond - str2=strTime[0:19]+' '+str(msec) - return datetime.strptime(str2,'%Y %m %d %H %M %S %f') +# def str2time(strTime): +# if len(strTime)>20:# +# msec=int(float('0.'+strTime[20:])*1000000) # microsecond +# str2=strTime[0:19]+' '+str(msec) +# return datetime.strptime(str2,'%Y %m %d %H %M %S %f') #datetime to mjd def time2mjd(dateT): @@ -515,14 +515,14 @@ def deg2HMS(ra0, dec0): return hhmmss+dms_d+dms_m+dms_s ################################################################################ -def cut_radius(rcutstar, mstar, mag): - return rcutstar * 10**(0.4*2*(mstar-mag)/3.125) +# def cut_radius(rcutstar, mstar, mag): +# return rcutstar * 10**(0.4*2*(mstar-mag)/3.125) -def core_radius(rcorestar, mstar, mag): - return rcorestar * 10**(0.4*(mstar-mag)/2) +# def core_radius(rcorestar, mstar, mag): +# return rcorestar * 10**(0.4*(mstar-mag)/2) -def v_disp(sigstar, mstar, mag): - return sigstar * 10**(0.4*(mstar-mag)/3.57) +# def v_disp(sigstar, mstar, mag): +# return sigstar * 10**(0.4*(mstar-mag)/3.57) ############################################################################### ##################################################################################### @@ -561,56 +561,56 @@ def distortField(ra, dec, ch): return distortRa, distortDec ##################################################################################### -def world_to_pixel(sra, - sdec, - rotsky, - tra, - tdec, - x_center_pixel, - y_center_pixel, - pixelsize=0.05): - """_summary_ - - Parameters - ---------- - sra : np.array - star ra such as np.array([22,23,24]),unit:deg - sdec : np.array - stat dec such as np.array([10,20,30]),unit:deg - rotsky : float - rotation angel of the telescope,unit:deg - tra : float - telescope pointing ra,such as 222.2 deg - rdec : float - telescope pointing dec,such as 33.3 deg - x_center_pixel : float - x refer point of MCI,usually the center of the CCD,which start from (0,0) - y_center_pixel : float - y refer point of MCI,usually the center of the CCD,which start from(0,0) - pixelsize : float - pixelsize for MCI ccd, default :0.05 arcsec / pixel - - Returns - ------- - pixel_position:list with two np.array - such as [array([124.16937605, 99. , 99. ]), - array([ 97.52378661, 99. , 100.50014483])] - """ - theta_r = rotsky / 180 * np.pi - w = WCS(naxis=2) - w.wcs.crpix = [x_center_pixel + 1, y_center_pixel + 1] # pixel from (1, 1) - w.wcs.cd = np.array([[ - -np.cos(-theta_r) * pixelsize / 3600, - np.sin(-theta_r) * pixelsize / 3600 - ], - [ - np.sin(-theta_r) * pixelsize / 3600, - np.cos(-theta_r) * pixelsize / 3600 - ]]) - w.wcs.crval = [tra, tdec] - w.wcs.ctype = ["RA---TAN", "DEC--TAN"] - pixel_position = w.all_world2pix(sra, sdec, 0) - return pixel_position +# def world_to_pixel(sra, +# sdec, +# rotsky, +# tra, +# tdec, +# x_center_pixel, +# y_center_pixel, +# pixelsize=0.05): +# """_summary_ + +# Parameters +# ---------- +# sra : np.array +# star ra such as np.array([22,23,24]),unit:deg +# sdec : np.array +# stat dec such as np.array([10,20,30]),unit:deg +# rotsky : float +# rotation angel of the telescope,unit:deg +# tra : float +# telescope pointing ra,such as 222.2 deg +# rdec : float +# telescope pointing dec,such as 33.3 deg +# x_center_pixel : float +# x refer point of MCI,usually the center of the CCD,which start from (0,0) +# y_center_pixel : float +# y refer point of MCI,usually the center of the CCD,which start from(0,0) +# pixelsize : float +# pixelsize for MCI ccd, default :0.05 arcsec / pixel + +# Returns +# ------- +# pixel_position:list with two np.array +# such as [array([124.16937605, 99. , 99. ]), +# array([ 97.52378661, 99. , 100.50014483])] +# """ +# theta_r = rotsky / 180 * np.pi +# w = WCS(naxis=2) +# w.wcs.crpix = [x_center_pixel + 1, y_center_pixel + 1] # pixel from (1, 1) +# w.wcs.cd = np.array([[ +# -np.cos(-theta_r) * pixelsize / 3600, +# np.sin(-theta_r) * pixelsize / 3600 +# ], +# [ +# np.sin(-theta_r) * pixelsize / 3600, +# np.cos(-theta_r) * pixelsize / 3600 +# ]]) +# w.wcs.crval = [tra, tdec] +# w.wcs.ctype = ["RA---TAN", "DEC--TAN"] +# pixel_position = w.all_world2pix(sra, sdec, 0) +# return pixel_position ############################################################################### @@ -635,39 +635,39 @@ def cal_pos(center_ra, center_dec, rotTel, rotSky, sRa, sDec): ##################################################################################### -def krebin(a, sample): - """ Fast Rebinning with flux conservation - New shape must be an integer divisor of the current shape. - This algorithm is much faster than rebin_array - Parameters - ---------- - a : array_like - input array - sample : rebinned sample 2 or 3 or 4 or other int value - """ - # Klaus P's fastrebin from web - size=a.shape - shape=np.zeros(2,dtype=int) - shape[0]=int(size[0]/sample) - shape[1]=int(size[1]/sample) - - sh = shape[0], a.shape[0] // shape[0], shape[1], a.shape[1] // shape[1] - - return a.reshape(sh).sum(-1).sum(1) +# def krebin(a, sample): +# """ Fast Rebinning with flux conservation +# New shape must be an integer divisor of the current shape. +# This algorithm is much faster than rebin_array +# Parameters +# ---------- +# a : array_like +# input array +# sample : rebinned sample 2 or 3 or 4 or other int value +# """ +# # Klaus P's fastrebin from web +# size=a.shape +# shape=np.zeros(2,dtype=int) +# shape[0]=int(size[0]/sample) +# shape[1]=int(size[1]/sample) + +# sh = shape[0], a.shape[0] // shape[0], shape[1], a.shape[1] // shape[1] + +# return a.reshape(sh).sum(-1).sum(1) ##################################################################### -def float2char(a, z): - ### a is input float value - ### transfer float a to chars and save z bis - b=str(a) - n=len(b) +# def float2char(a, z): +# ### a is input float value +# ### transfer float a to chars and save z bis +# b=str(a) +# n=len(b) - if z>=n: - return b - else: - return b[:z] +# if z>=n: +# return b +# else: +# return b[:z] ######################################################### def centroid(data): @@ -744,270 +744,270 @@ def dft2(ary, Q, samples, shift=None): return out ############################################################################## -def idft2(ary, Q, samples, shift=None): - """Compute the two dimensional inverse Discrete Fourier Transform of a matrix. +# def idft2(ary, Q, samples, shift=None): +# """Compute the two dimensional inverse Discrete Fourier Transform of a matrix. - Parameters - ---------- - ary : `numpy.ndarray` - an array, 2D, real or complex. Not fftshifted. - Q : `float` - oversampling / padding factor to mimic an FFT. If Q=2, Nyquist sampled - samples : `int` or `Iterable` - number of samples in the output plane. - If an int, used for both dimensions. If an iterable, used for each dim - shift : `float`, optional - shift of the output domain, as a frequency. Same broadcast - rules apply as with samples. +# Parameters +# ---------- +# ary : `numpy.ndarray` +# an array, 2D, real or complex. Not fftshifted. +# Q : `float` +# oversampling / padding factor to mimic an FFT. If Q=2, Nyquist sampled +# samples : `int` or `Iterable` +# number of samples in the output plane. +# If an int, used for both dimensions. If an iterable, used for each dim +# shift : `float`, optional +# shift of the output domain, as a frequency. Same broadcast +# rules apply as with samples. - Returns - ------- - `numpy.ndarray` - 2D array containing the shifted transform. - Equivalent to ifftshift(ifft2(fftshift(ary))) modulo output - sampling/grid differences +# Returns +# ------- +# `numpy.ndarray` +# 2D array containing the shifted transform. +# Equivalent to ifftshift(ifft2(fftshift(ary))) modulo output +# sampling/grid differences - """ - # this is for dtype stabilization - Q = float(Q) ## Q=lambda*f/D/pixelsize +# """ +# # this is for dtype stabilization +# Q = float(Q) ## Q=lambda*f/D/pixelsize - n, m = ary.shape ###ary maybe is the pupil function - N, M = samples +# n, m = ary.shape ###ary maybe is the pupil function +# N, M = samples - X, Y, U, V = (fftrange(n) for n in (m, n, M, N)) +# X, Y, U, V = (fftrange(n) for n in (m, n, M, N)) - ############################################################################### - nm = n*m - NM = N*M - r = NM/nm - a = 1 / Q +# ############################################################################### +# nm = n*m +# NM = N*M +# r = NM/nm +# a = 1 / Q - ################################################################### - # Eout_fwd = np.exp(-1j * 2 * np.pi * a / n * np.outer(Y, V).T) - # Ein_fwd = np.exp(-1j * 2 * np.pi * a / m * np.outer(X, U)) - Eout_rev = np.exp(1j * 2 * np.pi * a / n * np.outer(Y, V).T) * (1/r) - Ein_rev = np.exp(1j * 2 * np.pi * a / m * np.outer(X, U)) * (1/nm) +# ################################################################### +# # Eout_fwd = np.exp(-1j * 2 * np.pi * a / n * np.outer(Y, V).T) +# # Ein_fwd = np.exp(-1j * 2 * np.pi * a / m * np.outer(X, U)) +# Eout_rev = np.exp(1j * 2 * np.pi * a / n * np.outer(Y, V).T) * (1/r) +# Ein_rev = np.exp(1j * 2 * np.pi * a / m * np.outer(X, U)) * (1/nm) - ########################################################################### +# ########################################################################### - out = Eout_rev @ ary @ Ein_rev +# out = Eout_rev @ ary @ Ein_rev - return out +# return out ############################################################################### ############################################################################### -def opd2psf(opd,cwave,oversampling): - ''' - calculate psf from opd data with defined wavelength +# def opd2psf(opd,cwave,oversampling): +# ''' +# calculate psf from opd data with defined wavelength - Parameters - ---------- - opd : TYPE - the input opd data. - cwave : TYPE - the defined wavelength. +# Parameters +# ---------- +# opd : TYPE +# the input opd data. +# cwave : TYPE +# the defined wavelength. - ''' +# ''' - D=2; - focLength=41.253 # % MCI focal length in meters - pixel=10/oversampling # % pixel size in microns +# D=2; +# focLength=41.253 # % MCI focal length in meters +# pixel=10/oversampling # % pixel size in microns - opd=opd*1e9; ## convert opd unit of meter to nm +# opd=opd*1e9; ## convert opd unit of meter to nm - zoomopd=ndimage.zoom(opd, 2, order=1) +# zoomopd=ndimage.zoom(opd, 2, order=1) - m,n=np.shape(zoomopd); - clambda=cwave; +# m,n=np.shape(zoomopd); +# clambda=cwave; - wfe=zoomopd/clambda; #%% the main telescope aberration ,in wavelength; +# wfe=zoomopd/clambda; #%% the main telescope aberration ,in wavelength; - Q=clambda*1e-3*focLength/D/pixel +# Q=clambda*1e-3*focLength/D/pixel - pupil=np.zeros((m,m)); +# pupil=np.zeros((m,m)); - pupil[abs(zoomopd)>0]=1; - #idx=pupil>0; +# pupil[abs(zoomopd)>0]=1; +# #idx=pupil>0; - Nt=512 +# Nt=512 - phase=2*np.pi*wfe; #% wavefront phase in radians +# phase=2*np.pi*wfe; #% wavefront phase in radians - #%generalized pupil function of channel1; - pk=pupil*np.exp(cmath.sqrt(-1)*(phase)) +# #%generalized pupil function of channel1; +# pk=pupil*np.exp(cmath.sqrt(-1)*(phase)) - pf=dft2(pk, Q, Nt) +# pf=dft2(pk, Q, Nt) - pf2 = abs(pf)**2 - psfout=pf2/pf2.sum() +# pf2 = abs(pf)**2 +# psfout=pf2/pf2.sum() - cx,cy=centroid(psfout) +# cx,cy=centroid(psfout) - psf=ndimage.shift(psfout,[Nt/2-1-cy, Nt/2-1-cx],order=1, mode='nearest' ) +# psf=ndimage.shift(psfout,[Nt/2-1-cy, Nt/2-1-cx],order=1, mode='nearest' ) - return psf +# return psf ######################################################################## -def cal_PSF_new(channel,wfetimes,oversampling): - # filed postion: Ra=ys1, Dec=ys2 +# def cal_PSF_new(channel,wfetimes,oversampling): +# # filed postion: Ra=ys1, Dec=ys2 - ############ use opd cal PSF on the defined field; - #% psf interpolation test 20210207 - # xfmin=-0.07 # % the filed minmum value in the x direction in degrees - # xfmax= 0.07 #% the filed maxmum value in the x direction in degrees - # yfmin=-0.35 #% the filed minmum value in the y direction in degrees - # yfmax=-0.49 # % the filed maxmum value in the y direction in degrees +# ############ use opd cal PSF on the defined field; +# #% psf interpolation test 20210207 +# # xfmin=-0.07 # % the filed minmum value in the x direction in degrees +# # xfmax= 0.07 #% the filed maxmum value in the x direction in degrees +# # yfmin=-0.35 #% the filed minmum value in the y direction in degrees +# # yfmax=-0.49 # % the filed maxmum value in the y direction in degrees - cfx= 0.2*3600; # % field center in x derection in arcsec; - cfy= -0.45*3600; # % field center in y derection in arcsec; +# cfx= 0.2*3600; # % field center in x derection in arcsec; +# cfy= -0.45*3600; # % field center in y derection in arcsec; - #wavelength=[255,337,419,501,583,665,747,829,911,1000] +# #wavelength=[255,337,419,501,583,665,747,829,911,1000] - cwave=dict() +# cwave=dict() - if channel == 'g': - cwave[channel]=475 - waven=4 - if channel == 'r': - cwave[channel]=625 - waven=6 - if channel == 'i': - cwave[channel]=776 - waven=7 +# if channel == 'g': +# cwave[channel]=475 +# waven=4 +# if channel == 'r': +# cwave[channel]=625 +# waven=6 +# if channel == 'i': +# cwave[channel]=776 +# waven=7 - n=np.arange(1,101,1) +# n=np.arange(1,101,1) - ffx= 0.1+ 0.02222222*((n-1)%10) - ffy=-0.35-0.02222222*np.floor((n-0.1)/10) +# ffx= 0.1+ 0.02222222*((n-1)%10) +# ffy=-0.35-0.02222222*np.floor((n-0.1)/10) - psf=dict() - fx=dict() - fy=dict() - fx[channel]=ffx*3600-cfx +# psf=dict() +# fx=dict() +# fy=dict() +# fx[channel]=ffx*3600-cfx - fy[channel]=ffy*3600-cfy +# fy[channel]=ffy*3600-cfy - Nt=512 +# Nt=512 - psf[channel]=np.zeros((len(n),Nt,Nt)) +# psf[channel]=np.zeros((len(n),Nt,Nt)) - for i in range(len(n)): - #waven - # waven=4 - # fieldn=10 - file=self.information['dir_path']+'MCI_inputData/MCI_wavefront/wave_'+str(waven+1)+'/wavefront/opd_'+str(i+1)+'.mat' - data=sio.loadmat(file) - opd=data['opd'] ## opd data; - psf[channel][i,:,:]=opd2psf(wfetimes*opd, cwave[channel],oversampling) +# for i in range(len(n)): +# #waven +# # waven=4 +# # fieldn=10 +# file=self.information['dir_path']+'MCI_inputData/MCI_wavefront/wave_'+str(waven+1)+'/wavefront/opd_'+str(i+1)+'.mat' +# data=sio.loadmat(file) +# opd=data['opd'] ## opd data; +# psf[channel][i,:,:]=opd2psf(wfetimes*opd, cwave[channel],oversampling) - hdu1=fits.PrimaryHDU(psf[channel]) +# hdu1=fits.PrimaryHDU(psf[channel]) - hdu1.header.append(('channel', channel, 'MCI PSF channel')) +# hdu1.header.append(('channel', channel, 'MCI PSF channel')) - hdu1.header['pixsize']=(0.05/oversampling, 'PSF pixel size in arcsec') +# hdu1.header['pixsize']=(0.05/oversampling, 'PSF pixel size in arcsec') - hdu1.header['wfetimes']=(wfetimes, 'wavefront error magnification to CSST primary wavefront ') +# hdu1.header['wfetimes']=(wfetimes, 'wavefront error magnification to CSST primary wavefront ') - dtime=datetime.datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S') - hdu1.header.add_history('PSF is generated on :'+dtime) +# dtime=datetime.datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S') +# hdu1.header.add_history('PSF is generated on :'+dtime) - hdu2=fits.ImageHDU(fx[channel]) +# hdu2=fits.ImageHDU(fx[channel]) - hdu2.header.append(('field_X', 'arcsec')) +# hdu2.header.append(('field_X', 'arcsec')) - hdu3=fits.ImageHDU(fy[channel]) +# hdu3=fits.ImageHDU(fy[channel]) - hdu3.header.append(('field_Y', 'arcsec')) +# hdu3.header.append(('field_Y', 'arcsec')) - newd=fits.HDUList([hdu1,hdu2,hdu3]) +# newd=fits.HDUList([hdu1,hdu2,hdu3]) - PSFfilename=self.information['dir_path']+'MCI_inputData/PSF/PSF_'+channel+'.fits' +# PSFfilename=self.information['dir_path']+'MCI_inputData/PSF/PSF_'+channel+'.fits' - newd.writeto(PSFfilename,overwrite=True) +# newd.writeto(PSFfilename,overwrite=True) - return +# return ############################################################################# -def cal_Filter_PSF(wfetimes): - # filed postion: Ra=ys1, Dec=ys2 +# def cal_Filter_PSF(wfetimes): +# # filed postion: Ra=ys1, Dec=ys2 - ############ use opd cal PSF on the defined field; - #% psf interpolation test 20210207 - # xfmin=-0.07 # % the filed minmum value in the x direction in degrees - # xfmax= 0.07 #% the filed maxmum value in the x direction in degrees - # yfmin=-0.35 #% the filed minmum value in the y direction in degrees - # yfmax=-0.49 # % the filed maxmum value in the y direction in degrees - oversampling=2 +# ############ use opd cal PSF on the defined field; +# #% psf interpolation test 20210207 +# # xfmin=-0.07 # % the filed minmum value in the x direction in degrees +# # xfmax= 0.07 #% the filed maxmum value in the x direction in degrees +# # yfmin=-0.35 #% the filed minmum value in the y direction in degrees +# # yfmax=-0.49 # % the filed maxmum value in the y direction in degrees +# oversampling=2 - cfx= 0.2*3600; # % field center in x derection in arcsec; - cfy= -0.45*3600; # % field center in y derection in arcsec; +# cfx= 0.2*3600; # % field center in x derection in arcsec; +# cfy= -0.45*3600; # % field center in y derection in arcsec; - wavelist =np.array([255, 337,419,501,583,665,747,829,911,1000]) +# wavelist =np.array([255, 337,419,501,583,665,747,829,911,1000]) - filterP=np.load(self.information['dir_path']+'MCI_inputData/MCI_filters/mci_filterPWTC.npy',allow_pickle=True).item() +# filterP=np.load(self.information['dir_path']+'MCI_inputData/MCI_filters/mci_filterPWTC.npy',allow_pickle=True).item() - fn=np.arange(1,101,1) ### PSF field point - ffx= 0.1+ 0.02222222*((fn-1)%10) - ffy=-0.35-0.02222222*np.floor((fn-0.1)/10) +# fn=np.arange(1,101,1) ### PSF field point +# ffx= 0.1+ 0.02222222*((fn-1)%10) +# ffy=-0.35-0.02222222*np.floor((fn-0.1)/10) - psf=dict() - fx=dict() - fy=dict() - fx=ffx*3600-cfx - fy=ffy*3600-cfy - Nt=512 +# psf=dict() +# fx=dict() +# fy=dict() +# fx=ffx*3600-cfx +# fy=ffy*3600-cfy +# Nt=512 - for filterk in filterP.keys(): - filtername=filterk +# for filterk in filterP.keys(): +# filtername=filterk - ## print(filtername) +# ## print(filtername) - filterwaveC=filterP[filterk]['Clambda'] - filterFwhm=filterP[filterk]['FWHM'] +# filterwaveC=filterP[filterk]['Clambda'] +# filterFwhm=filterP[filterk]['FWHM'] - fdis=abs(filterwaveC-wavelist) +# fdis=abs(filterwaveC-wavelist) - findex=np.argsort(fdis) +# findex=np.argsort(fdis) - waven=findex[0] +# waven=findex[0] - psf[filtername]=dict() +# psf[filtername]=dict() - psf[filtername]['psf_field_X'] =fx - psf[filtername]['psf_field_Y'] =fy - psf[filtername]['psf_scale'] =0.05/oversampling - psf[filtername]['wfetimes'] =wfetimes - psf[filtername]['psf_mat'] =np.zeros( (Nt,Nt,7,len(fn)) ) - psf[filtername]['filter_name'] =filtername - psf[filtername]['filter_channel']=filterP[filterk]['channel'] - psf[filtername]['psf_iwave'] =np.zeros(7) +# psf[filtername]['psf_field_X'] =fx +# psf[filtername]['psf_field_Y'] =fy +# psf[filtername]['psf_scale'] =0.05/oversampling +# psf[filtername]['wfetimes'] =wfetimes +# psf[filtername]['psf_mat'] =np.zeros( (Nt,Nt,7,len(fn)) ) +# psf[filtername]['filter_name'] =filtername +# psf[filtername]['filter_channel']=filterP[filterk]['channel'] +# psf[filtername]['psf_iwave'] =np.zeros(7) - for ii in range(len(fn)): - #waven - # waven=4 - # fieldn=10 - file=self.information['dir_path']+'MCI_input/MCI_wavefront/wave_'+str(waven+1)+'/wavefront/opd_'+str(ii+1)+'.mat' - data=sio.loadmat(file) - opd=data['opd'] ## opd data; - for kk in range(7): - ### - wavek=filterwaveC+filterFwhm*(kk-3)*0.25 +# for ii in range(len(fn)): +# #waven +# # waven=4 +# # fieldn=10 +# file=self.information['dir_path']+'MCI_input/MCI_wavefront/wave_'+str(waven+1)+'/wavefront/opd_'+str(ii+1)+'.mat' +# data=sio.loadmat(file) +# opd=data['opd'] ## opd data; +# for kk in range(7): +# ### +# wavek=filterwaveC+filterFwhm*(kk-3)*0.25 - psfd=opd2psf(wfetimes*opd, wavek, oversampling) +# psfd=opd2psf(wfetimes*opd, wavek, oversampling) - psf[filtername]['psf_mat'][:,:,kk,ii]=psfd[:,:] +# psf[filtername]['psf_mat'][:,:,kk,ii]=psfd[:,:] - if ii==0: - psf[filtername]['psf_iwave'][kk]=wavek +# if ii==0: +# psf[filtername]['psf_iwave'][kk]=wavek - np.save(self.information['dir_path']+'mci_sim_result/'+filtername+'_PSF.npy', psf[filtername]) +# np.save(self.information['dir_path']+'mci_sim_result/'+filtername+'_PSF.npy', psf[filtername]) - return +# return @@ -1169,13 +1169,9 @@ class MCIsimulator(): # parallelTrapfile =self.information['dir_path']+'MCI_inputData/data/cdm_euclid_parallel.dat', # serialTrapfile =self.information['dir_path']+'MCI_inputData/data/cdm_euclid_serial.dat', mode='same')) - - - - #################################################### - - - + +#################################################### + ############################################################################### def readConfigs(self,simnumber,source): """ @@ -1192,7 +1188,7 @@ class MCIsimulator(): -################################################################################################3 +########################################################################### def processConfigs(self): """ @@ -1287,16 +1283,16 @@ class MCIsimulator(): return ######################################################################################### - def make_c_coor(self, bs, nc): - ''' - Draw the mesh grids for a bs*bs box with nc*nc pixels - ''' - - ds=bs/nc - xx01 = np.linspace(-bs/2.0,bs/2.0-ds,nc)+0.5*ds - xx02 = np.linspace(-bs/2.0,bs/2.0-ds,nc)+0.5*ds - xg2,xg1 = np.meshgrid(xx01,xx02) - return xg1,xg2 + # def make_c_coor(self, bs, nc): + # ''' + # Draw the mesh grids for a bs*bs box with nc*nc pixels + # ''' + + # ds=bs/nc + # xx01 = np.linspace(-bs/2.0,bs/2.0-ds,nc)+0.5*ds + # xx02 = np.linspace(-bs/2.0,bs/2.0-ds,nc)+0.5*ds + # xg2,xg1 = np.meshgrid(xx01,xx02) + # return xg1,xg2 ########################################################################################## def _createEmpty(self): @@ -1316,26 +1312,26 @@ class MCIsimulator(): ######################################################################################################### - def smoothingWithChargeDiffusion(self, image, sigma=(0.32, 0.32)): - """ - Smooths a given image with a gaussian kernel with widths given as sigmas. - This smoothing can be used to mimic charge diffusion within the CCD. + # def smoothingWithChargeDiffusion(self, image, sigma=(0.32, 0.32)): + # """ + # Smooths a given image with a gaussian kernel with widths given as sigmas. + # This smoothing can be used to mimic charge diffusion within the CCD. - The default values are from Table 8-2 of CCD_273_Euclid_secification_1.0.130812.pdf converted - to sigmas (FWHM / (2sqrt(2ln2)) and rounded up to the second decimal. + # The default values are from Table 8-2 of CCD_273_Euclid_secification_1.0.130812.pdf converted + # to sigmas (FWHM / (2sqrt(2ln2)) and rounded up to the second decimal. - .. Note:: This method should not be called for the full image if the charge spreading - has already been taken into account in the system PSF to avoid float counting. + # .. Note:: This method should not be called for the full image if the charge spreading + # has already been taken into account in the system PSF to avoid float counting. - :param image: image array which is smoothed with the kernel - :type image: ndarray - :param sigma: widths of the gaussian kernel that approximates the charge diffusion [0.32, 0.32]. - :param sigma: tuple + # :param image: image array which is smoothed with the kernel + # :type image: ndarray + # :param sigma: widths of the gaussian kernel that approximates the charge diffusion [0.32, 0.32]. + # :param sigma: tuple - :return: smoothed image array - :rtype: ndarray - """ - return ndimage.filters.gaussian_filter(image, sigma) + # :return: smoothed image array + # :rtype: ndarray + # """ + # return ndimage.filters.gaussian_filter(image, sigma) def _loadGhostModel(self): @@ -1398,7 +1394,7 @@ class MCIsimulator(): now=datetime.utcnow() #data_time=now.strftime("%Y-%m-%d-%H-%M-%S") result_day=now.strftime("%Y-%m-%d") - #self.result_path='../MCI_simData_'+result_day + @@ -1722,109 +1718,109 @@ class MCIsimulator(): #################### cal_PSF_array ################################# - # if self.save_starpsf==100: + if self.save_starpsf==100: - # self.log.info('calculate and save star PSF data......') + self.log.info('calculate and save star PSF data......') - # primary_g=fits.PrimaryHDU() - # PSF_g=fits.HDUList([primary_g]) + primary_g=fits.PrimaryHDU() + PSF_g=fits.HDUList([primary_g]) - # primary_r=fits.PrimaryHDU() - # PSF_r=fits.HDUList([primary_r]) + primary_r=fits.PrimaryHDU() + PSF_r=fits.HDUList([primary_r]) - # primary_i=fits.PrimaryHDU() - # PSF_i=fits.HDUList([primary_i]) + primary_i=fits.PrimaryHDU() + PSF_i=fits.HDUList([primary_i]) - # fov=0.05*min(self.information['xsize'],self.information['ysize']) + fov=0.05*min(self.information['xsize'],self.information['ysize']) - # dec_arr,ra_arr=make_c_coor(fov,25) + dec_arr,ra_arr=make_c_coor(fov,25) - # k=0; + k=0; - # ####################################################################### - # for ii in range(len(ra_arr)): - # for jj in range(len(ra_arr)): + ####################################################################### + for ii in range(len(ra_arr)): + for jj in range(len(ra_arr)): - # ################################################################## + ################################################################## - # galRa = ra_arr[ii,jj] +center_ra*3600 # ra of PFS, arcsecond - # galDec = dec_arr[ii,jj]+center_dec*3600 # dec of PSF, arcsecond + galRa = ra_arr[ii,jj] +center_ra*3600 # ra of PFS, arcsecond + galDec = dec_arr[ii,jj]+center_dec*3600 # dec of PSF, arcsecond - # fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, galRa, galDec) + fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, galRa, galDec) - # ############# do field distottion ########## + ############# do field distottion ########## - # if self.distortion: + if self.distortion: - # for i in range(3): - # ch=channel[i] - # fpx,fpy=distortField(fsx, fsy, ch) + for i in range(3): + ch=channel[i] + fpx,fpy=distortField(fsx, fsy, ch) - # else: + else: - # fpx=fsx - # fpy=fsy + fpx=fsx + fpy=fsy - # ################################################################### + ################################################################### - # psf=dict() + psf=dict() - # for i in range(3): + for i in range(3): - # ch=channel[i] + ch=channel[i] - # psfmat=self.get_PSF(fpx, fpy, ch) + psfmat=self.get_PSF(fpx, fpy, ch) - # temp=0 + temp=0 - # for iwave in range(7): - # temp=temp+1/7.0*psfmat[:,:,iwave] + for iwave in range(7): + temp=temp+1/7.0*psfmat[:,:,iwave] - # temp=temp/temp.sum() + temp=temp/temp.sum() - # ####rotate the PSF data - # if abs(theta.deg)>0: - # psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=False) # here we choose reshape=False, the rotated image will - # else: - # psf[ch]=temp - # ################################### + ####rotate the PSF data + if abs(theta.deg)>0: + psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=False) # here we choose reshape=False, the rotated image will + else: + psf[ch]=temp + ################################### - # hdu_g=fits.ImageHDU(psf['g']) - # hdu_g.header['ra'] = ra_arr[ii,jj] - # hdu_g.header['dec'] = dec_arr[ii,jj] - # hdu_g.header['rot_deg']=-theta.deg - # PSF_g.append(hdu_g) - # ################################### - # hdu_r=fits.ImageHDU(psf['r']) - # hdu_r.header['ra'] = ra_arr[ii,jj] - # hdu_r.header['dec'] = dec_arr[ii,jj] - # hdu_r.header['rot_deg']=-theta.deg - # PSF_r.append(hdu_r) - # ################################### - # hdu_i=fits.ImageHDU(psf['i']) - # hdu_i.header['ra'] = ra_arr[ii,jj] - # hdu_i.header['dec']= dec_arr[ii,jj] - # hdu_i.header['rot_deg']=-theta.deg - # PSF_i.append(hdu_i) - # ################################### - # del hdu_g - # del hdu_r - # del hdu_i - # k=k+1 + hdu_g=fits.ImageHDU(psf['g']) + hdu_g.header['ra'] = ra_arr[ii,jj] + hdu_g.header['dec'] = dec_arr[ii,jj] + hdu_g.header['rot_deg']=-theta.deg + PSF_g.append(hdu_g) + ################################### + hdu_r=fits.ImageHDU(psf['r']) + hdu_r.header['ra'] = ra_arr[ii,jj] + hdu_r.header['dec'] = dec_arr[ii,jj] + hdu_r.header['rot_deg']=-theta.deg + PSF_r.append(hdu_r) + ################################### + hdu_i=fits.ImageHDU(psf['i']) + hdu_i.header['ra'] = ra_arr[ii,jj] + hdu_i.header['dec']= dec_arr[ii,jj] + hdu_i.header['rot_deg']=-theta.deg + PSF_i.append(hdu_i) + ################################### + del hdu_g + del hdu_r + del hdu_i + k=k+1 - # ############################################ - # file_g=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C1.fits' - # PSF_g.writeto(file_g,overwrite=True) + ############################################ + file_g=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C1.fits' + PSF_g.writeto(file_g,overwrite=True) - # file_r=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C2.fits' - # PSF_r.writeto(file_r,overwrite=True) + file_r=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C2.fits' + PSF_r.writeto(file_r,overwrite=True) - # file_i=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C3.fits' - # PSF_i.writeto(file_i,overwrite=True) + file_i=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C3.fits' + PSF_i.writeto(file_i,overwrite=True) - # del PSF_g - # del PSF_r - # del PSF_i + del PSF_g + del PSF_r + del PSF_i ########## finish save PSF fits ########################################## ############################################################################ @@ -1919,9 +1915,6 @@ class MCIsimulator(): f2 = interp1d(earthshine_wave0, earthshine_flux0) earthshine_mci = f2(wave_mci) - # df = pd.DataFrame({'wave': wave_ifs, - # 'zodi': zodi_ifs, - # 'earthshine': earthshine_ifs}) self.zodiacal_wave = wave_mci # in A self.zodiacal_flux = zodi_mci @@ -1929,9 +1922,7 @@ class MCIsimulator(): self.earthshine_wave = wave_mci # A self.earthshine_flux = earthshine_mci - ######################################################################################## - - + ####################################################################################### ######################################################################################### self.cal_sky_noise() @@ -1952,11 +1943,8 @@ class MCIsimulator(): self.information['CD2_1']= np.sin(-theta)*self.information['pixel_size']/3600.0 #### self.information['CD2_2']= np.cos(-theta)*self.information['pixel_size']/3600.0 #### - ####################################################################### - - if self.TianceEffect: - - + ####################################################################### + if self.TianceEffect: ra_list = self.star['ra_gaia'].tolist() dec_list = self.star['dec_gaia'].tolist() pmra_list = self.star['pmra_gaia'].tolist() @@ -1993,8 +1981,6 @@ class MCIsimulator(): st_magi=[] st_magz=[] - - #################### generate star image ########## if self.debug: @@ -2116,12 +2102,8 @@ class MCIsimulator(): ugriz = np.array([umag, gmag, rmag, imag, zmag]) star_flux = sed.Model_Galaxy_SED(wave, ugriz, redshift, t, self.information['dir_path']) - - - else: - - + umag = binary_star[j-nsrcs, 2] gmag = binary_star[j-nsrcs, 3] rmag = binary_star[j-nsrcs, 4] @@ -2143,9 +2125,7 @@ class MCIsimulator(): ##cal_SED_photons(self, flux_arr): intscales,PSF_eff_weight=self.cal_SED_photons(star_flux) - - - + gx=dict() gy=dict() @@ -2180,7 +2160,6 @@ class MCIsimulator(): st_magi.append(imag) st_magz.append(zmag) - ## self.log.info('begin cal PSF') ###################################################################### psf=dict() @@ -2227,9 +2206,6 @@ class MCIsimulator(): cx0,cy0=centroid(stamp_img.array) - - - if self.appFatt : ### apply treering and bright fatter and diffusion; @@ -2266,8 +2242,6 @@ class MCIsimulator(): photons.addTo(final_image[ch]) - - ############################### ###### debug code # fi= galsim.ImageF(int(self.information['xsize']), int(self.information['ysize'])) @@ -2482,21 +2456,16 @@ class MCIsimulator(): self.information['ra_obj'] =self.information['star_ra'] self.information['dec_obj'] =self.information['star_dec'] - ##################################################################### - # self.earthshine(self.earthshine_theta) - - # self.zodiacal(self.information['star_ra'], self.information['star_dec'], self.dt.strftime("%Y-%m-%d")) + #################################################################### ################################################################################## - ra = self.information['ra_pnt0'] - dec = self.information['dec_pnt0'] - + time_jd=time2jd(self.dt) x_sat=float(self.orbit_pars[self.orbit_exp_num,1]) y_sat=float(self.orbit_pars[self.orbit_exp_num,2]) z_sat=float(self.orbit_pars[self.orbit_exp_num,3]) - wave0, zodi0 = self.zodiacal(ra, dec, self.TianCe_day) # erg/s/cm^2/A/arcsec^2 + wave0, zodi0 = self.zodiacal(self.information['ra_pnt0'], self.information['dec_pnt0'], self.TianCe_day) # erg/s/cm^2/A/arcsec^2 # EarthShine from straylight sl = StrayLight(self.information['dir_path'], jtime=time_jd, sat=np.array([x_sat, y_sat, z_sat]), @@ -2632,7 +2601,7 @@ class MCIsimulator(): print('k2=',k2) # - filename=self.information['dir_path']+'MCI_inputData/galaxy_Input/Lens_SED_IMG_0.025_230626/Lens_img_cut_IMG_'+str(k2+1)+'.fits' + filename=self.information['dir_path']+'MCI_inputData/galaxy_Input/noLens_SED_IMG_0.025_230626/Lens_img_cut_IMG_'+str(k2+1)+'.fits' self.log.info('galaxy_Input image path is: %s' %(filename)) @@ -2643,7 +2612,7 @@ class MCIsimulator(): srcs_cat=fits.open(filename) #### load galaxy SED fitsfile ### - filename=self.information['dir_path']+'MCI_inputData/galaxy_Input/Lens_SED_IMG_0.025_230626/Lens_img_cut_SED_'+str(k2+1)+'.fits' + filename=self.information['dir_path']+'MCI_inputData/galaxy_Input/noLens_SED_IMG_0.025_230626/Lens_img_cut_SED_'+str(k2+1)+'.fits' srcs_sed=fits.open(filename) self.log.info('galaxy_Input SED path is: %s' %(filename)) @@ -2655,10 +2624,10 @@ class MCIsimulator(): for kkk in range(1,len(srcs_cat)): t1=srcs_cat[kkk].header['new_ra'] -self.information['gal_ra'] +self.information['star_ra'] - ra_list.append(t1) + ra_list.append(float(t1)) t2=srcs_cat[kkk].header['new_dec'] -self.information['gal_dec'] +self.information['star_dec'] - dec_list.append(t2) + dec_list.append(float(t2)) if self.TianceEffect: @@ -2669,7 +2638,7 @@ class MCIsimulator(): rv_list = [0.0 for i in range(len(ra_list))] ################################################ - newRa, newDec = shao.onOrbitObsPosition(ra_list, dec_list, pmra_list, \ + newRa, newDec = shao.onOrbitObsPosition(self.information['dir_path'],ra_list, dec_list, pmra_list, \ pmdec_list, rv_list, parallax_list, len(ra_list), \ self.information['pos_x'], self.information['pos_y'], self.information['pos_z'], self.information['velocity_x'],self.information['velocity_y'],self.information['velocity_z'], "J2000", self.TianCe_day, self.TianCe_exp_start) else: diff --git a/csst_mci_sim/mci_data/mci_all_9K.config b/csst_mci_sim/mci_data/mci_all_9K.config index 885a24d..5b45b7e 100755 --- a/csst_mci_sim/mci_data/mci_all_9K.config +++ b/csst_mci_sim/mci_data/mci_all_9K.config @@ -62,7 +62,7 @@ TianceEffect = yes intscale = yes -ghosts = no +ghosts = yes shutterEffect = yes @@ -70,7 +70,7 @@ flatfieldM = yes PRNUeffect = yes -appFatt = no +appFatt = yes sky_shift_rot = yes @@ -80,9 +80,9 @@ sim_star = yes sim_galaxy = yes -save_starpsf = no +save_starpsf = yes -save_cosmicrays = no +save_cosmicrays = yes ############################################## ############################################## -- GitLab