From a10b95ced621ba3ff0da44b3fcc6553c424b3ea7 Mon Sep 17 00:00:00 2001 From: Hui-Hong Date: Wed, 2 Oct 2024 08:27:41 -0700 Subject: [PATCH] Hui-Hong HW1 --- .DS_Store | Bin 0 -> 6148 bytes Problem1.3/cal_Schwarzchild_radius.py | 15 + Problem1.3/matrix_multiply.py | 66 + concentration.pdf | Bin 0 -> 15679 bytes halomergerrate.pdf | Bin 0 -> 17335 bytes pbhmergers.py | 6 +- power_specs.py | 1902 +++++++++++++++++-------- volumemergerrate.pdf | Bin 0 -> 21316 bytes 8 files changed, 1353 insertions(+), 636 deletions(-) create mode 100644 .DS_Store create mode 100644 Problem1.3/cal_Schwarzchild_radius.py create mode 100644 Problem1.3/matrix_multiply.py create mode 100644 concentration.pdf create mode 100644 halomergerrate.pdf create mode 100644 volumemergerrate.pdf diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..28d6b3a2e4b0121ea441ce3e1e9197d2e489fa56 GIT binary patch literal 6148 zcmeHK%}T>S5T32I-6}#4iXH=Ai}pta@e*sji;;Rzsf{T$m}X0p+CwSis2+SB-@r%k zd7S;JD78>}5GgY-`_0bI>}J1&-3$P*YTqjX{xIRa=w z2THI5m5FF_{6_}lyIX-X=z#zizP&$p5Ov)^xH$gib^||(i^UI?T1qcxvOLRI`Hgcg z22R|J`%$&oyOO;NAp&=7H{H758Z`3jCnAiSe%R`%1mA06$n~WkdSXx&{m_e4uA>Kd zj^`Tr&EasrSg5c=Z#dSojW+>n%@D#t+-;bxqdNap%?`K>8(dmVO2-~>kc#q#6nM`7U7$64b zXTTgPe0BaQD3cf<1{Rb7c|KUE5^aNtMs;*xqe=k8G`f|bEo%vxV;Qs!CK_=Bh3Qm8 zoyts!!E`$MEfZ%OOf>3rV5azBW@Tn76sA_k`7IR=%+^RPF+dFbWMD~m^YZ>bonHU{ zOrjn!KnyGt11w#0YE|5lxmzbTC+}Jd^%hl%j7v0rrJ$j%VvLnn@fNBQ^jp$Ev<)U2 R(St%i0*VG|h=D(4;1e!jYpDPL literal 0 HcmV?d00001 diff --git a/Problem1.3/cal_Schwarzchild_radius.py b/Problem1.3/cal_Schwarzchild_radius.py new file mode 100644 index 0000000..116b64a --- /dev/null +++ b/Problem1.3/cal_Schwarzchild_radius.py @@ -0,0 +1,15 @@ +"Problem1.3(a) Write a python script, using pint, which finds the Schwarzchild radius of the Sun, in m." +import numpy as np +import pint +from astropy import constants as const + +units=pint.UnitRegistry() +def Cal_R_Sch(m_in_solarms): + const_G=const.G.value*units.meter**3/(units.kilogram*units.second**2) + const_solarmass=const.M_sun.value*units.kilogram + const_sol=const.c.value*units.meter/units.second + + return 2*const_G*const_solarmass*m_in_solarms/const_sol**2 + +radius=Cal_R_Sch(1).magnitude +print("The Schwarzchild radius of the Sun is "+"%.2f"%radius+"m.") \ No newline at end of file diff --git a/Problem1.3/matrix_multiply.py b/Problem1.3/matrix_multiply.py new file mode 100644 index 0000000..2c38d66 --- /dev/null +++ b/Problem1.3/matrix_multiply.py @@ -0,0 +1,66 @@ +"Write a simple routine to multiply two matrices together in python" +import numpy as np +import time + +"routine using nested for loops" +def matrix_mul_1(Matrix_A,Matrix_B): + rowA=len(Matrix_A);colA=len(Matrix_A[0]) + rowB=len(Matrix_B);colB=len(Matrix_B[0]) + + Matrix_C=[] + for i in range(rowA): + row_C=[] + for j in range(colB): + sum=0 + for m in range(colA): + sum+=Matrix_A[i][m]*Matrix_B[m][j] + row_C.append(sum) + Matrix_C.append(row_C) + + return np.array(Matrix_C) + +"routine using a list comprehension" +def matrix_mul_2(Matrix_A,Matrix_B): + return np.array([[np.sum(np.array([Matrix_A[j][m]*Matrix_B[m][i] for m in range(len(Matrix_A[0]))])) for i in range(len(Matrix_B[0]))] for j in range(len(Matrix_A))]) + +"routine using the built-in numpy matrix multiplication" +def matrix_mul_3(Matrix_A,Matrix_B): + return np.dot(Matrix_A,Matrix_B) + + +np.set_printoptions(threshold=10) +"test" +"generate dimension=dim diagonal matrix so that it is inversible." +dim=100 +mA=np.diag([i+1 for i in range(dim)]) +mB=np.linalg.inv(mA) +print("Matrix A is") +print(mA) +print("Matrix B is") +print(mB) + +"test1" +print("test1 for routine using nested for loops") +print("Multiplication of A and B is") +stime=time.time() +print(matrix_mul_1(mA,mB)) +etime=time.time() +print("time="+"%.5f"%(etime-stime)+" s") +print('\n') + +"test2" +print("test2 using a list comprehension") +print("Multiplication of A and B is") +stime=time.time() +print(matrix_mul_2(mA,mB)) +etime=time.time() +print("time="+"%.5f"%(etime-stime)+" s") +print('\n') + +"test3" +print("test3 using the built-in numpy matrix multiplication") +print("Multiplication of A and B is") +stime=time.time() +print(matrix_mul_3(mA,mB)) +etime=time.time() +print("time="+"%.5f"%(etime-stime)+" s") \ No newline at end of file diff --git a/concentration.pdf b/concentration.pdf new file mode 100644 index 0000000000000000000000000000000000000000..cf5c579162838899a300caf97db54d24f92bd9e9 GIT binary patch literal 15679 zcmeHu2{hHw_phnqhRReDzDR~}=gZBIdB}W8hA7iTnc|X^GG!=9<~c-CNM%YYDIp3; z2vLMYDhZ{$^Sx5Y>)-$P{%gIp-g>OF&iC%`8TQ_1-?R5V=W~S(RMjQXQaFTgem8Wb z0)d85kekf`gq$3NG_ZEC_k=KjLV}RSN8Ifpq>44g+S$zkQcytHyV}7GWfpZ%^Py-M zQLHKU5N^Jcx|=Hn!p=V%SW_tW9j1ctU8@JO?U(8+%tk zAkKf*wYPJ!R(A6NJisk*1eBBvgvH^d2#_UU(AonuUJ##T@9E~{VGHyD3;atTu!Rdj zF`3+3Sw7Q$fAX@OU6|jr~>7bUQ>OiC9&IiV6+jo1dvhbTLx+H3; zG8}l~W~!BbPn?!b>C3XMZ-?)PrFknWC*{Z#&xlM4j^F;E&_(j#yDqPo>-T{~78t6j z6`9N&^J6UR=u`0uVi)=nwo&pB11T@WIQ`jaGZj^o9g}#GctwUrM-ln^tzr zO4PI?r}$~ax~U#dVugqolDCJ4U8zE}+_$&tp4KFhp+=?a&jbrk;?;Mam%NKB8*yLf ztGG^WHm&4V+Kqw#{@|XID96}xEF5dL68CFc(X(%oEtUOlw+t#B zeXhOgJ^RIMcrP5hHvM5TyFU6f=ee}UCtO=cT@)JIRASrFN2|#sHDsK~sfm!nPxl6u z#Rk;9^4gr}DYu5~WHNuyPADBvT#e3?(cyWpdo+yMCSdQbI1ZC*LYOeMOvI5dqTZ^j zuQ5zk?Ijl2XxMSwk0y~~FLaz@wn-vic~O6*V*~NjkqxuER617Mtv7n4?-bVM&ukO! zpw07|00nJ|FJ}uaXQ_isvK|+*9O8*m%6~dQ3V|ulW8JFPm(g<9h^92dgk~D^Z$mFu z=1A~4Kq z>Qs{0aoo$YtiHySB4Di?CL~Z{*E%KqCZmgsce^fbwV8wy^1=OYc7FLc9Sa67D!CbY zeMo%}%^@X2N|9-M^@`!1NZ)1``+@UXAuz^=0jxoGy>i#apuY-Slv^R4Z5+DImc zSlT98t2_?^zKy;WKa^&^)4Asb@?}m(Ll7N#wX8bo+p5$sro1kdMmA4ww#i6WmQ5X# z>v8b`CN2FL3X?hK8#dQ45AOTWVP;?$7%zTgw`4voTUPt}wT2U7wb|yu!+jR*bdu6% zLmS-`)DM&o=T&^oHqkEZJC-l({e3grd5a+TllZtHL>p~k_yy6+Bbg<2DXosmec}Tu z^r7l+Z?WyLl02`XqF)27z4L+R$o070 z>r#xi=d?xcSvxA{i)b@2Y@S?SvSo8(3k_Kub@)i^>O%sh1y2t}xW-x;$#eL( zW{t_@X%{QLYfW^nxcoJgWnVqcnDpY%uGJ|PtbL90Z0m990p%?*g2S2TgY2%p%JuaM zs5ZYa8NJ0Lr{B2f%7%SKOx=%*JKY+CPx&PmB-cH>)WRgk!C5L2JhZ2DdN z`8r`c&E~>By&Dp@vjqh_KU=G&)3rYt7Tf!NVuv45{z>H=+1KUA#?7f>9||&S^L%)$ z4IN%`2Yj=+;*LMLm*SMQJ*{|jkcwHUSYoUuU1Fvr=Ck_B2BxH+S4
    -()WtN#vaBUrP{wqQay?TE zd&zQ?i?X9EF0oW15&tr8va{G)>2nohPddxm>!DQ#ZntT0Jxd^$V;as%Xto_XlqP1I zp)_1*S1omP^Qas9-t9*>GaH1P^u({GcOtJ%G4s?8Vi|mQkKZY3%ctREBUM4TE#It$ zg=IcL#r6AH{7!KD7+BxVb~alt%)*)0slgS@<{?DD6z$IrVZRl*^Lj?1ud849%{PTT zYChk##SWDr^>_B3IC847|Ek5FfJwpWsL0d&)w6dd5Qxx;C3=6kWyp4^AFCcXTpx=Zf7K7QP%rUyx9Tzlds zb7E^kK=_VNH3g-u64pv#L=Q>z14MNP5f9oS@@tJv%^k%D#Cb|jP7ujw3DxHRW}kM`Rt@l$-L z6$DxwqW&!^0|GB))Bqh#k7j)L9g$AMSoA&p_RX7uQ*(ssC)pjGK~WK+!Iz+!2Jxhv zPqspyh>TA=ZGxU13_e~|H!au36|AjuK_Ps7l%IMdN8@8Xrsoyq6DZ+V&)lWx*z#t* zv5vV>22YT=KVD;f&lXmuiPyyJ@Q>iJG_LV>4-zX4dl7LY-R00dF1_pRA7%<7wkLS4 zdg09O@WXE7Nc0eGRK1;?RJ`#?%YlkFWjD&~+#lt63B}Yy%w6pfdpC}u+8p}iRQ4Mh z9ZF5H?t9uFmeD#i+ceE{*Nw=jSGa?hCSFF? zuAY+mHqK6myx=sz(b@wJ;6ZHtGagrS1pzJ`pTmiRsuP)P529x{T(*Fur6Hs)2*Pm? z(h-j2=QD|49TXHmjO^s$Nr4B2PzWSEyYhTkE{R5g2_&(2G?<8h0)r4_;3429s`j3? z9!~BQHxCF+O*iIy>RMAgoTwqa6bc0=1K|JH6X>dA{c|P^9*;okS-Svn{vxT7mkouQ z`M_Lb7nq+{&>kl{iX%_~HNinF%0Xau7WAp)X$z+jL<|wmp4>I&epf5ABdsRc4MCtD>~2WNYDI$diYD!+hTJosotv3D_n>GP|xOrc2A zUwkeS3zlfrv7|#=&PS2Ev2MrQw3d zKo}w(0SZLG6CgYu3_`#{(m2po1~7mE&u|UD$ASX#R%y^59+N0T1ob=(uBpOcr4euq zGXV<)7XmOmxGx?EK0gL%3+AOhV{t$cGSq^XCQ@snG~gG|31|!uaDHAKm=})$iXp=N za1b8UGC)~yO)VI>9f7I@xQ6GaE*h*rzynY$0-h5O`cRd!z#LptEe)Ot2Y*1pDuUUD zMT2WlsOthxH(wBVAR@IVtl!@cfWi48F!QtHVXgq@3udseTnjZ=GQ2G3f$aibRk;2g48VS1Vb24*{rtHM*zkc8;q#K;+X!Z18-YE= zU)?d7V{zPoH`H>aK{MVO;{2hQQzSB;1k#i5b0r~$NP7B?j;-VO)VC@;h<%XGVDr`Drn`>5Rqxr6fbo);Gr}N^EG@kcJl<9Ro+trihFsilJ+e+WXY2PL< z4(V0rSoviu*UShGN=WTE=lUf@EHo}#Kt|G^KRna8iS_D+ZO6%5eIBLjjbFf3eo<BSYE6@UFWGayuSYqi z{3%TX=laGKm{=Sp{EZ1T_@*!BLOx2H?LNe;GL@>VK{Mw+QQZ1Au=|8XLlPy3!?ppb zx`#3v7h3x;_RN%{t$wiN-ZzxZdPdx5Pv0V2W;9r1THPv>Z&n#wpKBY<`nstyb0VJ{ zdUZDz9W0&1wt2t!VUKl_lescuUC~j_d#sz&HfP^JHJ+Q5UvK2u>=|}$R%oBcLF?wE z!HpXwuW%RG>$HqoUv`#c9yz0Q_|fDu-@*yhhVcpdUL^&)6+~Wa)qfKii;`X{aB>=lIV;uO!^3Md!Ub&;HV*YGF}>YY z?PwkrcUX;H>HVkTM@O07cfCANcYRG>;uBdha?v{qNy>F&r~qHxr%NB+T;v zdH*)^`blRs6*13qZ6o{cKUsS)-kHm1*Vk2i7kaj7vzl~vUD6vbkn3!o;x})n@UNiz ze@!niU~gP(ZtyNC4 zt<_!~C*>bq)N@3QFXMl;T+GC7pu-gEsp(F8W6Z|n2~j5n;Huz${Tfa6-CFnwY)g8*~G zk&_Cj-J7};@8;%~5;|s1R75A~MST#cZ)N9A@<|dBDk9f*ncG3+J|mr)yI+dev0t9* z;}eWMz1BM9ocr_qe%JKKjdR{{O{Qzc3Dc*{tk!o_>pb7Xk z8Lxs|BU24jCUY>Q^nNVlX`KDJJ`yuHEOSmI{q*>kE(veqan^Th31%u6b7YXOsckrfS zcUnu^MpmhzSEs(ET$Eq4jrgq^5$wrhKWd^W5`I}Qas(SGmM(tuybqVeSRKJaBwBs| zdPO3I*XNv%eYDf%Q7mEyxrIT3?Pi}HyRJ;bEiF<6cPwUv@1nP$U%p)3g&o6ox=-?j zk~&-UyOO52xR>pFVW$^YOW|vrPAz^CMC7cP>JKNZ(LM8Gcyx;0NgPSK#cQ5++>0%) zXVcIowgiV0+`V>lGi!c?nhXA^K@foacKUzY9+zuS*HG&eJ@|?aaeh@4u!^b7RY_7A znrnZ=dXmV%bTv9AKGAgTI&?63rU`L*6^P!XTWCC+WKi>ZeALU6JXb4_S(W~o?Tsj(B}SY z^A@=?-LBu%$@bfVSQ`~WMgC{ zx0f7!Fyb2Ub!P+V!}t2BY2K^r9?$MgT|tA3udDv1LFuJ?92Qe%0317xmMw&98)A^W zn(El^hD(MCw{$&)J7eNHH~yF@$Lx0^m}ouX>K~Oh^#Q}zs7`UkV7rr zy40d%+2{qBRgUUq3QY^qD>RWYVuUF-+YV<^RjVl!bc;Rk*^9vN;lR*R;ht6F1q=?h zvu_Emc(wrsFMpHle#cmn-_07QILb#yg6HhAr z`{cgKyEZ+2(nof74<)|YrFH?k8dB{ZS7>H-N00qc8jdY)NMq3 zaYcS7KVrw16=eO_SPqL>y5@MpTP`&y=285rjVS22@Z)|+^IX}96xkPfnyr=myR9%e zd!@YU#m%fg<*0a_x_B&|d~ud_GO9n3uT!AH=z09nOT=WCN5G1KxZ)81Qeu3_>( zn++SWs+|%mt2)Ll(Rt21^__3AwIvFg}%e0H>gc|x#wnWQ-;#ZX{?;X z;Yb74WCO;}J#W3M?l2G#UoIK(KKfuY+Gyg|Fu&hX(eC3uxC!Ob%a zXyyjXlbcPee8=(}xK3+docpXf^F^?QON}IsPl^^pdc_Kx8NJ+bEdBxOk&rj8b_8K@ zHDpYY@KiJ6qhI(S$5c;S{0$Xcf+O+xTCDhG9^O`MLbAI!dqqa~ee72j7Ye;eDBm>; zs5yIodnx2myKWQny(IS2g3P;gSI-rCoIP)Knn&UOM-`59(X4~`In(c-JXx+9ofSFf z>UQU~Rcan9{e~^k$(2Hurq`s{HaOvxbepQoCH-%eS^eNC7OG%Y{1X1)P%Yrmsi#XkLX=EQr~ac@7`{WU0R z9=-G&9%TRjW}IknI=Waq{MiN;AfC#2a^PBq=BDd%`lm8NDRT#fb5}7ZqUuFDxq{9` zsOan;O|*R_Y^-0PuvMt6LW{UQ&p?_%-%*Ddo0_fN7x(s#yo>|d+a9CnjwCX*sXkLR{+b~&CO&v|!Xw6c-AfaT_8&s_%B1x7DehaX&K~{&uzTtp3Q}*bllvSEJEiHyKyT21RI1Td8%MqsM-zqONg& zVRBR3wkkOMT3+X!ru`v}Q;2q%+6fk&Kn^On1zh4Y#_Eh`MGvThB`vVsIQ^MDjU!EOHZY0%Aaoyl(lDw+X zw4qEA!C=TPZ<>VU6nxh7d3eVq_g&`$_Ar)l@KgqA(F7fsUc;Azg$#xiAKtRTkmv)q zc6d4+Vhs1#7JeeTGQC%%AUpH3yjfC1Z;WGpISj3*SH!krxf3jeY z%jPZU=k7$kY(N>=WWFj|^=#0_KJfjUe*D$yurb3J%lq!{;~yVpY;pX-!}X}omNVL4 zBlj8`i<3-MSXH`}OpVn9JCU_2%r-4qg4xz_Tdr?pcvW-yjU4NTwS4B26YGrM9t&Mz zA^tVk!sC`09F{=XgQ8Kk)D~K03w3qZkpya@8tR@gG&a5oWEylg+Z43>jXSrZQ)1SZ zLgqpo&9`9I74li!zlK)Jd<~XRY=H8=&!GZF<{hLbtsvJwM^D&goAg<2Rt5n^)v461 z$f{cYp+lIrCq3uJ{L@_7ZRJlI2iy$Il$Sb4pT@~Az+L*WEH6G%%}l}Amt?Z8zArTY zVS_+-hQx<5mic!PSgv1>Gs=l9I(-3k_DiudC*Qs-3qh~wok0)8 zk=^|jymdLxZ@tVqedk0rEt+WsY5qAr0$;4f`vRUsbOlZsT-W(Q{;pl)oghA%y~~o< z$w7aJLW-3&H}-A_J!xDrN|`B+d%-3>aM=Fj1E%^bY-6Ml%}bRg0&U|-n5$)Y2ZlOQ zPW3?cz@hrUo%eHc*VfLIUc36Hj`u8|DW_1TY89sD!SLqzC+~~M7RTeAs?PcZr*3~P zS!=qA#M2`)c4{i~h~)LvVYidwPa$|s!nuwseaz4}n7_`5 zI1fqeseE(pisSuWEIP1vRp{~U?JlMRipG$ZtpPsM`=jq?zK%nQOpP(mit-Xw zx`%6Cp}8H(+U=jMAGN$HpQKcuzxRPK`$nBMYas{Lok?aILY!QX(X9yk~-*S*=DIvShTtbg}O1kW8(v>;bJZhT%X&{!*@{p9s zMr8h$&yxiDjfPJ7pAL^>!UsNgsZjx=;QiN)yDg4QrcY=32m=?E5rJ~ zMO+;knwh2Fx8EXc1s(otNQDL$Hj52hH%g0D1tdFDsnr`V1!Jp*Hgv{i~?HyjbUJ!MOO)eXq@b&JegeYwH|v{Fx&q-lb8Z6`nuRb$X}rRw#a za#>aDYWFg-*(LDbse%M8&(4@kn8kcp+6Iuq}J&r+^Vb5f)f1v`l)${>)xu4UETJX{toF&8kQi*=?(B1VE>Xq37s9NFco z@2__1@myCLw{XXu{m?R)My)0bS&Tm1qN@Z@Nu zssm|vwCAVIZd;k~v{O6YoQYE*igf3CJEP6AHWl$4BO3&)H9Zy7*`M?QlN+35xG`4Y z0|A<~u6?A!fvY`5Mk`L4yQ?>eKr$BZsx0EBC3I#l=Q#3*`)o3f=r`R@tN+b(U>dDU z^H5Rop>ha&tY}Lno%+n}^;S=s_*qPIR>&^E&DX!hxZp~4@v34=a3qugiU964$m)oT z+XffiYuF*kz#E$-i~1x^30%BRC{=F|K`g(L$?QQXQ7FUVjEw< zO#d}31-G<|S;B*a#~;P4dQ|__hSzsH*9agw?$&6z+xGUvSdm^UzS*%M@?wC9(qxhq z)cKayMD4YqR^y!R?L0N3ri)`&bQwgjk&|9qE+x16Hk*sO#BjZI_oGSW9s83(kj(s> z&M@HKZwZHnL{@NCGCBqB&$dP&ee|*j8<(5!-oTtXl<#l=ps$quj?$x{VL#MBXLo*A?2Um3 zZr48FSiR@#58Cq@#_@mZ2W%+6_46CeqO5CZxCG6z7{vlkTBxA(U!W8@Im9v)3%F!J zdb(TN0&F`#O}NA0aA$imfUHv=J*cP-a43fKbhh?%gpj`W9&QlQ%@v?aD2Krl#nHpw z9>QZ0NV1bRAbJ9f3cv&)U7bJ^0(ek?qyY>>2LcENqya<(6l;hGFeA1Q(hfq}gDY!1 zgmeOx41{!s&}aaZaRn2ACxEK~5Ey`yK>}{PAfz{hJPZTa|G?P*Y{_q&%{++azuJeF z;%q>ELRAxt7WwbsY@nYwn?KfI(MJ9kd=31KTylfF_#5!w!q;FiKnF|UMJu}>&}1H8 zvjWHedt42O@c$!R&BBgj@e3ERfRFeO>;I3d`F|W&11#--7FPpMBrxLUKf~1k$PhTJ z`ai(cV18k0U~J9*U0es%I|x9gs5ElIw0!R9cbzVKveS$oFU-i zs7+k~-0SBP?41Dghe~z>=rZ^TL3L{Y07j+$_JjZ@_nR{B5UFRFU=|MmFq?k~`#tJ` z7#L>WRpNfm|I@kt*bd>5e2;KsjzyYs&a{5>5=^nY=NQ~|$cqd2*_s=_}46jj}U zL1AzJY==R~VDM<1BuZKYg%S}%%(nx`T~QE5xI4R1oSkeSZ>$tv3L^nIQYh}8JCI11 zU$3OxJRCqLfH=1EvIW1XTinatjttpY+a3gq_HzWWMOyRmuqPwnjST@1g+KoQAQKNP znGF4+!60V%4|4rUL*uaUaosW+iU3Sv84bW1|DeIQ-OJly04C}W8WE^&$#~!c0W9<~ z8W#QrE~R1cU>TRuU@b4FVekM1x~v_R0PhY<$HU?P&~G^n1tQXAG#mlMr^{#<0{RbL zz&^Nq4m5@Uvbd#v(HJn`avGWlK())-VWr`H4F3DIMkt&#yg4qV;Zg9$wu}Z}8ZW2e zP|N2)0ag7$BcT3VBLWtHE0>K&kcJVDOXrnF!JFPP8XDY_EvI2of5;(?hqumUeTne4 zyp$$`f&IWT8WsS4m(vJ;&LJZM`-`Q0i75OZeTg^#@Lk?+1-l{={#Z*CQ3jl0E$xd& zVg9rY6bx!z-cAO;TnA`041!+P4vmLB5B&GbHb6E6z_-h2INV|#P&}-iob5fRKXle| uCA$F(GIgIrlHARrJf4UU7FN?&NBl3kyI&~) literal 0 HcmV?d00001 diff --git a/halomergerrate.pdf b/halomergerrate.pdf new file mode 100644 index 0000000000000000000000000000000000000000..a379ab4a0372aed5d991d04d8f08e9fb3f68f23d GIT binary patch literal 17335 zcmeHv2{hHw_b(ybM445%<_LFkFSn3+p2--Qr%Z84nvgN`JS0=*A&F8#${a$Iga)%B zl}bu^=X<4)*T4Vo{nvVHz4f%tI_IA6KIiPS&;Fjh_c@<$d^!qBA}CQY7+={CytEdE zf+OMX*3K|VNjO5s%E8VHjs`V)aD>4L4?8$Q-il!5;_d*KmWJ86*+PU8YZ@r}5|s4` zRs=h^*h(WMcQ*nYvr^WvA`t96-QZa0Ck&zE>29NMM}V7waTJt5|8~9vI6~DGj4Ags zm;afo!p-3beMc)>CpQN;Zl$@gwzadJ4FN7r?EE)7sI?Q2b#rqkc)?M~6%OP9H+F8I z0>AQG)6Uk(O3vLE@BmSW;gF&daEusM6bCm43|e^t;??Qv*?GBpd)febKokBY4@g97 zgaLKDDm^797cgHqLdgXPM8VF+-4@7C!_Lit;0Q;p2o4)qYe&XNn-u7v*GZ6??6in8cpQ*1D*`GC0 zn)2`>a!H+BB(8RYPV+IadmD;a7H5KE?%1Y0ka(U# zIbes7oD(tSs(Gha@92OCbZ;HkD^W8QxizEck>}R$i9ItAX=)g@H)JMoZ`r;tRVUlN z?*8OEQ0{o{eI};%ZFG0v8QtR#Jd_J#k{WwH)dt*?yq&gC)1ZVMNSQDbpG(FlYcpfc zn0I^h6}AkuoOdvPs}WK8b;`W_&F+eCyYJi>lM~fSmaeGVg>p~NXLc}qYV)xv>eEI0 z1?q%l5cawru10UpPtY35Tmw^GLyXC7C~Cv z@f1G77*A?5)}@ z89^+=xvKmbwOf{Z578C$P0kmHk*f*rHgom)rY4`AD0FO_^DQn`wHdcl&aUKYs)^~H z{FIG$F@f=Wv5q`#yOpmbDwH1D<9;ScT{Y1ZZloX9M8?#r;E+EpNEJhp755_uT(U~@e7ic?l=0#Jj+3qgf5DeE#SiCfb9x328El1s#nN4~o zs_T5E!!+QVl$%=ans~t*2O;XHE`*``eqJ?^hfSxB5yw~A z-8fu)V|%zBgHAZM6ch|gTl)kK4TUHaW2Z@R`yh_ug0PHCHC%npNm{n`fFm00y(hLL zHZCQ+CGS;xZt)HGh`i)C89SS>wa&r$@ct(W*LK%dN6Yvz=Dn%SzwP>=&}@FlzVOAh zdoAzZKcM7reUsI_+=FxTGc!A6cI3*`$o=IXb{uWM>3cU;S~0CXP#$@J zPfksu5cVSUQhs{gz3DSd%>r8!&)j~_#ARM@_`!xY*l?dfaCJcPj=&ub6KWFCZ}cej zwKCt5$ySGGULe{1-e2pE+v|+6ALk{f^i;e8jc(vizLLItj<1+%py){m-blXoYEIZR za?Co)sr55|os5oFyF%D7e`HPHWgi74lN8=;ny;4`KPq&wwY-))R@LQA4i5{}*)454 zmEpwzvmr(8Wa;v}Q=q`i7kI=gal-#t{P?xF{r7ifjOPYibW*N=KLWFsdm>JeWR+iB zd;hsa(<99VW;F}XkEHm)tz`__QfZu*Y^5X5Cug}{Zs=4MGrMT(LfB&$VR7nHhNZQU z%Q==26^+YpD9@v+@9rew?6{T(&ZPj_3CdK zzAO9qwDH5c$30{ET(%t$KsiZqD;JO3mzzCJj4CcoQHxVVdM3-DzEx+2ag@wb`D$LI zJY=Rmqx;goC_1P+`V3PlPwnD^q|}(O5gO67dlx#q@eMK}8HXKlVR|iEaR=D0wZEgn z#Gmv;8-KQ=zFPE{KQi2wRP5_8Rhy-?`rQcLEALu|*kpISjt^($4UucUd}es7a$EuP zHzq_I7FINKI=A_D$A`s*GVl| zGR9ep-ac@xK!3Ahy`c9ID^|n$Fvz57=fE{iF5CP;EkiP7&=$Pg6}7~t3Q7!=bBu** z$Z$;Lu~VhI$`g`(i9gPstrb5@^HRie!GmAm;f!X+GR1xA%wb+#%fnkt19q`3a(~b6 z2!2EfuTYUMyUYQn3ykQ_spXltb1h+xNgS@NM)u?d##V+yIpc6u{b$NCrUd^23um$q zxu?2=1sADKG&3C25NJVHAkcoPA?QK|T-!R!{HQ~hjSI;hIcb|rnXzY$bg2Rlspd}}u*npEp3HdJ z$|8=d&3o*iCC46Gjalkbjdk1)wcdW-gVdI57keF@I$7!;KihVuQS$Ck%vW+hOY+Xj z0k)4y6ChLbW<=dU!6 zySdBP3z@gHK3|e%`<&iol$|L#H+B4LILRT^_ygf(k$MR4B-5M7u$GJG6Onm-B|6bT z?>(t2OLSoQ$Sfb}LRg~6g(1;#o742HNIZi5bfQE4056+k^kY?p`csQjS}lTM!c*SD znmQuXFRS&PT7Q&VsP;dl8a;Kr{DSXZP18MzXDW4*_!nh5T3S>NrE5N8*63$f*H;T@ zI4-2I6klpwV~mGY)wsSI8#hjkj~@u^TQo?OGps+{eGyhg|A^J4MvtX{QR3!Mxj`^> zeeSU!f^esdan6Aj&f>7Tk*f0DY!WBF7RzW?^u{%ygPwj*;cT;&Je}gZv`Ce8`(sqi zlf`%L20p>x^}e4zI`o5I_@LgR))~>J7q{Hrei%)oCRn&Rszu3&^2xirMptaFtnTVP z>kRkm7G-{Q+06EbW4yGfsnO$7)d(|_*#lJjzSGO5ezGt3KIHmDsMa*<2&a{R>{gzw zx6PBk`|Ya|t&>ifo&SEU=6wtQQ*+&YZ?&?xefSPi)K-0GEYi;BZFbivm>=&TaC-;Y zTNS|DJ@=09^;f<~v8Q^r<3{GYuROe;POH7~IcHXUHt&dVw@qwKC;OeI#nvP44}vta z9=g|^UK*J^E>1c~`9kEN`;BI~No~uVh?91U1qC$AA+}yR>b0Q}$Ladb5q{Lk{o`IK z!%Me3=RaL)i>XdoC|D-lHQQ1Tg&IHoD@WoD{ zwEPb?JYKNu4+pK!Jvth6=1Th)$;Zq=>Kd1%!*|4eixuQhQ?_$@l7!hv*ji z;w2x9V_~e$eMI4p5sY8V<+_VPEeZSmBX^q{wD@?O6TBYwEb>IA>#;k`T37quEtW^_ zOYx?7=EC6c!xXd@vf8{cv)r=gaQf?sz7x%bN%Y3L94!QqIcDju=`Zl;w?#6%6+&*g0Dn zdh1)cd5LISyEq;5hC)?GD^DoM1o7X`m{HLU1YuAN2?e(bPWJY8Ai{$}Iy1PqI2@q~ z0!1-6!V!u#SEAit4Wy+(9OvZeMSwbmBVh<=bh(w#Q3Qnq1BhU-C@>HX2|B?^KwUs2 z3U*#Lo=zSFcTYHq7!t3v)U+aaIuXN8Q6v%ySi!$vCD2sE>gPykEEa~)vT_B}`O8f8 zy{!qv=p5o2v&#I6fQ+1M364Mn#2_5DCQgRfS(T@(mkku^;?a00!uL?Ib8>JbfaeH^ zf9Qz_8iON7vJN14#{Kj3{|nAx2vu;J;AA7~=HOxn4X0`4OXL@@iv>UR33jfAQ2okk zY!E2I_!pmRrUlJ&We>n1k?{Y41p80>z@Q}HVqkT_u>xEIBLT-@a4;L#^Z(aftqJ3AXwi18AWH^|6|;7@&+uNB|OIAc2nX5C+2u@z%QuBq0nH0E8~iR zaj|G17(CQY42}i41P~UK6B8Pu!x4o5<c# zK|hd?h#$uGJhY8Elmi zB3EMKD9~VK9U(nf$0{@@qA~%6f?ECR1!4+`gb|

    ;R9$gO_wAL5i@V1CSo9PefS* zy?}ILO(N<8q!X)&s4G8p0@8=I3F3?R8}tF`#%h9Qjzda-H3YmuIRhQNp)V*W70bgD1&8O|#Ez!r@w@{5K9!;JsSMVyrrS6P!kVAwy1?WZ8ed z`r&lo&{?z2G{QMXn@)s+5n(1Nw5>HUYQfP)J4od4IANETK3h!0Eqn9qPAfEUE7NZ_ z8d$|Wp2_{nT9-3lW*=H{2!jd|Pov*uDRkVE?$f72iP^#USQjI!-iy2Ps*qiAOHw=Z zy?VXE;+A-hcpkOtO`Bk05h-OWx6|mGu_|#9p?Mu8d;H#~M}C**kvr$+$zI4x+io)D zwL0{-DPxf0>!%#A?nZ_JQ1UAag!%Z*A)Hvx5SB)93dLFe^GI*J_yPbHs`z z1f89J=TUj_!vU8bd;6EBxhdh}iIcrd)^Br94Rk4)n+UYL4zk~0Go_1Xp~PvHWR(i% z4~zQJ_VZ_&ecpMnt~ZU$@t-wr-^Zof)doO=I&p#7CM1U*qVmQXH^ z_76qx#xFfRD%bPn%09l2Do;Z&o&s+jTPGg1kHr3X9!0JmbJ>w={$+*R6H&`L{r!s8 zwpyNLP7#N+7IYXN-fXI;+40FmQC`q1?(yp*P4~AQO?F}SRr^Z8acP)eoz`$*ut;mJ zTymgyfzz~~z`2R+|J9jCgI#Z}x?#2P@#I`IohO2&k%w3ZWo{Q1UdKIMJX;%|rWN~+ ztK%U9dzx<=56>0*wn0-{xcy6gjh5S=tLYd@7G832Cq`_u3W@W0Qa0+A8O^fnlhk8O zIfwgl-o#?Z)4Li^jF=rn2PtfvGMe!nZjt5;l~La1h5E+I$Q+s=!u6&^xR*iR$Bt<9 z#FRvwKxV|;=Rsi~{2AIe+i+1_Op{-P`~q?3a%WD=ZQ121yKFxEg+uYy6d!yyHS-f) z>$e%>K}YI+2pre;(tF9J%Zf5^bx#Ww>`2!stGbk}r|VHHT{u;ib$0y7VUbPb_0K_t z7=E3+q;w8MD{#@g!3QJW@qXa(;7}?Rv7dW!nJisUN8`+O3~%qw0A(6fV#O$5dG_Ui z+Yxf_-A3M=y9+Os*Mr~pqYER<6yV|*_OoYT4qA!ZL^mp;wq{Hd3Zy9c@!#&kg1mOy z%^0c(gqLtfzs5uhW(u9m@?{pDZO54j#7m9AU+L+EcjRX!-aF`eFA=uizK>j({^mDXC90fkyr^_VUFYTYQ)x2NElQ!^Bdob$M&uq1l+6Ww7Uoxui zpTjfNE{ulbC^e&=P0cJYI0+&2Zn2wQJmXEDG|W25N}u9zmhFY@@*?GrP*d(rB>1lv zI_M1&)YMfyPllwqHOcfweKI61A2y58K_7v-6!)&QW@k%KG=4V0Y0H^?zL zq&BBj(Km}GlDHVOXuS8(n&#cTan!nL#I#Q`YRK)IlD!>&(DoQxxgl2nF!-A-S6PJ< zaX1FF$$hff?H$`qE&Gdi>qZ3bOMF62$Te&(>3Uf;@x(`5UgX2eqSNocIvS=}2KT9d zGWISh--9(TW_wsi6B)Gx5T zcjYD@4{v{>_GxbUvVtWsPhReVFRIiOK1nfcN7c)sX5xw%v!1;lt^BA|Fu-K`>1|(r zcdl;ty>@*2j+{)U+N&oY(GI-xJFHSWm)Om zt*#suJvVEIezg2Gl=<-7=>_XHQGe#L$%Z?$Bt?Y`)fDY)&)(t$9Z3$=MyTGs^uf=~ zC)efDq3OjvBp6zxXJZ(j8GA$DbCTrK=2MDM=fCfh(KV4-YJco==Q*lF#e;cbYjnMb zWJ{Wdx?qhnTt zrkENL%=WX&Y;Rk?qHXFP8^Z;r$Y|oiOZ)P^`n0f>d>nghefd)EVZ?rJ`mXrVQW>68 z`>(boT)&R_PO#vd;hJMu7806tvmzBOc`LK1$NBzI%z*Ce^qd%>eMdQ*pBV+2?^GL< z7A-_;4cC*IzN$|um3lUJki%#HhfOB?udb&!dff|yLD%R2*Hes?K7^ScHetV&=z2C3 zhu@SFWl^DNzDHU`<6q5}5g6v+p9*)~(E?{ys7woIQ((^^jq7}=Mn=Ln$lS3BKmY25 zimmr|HsvRIhoLBohA?HU7GVEH{^7ly#h_Xb)K#YnK2m1vI zjlmF4TK;m=j={qazs^UL+&x`aPZL%F1Z(|-QM4FRABn>rY9B|7;90TC!#d7G z7U=xLqTU@sCRQKv<-N}roXWH>Sfc$DJDSQdz*VdNB>5yOb8`3zHTKGM;RTOvxCWAR zdarHo7pUIt^Ve-u$aZ_bKv!!*FSCnriE%2uG}+jC%l3%i@&Tq1!(l<=8>L@-cyCft zd!kIH^l25(3)VEAViO*SGtGG8X#!hzJH2I>Rj|ZNRfNLJ6OCjVygzLHU3#R=W?sKm z#QB6sO!xXn_+RMeC6(2eD3!ip_09R##AMoD0SWw*id7{gT zMv7waD{K77xw;i_yps^@Z;nvd9f5n{F!YSp_1L~57lw$9zqPY*+Wg3uiIS-WQL&} zm1gkXJj?#;aL+b6R+>9$3=!NkYMNV@FMGyhSw!rXZh9}z7#B}Ffn7HK{=sW|g?@}c zoSS>Yh($&*E!oaJ@#%Fu=Ejwx^gErfvYI`OrXv2gYAk;2uI8zwk@+0nd@NlrI_LIL zkNh6QxPboCO|e^=*{Rb6 z-7rn^TQg~dL(Go_d`w1eJFKh|eR1ZJ)#HZO~3M!JFJ>I9j;*YzTE;xDveWnZgEgXNM8MZR{kD!e2Sy`NX~~uv*der3p{xDj@ohg zP^~1V4{p8P6_r0O(QNTq)-l*o&R40{b@+i>rI;#Wph+Ss<>NhiCYOBQjXr%S z7%(k`zGGiik4KhD-KUM!y3tzk`r&~jUb67o*m=2p*c*;Ve7Q0E22=GeYKe@qresZ> zJQYnM)$%-QsIRS)B~s~{ND;WIm!Ly3==UWZD~$_CWY{Y2pll?zEx-MqpB(mX6h0@d zzRAOqtg|X8FTY%`j?RvBYK!@^{i9-Pfzc#)U+j(i^mS{=l7;x6aufSkCkq9R;nr!r zHlDT#Ni+Bzb_VY6ZG8QC0{2wxkF(`d+3e&yvwhpIUb3N}MFps~?Fmn*=3Eksk`X?a z&f_Q7f)w_^sRo3~8mbyen`LNlU8As2`SOHOa<`^!U}jqw`;CN-DY7{|vr6wN(%Bt{ zc+IP5Y%0rZJ*2G2u`PZ~6Q7Va;ql5@ht%JYQZRgZuz2>Z+nmp7Qp*-3@w8_BZ4dCd z{_Ukkfiw8E)5FGeZU@topGyziM&8?VRZ{zWb|_)_C|@B3O)9cOV1W5tT%^2)`pcy020O_bUwru^ii@%B!zQj`3T)kyA(skg zn6$nkyb_9%*z9QA%k+n}+V2H2rLsKoe9?2G+%t!Pm{6E{CH z6~pw)g0iuv)yr9p-p;okE+i1AA38RBYl^u5+B@NjI)CPfarrs?vK^?bLGj_JHZ* zc7-U9sHO1ftp-cluMa1_(>zxZkNV1LP<0yQbC@Y0JG9)ogRs#q6|HHHG+;H>WG*K87^LOG+!FI{asF z;jXz9{F5B{d|l0q`>(WZsk2Xid}zw@E=QEvW?D$J zSj;ovX&b|AE~X39*f{bR`bCE~_GVV)TeY@vn0}h4GnhUVy2(QPt5b(v?@(iKI3py9 zoVhv=g$;ags9le%CAPEu5qVeFcp&wJhY9PsL*pK7GES+vdoI&l79;r0&WK8NQoXLyU4`J#diIDO>uIxqKUMDpPnII$IH;N)=sXWY~go- zszM=LIuWV6;zaD+YWxRr+pk_t($9~+5^)I`^SRoEiQ`DF8Mk-*!E?<*h@sJ1_V)ZI zx1Q%lG@Q*NMNw}un}2m?Fu(w;oefqG?*_~vxD-AO_jend8vuUg((UcV1B_(H2zrSJ zObvWGLxT;j%@7u=lb+Fwj~%xQZl>-irJvOcQ7Ni3$!V%Gmcrq|KN@Z}+go z?D>V96Czi)hSjGfpNFv- zCvE4piOuCMO+Kcf7qg3sZ>21Mh7ZY2(>3(1kx$s(jcVhXr_F3)qgMMW$#^@4kAXe0 z{RQpzoU=zbE=6YGo^P4RJ?P3rNjp;2WN|ZqNsZ-b`~qpQ##F%ca>e^=4gLL`?6dx@ zn=H=1x>4v2MoC0h)ssw?>Q>Lbk3j%{XOprobZohw>K^6R?GZ8r(~<1REV!^y-FRH7 zW77)^D)0qG=$U=}uEt|B25?my9c+%zd%uqyPmiTicO^WzLStV>w-hNjG}ZD7#pY1c zZ}(`&jCqArnruhe;buMt7LCVNJPx!6)6k;TUQ*)lNhRYHt5`R~T^Rakv0xRh;F1e@ z_sM1HBQGmO;XE0?^bS}WgL_Cy_8#LE zCD=|yM(?|@XSA^jZZpy+vD|S(x)~)qoGW6S;?p+JMelVNjoNZusbfGgw~?;xFcrOR z3THzjoZCER(Qw`*;bTwDf#I0iYWpKnxsBm7R4rFXjr_mYwgsr)U9@;OmLkMw!w~u* z;5D}%?m*6tYX&o(CFiS#FYNPR`nc35|8|z*8js|QKc@|xu;2d7|50o!@!tq;74`7Z zz!{n$ZVD&V{IvZQ&v!2mUS#8Ys#dLfpU=F{tGVobgJ*ilkwX?C=ZBYH6OI@S>^UTU z+%U$O(;}^BL({#S+j8jGv!8CE(<#61E8I+Gf9pM^+Gi0bFI`MEed>dQ7LJM& z@3vP=J*X?#=JND%+LqGq!R~ZBuM2Akk1OjPSGVnnb%BkMf)}7JxOug_^i} z!Xvt9V&zz;7ECj5aKAw9s$XCVFrR{H91n4MAB1RbEP49?07EF-sVd$3jqbRRUAtPm z^%*SS{>f+s2fah_ULSV3^HO6;7xs@wCCTFjh6;ULP$s#oS9YJW*9q8WeE!_PXxck; zVNjYbOQQ5U9DK>D?VdcNSlfAfl+1#uhfC@t3Su(;pDx!pK_lc)@A62EogZ#EHJ+^>BTFxSNN zkUB^lxN)>sf)(&jthsTth1mmZY?3Z$Y)rz&Z$F{Kl^U`~)JEScn^Z=zwZL*$$-&Tl z=k{$UNs6VP+eA(O)hz{BfVC`PL7?c5q)|Aj{c7jP?Ey+I*wfoBsvb5kh7&CGMr6kK z?~E=85Rm7wS12L zC(HDmh6;``G_zwVUJJg?xKy{}W>&GoQ?dAAv8s5l8{2Q6?9uOEE-bC*TJBrzPW}nOR`^A;TMq@YT-U+x_EqcVwR-$DP9=)rnkbAM_H6xW@IR z{nA7;YKQu}oW&-St3|y0Z(OfKl*c!Q+lQUIVKCcrOn?*U_w7C#dfk%GbnOPE=CtTW_V*N_5e-)sC8$VGg;k%??^2Y5{6@)TWqt` z6063V-77yDGA$v z00ZC!gqMev4fu8mAdWnskAGb3?7=5V#G)q=5e1Go5ne7W9O3Q;kXeM| zphR%=w6lXD?7>GuaD=^+52*J75HA3!LAW^qDsh0H(tt2#0R4l|f#X5F6$J5sBW&Ra zJ5U0+Ee8mEhHwHuz`pJRM*$240(9*GN&uGzAY=e%h5#IS!x27k#Bp%`2?$O=-zxn9 zR|62O-*B}REZTo%Xx78kKmbOR7zDcd@4(gIKjCVBEXNv4^%s~LULpa} zA`x%l0ks2wBd%0C!hyjf5;y_iCD~adOp!r&p`Exu^I8bC_L&z+O{Ty+{n2Hfc{bhcMbA$myC-##8qHTo(<4bH0 znOlHHB$mD50Czzw`2c+8O6fSz?Um9ApkGi4wqmWpSr$YJT&ABN2dr}+QE;MT0cNR$ zhJ$5ZPrBA8lSL5$X)z>r{9 zHB{c+8z95M+|-?HAqYM2Rv?-E0?)5}J^q^&>>92S^5Cmyf|I+O0`#%?UWNT=Bw7r3 zLujM~I4;47AjJicNC83E3LQYt?*%@Hhl@MG#mO4(gAv7wqJ`m(1cHaxegwkx*CSDP zPY2Klp#E*WZNRtJYg>8P+QY4_Y>tAp|Jj4!nu+;(+S$Xvu84yHRPWC}0E)%}`PjpM z)j=N;K!0$zpLHlK9>R`psDpNdKkA_AZ=e$gUf{+$2-y5b9rUKHYYTqhHui@^VxZS| zJsk!@ns2Csq`a{XgNELY^=-v4z$x8Whm`oUEd(UrKnDc3QC3I{1SsA>hk-8h*Vlny z4FFm<)M4=8HgaQ~*q?M*G_e6MQQ24r1pP;y1hgrxZ!3Y6_=9H&G-MJs&|x-d3*M?f= value-eps)*(array < value+eps)) + # Floating point inaccuracy. + eps = 1e-7 + return np.where((array > value - eps) * (array < value + eps)) + """Just rebins the data""" -def rebin(data, xaxis,newx): - if newx[0] < xaxis[0] or newx[-1]> xaxis[-1]: - raise ValueError("A value in newx is beyond the interpolation range") - intp=scipy.interpolate.InterpolatedUnivariateSpline(np.log(xaxis),data) - newdata=intp(np.log(newx)) - return newdata + + +def rebin(data, xaxis, newx): + if newx[0] < xaxis[0] or newx[-1] > xaxis[-1]: + raise ValueError("A value in newx is beyond the interpolation range") + intp = scipy.interpolate.InterpolatedUnivariateSpline(np.log(xaxis), data) + newdata = intp(np.log(newx)) + return newdata + """Saves the figure, automatically determining file extension""" + + def save_figure(path): - bk=matplotlib.backends.backend - if path == "": - return - elif bk == 'TkAgg' or bk == 'Agg' or bk == 'GTKAgg': - path = path+".png" - elif bk == 'PDF' or bk == 'pdf': - path = path+".pdf" - elif bk == 'PS' or bk == 'ps': - path = path+".ps" - return plt.savefig(path) + bk = matplotlib.backends.backend + if path == "": + return + elif bk == "TkAgg" or bk == "Agg" or bk == "GTKAgg": + path = path + ".png" + elif bk == "PDF" or bk == "pdf": + path = path + ".pdf" + elif bk == "PS" or bk == "ps": + path = path + ".ps" + return plt.savefig(path) + """Little function to adjust a table so it has a different central value""" -def corr_table(table, dvecs,table_name): - new=np.array(table) - new[12:,:] = table[12:,:]+2*table[0:12,:]*dvecs - pkd="/home/spb41/cosmomc-src/cosmomc/data/lya-interp/" - np.savetxt(pkd+table_name,new,("%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g")) - return new + + +def corr_table(table, dvecs, table_name): + new = np.array(table) + new[12:, :] = table[12:, :] + 2 * table[0:12, :] * dvecs + pkd = "/home/spb41/cosmomc-src/cosmomc/data/lya-interp/" + np.savetxt( + pkd + table_name, + new, + ( + "%1.3g", + "%1.3g", + "%1.3g", + "%1.3g", + "%1.3g", + "%1.3g", + "%1.3g", + "%1.3g", + "%1.3g", + "%1.3g", + "%1.3g", + ), + ) + return new + """ A class to be derived from by flux and matter power_spec classes. Stores various helper methods.""" + + class power_spec: - # Snapshots - Snaps=() - #SDSS redshift bins. - Zz=np.array([1.0]) - #SDSS kbins, in s/km units. - sdsskbins=np.array([1.0]) - # Omega_matter - om=0.267 - #Hubble constant - H0=0.71 - #Boxsize in Mpc/h - box=60.0 - #Size of the best-fit box, for testing varying boxsize - bfbox=60.0 - #Some paths - knotpos=np.array([1.0]) - base="" - pre="" - suf="" - ext="" - #For plotting - ymin=0.8 - ymax=1.2 - figprefix="/figure" - def __init__(self, Snaps=("snapshot_000", "snapshot_001","snapshot_002","snapshot_003","snapshot_004","snapshot_005","snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"), - Zz=np.array([4.2,4.0,3.8,3.6,3.4,3.2,3.0,2.8,2.6,2.4,2.2,2.0]), - sdsskbins=np.array([0.00141,0.00178,0.00224,0.00282,0.00355,0.00447,0.00562,0.00708,0.00891,0.01122,0.01413,0.01778]), - knotpos=np.array([0.07,0.15,0.475, 0.75, 1.19, 1.89,4,25]), om=0.267, H0=0.71,box=60.0, - base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/",suf="/", ext=".txt"): - if len(Snaps) != np.size(Zz): - raise DataError("There are "+str(len(Snaps))+" snapshots, but "+str(np.size(Zz))+"redshifts given.") - self.Snaps=Snaps - self.Zz=Zz - self.sdsskbins=sdsskbins - self.om=om - self.H0=H0 - self.box=box - self.knotpos=knotpos*H0 - self.base=base - self.suf=suf - self.ext=ext - return - - """ Get redshift associated with a snapshot """ - def GetZ(self, snap): - ind=np.where(np.array(self.Snaps) == snap) - if np.size(ind): - return self.Zz[ind] - else: - raise DataError(str(snap)+" does not exist!") - - """ Get snapshot associated with a redshift """ - def GetSnap(self, redshift): - ind=wheref(self.Zz, redshift) - if np.size(ind): - return str(np.asarray(self.Snaps)[ind][0]) - else: - raise DataError("No snapshot at redshift "+str(redshift)) - - """Get the k bins at a given redshift in h/Mpc units""" - def GetSDSSkbins(self, redshift): - return self.sdsskbins*self.Hubble(redshift)/(1.0+redshift) -#Corr is /sqrt(self.H0)... - - """ Hubble parameter. Hubble(Redshift) """ - def Hubble(self, zz): - return 100*self.H0*math.sqrt(self.om*(1+zz)**3+(1-self.om)) - #Conversion factor between s/km and h/Mpc is (1+z)/H(z) - - """ Do correct units conversion to return k and one-d power """ - def loaddata(self, file, box): - #Adjust Fourier convention. - flux_power=np.loadtxt(file) - scale=self.H0/box - k=flux_power[1:,0]*scale*2.0*math.pi - PF=flux_power[1:,1]/scale - return (k, PF) - - """ Plot comparisons between a bunch of sims on one graph + # Snapshots + Snaps = () + # SDSS redshift bins. + Zz = np.array([1.0]) + # SDSS kbins, in s/km units. + sdsskbins = np.array([1.0]) + # Omega_matter + om = 0.267 + # Hubble constant + H0 = 0.71 + # Boxsize in Mpc/h + box = 60.0 + # Size of the best-fit box, for testing varying boxsize + bfbox = 60.0 + # Some paths + knotpos = np.array([1.0]) + base = "" + pre = "" + suf = "" + ext = "" + # For plotting + ymin = 0.8 + ymax = 1.2 + figprefix = "/figure" + + def __init__( + self, + Snaps=( + "snapshot_000", + "snapshot_001", + "snapshot_002", + "snapshot_003", + "snapshot_004", + "snapshot_005", + "snapshot_006", + "snapshot_007", + "snapshot_008", + "snapshot_009", + "snapshot_010", + "snapshot_011", + ), + Zz=np.array([4.2, 4.0, 3.8, 3.6, 3.4, 3.2, 3.0, 2.8, 2.6, 2.4, 2.2, 2.0]), + sdsskbins=np.array( + [ + 0.00141, + 0.00178, + 0.00224, + 0.00282, + 0.00355, + 0.00447, + 0.00562, + 0.00708, + 0.00891, + 0.01122, + 0.01413, + 0.01778, + ] + ), + knotpos=np.array([0.07, 0.15, 0.475, 0.75, 1.19, 1.89, 4, 25]), + om=0.267, + H0=0.71, + box=60.0, + base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/", + suf="/", + ext=".txt", + ): + if len(Snaps) != np.size(Zz): + raise DataError( + "There are " + + str(len(Snaps)) + + " snapshots, but " + + str(np.size(Zz)) + + "redshifts given." + ) + self.Snaps = Snaps + self.Zz = Zz + self.sdsskbins = sdsskbins + self.om = om + self.H0 = H0 + self.box = box + self.knotpos = knotpos * H0 + self.base = base + self.suf = suf + self.ext = ext + return + + """ Get redshift associated with a snapshot """ + + def GetZ(self, snap): + ind = np.where(np.array(self.Snaps) == snap) + if np.size(ind): + return self.Zz[ind] + else: + raise DataError(str(snap) + " does not exist!") + + """ Get snapshot associated with a redshift """ + + def GetSnap(self, redshift): + ind = wheref(self.Zz, redshift) + if np.size(ind): + return str(np.asarray(self.Snaps)[ind][0]) + else: + raise DataError("No snapshot at redshift " + str(redshift)) + + """Get the k bins at a given redshift in h/Mpc units""" + + def GetSDSSkbins(self, redshift): + return self.sdsskbins * self.Hubble(redshift) / (1.0 + redshift) + + # Corr is /sqrt(self.H0)... + + """ Hubble parameter. Hubble(Redshift) """ + + def Hubble(self, zz): + return 100 * self.H0 * math.sqrt(self.om * (1 + zz) ** 3 + (1 - self.om)) + + # Conversion factor between s/km and h/Mpc is (1+z)/H(z) + + """ Do correct units conversion to return k and one-d power """ + + def loaddata(self, file, box): + # Adjust Fourier convention. + flux_power = np.loadtxt(file) + scale = self.H0 / box + k = flux_power[1:, 0] * scale * 2.0 * math.pi + PF = flux_power[1:, 1] / scale + return (k, PF) + + """ Plot comparisons between a bunch of sims on one graph plot_z(Redshift, Sims to use ( eg, A1.14). Note this will clear current figures.""" - def plot_z(self,Knot,redshift,title="",ylabel="", legend=True): - #Load best-fit - (simk,BFPk)=self.loadpk(Knot.bstft+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.bfbox) - #Setup figure plot. - ind=wheref(self.Zz, redshift) - plt.figure(ind[0][0]) - plt.clf() - if title != '': - plt.title(title+" at z="+str(redshift)) - plt.ylabel(ylabel) - plt.xlabel(r"$k\; (\mathrm{Mpc}^{-1})$") - line=np.array([]) - legname=np.array([]) - for sim in Knot.names: - (k,Pk)=self.loadpk(sim+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.box) - oi = np.where(simk <= k[-1]) - ti = np.where(simk[oi] >= k[0]) - relP=rebin(Pk, k, simk[oi][ti]) - relP=relP/rebin(BFPk, simk, simk[oi][ti]) - line=np.append(line, plt.semilogx(simk[oi][ti]/self.H0,relP,linestyle="-")) - legname=np.append(legname,sim) - if legend: - plt.legend(line,legname) - plt.semilogx(self.knotpos,np.ones(len(self.knotpos)),"ro") - plt.ylim(self.ymin,self.ymax) - plt.xlim(simk[0]*0.8, 10) - return - - """ Plot a whole suite of snapshots: plot_all(Knot, outdir) """ - def plot_all(self, Knot,zzz=np.array([]), out=""): - if np.size(zzz) == 0: - zzz=self.Zz #lolz - for z in zzz: - self.plot_z(Knot,z) - if out != "": - save_figure(out+self.figprefix+str(z)) - return - - """ Plot absolute power spectrum, not relative""" - def plot_power(self,path, redshift,colour="black"): - (k_g,Pk_g)=self.loadpk(path+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.box) - plt.loglog(k_g,Pk_g, color=colour) - plt.xlim(0.01,k_g[-1]*1.1) - plt.ylabel("P(k) /(h-3 Mpc3)") - plt.xlabel("k /(h MPc-1)") - plt.title("Power spectrum at z="+str(redshift)) - return(k_g, Pk_g) - - """ Plot absolute power for all redshifts """ - def plot_power_all(self, Knot,zzz=np.array([]), out=""): - if np.size(zzz) == 0: - zzz=self.Zz #lolz - for z in zzz: - ind=wheref(self.Zz, z) - plt.figure(ind[0][0]) - for sim in Knot.names: - self.plot_power(sim,z) - if out != "": - save_figure(out+self.figprefix+str(z)) - return - - """ Compare two power spectra directly. Smooths result. + + def plot_z(self, Knot, redshift, title="", ylabel="", legend=True): + # Load best-fit + (simk, BFPk) = self.loadpk( + Knot.bstft + self.suf + self.pre + self.GetSnap(redshift) + self.ext, + self.bfbox, + ) + # Setup figure plot. + ind = wheref(self.Zz, redshift) + plt.figure(ind[0][0]) + plt.clf() + if title != "": + plt.title(title + " at z=" + str(redshift)) + plt.ylabel(ylabel) + plt.xlabel(r"$k\; (\mathrm{Mpc}^{-1})$") + line = np.array([]) + legname = np.array([]) + for sim in Knot.names: + (k, Pk) = self.loadpk( + sim + self.suf + self.pre + self.GetSnap(redshift) + self.ext, self.box + ) + oi = np.where(simk <= k[-1]) + ti = np.where(simk[oi] >= k[0]) + relP = rebin(Pk, k, simk[oi][ti]) + relP = relP / rebin(BFPk, simk, simk[oi][ti]) + line = np.append( + line, plt.semilogx(simk[oi][ti] / self.H0, relP, linestyle="-") + ) + legname = np.append(legname, sim) + if legend: + plt.legend(line, legname) + plt.semilogx(self.knotpos, np.ones(len(self.knotpos)), "ro") + plt.ylim(self.ymin, self.ymax) + plt.xlim(simk[0] * 0.8, 10) + return + + """ Plot a whole suite of snapshots: plot_all(Knot, outdir) """ + + def plot_all(self, Knot, zzz=np.array([]), out=""): + if np.size(zzz) == 0: + zzz = self.Zz # lolz + for z in zzz: + self.plot_z(Knot, z) + if out != "": + save_figure(out + self.figprefix + str(z)) + return + + """ Plot absolute power spectrum, not relative""" + + def plot_power(self, path, redshift, colour="black"): + (k_g, Pk_g) = self.loadpk( + path + self.suf + self.pre + self.GetSnap(redshift) + self.ext, self.box + ) + plt.loglog(k_g, Pk_g, color=colour) + plt.xlim(0.01, k_g[-1] * 1.1) + plt.ylabel("P(k) /(h-3 Mpc3)") + plt.xlabel("k /(h MPc-1)") + plt.title("Power spectrum at z=" + str(redshift)) + return (k_g, Pk_g) + + """ Plot absolute power for all redshifts """ + + def plot_power_all(self, Knot, zzz=np.array([]), out=""): + if np.size(zzz) == 0: + zzz = self.Zz # lolz + for z in zzz: + ind = wheref(self.Zz, z) + plt.figure(ind[0][0]) + for sim in Knot.names: + self.plot_power(sim, z) + if out != "": + save_figure(out + self.figprefix + str(z)) + return + + """ Compare two power spectra directly. Smooths result. plot_compare_two(first P(k), second P(k))""" - def plot_compare_two(self, one, onebox, two,twobox,colour=""): - (onek,oneP)=self.loadpk(one,onebox) - (twok,twoP)=self.loadpk(two,twobox) - onei = np.where(onek <= twok[-1]) - twoi= np.where (onek[onei] >= twok[0]) - relP=rebin(twoP, twok, onek[onei][twoi]) - relP=relP/rebin(oneP, onek, onek[onei][twoi]) - onek=onek[onei][twoi] - plt.title("Relative Power spectra "+one+" and "+two) - plt.ylabel(r"$P_2(k)/P_1(k)$") - plt.xlabel(r"$k\; (h\,\mathrm{Mpc}^{-1})$") - if colour == "": - line=plt.semilogx(onek,relP) - else: - line=plt.semilogx(onek,relP,color=colour) - plt.semilogx(self.knotpos,np.ones(len(self.knotpos)),"ro") - ind=np.where(onek < 10) - plt.ylim(min(relP[ind])*0.98,max(relP[ind])*1.01) - plt.xlim(onek[0]*0.8, 10) - return line - - """Get the difference between two simulations on scales probed by + + def plot_compare_two(self, one, onebox, two, twobox, colour=""): + (onek, oneP) = self.loadpk(one, onebox) + (twok, twoP) = self.loadpk(two, twobox) + onei = np.where(onek <= twok[-1]) + twoi = np.where(onek[onei] >= twok[0]) + relP = rebin(twoP, twok, onek[onei][twoi]) + relP = relP / rebin(oneP, onek, onek[onei][twoi]) + onek = onek[onei][twoi] + plt.title("Relative Power spectra " + one + " and " + two) + plt.ylabel(r"$P_2(k)/P_1(k)$") + plt.xlabel(r"$k\; (h\,\mathrm{Mpc}^{-1})$") + if colour == "": + line = plt.semilogx(onek, relP) + else: + line = plt.semilogx(onek, relP, color=colour) + plt.semilogx(self.knotpos, np.ones(len(self.knotpos)), "ro") + ind = np.where(onek < 10) + plt.ylim(min(relP[ind]) * 0.98, max(relP[ind]) * 1.01) + plt.xlim(onek[0] * 0.8, 10) + return line + + """Get the difference between two simulations on scales probed by the SDSS power spectrum""" - def compare_two(self, one, two,redshift): - (onek,oneP)=self.loadpk(one,self.bfbox) - (twok,twoP)=self.loadpk(two,self.box) - onei = np.where(onek <= twok[-1]) - twoi= np.where (onek[onei] >= twok[0]) - relP=rebin(twoP, twok, onek[onei][twoi]) - relP=relP/rebin(oneP, onek, onek[onei][twoi]) - onek=onek[onei][twoi] - sdss=self.GetSDSSkbins(redshift) - relP_r=np.ones(np.size(sdss)) - ind = np.where(sdss > onek[0]) - relP_r[ind]=rebin(relP,onek,sdss[ind]) - return relP_r - - """Do the above 12 times to get a correction table""" - def compare_two_table(self,onedir, twodir): - nk=np.size(self.sdsskbins) - nz=np.size(self.Zz)-1 - table=np.empty([nz,nk]) - for i in np.arange(0,nz): - sim=self.pre+self.GetSnap(self.Zz[i])+self.ext - table[-1-i,:]=self.compare_two(onedir+sim,twodir+sim,self.Zz[i]) - return table - - """ Plot a whole redshift range of relative power spectra on the same figure. + + def compare_two(self, one, two, redshift): + (onek, oneP) = self.loadpk(one, self.bfbox) + (twok, twoP) = self.loadpk(two, self.box) + onei = np.where(onek <= twok[-1]) + twoi = np.where(onek[onei] >= twok[0]) + relP = rebin(twoP, twok, onek[onei][twoi]) + relP = relP / rebin(oneP, onek, onek[onei][twoi]) + onek = onek[onei][twoi] + sdss = self.GetSDSSkbins(redshift) + relP_r = np.ones(np.size(sdss)) + ind = np.where(sdss > onek[0]) + relP_r[ind] = rebin(relP, onek, sdss[ind]) + return relP_r + + """Do the above 12 times to get a correction table""" + + def compare_two_table(self, onedir, twodir): + nk = np.size(self.sdsskbins) + nz = np.size(self.Zz) - 1 + table = np.empty([nz, nk]) + for i in np.arange(0, nz): + sim = self.pre + self.GetSnap(self.Zz[i]) + self.ext + table[-1 - i, :] = self.compare_two(onedir + sim, twodir + sim, self.Zz[i]) + return table + + """ Plot a whole redshift range of relative power spectra on the same figure. plot_all(onedir, twodir) Pass onedir and twodir as relative to basedir. ie, for default settings something like best-fit/flux-power/""" - def plot_compare_two_sdss(self, onedir, twodir, zzz=np.array([]), out="", title="", ylabel="", ymax=0,ymin=0, colour="",legend=False): - if np.size(zzz) == 0: - zzz=self.Zz #lolz - line=np.array([]) - legname=np.array([]) - sdss=self.sdsskbins - plt.xlabel(r"$k_v\; (s\,\mathrm{km}^{-1})$") - for z in zzz: - sim=self.pre+self.GetSnap(z)+self.ext - Pk=self.compare_two(onedir+sim,twodir+sim,z) - if colour == "": - line=np.append(line, plt.semilogx(sdss,Pk)) - else: - line=np.append(line, plt.semilogx(sdss,Pk,color=colour)) - legname=np.append(legname,"z="+str(z)) - plt.title(title) - if ylabel != "": - plt.ylabel(ylabel) - if legend: - plt.legend(line, legname,bbox_to_anchor=(0., 0, 1., .25), loc=3,ncol=3, mode="expand", borderaxespad=0.) - if ymax != 0 and ymin !=0: - plt.ylim(ymin,ymax) - plt.xlim(sdss[0],sdss[-1]) - plt.xticks(np.array([sdss[0],3e-3,5e-3,0.01,sdss[-1]]),("0.0014","0.003","0.005","0.01","0.0178")) - if out != "": - save_figure(out) - return plt.gcf() - - """ Plot a whole redshift range of relative power spectra on the same figure. + + def plot_compare_two_sdss( + self, + onedir, + twodir, + zzz=np.array([]), + out="", + title="", + ylabel="", + ymax=0, + ymin=0, + colour="", + legend=False, + ): + if np.size(zzz) == 0: + zzz = self.Zz # lolz + line = np.array([]) + legname = np.array([]) + sdss = self.sdsskbins + plt.xlabel(r"$k_v\; (s\,\mathrm{km}^{-1})$") + for z in zzz: + sim = self.pre + self.GetSnap(z) + self.ext + Pk = self.compare_two(onedir + sim, twodir + sim, z) + if colour == "": + line = np.append(line, plt.semilogx(sdss, Pk)) + else: + line = np.append(line, plt.semilogx(sdss, Pk, color=colour)) + legname = np.append(legname, "z=" + str(z)) + plt.title(title) + if ylabel != "": + plt.ylabel(ylabel) + if legend: + plt.legend( + line, + legname, + bbox_to_anchor=(0.0, 0, 1.0, 0.25), + loc=3, + ncol=3, + mode="expand", + borderaxespad=0.0, + ) + if ymax != 0 and ymin != 0: + plt.ylim(ymin, ymax) + plt.xlim(sdss[0], sdss[-1]) + plt.xticks( + np.array([sdss[0], 3e-3, 5e-3, 0.01, sdss[-1]]), + ("0.0014", "0.003", "0.005", "0.01", "0.0178"), + ) + if out != "": + save_figure(out) + return plt.gcf() + + """ Plot a whole redshift range of relative power spectra on the same figure. plot_all(onedir, twodir) Pass onedir and twodir as relative to basedir. ie, for default settings something like best-fit/flux-power/ onedir uses bfbox, twodir uses box""" - def plot_compare_two_all(self, onedir, twodir, zzz=np.array([]), out="", title="", ylabel="", ymax=0,ymin=0, colour="",legend=False): - if np.size(zzz) == 0: - zzz=self.Zz #lolz - line=np.array([]) - legname=np.array([]) - for z in zzz: - line=np.append(line, self.plot_compare_two(onedir+self.pre+self.GetSnap(z)+self.ext,self.bfbox,twodir+self.pre+self.GetSnap(z)+self.ext,self.box,colour)) - legname=np.append(legname,"z="+str(z)) - if title == "": - plt.title("Relative Power spectra "+onedir+" and "+twodir) - else: - plt.title(title) - if ylabel != "": - plt.ylabel(ylabel) - if legend: - plt.legend(line, legname,bbox_to_anchor=(0., 0, 1., .25), loc=3,ncol=3, mode="expand", borderaxespad=0.) - if ymax != 0 and ymin !=0: - plt.ylim(ymin,ymax) - if out != "": - save_figure(out) - return plt.gcf() - - """Get a power spectrum in the flat format we use""" - """for inputting some cosmomc tables""" - def GetFlat(self,dir, si=0.0): - Pk_sdss=np.empty([11, 12]) - #Note this omits the z=2.0 bin - #SiIII corr now done on the fly in lya_sdss_viel.f90 - for i in np.arange(0,np.size(self.Snaps)-1): - scale=self.Hubble(self.Zz[i])/(1.0+self.Zz[i]) - (k,Pk)=self.loadpk(dir+self.Snaps[i]+self.ext, self.box) - Fbar=math.exp(-0.0023*(1+self.Zz[i])**3.65) - a=si/(1-Fbar) - #The SiIII correction is kind of oscillatory, so we want - #to average over the whole interval being probed. - sdss=self.sdsskbins - kmids=np.zeros(np.size(sdss)+1) - sicorr=np.empty(np.size(sdss)) - for j in np.arange(0,np.size(sdss)-1): - kmids[j+1]=math.exp((math.log(sdss[j+1])+math.log(sdss[j]))/2.) - #Final segment should make no difference - kmids[-1]=2*math.pi/2271+kmids[-2] - #Near zero bad things will happen - sicorr[0]=1+a**2+2*a*math.cos(2271*sdss[0]) - for j in np.arange(1,np.size(sdss)): - sicorr[j]=1+a**2+2*a*(math.sin(2271*kmids[j+1])-math.sin(2271*kmids[j]))/(kmids[j+1]-kmids[j])/2271 - sdss=self.GetSDSSkbins(self.Zz[i]) - Pk_sdss[-i-1,:]=rebin(Pk, k, sdss)*scale*sicorr - return Pk_sdss - - """ Calculate the flux derivatives for a single redshift + + def plot_compare_two_all( + self, + onedir, + twodir, + zzz=np.array([]), + out="", + title="", + ylabel="", + ymax=0, + ymin=0, + colour="", + legend=False, + ): + if np.size(zzz) == 0: + zzz = self.Zz # lolz + line = np.array([]) + legname = np.array([]) + for z in zzz: + line = np.append( + line, + self.plot_compare_two( + onedir + self.pre + self.GetSnap(z) + self.ext, + self.bfbox, + twodir + self.pre + self.GetSnap(z) + self.ext, + self.box, + colour, + ), + ) + legname = np.append(legname, "z=" + str(z)) + if title == "": + plt.title("Relative Power spectra " + onedir + " and " + twodir) + else: + plt.title(title) + if ylabel != "": + plt.ylabel(ylabel) + if legend: + plt.legend( + line, + legname, + bbox_to_anchor=(0.0, 0, 1.0, 0.25), + loc=3, + ncol=3, + mode="expand", + borderaxespad=0.0, + ) + if ymax != 0 and ymin != 0: + plt.ylim(ymin, ymax) + if out != "": + save_figure(out) + return plt.gcf() + + """Get a power spectrum in the flat format we use""" + """for inputting some cosmomc tables""" + + def GetFlat(self, dir, si=0.0): + Pk_sdss = np.empty([11, 12]) + # Note this omits the z=2.0 bin + # SiIII corr now done on the fly in lya_sdss_viel.f90 + for i in np.arange(0, np.size(self.Snaps) - 1): + scale = self.Hubble(self.Zz[i]) / (1.0 + self.Zz[i]) + (k, Pk) = self.loadpk(dir + self.Snaps[i] + self.ext, self.box) + Fbar = math.exp(-0.0023 * (1 + self.Zz[i]) ** 3.65) + a = si / (1 - Fbar) + # The SiIII correction is kind of oscillatory, so we want + # to average over the whole interval being probed. + sdss = self.sdsskbins + kmids = np.zeros(np.size(sdss) + 1) + sicorr = np.empty(np.size(sdss)) + for j in np.arange(0, np.size(sdss) - 1): + kmids[j + 1] = math.exp( + (math.log(sdss[j + 1]) + math.log(sdss[j])) / 2.0 + ) + # Final segment should make no difference + kmids[-1] = 2 * math.pi / 2271 + kmids[-2] + # Near zero bad things will happen + sicorr[0] = 1 + a ** 2 + 2 * a * math.cos(2271 * sdss[0]) + for j in np.arange(1, np.size(sdss)): + sicorr[j] = ( + 1 + + a ** 2 + + 2 + * a + * (math.sin(2271 * kmids[j + 1]) - math.sin(2271 * kmids[j])) + / (kmids[j + 1] - kmids[j]) + / 2271 + ) + sdss = self.GetSDSSkbins(self.Zz[i]) + Pk_sdss[-i - 1, :] = rebin(Pk, k, sdss) * scale * sicorr + return Pk_sdss + + """ Calculate the flux derivatives for a single redshift Output: (kbins d2P...kbins dP (flat vector of length 2xkbins))""" - def calc_z(self, redshift,s_knot, kbins): - #Array to store answers. - #Format is: k x (dP, d²P, χ²) - kbins=np.array(kbins) - nk=np.size(kbins) - snap=self.GetSnap(redshift) - results=np.zeros(2*nk) - if np.size(s_knot.qvals) > 0: - results = np.zeros(4*nk) - pdifs=s_knot.pvals-s_knot.p0 - qdifs=np.array([]) - #If we have somethign where the parameter is redshift-dependent, eg, gamma. - if np.size(s_knot.p0) > 1: - i=wheref(redshift, self.Zz) - pdifs = s_knot.pvals[:,i]-s_knot.p0[i] - if np.size(s_knot.qvals) > 0: - qdifs=s_knot.qvals[:,i] - s_knot.q0[i] - npvals=np.size(pdifs) - #Load the data - (k,PFp0)=self.loadpk(s_knot.bstft+self.suf+snap+self.ext,s_knot.bfbox) - #This is to rescale by the mean flux, for generating mean flux tables. - ### - #tau_eff=0.0023*(1+redshift)**3.65 - #tmin=0.2*((1+redshift)/4.)**2 - #tmax=0.5*((1+redshift)/4.)**4 - #teffs=tmin+s_knot.pvals*(tmax-tmin)/30. - #pdifs=teffs/tau_eff-1. - ### - PowerFluxes=np.zeros((npvals,np.size(k))) - for i in np.arange(0,np.size(s_knot.names)): - (k,PowerFluxes[i,:])=self.loadpk(s_knot.names[i]+self.suf+snap+self.ext, s_knot.bfbox) - #So now we have an array of data values, which we want to rebin. - ind = np.where(kbins >= k[0]) - difPF_rebin=np.ones((npvals,np.size(kbins))) - for i in np.arange(0, npvals): - difPF_rebin[i,ind]=rebin(PowerFluxes[i,:]/PFp0,k,kbins[ind]) - #Set the things beyond the range of the interpolator - #equal to the final value. - if ind[0][0] > 0: - difPF_rebin[i,0:ind[0][0]]=difPF_rebin[i,ind[0][0]] - #So now we have an array of data values. - #Pass each k value to flux_deriv in turn. - for k in np.arange(0,np.size(kbins)): - #Format of returned data is: - # y = ax**2 + bx + cz**2 +dz +e xz - # derives = (a,b,c,d,e) - derivs=self.flux_deriv(difPF_rebin[:,k], pdifs,qdifs) - results[k]=derivs[0] - results[nk+k]=derivs[1] - if np.size(derivs) > 2: - results[2*nk+k]=derivs[2] - results[3*nk+k]=derivs[3] - return results - - """ Calculate the flux derivatives for all redshifts + + def calc_z(self, redshift, s_knot, kbins): + # Array to store answers. + # Format is: k x (dP, d²P, χ²) + kbins = np.array(kbins) + nk = np.size(kbins) + snap = self.GetSnap(redshift) + results = np.zeros(2 * nk) + if np.size(s_knot.qvals) > 0: + results = np.zeros(4 * nk) + pdifs = s_knot.pvals - s_knot.p0 + qdifs = np.array([]) + # If we have somethign where the parameter is redshift-dependent, eg, gamma. + if np.size(s_knot.p0) > 1: + i = wheref(redshift, self.Zz) + pdifs = s_knot.pvals[:, i] - s_knot.p0[i] + if np.size(s_knot.qvals) > 0: + qdifs = s_knot.qvals[:, i] - s_knot.q0[i] + npvals = np.size(pdifs) + # Load the data + (k, PFp0) = self.loadpk(s_knot.bstft + self.suf + snap + self.ext, s_knot.bfbox) + # This is to rescale by the mean flux, for generating mean flux tables. + ### + # tau_eff=0.0023*(1+redshift)**3.65 + # tmin=0.2*((1+redshift)/4.)**2 + # tmax=0.5*((1+redshift)/4.)**4 + # teffs=tmin+s_knot.pvals*(tmax-tmin)/30. + # pdifs=teffs/tau_eff-1. + ### + PowerFluxes = np.zeros((npvals, np.size(k))) + for i in np.arange(0, np.size(s_knot.names)): + (k, PowerFluxes[i, :]) = self.loadpk( + s_knot.names[i] + self.suf + snap + self.ext, s_knot.bfbox + ) + # So now we have an array of data values, which we want to rebin. + ind = np.where(kbins >= k[0]) + difPF_rebin = np.ones((npvals, np.size(kbins))) + for i in np.arange(0, npvals): + difPF_rebin[i, ind] = rebin(PowerFluxes[i, :] / PFp0, k, kbins[ind]) + # Set the things beyond the range of the interpolator + # equal to the final value. + if ind[0][0] > 0: + difPF_rebin[i, 0 : ind[0][0]] = difPF_rebin[i, ind[0][0]] + # So now we have an array of data values. + # Pass each k value to flux_deriv in turn. + for k in np.arange(0, np.size(kbins)): + # Format of returned data is: + # y = ax**2 + bx + cz**2 +dz +e xz + # derives = (a,b,c,d,e) + derivs = self.flux_deriv(difPF_rebin[:, k], pdifs, qdifs) + results[k] = derivs[0] + results[nk + k] = derivs[1] + if np.size(derivs) > 2: + results[2 * nk + k] = derivs[2] + results[3 * nk + k] = derivs[3] + return results + + """ Calculate the flux derivatives for all redshifts Input: Sims to load, parameter values, mean parameter value Output: (2*kbins) x (zbins)""" - def calc_all(self, s_knot,kbins): - flux_derivatives=np.zeros((2*np.size(kbins),np.size(self.Zz))) - if np.size(s_knot.qvals) > 1: - flux_derivatives=np.zeros((4*np.size(kbins),np.size(self.Zz))) - #Call flux_deriv_const_z for each redshift. - for i in np.arange(0,np.size(self.Zz)): - flux_derivatives[:,i]=self.calc_z(self.Zz[i], s_knot,kbins) - return flux_derivatives - - """Calculate the flux-derivative for a single redshift and k bin""" - def flux_deriv(self, PFdif, pdif, qdif=np.array([])): - pdif=np.ravel(pdif) - if np.size(pdif) != np.size(PFdif): - raise DataError(str(np.size(pdif))+" parameter values, but "+str(np.size(PFdif))+" P_F values") - if np.size(pdif) < 2: - raise DataError(str(np.size(pdif))+" pvals given. Need at least 2.") - PFdif=PFdif-1.0 - if np.size(qdif) > 2: - qdif=np.ravel(qdif) - mat=np.vstack([pdif**2, pdif, qdif**2, qdif] ).T - else: - mat=np.vstack([pdif**2, pdif] ).T - (derivs, residues,rank, sing)=np.linalg.lstsq(mat, PFdif) - return derivs - - """ Get the error on one test simulation at single redshift """ - def Get_Error_z(self, Sim, bstft,box, derivs, params, redshift,qarams=np.empty([])): - #Need to load and rebin the sim. - (k, test)=self.loadpk(Sim+self.suf+self.GetSnap(redshift)+self.ext, box) - (k,bf)=self.loadpk(bstft+self.suf+self.GetSnap(redshift)+self.ext,box) - kbins=self.Getkbins() - ind = np.where(kbins >= 1.0*2.0*math.pi*self.H0/box) - test2=np.ones(np.size(kbins)) - test2[ind]=rebin(test/bf,k,kbins[ind])# - if ind[0][0] > 0: - test2[0:ind[0][0]]=test2[ind[0][0]] - if np.size(qarams) > 0: - guess=derivs.GetPF(params, redshift,qarams)+1.0 - else: - guess=derivs.GetPF(params, redshift)+1.0 - return np.array((test2/guess))[0][0] + + def calc_all(self, s_knot, kbins): + flux_derivatives = np.zeros((2 * np.size(kbins), np.size(self.Zz))) + if np.size(s_knot.qvals) > 1: + flux_derivatives = np.zeros((4 * np.size(kbins), np.size(self.Zz))) + # Call flux_deriv_const_z for each redshift. + for i in np.arange(0, np.size(self.Zz)): + flux_derivatives[:, i] = self.calc_z(self.Zz[i], s_knot, kbins) + return flux_derivatives + + """Calculate the flux-derivative for a single redshift and k bin""" + + def flux_deriv(self, PFdif, pdif, qdif=np.array([])): + pdif = np.ravel(pdif) + if np.size(pdif) != np.size(PFdif): + raise DataError( + str(np.size(pdif)) + + " parameter values, but " + + str(np.size(PFdif)) + + " P_F values" + ) + if np.size(pdif) < 2: + raise DataError(str(np.size(pdif)) + " pvals given. Need at least 2.") + PFdif = PFdif - 1.0 + if np.size(qdif) > 2: + qdif = np.ravel(qdif) + mat = np.vstack([pdif ** 2, pdif, qdif ** 2, qdif]).T + else: + mat = np.vstack([pdif ** 2, pdif]).T + (derivs, residues, rank, sing) = np.linalg.lstsq(mat, PFdif) + return derivs + + """ Get the error on one test simulation at single redshift """ + + def Get_Error_z( + self, Sim, bstft, box, derivs, params, redshift, qarams=np.empty([]) + ): + # Need to load and rebin the sim. + (k, test) = self.loadpk(Sim + self.suf + self.GetSnap(redshift) + self.ext, box) + (k, bf) = self.loadpk(bstft + self.suf + self.GetSnap(redshift) + self.ext, box) + kbins = self.Getkbins() + ind = np.where(kbins >= 1.0 * 2.0 * math.pi * self.H0 / box) + test2 = np.ones(np.size(kbins)) + test2[ind] = rebin(test / bf, k, kbins[ind]) # + if ind[0][0] > 0: + test2[0 : ind[0][0]] = test2[ind[0][0]] + if np.size(qarams) > 0: + guess = derivs.GetPF(params, redshift, qarams) + 1.0 + else: + guess = derivs.GetPF(params, redshift) + 1.0 + return np.array((test2 / guess))[0][0] + """ A class written to store the various methods related to calculating of the flux derivatives and plotting of the flux power spectra""" + + class flux_pow(power_spec): - figprefix="/flux-figure" - kbins=np.array([]) - def __init__(self, Snaps=("snapshot_000", "snapshot_001","snapshot_002","snapshot_003","snapshot_004","snapshot_005","snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"), - Zz=np.array([4.2,4.0,3.8,3.6,3.4,3.2,3.0,2.8,2.6,2.4,2.2,2.0]), - sdsskbins=np.array([0.00141,0.00178,0.00224,0.00282,0.00355,0.00447,0.00562,0.00708,0.00891,0.01122,0.01413,0.01778]), - knotpos=np.array([0.07,0.15,0.475, 0.75, 1.19, 1.89,4,25]), om=0.266, H0=0.71,box=60.0,kmax=4.0, - base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/",bf="best-fit/",suf="flux-power/", ext="_flux_power.txt"): - power_spec.__init__(self, Snaps,Zz,sdsskbins,knotpos, om, H0,box,base,suf, ext) - (k_bf,Pk_bf)= self.loadpk(bf+suf+"snapshot_000"+self.ext,self.bfbox) - ind=np.where(k_bf <= kmax) - self.kbins=k_bf[ind] - - def plot_z(self,Sims,redshift,title="Relative Flux Power",ylabel=r"$\mathrm{P}_\mathrm{F}(k,p)\,/\,\mathrm{P}_\mathrm{F}(k,p_0)$", legend=True): - power_spec.plot_z(self,Sims,redshift,title,ylabel,legend) - if legend: - kbins=self.GetSDSSkbins(redshift) - plt.axvspan(kbins[0], kbins[-1], color="#B0B0B0") - plt.ylim(self.ymin,self.ymax) - plt.xlim(self.kbins[0]*0.8, 10) - - """Get the kbins to interpolate onto""" - def Getkbins(self): - return self.kbins - - """Load the SDSS power spectrum""" - def MacDonaldPF(self,sdss, zz): - psdss=sdss[np.where(sdss[:,0] == zz)][:,1:3] - fbar=math.exp(-0.0023*(1+zz)**3.65) - #multiply by the hubble parameter to be in 1/(km/s) - scale=self.Hubble(zz)/(1.0+zz) - PF=psdss[:,1]*fbar**2/scale - k=psdss[:,0]*scale - return (k, PF) - - """Load a Pk. Different function due to needing to be different for each class""" - def loadpk(self, path,box): - #Adjust Fourier convention. - flux_power=np.loadtxt(self.base+path) - scale=self.H0/box - k=(flux_power[1:,0]-0.5)*scale*2.0*math.pi - PF=flux_power[1:,1]/scale - return (k, PF) + figprefix = "/flux-figure" + kbins = np.array([]) + + def __init__( + self, + Snaps=( + "snapshot_000", + "snapshot_001", + "snapshot_002", + "snapshot_003", + "snapshot_004", + "snapshot_005", + "snapshot_006", + "snapshot_007", + "snapshot_008", + "snapshot_009", + "snapshot_010", + "snapshot_011", + ), + Zz=np.array([4.2, 4.0, 3.8, 3.6, 3.4, 3.2, 3.0, 2.8, 2.6, 2.4, 2.2, 2.0]), + sdsskbins=np.array( + [ + 0.00141, + 0.00178, + 0.00224, + 0.00282, + 0.00355, + 0.00447, + 0.00562, + 0.00708, + 0.00891, + 0.01122, + 0.01413, + 0.01778, + ] + ), + knotpos=np.array([0.07, 0.15, 0.475, 0.75, 1.19, 1.89, 4, 25]), + om=0.266, + H0=0.71, + box=60.0, + kmax=4.0, + base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/", + bf="best-fit/", + suf="flux-power/", + ext="_flux_power.txt", + ): + power_spec.__init__( + self, Snaps, Zz, sdsskbins, knotpos, om, H0, box, base, suf, ext + ) + (k_bf, Pk_bf) = self.loadpk(bf + suf + "snapshot_000" + self.ext, self.bfbox) + ind = np.where(k_bf <= kmax) + self.kbins = k_bf[ind] + + def plot_z( + self, + Sims, + redshift, + title="Relative Flux Power", + ylabel=r"$\mathrm{P}_\mathrm{F}(k,p)\,/\,\mathrm{P}_\mathrm{F}(k,p_0)$", + legend=True, + ): + power_spec.plot_z(self, Sims, redshift, title, ylabel, legend) + if legend: + kbins = self.GetSDSSkbins(redshift) + plt.axvspan(kbins[0], kbins[-1], color="#B0B0B0") + plt.ylim(self.ymin, self.ymax) + plt.xlim(self.kbins[0] * 0.8, 10) + + """Get the kbins to interpolate onto""" + + def Getkbins(self): + return self.kbins + + """Load the SDSS power spectrum""" + + def MacDonaldPF(self, sdss, zz): + psdss = sdss[np.where(sdss[:, 0] == zz)][:, 1:3] + fbar = math.exp(-0.0023 * (1 + zz) ** 3.65) + # multiply by the hubble parameter to be in 1/(km/s) + scale = self.Hubble(zz) / (1.0 + zz) + PF = psdss[:, 1] * fbar ** 2 / scale + k = psdss[:, 0] * scale + return (k, PF) + + """Load a Pk. Different function due to needing to be different for each class""" + + def loadpk(self, path, box): + # Adjust Fourier convention. + flux_power = np.loadtxt(self.base + path) + scale = self.H0 / box + k = (flux_power[1:, 0] - 0.5) * scale * 2.0 * math.pi + PF = flux_power[1:, 1] / scale + return (k, PF) """ A class to plot matter power spectra """ + + class matter_pow(power_spec): - ob=0.0 - #For plotting - ymin=0.4 - ymax=1.6 - figprefix="/matter-figure" - def __init__(self, Snaps=("snapshot_000", "snapshot_001","snapshot_002","snapshot_003","snapshot_004","snapshot_005","snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"), - Zz=np.array([4.2,4.0,3.8,3.6,3.4,3.2,3.0,2.8,2.6,2.4,2.2,2.0]), - sdsskbins=np.array([0.00141,0.00178,0.00224,0.00282,0.00355,0.00447,0.00562,0.00708,0.00891,0.01122,0.01413,0.01778]), - knotpos=np.array([0.07,0.15,0.475, 0.75, 1.19, 1.89,4,25]), om=0.266,ob=0.0449, H0=0.71,box=60.0, - base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/",suf="matter-power/", ext=".0", matpre="PK-by-"): - power_spec.__init__(self, Snaps,Zz,sdsskbins,knotpos, om, H0,box,base,suf,ext) - self.ob=ob - self.pre=matpre - - def plot_z(self,Sims,redshift,title="Relative Matter Power",ylabel=r"$\mathrm{P}(k,p)\,/\,\mathrm{P}(k,p_0)$"): - power_spec.plot_z(self,Sims,redshift,title,ylabel) - - """Load a Pk. Different function due to needing to be different for each class""" - def loadpk(self, path,box): - #Load baryon P(k) - matter_power=np.loadtxt(self.base+path) - scale=self.H0/box - #Adjust Fourier convention. - simk=matter_power[1:,0]*scale*2.0*math.pi - Pkbar=matter_power[1:,1]/scale**3 - #Load DM P(k) - matter_power=np.loadtxt(self.base+re.sub("by","DM",path)) - PkDM=matter_power[1:,1]/scale**3 - Pk=(Pkbar*self.ob+PkDM*(self.om-self.ob))/self.om - return (simk,Pk) - - """ Plot absolute power spectrum, not relative""" - def plot_power(self,path, redshift,camb_filename=""): - (k_g, Pk_g)=power_spec.plot_power(self,path,redshift) - sigma=2.0 - pkg=np.loadtxt(self.base+path+self.suf+self.pre+self.GetSnap(redshift)+self.ext) - samp_err=pkg[1:,2] - sqrt_err=np.array(np.sqrt(samp_err)) - plt.loglog(k_g,Pk_g*(1+sigma*(2.0/sqrt_err+1.0/samp_err)),linestyle="-.",color="black") - plt.loglog(k_g,Pk_g*(1-sigma*(2.0/sqrt_err+1.0/samp_err)),linestyle="-.",color="black") - if camb_filename != "": - camb=np.loadtxt(camb_filename) - #Adjust Fourier convention. - k=camb[:,0]*self.H0 - #NOW THERE IS NO h in the T anywhere. - Pk=camb[:,1] - plt.loglog(k/self.H0, Pk, linestyle="--") - plt.xlim(0.01,k_g[-1]*1.1) - return(k_g, Pk_g) + ob = 0.0 + # For plotting + ymin = 0.4 + ymax = 1.6 + figprefix = "/matter-figure" + + def __init__( + self, + Snaps=( + "snapshot_000", + "snapshot_001", + "snapshot_002", + "snapshot_003", + "snapshot_004", + "snapshot_005", + "snapshot_006", + "snapshot_007", + "snapshot_008", + "snapshot_009", + "snapshot_010", + "snapshot_011", + ), + Zz=np.array([4.2, 4.0, 3.8, 3.6, 3.4, 3.2, 3.0, 2.8, 2.6, 2.4, 2.2, 2.0]), + sdsskbins=np.array( + [ + 0.00141, + 0.00178, + 0.00224, + 0.00282, + 0.00355, + 0.00447, + 0.00562, + 0.00708, + 0.00891, + 0.01122, + 0.01413, + 0.01778, + ] + ), + knotpos=np.array([0.07, 0.15, 0.475, 0.75, 1.19, 1.89, 4, 25]), + om=0.266, + ob=0.0449, + H0=0.71, + box=60.0, + base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/", + suf="matter-power/", + ext=".0", + matpre="PK-by-", + ): + power_spec.__init__( + self, Snaps, Zz, sdsskbins, knotpos, om, H0, box, base, suf, ext + ) + self.ob = ob + self.pre = matpre + + def plot_z( + self, + Sims, + redshift, + title="Relative Matter Power", + ylabel=r"$\mathrm{P}(k,p)\,/\,\mathrm{P}(k,p_0)$", + ): + power_spec.plot_z(self, Sims, redshift, title, ylabel) + + """Load a Pk. Different function due to needing to be different for each class""" + + def loadpk(self, path, box): + # Load baryon P(k) + matter_power = np.loadtxt(self.base + path) + scale = self.H0 / box + # Adjust Fourier convention. + simk = matter_power[1:, 0] * scale * 2.0 * math.pi + Pkbar = matter_power[1:, 1] / scale ** 3 + # Load DM P(k) + matter_power = np.loadtxt(self.base + re.sub("by", "DM", path)) + PkDM = matter_power[1:, 1] / scale ** 3 + Pk = (Pkbar * self.ob + PkDM * (self.om - self.ob)) / self.om + return (simk, Pk) + + """ Plot absolute power spectrum, not relative""" + + def plot_power(self, path, redshift, camb_filename=""): + (k_g, Pk_g) = power_spec.plot_power(self, path, redshift) + sigma = 2.0 + pkg = np.loadtxt( + self.base + path + self.suf + self.pre + self.GetSnap(redshift) + self.ext + ) + samp_err = pkg[1:, 2] + sqrt_err = np.array(np.sqrt(samp_err)) + plt.loglog( + k_g, + Pk_g * (1 + sigma * (2.0 / sqrt_err + 1.0 / samp_err)), + linestyle="-.", + color="black", + ) + plt.loglog( + k_g, + Pk_g * (1 - sigma * (2.0 / sqrt_err + 1.0 / samp_err)), + linestyle="-.", + color="black", + ) + if camb_filename != "": + camb = np.loadtxt(camb_filename) + # Adjust Fourier convention. + k = camb[:, 0] * self.H0 + # NOW THERE IS NO h in the T anywhere. + Pk = camb[:, 1] + plt.loglog(k / self.H0, Pk, linestyle="--") + plt.xlim(0.01, k_g[-1] * 1.1) + return (k_g, Pk_g) + """The PDF is an instance of the power_spec class. Perhaps poor naming""" + + class flux_pdf(power_spec): - def __init__(self, Snaps=("snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"), - Zz=np.array([3.0,2.8,2.6,2.4,2.2,2.0]), - sdsskbins=np.arange(0,20), - knotpos=np.array([]), om=0.266,ob=0.0449, H0=0.71,box=48.0, - base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/",suf="flux-pdf/", ext="_flux_pdf.txt",): - power_spec.__init__(self, Snaps,Zz,sdsskbins,knotpos, om, H0,box,base,suf,ext) - - def loadpk(self, path, box): - flux_pdf = np.loadtxt(self.base+path) - return(flux_pdf[:,0], flux_pdf[:,1]) - - """ Compare two power spectra directly. Smooths result. + def __init__( + self, + Snaps=( + "snapshot_006", + "snapshot_007", + "snapshot_008", + "snapshot_009", + "snapshot_010", + "snapshot_011", + ), + Zz=np.array([3.0, 2.8, 2.6, 2.4, 2.2, 2.0]), + sdsskbins=np.arange(0, 20), + knotpos=np.array([]), + om=0.266, + ob=0.0449, + H0=0.71, + box=48.0, + base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/", + suf="flux-pdf/", + ext="_flux_pdf.txt", + ): + power_spec.__init__( + self, Snaps, Zz, sdsskbins, knotpos, om, H0, box, base, suf, ext + ) + + def loadpk(self, path, box): + flux_pdf = np.loadtxt(self.base + path) + return (flux_pdf[:, 0], flux_pdf[:, 1]) + + """ Compare two power spectra directly. Smooths result. plot_compare_two(first P(k), second P(k))""" - def plot_compare_two(self, one, onebox, two,twobox,colour=""): - (onek,oneP)=self.loadpk(one,onebox) - (twok,twoP)=self.loadpk(two,twobox) - relP=oneP/twoP - plt.title("Relative flux PDF "+one+" and "+two) - plt.ylabel(r"$F_2(k)/F_1(k)$") - plt.xlabel(r"$Flux$") - line=plt.semilogy(onek,relP) - return line - - """ Plot absolute power spectrum, not relative""" - def plot_power(self,path, redshift): - (k,Pdf)=self.loadpk(path+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.box) - plt.semilogy(k,Pdf, color="black", linewidth="1.5") - plt.ylabel("P(k) /(h-3 Mpc3)") - plt.xlabel("k /(h MPc-1)") - plt.title("PDF at z="+str(redshift)) - return(k, Pdf) - - """ Calculate the flux derivatives for a single redshift + + def plot_compare_two(self, one, onebox, two, twobox, colour=""): + (onek, oneP) = self.loadpk(one, onebox) + (twok, twoP) = self.loadpk(two, twobox) + relP = oneP / twoP + plt.title("Relative flux PDF " + one + " and " + two) + plt.ylabel(r"$F_2(k)/F_1(k)$") + plt.xlabel(r"$Flux$") + line = plt.semilogy(onek, relP) + return line + + """ Plot absolute power spectrum, not relative""" + + def plot_power(self, path, redshift): + (k, Pdf) = self.loadpk( + path + self.suf + self.pre + self.GetSnap(redshift) + self.ext, self.box + ) + plt.semilogy(k, Pdf, color="black", linewidth="1.5") + plt.ylabel("P(k) /(h-3 Mpc3)") + plt.xlabel("k /(h MPc-1)") + plt.title("PDF at z=" + str(redshift)) + return (k, Pdf) + + """ Calculate the flux derivatives for a single redshift Output: (kbins d2P...kbins dP (flat vector of length 2x21))""" - def calc_z(self, redshift,s_knot): - #Array to store answers. - #Format is: k x (dP, d²P, χ²) - npvals=np.size(s_knot.pvals) - nk=21 - results=np.zeros(2*nk) - pdifs=s_knot.pvals-s_knot.p0 - #This is to rescale by the mean flux, for generating mean flux tables. - ### - #tau_eff=0.0023*(1+redshift)**3.65 - #tmin=0.2*((1+redshift)/4.)**2 - #tmax=0.5*((1+redshift)/4.)**4 - #teffs=tmin+s_knot.pvals*(tmax-tmin)/30. - #pdifs=teffs/tau_eff-1. - ### - ured=np.ceil(redshift*5)/5. - lred=np.floor(redshift*5)/5. - usnap=self.GetSnap(ured) - lsnap=self.GetSnap(lred) - #Load the data - (k,uPFp0)=self.loadpk(s_knot.bstft+self.suf+usnap+self.ext,s_knot.bfbox) - uPower=np.zeros((npvals,np.size(k))) - for i in np.arange(0,np.size(s_knot.names)): - (k,uPower[i,:])=self.loadpk(s_knot.names[i]+self.suf+usnap+self.ext, s_knot.bfbox) - (k,lPFp0)=self.loadpk(s_knot.bstft+self.suf+lsnap+self.ext,s_knot.bfbox) - lPower=np.zeros((npvals,np.size(k))) - for i in np.arange(0,np.size(s_knot.names)): - (k,lPower[i,:])=self.loadpk(s_knot.names[i]+self.suf+lsnap+self.ext, s_knot.bfbox) - PowerFluxes=5*((redshift-lred)*uPower+(ured-redshift)*lPower) - PFp0=5*((redshift-lred)*uPFp0+(ured-redshift)*lPFp0) - #So now we have an array of data values. - #Pass each k value to flux_deriv in turn. - for k in np.arange(0,nk): - (dPF, d2PF,chi2)=self.flux_deriv(PowerFluxes[:,k]/PFp0[k], pdifs) - results[k]=d2PF - results[nk+k]=dPF - return results - - """Get the kbins to interpolate onto""" - def Getkbins(self): - return np.arange(0,20,1)+0.5 - """ Plot comparisons between a bunch of sims on one graph + + def calc_z(self, redshift, s_knot): + # Array to store answers. + # Format is: k x (dP, d²P, χ²) + npvals = np.size(s_knot.pvals) + nk = 21 + results = np.zeros(2 * nk) + pdifs = s_knot.pvals - s_knot.p0 + # This is to rescale by the mean flux, for generating mean flux tables. + ### + # tau_eff=0.0023*(1+redshift)**3.65 + # tmin=0.2*((1+redshift)/4.)**2 + # tmax=0.5*((1+redshift)/4.)**4 + # teffs=tmin+s_knot.pvals*(tmax-tmin)/30. + # pdifs=teffs/tau_eff-1. + ### + ured = np.ceil(redshift * 5) / 5.0 + lred = np.floor(redshift * 5) / 5.0 + usnap = self.GetSnap(ured) + lsnap = self.GetSnap(lred) + # Load the data + (k, uPFp0) = self.loadpk( + s_knot.bstft + self.suf + usnap + self.ext, s_knot.bfbox + ) + uPower = np.zeros((npvals, np.size(k))) + for i in np.arange(0, np.size(s_knot.names)): + (k, uPower[i, :]) = self.loadpk( + s_knot.names[i] + self.suf + usnap + self.ext, s_knot.bfbox + ) + (k, lPFp0) = self.loadpk( + s_knot.bstft + self.suf + lsnap + self.ext, s_knot.bfbox + ) + lPower = np.zeros((npvals, np.size(k))) + for i in np.arange(0, np.size(s_knot.names)): + (k, lPower[i, :]) = self.loadpk( + s_knot.names[i] + self.suf + lsnap + self.ext, s_knot.bfbox + ) + PowerFluxes = 5 * ((redshift - lred) * uPower + (ured - redshift) * lPower) + PFp0 = 5 * ((redshift - lred) * uPFp0 + (ured - redshift) * lPFp0) + # So now we have an array of data values. + # Pass each k value to flux_deriv in turn. + for k in np.arange(0, nk): + (dPF, d2PF, chi2) = self.flux_deriv(PowerFluxes[:, k] / PFp0[k], pdifs) + results[k] = d2PF + results[nk + k] = dPF + return results + + """Get the kbins to interpolate onto""" + + def Getkbins(self): + return np.arange(0, 20, 1) + 0.5 + + """ Plot comparisons between a bunch of sims on one graph plot_z(Redshift, Sims to use ( eg, A1.14). Note this will clear current figures.""" - def plot_z(self,Knot,redshift,title="",ylabel=""): - #Load best-fit - (simk,BFPk)=self.loadpk(Knot.bstft+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.bfbox) - #Setup figure plot. - ind=wheref(self.Zz, redshift) - plt.figure(ind[0][0]) - plt.clf() - if title != '': - plt.title(title+" at z="+str(redshift),) - plt.ylabel(ylabel) - plt.xlabel(r"$\mathcal{F}$") - line=np.array([]) - legname=np.array([]) - for sim in Knot.names: - (k,Pk)=self.loadpk(sim+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.box) - line=np.append(line, plt.semilogy(simk,Pk/BFPk,linestyle="-", linewidth=1.5)) - legname=np.append(legname,sim) - plt.legend(line,legname) - return - - """Get a power spectrum in the flat format we use""" - """for inputting some cosmomc tables""" - def GetFlat(self,dir): - Pk_sdss=np.empty([11, 12]) - #For z=2.07 we need to average snap_011 and snap_010 - z=2.07 - (k,PF_a)=self.loadpk(dir+self.suf+"snapshot_011"+self.ext, self.box) - (k,PF_b)=self.loadpk(dir+self.suf+"snapshot_010"+self.ext, self.box) - PF1=(z-2.0)*5*(PF_b-PF_a)+PF_a - z=2.52 - (k,PF_a)=self.loadpk(dir+self.suf+"snapshot_009"+self.ext, self.box) - (k,PF_b)=self.loadpk(dir+self.suf+"snapshot_008"+self.ext, self.box) - PF2=(z-2.4)*5*(PF_b-PF_a)+PF_a - z=2.94 - (k,PF_a)=self.loadpk(dir+self.suf+"snapshot_007"+self.ext, self.box) - (k,PF_b)=self.loadpk(dir+self.suf+"snapshot_006"+self.ext, self.box) - PF3=(z-2.8)*5*(PF_b-PF_a)+PF_a - PDF = np.array([PF1,PF2,PF3]) - np.savetxt(sys.stdout,PDF.T("%1.8f","%1.8f","%1.8f")) - return (PF1, PF2, PF3) - -if __name__=='__main__': - flux=flux_pow() - matter=matter_pow() - fpdf=flux_pdf() - A_knot=knot(("A0.54/","A0.74/","A0.84/","A1.04/", "A1.14/","A1.34/"), (0.54,0.74,0.84,1.04,1.14,1.34),0.94,"best-fit/", 60) - AA_knot=knot(("AA0.54/","AA0.74/","AA1.14/","AA1.34/"), (0.54,0.74,1.14,1.34),0.94,"boxcorr400/", 120) - B_knot=knot(("B0.33/","B0.53/","B0.73/", "B0.83/", "B1.03/","B1.13/", "B1.33/"), (0.33,0.53,0.73,0.83,1.03, 1.13,1.33),0.93,"best-fit/", 60) - C_knot=knot(("C0.11/", "C0.31/","C0.51/","C0.71/","C1.11/","C1.31/","C1.51/"),(0.11, 0.31,0.51,0.71, 1.11,1.31,1.51),0.91,"bf2/", 60) - D_knot=knot(("D0.50/","D0.70/","D1.10/","D1.20/", "D1.30/", "D1.50/", "D1.70/"),(0.50, 0.70,1.10,1.20,1.30, 1.50, 1.70),0.90,"bfD/", 48) - interp=flux_interp(flux, (AA_knot, B_knot, C_knot, D_knot)) -#Thermal parameters for models - -#Models where gamma is varied - bf2zg=np.array([1.5704182,1.5754168,1.5797292,1.5836439,1.5870027,1.5915027,1.5943816,1.5986097,1.6027860,1.6073484,1.6116327,1.6158375]) - bf2zt=np.array([21221.521,21915.172,22303.908,22432.774,22889.881,22631.245,22610.916,22458.702,22020.437,21752.465,21278.358,20720.927])/1e3 - bf2g=np.array([1.4260277,1.4332204,1.4354097,1.4422902,1.4455976,1.4502961,1.4547729,1.4593188,1.4637958,1.4682989,1.4728728,1.4769461 ]) - bf2t=np.array([21453.313,21917.668,22655.974, 22600.650, 23061.612, 22798.62,22694.608,22721.762,22208.044, 21826.779, 21552.183,20908.159])/1e3 - bf2bg=np.array([1.1853212,1.1950833,1.2041183,1.2128127,1.2183342,1.2275667,1.2328254,1.2391863,1.2463307,1.2518918,1.2576843,1.263106]) - bf2bt=np.array([21464.978,22181.427,22594.685,22749.661,23218.541,22985.179,22974.422,22822.704,22387.685,22101.844,21612.527,21031.671])/1e3 - bf2cg=np.array([1.0216708,1.0334782,1.0445199,1.0552473,1.0616736,1.0729254,1.0792092,1.0865544,1.0950485,1.1010989,1.1077290,1.1138842]) - bf2ct=np.array([21561.621,22289.611,22714.780,22882.540,23358.805,23138.204,23134.987,22988.727,22560.401,22273.657,21786.965,21207.151])/1e3 - bf2dg=np.array([0.69947395,0.71491958,0.72959149,0.74404393,0.75198024,0.76702529,0.77516276,0.78435663,0.79540249,0.80256333,0.81069041,0.81829536]) - bf2dt=np.array([21769.984,22520.422,22969.096,23161.952,23655.146,23460.672,23475.417,23345.052,22934.061,22657.571,22182.125,21612.365])/1e3 - -#Models where T0 is varied; gamma changes a bit also - T15g=np.array([1.3485881,1.3577712,1.3659720,1.3736938,1.3795726,1.3878026,1.3931488,1.3998762,1.4070048,1.4137031,1.4203951,1.4270162]) - T15t=np.array([13985.271,14489.076,14788.379,14914.136,15257.258,15123.863,15145.278,15078.820,14821.798,14677.305,14395.929,14059.430])/1e3 - T35g=np.array([1.3398174,1.3480600,1.3555022,1.3624323,1.3671665,1.3741617,1.3781567,1.3829295,1.3879474,1.3919720,1.3960214,1.3998402]) - T35t=np.array([32144.682,33160.696,33728.095,33910.458,34564.844,34159.548,34094.253,33809.077,33095.110,32603.121,31807.908,30881.237])/1e3 - T45g=np.array([1.3248931,1.3350548,1.3440339,1.3521382,1.3579792,1.3654530,1.3699611,1.3748911,1.3798712,1.3838409,1.3877698,1.3915019]) - T45t=np.array([40273.406,41586.734,42333.225,42588.319,43436.152,42930.669,42854.484,42486.193,41569.709,40929.328,39905.296,38717.452])/1e3 - -#Check knot - bf2ag=np.array([1.0184905 ,1.0332849 ,1.0404092 ,1.0538235 ,1.0598905 ,1.0695964 ,1.0771414 ,1.0827464 ,1.0904784 ,1.0964669 ,1.1009428,1.1068357]) - bf2at=np.array([26914.856, 27497.555, 28407.030, 28357.785, 28919.376, 28610.038, 28478.749, 28484.856, 27841.160, 27335.328, 26933.453,26099.643])/1e3 - - G_knot=knot(("bf2z/","bf2b/","bf2c/","bf2d/", "bf2T15/","bf2T35/","bf2T45/","bf2a/"),(bf2zg,bf2bg,bf2cg,bf2dg,T15g,T35g,T45g,bf2ag),bf2g,"bf2/",60, (bf2zt,bf2bt,bf2ct,bf2dt,T15t,T35t,T45t,bf2at),bf2t) - g_int=flux_interp(flux,(G_knot)) + def plot_z(self, Knot, redshift, title="", ylabel=""): + # Load best-fit + (simk, BFPk) = self.loadpk( + Knot.bstft + self.suf + self.pre + self.GetSnap(redshift) + self.ext, + self.bfbox, + ) + # Setup figure plot. + ind = wheref(self.Zz, redshift) + plt.figure(ind[0][0]) + plt.clf() + if title != "": + plt.title(title + " at z=" + str(redshift),) + plt.ylabel(ylabel) + plt.xlabel(r"$\mathcal{F}$") + line = np.array([]) + legname = np.array([]) + for sim in Knot.names: + (k, Pk) = self.loadpk( + sim + self.suf + self.pre + self.GetSnap(redshift) + self.ext, self.box + ) + line = np.append( + line, plt.semilogy(simk, Pk / BFPk, linestyle="-", linewidth=1.5) + ) + legname = np.append(legname, sim) + plt.legend(line, legname) + return + + """Get a power spectrum in the flat format we use""" + """for inputting some cosmomc tables""" + + def GetFlat(self, dir): + # Pk_sdss=np.empty([11, 12]) + # For z=2.07 we need to average snap_011 and snap_010 + z = 2.07 + (k, PF_a) = self.loadpk(dir + self.suf + "snapshot_011" + self.ext, self.box) + (k, PF_b) = self.loadpk(dir + self.suf + "snapshot_010" + self.ext, self.box) + PF1 = (z - 2.0) * 5 * (PF_b - PF_a) + PF_a + z = 2.52 + (k, PF_a) = self.loadpk(dir + self.suf + "snapshot_009" + self.ext, self.box) + (k, PF_b) = self.loadpk(dir + self.suf + "snapshot_008" + self.ext, self.box) + PF2 = (z - 2.4) * 5 * (PF_b - PF_a) + PF_a + z = 2.94 + (k, PF_a) = self.loadpk(dir + self.suf + "snapshot_007" + self.ext, self.box) + (k, PF_b) = self.loadpk(dir + self.suf + "snapshot_006" + self.ext, self.box) + PF3 = (z - 2.8) * 5 * (PF_b - PF_a) + PF_a + PDF = np.array([PF1, PF2, PF3]) + np.savetxt(sys.stdout, PDF.T("%1.8f", "%1.8f", "%1.8f")) + return (PF1, PF2, PF3) + + +if __name__ == "__main__": + flux = flux_pow() + matter = matter_pow() + fpdf = flux_pdf() + # A_knot=knot(("A0.54/","A0.74/","A0.84/","A1.04/", "A1.14/","A1.34/"), (0.54,0.74,0.84,1.04,1.14,1.34),0.94,"best-fit/", 60) + # AA_knot=knot(("AA0.54/","AA0.74/","AA1.14/","AA1.34/"), (0.54,0.74,1.14,1.34),0.94,"boxcorr400/", 120) + # B_knot=knot(("B0.33/","B0.53/","B0.73/", "B0.83/", "B1.03/","B1.13/", "B1.33/"), (0.33,0.53,0.73,0.83,1.03, 1.13,1.33),0.93,"best-fit/", 60) + # C_knot=knot(("C0.11/", "C0.31/","C0.51/","C0.71/","C1.11/","C1.31/","C1.51/"),(0.11, 0.31,0.51,0.71, 1.11,1.31,1.51),0.91,"bf2/", 60) + # D_knot=knot(("D0.50/","D0.70/","D1.10/","D1.20/", "D1.30/", "D1.50/", "D1.70/"),(0.50, 0.70,1.10,1.20,1.30, 1.50, 1.70),0.90,"bfD/", 48) + # interp=flux_interp(flux, (AA_knot, B_knot, C_knot, D_knot)) + # Thermal parameters for models + + # Models where gamma is varied + bf2zg = np.array( + [ + 1.5704182, + 1.5754168, + 1.5797292, + 1.5836439, + 1.5870027, + 1.5915027, + 1.5943816, + 1.5986097, + 1.6027860, + 1.6073484, + 1.6116327, + 1.6158375, + ] + ) + bf2zt = ( + np.array( + [ + 21221.521, + 21915.172, + 22303.908, + 22432.774, + 22889.881, + 22631.245, + 22610.916, + 22458.702, + 22020.437, + 21752.465, + 21278.358, + 20720.927, + ] + ) + / 1e3 + ) + bf2g = np.array( + [ + 1.4260277, + 1.4332204, + 1.4354097, + 1.4422902, + 1.4455976, + 1.4502961, + 1.4547729, + 1.4593188, + 1.4637958, + 1.4682989, + 1.4728728, + 1.4769461, + ] + ) + bf2t = ( + np.array( + [ + 21453.313, + 21917.668, + 22655.974, + 22600.650, + 23061.612, + 22798.62, + 22694.608, + 22721.762, + 22208.044, + 21826.779, + 21552.183, + 20908.159, + ] + ) + / 1e3 + ) + bf2bg = np.array( + [ + 1.1853212, + 1.1950833, + 1.2041183, + 1.2128127, + 1.2183342, + 1.2275667, + 1.2328254, + 1.2391863, + 1.2463307, + 1.2518918, + 1.2576843, + 1.263106, + ] + ) + bf2bt = ( + np.array( + [ + 21464.978, + 22181.427, + 22594.685, + 22749.661, + 23218.541, + 22985.179, + 22974.422, + 22822.704, + 22387.685, + 22101.844, + 21612.527, + 21031.671, + ] + ) + / 1e3 + ) + bf2cg = np.array( + [ + 1.0216708, + 1.0334782, + 1.0445199, + 1.0552473, + 1.0616736, + 1.0729254, + 1.0792092, + 1.0865544, + 1.0950485, + 1.1010989, + 1.1077290, + 1.1138842, + ] + ) + bf2ct = ( + np.array( + [ + 21561.621, + 22289.611, + 22714.780, + 22882.540, + 23358.805, + 23138.204, + 23134.987, + 22988.727, + 22560.401, + 22273.657, + 21786.965, + 21207.151, + ] + ) + / 1e3 + ) + bf2dg = np.array( + [ + 0.69947395, + 0.71491958, + 0.72959149, + 0.74404393, + 0.75198024, + 0.76702529, + 0.77516276, + 0.78435663, + 0.79540249, + 0.80256333, + 0.81069041, + 0.81829536, + ] + ) + bf2dt = ( + np.array( + [ + 21769.984, + 22520.422, + 22969.096, + 23161.952, + 23655.146, + 23460.672, + 23475.417, + 23345.052, + 22934.061, + 22657.571, + 22182.125, + 21612.365, + ] + ) + / 1e3 + ) + + # Models where T0 is varied; gamma changes a bit also + T15g = np.array( + [ + 1.3485881, + 1.3577712, + 1.3659720, + 1.3736938, + 1.3795726, + 1.3878026, + 1.3931488, + 1.3998762, + 1.4070048, + 1.4137031, + 1.4203951, + 1.4270162, + ] + ) + T15t = ( + np.array( + [ + 13985.271, + 14489.076, + 14788.379, + 14914.136, + 15257.258, + 15123.863, + 15145.278, + 15078.820, + 14821.798, + 14677.305, + 14395.929, + 14059.430, + ] + ) + / 1e3 + ) + T35g = np.array( + [ + 1.3398174, + 1.3480600, + 1.3555022, + 1.3624323, + 1.3671665, + 1.3741617, + 1.3781567, + 1.3829295, + 1.3879474, + 1.3919720, + 1.3960214, + 1.3998402, + ] + ) + T35t = ( + np.array( + [ + 32144.682, + 33160.696, + 33728.095, + 33910.458, + 34564.844, + 34159.548, + 34094.253, + 33809.077, + 33095.110, + 32603.121, + 31807.908, + 30881.237, + ] + ) + / 1e3 + ) + T45g = np.array( + [ + 1.3248931, + 1.3350548, + 1.3440339, + 1.3521382, + 1.3579792, + 1.3654530, + 1.3699611, + 1.3748911, + 1.3798712, + 1.3838409, + 1.3877698, + 1.3915019, + ] + ) + T45t = ( + np.array( + [ + 40273.406, + 41586.734, + 42333.225, + 42588.319, + 43436.152, + 42930.669, + 42854.484, + 42486.193, + 41569.709, + 40929.328, + 39905.296, + 38717.452, + ] + ) + / 1e3 + ) + + # Check knot + bf2ag = np.array( + [ + 1.0184905, + 1.0332849, + 1.0404092, + 1.0538235, + 1.0598905, + 1.0695964, + 1.0771414, + 1.0827464, + 1.0904784, + 1.0964669, + 1.1009428, + 1.1068357, + ] + ) + bf2at = ( + np.array( + [ + 26914.856, + 27497.555, + 28407.030, + 28357.785, + 28919.376, + 28610.038, + 28478.749, + 28484.856, + 27841.160, + 27335.328, + 26933.453, + 26099.643, + ] + ) + / 1e3 + ) + + # G_knot=knot(("bf2z/","bf2b/","bf2c/","bf2d/", "bf2T15/","bf2T35/","bf2T45/","bf2a/"),(bf2zg,bf2bg,bf2cg,bf2dg,T15g,T35g,T45g,bf2ag),bf2g,"bf2/",60, (bf2zt,bf2bt,bf2ct,bf2dt,T15t,T35t,T45t,bf2at),bf2t) + # g_int=flux_interp(flux,(G_knot)) diff --git a/volumemergerrate.pdf b/volumemergerrate.pdf new file mode 100644 index 0000000000000000000000000000000000000000..4b135571804f5b221c43183d119888f295d375d3 GIT binary patch literal 21316 zcmeIa2{e^m_&Ay=Gi6qWgOK@54#$*vmN6M3^N{&Cgi_{t9y26lo=PMcLgtxLk~AO` zQmIhlzVA`rhurVq|9Agu-L>wzTD!f^e)sU~XFu;V?ftymI*Q6dFkvJacf|`xc>@^? z0)@C)T_BT^f{5r?+S__U;6RNYM8v??-4-IEV2QPKcC&}b$dK8(+JFtkcXd$m!K&zE zEwQ!`aXA0BCGo zfeQ4_Z!KFJM@xA(9{>ih1rh}n7Kb2^qQWSM833Rq259_y`FggVZeAE`03YCj|G@`{ zqFu-U2mXycWk+XVy$}&)X8<9Jw$^So0RA*>UG1?B5ZDe<0kA!>7+XshGN0_`Obt)E z-lG94KLx&@WF0G8o2h5YB6%1GjkfZ)@YkWh4+;IqY4mXW<#|J|VT17jlECKYl(~|{ z4P|hF)&7AT^i{vFKfcbrpSRcUI@QcFbn#+Vqg?#T#;Ui$%%!c%tY6n-PPq|xPKCD- ziE>xD&Iv4T($(uvmS1fDXeYiMsq%qo&H3ZeYOUl)ayjF!QFixRCSD4NIVK+&dowxd zi_*8h)PMRJ{o<9>`%3T-sBzZyxJIR8iCLG8?R|pu&WFERH;uKA@--@gK)$+t|GWig-p#J8i_?j93)Nx}Oa+vIcVwk&iBUa*jnz`=WAt97HmKog#R zJ&BKTm9BcgqQCXjc?yxR{NCK=>nwBR-o^3)l2yUpRdMj%hV>GY5OeaGv01G+{*_t> z+wQwi!*Si-AFtWZM89ydq};MU_t~{r3Fr{|1n)+wE)UOYiOCqVJCR#(9o!Rzk6J+{`$~M{o^&C^dnTSIeHu< zeOLOGx4V5F+1=?%3y$uVg0h~y-kq##$mXG+Isa<>NBh_1^hN>`CwSIHSKUkZrOWOa z3T7=?^RaD-cizmQD>0_K;TA8zHeGP{1Kl>qtdye18l=+~HEzsSRbTyit-`w*T~^yR zcy}<~>-ikDbH_}Qvy@g7VsuS9B1cI2f+{R{uBR+`2T#Jo zcqqDu}HFdhIVO+et z{^X~HXGIpRN|yABR$;YbbjeR6Jdn!dIj+Qr1U_;6@7jshqKg}+Qp=>Avee_c(|OM7C)G0X@vHU!K2g*2oZEc;anSajQvN zY8l~;LvczuJjV%^#W1Er_tzsV_2f@U%&kp}B!()*J%I$>aqiaWeru#LYeY#fUDQ+f zf{WVFx2xTqYM?c4l*>2chF4LF)zfQjmhOxR5iB*nK0JADZ)?TrKG8%^XrhE`*zZ(o zEb%g;Q3F)G9Af$=zNXoI>amqtHQ(W6TySBbclr$Lq#dsoZZrns&WgFjmNXQ zB7G~GAo9lN)f5v}jS^=*g@p7MZl0;)4;_V;U-avsFF0HlZB?VGSLZ#@p{ZANEonV- zfs&}K!#{;UFxdhY?Se+33Emn+Ue2dV3Cu(u&Of2>t?ANQNN?L&n%<8XAD9f`ReJwh z%Gt3A*)P_$OMdu`L8%cpk*YC4mKDWVs*20(hay+`I&yo#STNWW>a;$i zfofT`R7esgZAGfGOB@aw&CO`GhM?Nq)GwnbVy4s0RmNc>jP)t2uFO4EtxhG-3rtKQ zj!&db4*Ru*iFKNurLgIG9kCP_=;J@clDWvt#2du)V~qky=cTCmf%#;NmE2>*$fME2 z*IkBXrN+ZA)X)>CeLDH!8`HJ7o8mX(0#9d7V9Ept<*R4y=W0np3=QiZb3kuK@UCR~ z6WWo(s3AvYh6jwfplUpmD!G`UGX<4LRq1En-x0#EzK9><%N`ie`FxBkK|`vb>_7Rs`5b1 z-&arto|`71zaWI)srPLB8%fGIvNuJsE-d%6A^o_qLW1X7r?n}G-={QREP)&;^IP@o z+9QNg=8qHz zUH{N=h%&m6+Y!syq&V_W*uu@V>w^iul83!v= zkZhDYRjSLoGxdC9rD{7qWFqBs=BSm1)XT5*9QJsJgW~*^R3E)vAQLl08YG>oSrL#g z-F&FmdlxyTSG^%|RY5_@tnbxr6$Mg@rco?1#Lob_x@CIIFH9$CHZIHu_N}7+f(KKl zM4n-k5EVN6B8VaQ$4mO_J3Wk)_@r$~rN^V&3GX6o zP7iuRu^;0*&4?8sN~R*30ZPg%L%Y%NIKIFqi1gjYkwfJXWH%SZPUhSfoRy#Uy#i zY(EAYVBW?>K38|e)TJy-@V$%XpKYwtWLm4w1e3p9-ts@+HIC(|Qn z)&gT_2aZL(ZT2QfHKC-_5PhrQ(sS%&l$74s*ooVr`LC4&?^^XUM4fRpIL`&AM17IDlSwzXRpW05yOOnA88xUOFNU9>ncvYa(!+#1Sb`6jPqEj zr;l=AuDjks-?e<%bmGoh%6$yLpHuc#sRna+pzVGa*nGZ_MWv}p>j%jcC zk%-#yqpIrwCGUOWB9ty3>^k`_gB##d6~njyzpW{7mH1B_ZeL|rHH zo=zn-@u-%2vi;;y&Xf}|T&;|nU998PafQnbm;EB?xkK9g;F%)B!^zSz{G)i8&?|bpC@7l)u$Lp4Uo)TRfH{T@$Ql-Vt1owB% z*kK-iq5J-NBJym}h0a$lf+~?l>@u!Z^JnjEH#}%=?B)~++s1zuqtFXlV86_VU$z*| ztN`r~1PqS+ZJ^<-AF$FCuxfrc1kI+uF;vhf(7O4@NRf_1(*ztMA-H7TMmw7@b1Kd{ zaNA8uD)V--S^%_-Xu?#U;ySGu{v8%6xt?$`tffsP54nHytw-h;Zk;qj*hwHVHE$sb z?(tSqS4(8-(N=3-ktx}BZ^QdR#2>{5Jzzn*t;_#1sNk@Dfa<`K2w)_OKgn_qk=*{7 z(|rFv=h`;v;ZVUd`rw3^lOZLLjeh=&!q3)Ro@BY7PgwQ|xMZ^qn(@o7poy%7$@2|7a}#aq(|zxb-A>-?R`5&<#C%?vud z-Glk+BA&~a&M_`SpV+^aQn1j~_sGt&eEn!9BKPt9=J5LQ7B@8g>8b+|1X|O- zH&Wix(-voOi>N4^HBvMbP_(^ZY3QYI>FOz@ZRPCf;RRao4we|u_6E%2U#7H@D`1;~ zrZs5uD>~ZQ*#d?yXz7|l#Ka&XT7dnGgorqRM)-~s;GYgMGJtvOi1EaNgF>KWBH--u zI~KJN3<^vjgb;-R6QQ8MASiKg2(XEwt*14{(H-lCfxvK{f}NgPmRO7<&KehnLP0wo z`1el<=&EV?YbLm;D4EDyJy;CvHgyW9sL zIKS~J=V=Xk0MKwW==g9~wRN<2zyi-jK={EYLU05MXW-idJ_XcYFaLi5IhlwW;AC;M zmUFdtwgsotvh=~h3xF#M{M5(Vx)_4>JG-$Dp&~~Az`1K#;5v7B0166){2xeU|I;xL zFmVVH*j>Qe0ue`uLr@4585&rYm^c&y7ZZgb&`?|w1>gjqgY6M0*iMcR1Ka;sBK!Ru z0kj1ZIL7ZJiU5ca3L-!p2_#@3v^YTHz;6T`A_`0{1}0HB1dbLZ0}=!!ih_uW0)wCs z5HTdsRvZ8T3Y5Vde2xGT;BpcJ`h#Pl#nC{%QwQd_#URAUz#IeuTojm40O5gsMS+Fy zi~+O-=Eao}NB|MyxFjlu#^q=+0A8RT1%m?%+?f{%%qt295Q7H$At9nbE)F0I%y9`0 zwnO2N0OsKQxJ?5Q2*3a!mJFOz6zGFP%5UJn97k#3Oi1tt6hsjSHn=n}2NG^y!0C1p znJBPCTu%_c{}uoMcZPuP%q|K-1wj6r41RCd?;O}NvfqH0}NM=QT* z0;GrC34{yx8yEwmjo%4eISMKc>>&UvNCO8XG%yE_E?^+0{!PDW;&%?v0stS943~rD zT?x0>JI}#3|0Iwuz+DCNf0Y4H4*ceM0Jq<9sR9fiKqBbs{FjZuf3p!#Qv9blMvCpx z8?s1 z;7ah*?Db>S)2E&*7o?stFon>eIP=hgmT&D|F385as2bj~V3nef z&uTTtG9UBV^*EPwyYL&r$;8THqJ9`Cn^1p74Z#DFXF}j6+TImB`I}?B=?RRq9uI{j zB^j(}$&ZKEvZu;aZ(d#;T+ec>EKKl_<3x5AzOrxXp7+^4l&XCWuu1moXYk$9<-hO; z1NMEd-lai9z(GQ^mOX9O%HIaToL=R4KTe}6GY&Hd&8*mr9h{{=hnj_w8>Au4(ns&l*c5KU2$GTEGR8W~^MIa!&A z)P7qBM=HCw0mcknvWHcfvn>x5v}Jl)_+jed^^Xqq(;pc)0EpeB`d>f*%xE_d(sBKi ztxyt$70ou1H~RQL49-yZZjFas_G_ay?WreapkTlGJT#TSEvTB;%fif@fL@+~*=BZ9 zk}+D{j-e%!rI|5mex0X%gNhzS}r?GOru zwQPBnL)!tXwU6Hgy|`-HpMedgv+fsBJd0gQ4e#zsiCuHB)(#Om_ZEBXv_5lubfcYF zZoeh`v0GE-{Z<3Z#3xJn-=1N=5(Iv|=XCaFBv2qvtzK;1-twN8k ztP;GEld(Bq$$zV)BanMZLsU=lQI;IZVvGI)eva|We&6igf_2Yr6uuf7k+E8=ej~APzVsCUeCYXFQC}6Kgxy<}Ds^K9^Qx*)J=zhS?J?nE}^EkdV zqvOLEO=YSR1RgC<=ZF20`NCfXAYD$gbE{G36E&K)$BB@`RsHnZldf;xt|BsI%p>)* z);hseY*;k6@eFaUi`M!?NhZ;WceHDzK}Q|CtmiAZso^ZzU#B~7-TP!R2m5-bOtuTv z#^TJ2ctfW|IL(-HLd=6vi?wrfz{s<4L`QG-L1TDwviVwia!IQIM|P*lLzsfL@tjSu z9H~KN_8~;Y_yg$a$I7G^FElg@eEYgOa3W@s#oJK5s0~f^!8REkWOGZB$NAyRk=2*d zL5bJITB3!51k%ytHs7n+m5#p)p8lR)&R0Xvy{V8YRXu%d`^Kvu8>^dJ#4#tYa2>Ea ze|si{fbQkl(n%Vw1TcUgWkQ9Iw0@t+SuHr(aP8&QNnY-EQ7jV-N*B@^G>!Gy5|GIu zSKobfufG0S(s|gO*J9?^VB$sVWRG3r@d-kR+AJV$)$JIxdHe=wgFKC)M>y=+(mois8v#>y{ z63JP002xeuRf`!11wn0$Gm=MHh{wINA9yW5LcJz-bIa1OSf z&m|w<-hAdHKfHchfcuN;voHjP_k%;%+*r?a!q&@JLXG$;2aeU(w_KkJTQuoQC|TK@ z##A^)pECKB=uMt0_L#Gqg4&@f4=Z1*(FLZzG4~GYKi*|~iVDHrlx{dw(3MQlqlA){ zItqWreMADIiKWYCuvaRjp2-VX4 zW&7zf91dnv#V*z*%Ie*ZTPU@(l1MfbN$u;99~ic-4mUp-|3Q74!{F!v2-|&n``-xL z$4a!(NrW5}{l1|x&@)VvvMt5Mb*N_>R~wQtPA7ci=zUDfn&Fed#dX`Rd(y-PV)t5K zv$N$}Ej4Y~+G{q>lxT9xutfK#6*I0mag5vEsl!IZE2#A-WAh`=9%?>4%U~}&No4Jq z-GT0PjWMIGj`gZ6);Cgt=23)T>)KNwo(2V9dgI`83*w2qIngWMCI!9GSEv@rQL!BK z^XvS6L8##TCEpeNV;H$@v#E79rN#wr$Z=AJr@B@jvfqZBZ}!Hr-91U;DW54POv^F& zj=$(grcOo8&0IZQ_fnbSg^C+j-=04wbO64-1vqLrA%T9%RUm zIs&^?gi559V9iFHvJ|t9YgK|B&VGk2k|yRs|9D6i;(6S5$xxLyvWzqCJtB@ThyT(I z9|pnY9+W9>lJqR(jh| zS4j2TlvuFQ8mi#Rn0Tx`nX%5|e((Hqo71V?Shj)n?AoDVG=0O`OeBg}EB5)q(i*KJ zzldHVtI72%UNoswO!G`M>GoHdU)gML5O0N>a2|lczg1}wfLgx$6_wW0aP5>OnVhFP zm;2~r*08VMnW?59ZcYk%rTk5iur8>7U&8&j8fG*ep}qIat4yVeE(@EQWcWrjM#pl? z)V)amz)@$06|NR%RryQ~zs|2%GFUd6w=SO+N}O=|V;=cKfArRjNI={8ZOZyL(r~%k zC+7I%yIggdGil_GTFxiZmlWvSFn;(-_ru~EoalmEOd+{7%NNpGB6%0p4*d7BhZL7m z3@Q`i>`&0P&`;T0u~^KH-odM4Hc(8vwA39@I#d$Vq}U;&t@Zt6nx>Fpm#aSW`OTD= zvCCK@{!--Ul+Pn)YDUP7B7JGTp*L06=bt)rgh{wQ{T%t09nlG=8M4qoTjZtBsJxyu z_BvX{dcS@M8_1se)LpZ7@k-BE?Jex2$)Hygh9qD1t(M%!Mye|1>W&V$%=K74zkI>l z$1!n^lWD_)KY@EQ`V_=et)g~<90 z8=j$_+&-`I%=AnUn!yJ?Ps#j>+xA}T+?>H0o7-T<-7RV zIzGT7<$_?bgO9n!kXAdR&_!CI!pihyX1Dqe>1!&fE{T*%!$RcplbBiU6>?=IOu7a+oJK6)R*Y%x{nPQH7!j zNpB@3r(H85r-p^F3repia(IgUJbopAqS4qdsdr-g+B~Pe?z&!yu$8=7#PQ+0qaUU7 zoM?=zn?IoiZ)SKtAZ8D(Cq(4oS^KNO@0ZM_a?qZp@M57dDR@}#WYs=xGAb4O!u5x; zoh?rYWxS323#Nbw$Rhq7CczXGiWY7>M6gWROKD;;UV2n-ZR&2!lYVDPUE2jm`#P>O9k+6o~d+%}LSB_n&T&65ns;n*Ya%NFxg zA3})8>wJ~jiBfsY>5YihtzYQG%QnBk2=%@oqIT!K`h3J7JMHmB$3(R=d@l-Ue7C9W zc@j!|N3E}~kz4=Zk~#p&zcn8AYtl*pGg5?+Ot$E|2dchFPgt!j164KCEX^l_uSqSO zefsr_Z2E5B=jJ*wo}%r~sZde0iNu`XELvA1CGt*F6tJ z)~X`Onx7^pUsJ!t%}I}Nb1#KBhvIR=lR^cd@{+!%0pD$BadbOK>EaJVFOeKgwk@nK>s`J?+9(egEvM5A1KUMb3R&WYky+ z={R6_{?^Wcq4p6`LbnGIph>Tk;d4U+$i7x|uTn2=-3t5iSk8}0p2AA|;q7y0JuI>9 z@%Mxflvh=lKXiSETQ@y1g77X7P$Wi{j~0CQ?qn|eGW*1;>Sq2q5eZJ3fuxyoS+2_x zce<17>JUG%=Il!xE417E{PV7sc*12LWH3%W^+_4w zV)*H)W&)Ep&8g+m&sR>dc}sjgV6lH|Hi+(BD*}Zufa*Z24WMw8E)<2tx1tTDQxMK2 z0Sta3+J3wm*)mLhZl9+Gkb!nS5fc^s&)JcGt2m48)teFUdpdx#0f|Qw#=t{1XLlIq zYj53hzIpq4v3D8fGA8pP#DCxk(&du=vac99~mYf%5oEX~r) zjQ4z?g`4kD)rliLi5yIX zT%KGI@)OY(ayAIy#iU;)yLKbJ~n~`H_`@;Jr9ap2nw&`o+4VQB(%5Qe1~!srI{u43NS2~aY)hgzu)PYiv8q( z@)RJiq-C|E@=LkNoO?uq{^fH1ET-tg(1wEh*YjUntftUek%7H^tR#c)^1eFcPkQfqqn#}23w`wSa_Lj-WD?&FSnJ>F3<&sME+nek=+X(L zSP~^-ghH-x56(bT6YpKkl6+pO`nZYxj5)mUoUm6fzp>@#LItm=qRTmUMVnM#6K1Zl zO>i{mKTW&D#E=&0tIk@TDY)iNj(Tt?lg5+$(;CS~fdJiB#a!1dV)?PuwRa#lVag^r@7h`0lk{W)_@CG|67X7?W+g`r{w7 z3>9ChiB^2=+e)B$V#_ALd057D>HT{ply{i;yODtCfULn2cyjvUc%(DQ3fSF ze$k91J`t}JZe#tI5|`6DsC>iTy4s+)`ISVHZ*#AWkbUxxoTFQtdXiS7fJ}ElUm-{E zmmOz)tbxjO=cjGRebJ8iPU(UrFbrp_Hie87Sd7;}Fx}Kl6m1!_(VP_OT8Fo*F!49c zqmRq9e^Q`JOrn|--8TC9*^}~?emrlYtJ{NV^XyV8f}<=+nN3_~M%BVJM;%4ww1!(v zgaR7xnQtAhA^6+H~YzdfqANbF8Epe6EOgg~K@%f?;MJfnU;~Xe+Nis`ULJbd|Mb_am=x4HxzvDZtXrPJ#I?B0(?;7G7_fu)g+I9lq?&#MZ3 zq)F)x8q;~7Q%uwPY}dHZZoLg}6~a%aI(~K?8N5}!b&~?$dt1lgNF{f*`oOtIuGL62 zk%@Nk*mS46ar32@LX2-}TN)Dr}e4GF9A zg=C9tW8C@i0+ZMDuAdfq%and&;nL-}L(-iuV_%GR_cO*Q-xVquy^tzS4W-o|P4y^G z3{0UttYEKl7D-;%)8{8I`Y;xqm(kqrZb8sr6H-uEsnA7>$nurP(sMgIAnO@7jiHwyM49?{8lkbEIx}($r!{rRs&dQi(YjWHrGFM%HN+)$( z%O)tNJA$=7xp#qJMbEU_YXNWh$e9ynH5As>l@0FFmV~06e)Mx+q1KT}DmTt(EaDN- zu0Pti`oVR@+aJ%O6N)=++;jF~sqY@xO5O88 zSB?=~Pw1I4c0npK4@!mIEMB5N{awT}xg?d@ma3;hf9iBkUl9E@#xcyR;rdEU9xYp) zKzqJqX9OmT*>U~nN`aW!`4>LWCkJfb-)d$6B6jZ{0#yqru}~zGI7vB7HZ3tyK8k}; z^-#ya+aS`p0OP~5$MVIAt9Zf%f(0lKy{V75#VHzHgTjKUK+phmT}wV9@IGqjQjXSFN1aFygczsGg+Q#*oCSQ*0d}c~VtV z)jCPFjxCrg)2PPmT806CjozDD$Igc>6DCi36l2|EHzVI2HrUjDe=g;tR`9JP*mov_ zCduF!wRLl)7bdXfEk$TG^EXmAB_5)X$m-IG2g4R&18ZdC;@zv1nn7ZV2P|~A?(|;{ zPy}!azI&lUdK#Ec;MSEVnOQDkop@I29k#)zdAlwtSxK47;tT9;Sfu)8W5~UmPdo0# z4#);v%)cx3Lgg zj+sv8=a%%{mh`hmx4k&ZF8Syew3nWEHWR%q)9EQhD0EqQB;fq)jazZAHNVoih2BRDdA>RM*igBa!5gJ|;heKu z+)pF*kx?;+Ie79Xfq zXFT*Z6R+8#(6_tm%R7w8A1ftUskYY`89`U&G^ZImV6$kO3POgvPHUfja2B!1fu<)d zwK6%|=17-eMLJ1yUSC%~(~o)U$L9~9B(H4Y33^`kJYZY?R!u~}_tKvz;K(2X92vuW z^|FN@wTD*A%Ay|tURBy=dQGUs44;VcBaH`R+>?dvv!?88+|}|g#&pcZp_do}o#O-eUCcDfhXBI3TvW4@-oF%U*VRQa$ zN}F21S~ac&9@u;)VL-J-CTtwE8$rn9(9GwaMJ!6CeI40z7%oKd)tvj=<duR!|&DFa2zn8olXH5nRSjp?|^mvt+oRL?55rCy;8ge6(I*;erD&# zwRW@5d%)j?dTuWVTz93lmJT%tydRV&E$l>aot~Y8x$e`w(zG}wV;KWKJwxi=*Wndi z{TwfH1wRI<@rQBfL|;3WC}it;FJKZydFM{5eqr40=$p{^Z?(?!Z0GY$IlYoj1$Xd^ zyqIZV?J0cP_%c8G!PNph80i74`CA1E0XRl>uSQf4?FuMJz%^HYNPz3)$^>9(Y_?FA zPS6o}VD(ZYO$@yI!$S@3E@3xnQ=ij_&0e$(?I7(fr&-nuQ!QySk}dF5 zsN0B0@A#DXa@1O$uh5;t5zQHCQDm%!kqlSlKIN*!`hR0GFoKLOv#7bpMCv34cn3aL zr9fX2^wwqnf;b7IU3h+6>4jV%chWU@e?mcnw2e*bBCdRn4$IHbHT1cskW4uU>*iRc%4ugNRsSwU z*AmZ7%bL>jii$Gt>Up-CG1;h>_;dNETa6_ZUWazP@_@++y-me4;f2B5Uul4EoBW3D=n1@#9o1$ zvk-SFsS%G8!dRPym^gtfmYLQXi1qZS_;zo-Ob1MEDqqMb-Mf2YfX4G79EM+~+&dwa z-%8zmj)cZ0o&7;8gwrg3!*JC&`OEM<$*K5;TD$Yo`K^&lB%QbM&IbH!=nm9)xMBWy zHl3f_nl}7%;CoIzlw{tKy9P^`vZ$J=ECF}=FPozZAC`&ka!I{9VA;Dh!2d^{9r6DS z*N&1NIu1~emvBtMSKUuLKq+c)d-6In_cQfcwIOb^QO}NwPY*DeW#`YBheb_ozsH_G zJHc{B?4n`(1$N^=&-mKPRo_naGl4FrL$bVc(%lU3UUbD(z6or;zGM{F5x0CSL&llt z2rKK&y3kW*T$me=gIEn@K80s`Ug*xFX3c$uPpwl~=OfrbVAuGGSpA!jqo*#Smgd@{ zQ|1mza~~;hEj(%}B6ogPm4RRWGt`axXq}*@;9C{FiyAh=2@W~kVf;)pgt3@;BVMw0 z2!A`raAfoVMS?uj!kS4=J?AUfvF0`UK(hrh&5L2qpF%`BTFX8>0wD(n$trafSsW`n zwPP)$CeW-LK&F4IL<6^*cY`DfxDEoK6pEL$-yEH8nIPsMd)Csa=5GCJD%o6bTK28P z(YT^OUb(Ls=8%bZc!o;Tx@rU4TFs?O`ol%b)A@+oz-gVXQyK|2mpAoq5H(*E-a24a zf2*6p_uU2+)kRYQzSN~P`YqWLF(N+N`6marS_WfA4cXexh1I5TWPH|3F5!5qFXNXX zoZawZeBc6dKCM>BGoZe%Nu)BVHPS93xXtleJzl01p2GNOiLL@8+sIg<{>j@f2xX$~ zxl&JA-wjN!H`NNJTgr8eSkz>Ij^04o69{4f0tZBP;zdMU9bIi9z$J7M z6(EKHBBBlv1DzxDH^h;qVIvo&4{F>p*^m>-@2KyX4w`gz13#@t%w5yLGbvtuP$~9KgYz zD5ZnXS%4-xVLAsu{C^*wg9K>ne-fVaTmLuxC%oc!n8p7v|Np~t{y!d`15n!k+3*}7 z7zGU6`JaX70Fg&P#LNGK@ErI*p*diP&i}jNIimj^o+A!~Oo3rryFzoo^6#*m-}(Qs z@En{-`+dh#1%ivg(F4T@E}8@=X+nTa1NkE$QQ)9X0hDwAcy>zq5FqjjS7!-?E`cR7 zK*Hc|CIjtkfIzIBr?wD4Fyh+nh%LZxg8O9$d=cqa2^4NX@DOeYoV)|B&*9qGLx4a# zT*(0f(7{e=M|=V&I=G*XptQx6oFKT{6aW~`5P*ttbuJL(ujyPNfb7N9;lw8L*9>kD z+z}uEF(7#%K-3kki8}X+gwrTcJJSH_^{)_Zu*2_NSB61w$`CNIG8_UB8#uzhO?RvBaxh#N zuYw%V=C}O+=c`;I)|L?5S%Zii{C6`zSNNqS9MlZ)yBYXttOFZ@4}dt}f1v@P`vv&F zmEh|BsNCZngbD46R(5f8_40%uf1yfw&oeJqM>rhNlTbU``0pIx@81|xaPtDCwNFeMG4iyJnJZK@P z7%vpc%SR@HYX`(Ip9C}^cV{=Ov!fNn8zC$z3>SnrV6pC=5+Wil|2z_Q!`K6zfY@sr zFKgh#V7q&{+t@*@EUle@z5g`=-!6>$U~KKk0FEP02EMezK$!pjI#KXN-q%j_fbm2DC3bH+;1@{Rd-_75-~r6OIv_Oc zk2(;`f3^de(cbZpq5vP>UxxsQ8!ZZlL(UmGvF{YysvAGGtGzB?gOQ S$1t?G7>bOWTS;4)?EeCc38bX} literal 0 HcmV?d00001