From 5d7321573cfa26d329ab7f6f97a01be170fcf7a0 Mon Sep 17 00:00:00 2001 From: Brett Morris Date: Thu, 26 Mar 2020 11:23:25 +0100 Subject: [PATCH] Adding mcmc example --- docs/index.rst | 8 ++-- example/fit_spectrum.py | 26 +++++++++++ example/generate_spectrum.py | 26 +++++++++++ example/mcmc.py | 42 +++++++++++++++++ retrieval/core.py | 67 +++++++++++++++------------- retrieval/data/example_spectrum.npy | Bin 0 -> 52928 bytes retrieval/spectrum.py | 6 ++- retrieval/tests/test_core.py | 11 +++++ 8 files changed, 150 insertions(+), 36 deletions(-) create mode 100644 example/fit_spectrum.py create mode 100644 example/generate_spectrum.py create mode 100644 example/mcmc.py create mode 100644 retrieval/data/example_spectrum.npy create mode 100644 retrieval/tests/test_core.py diff --git a/docs/index.rst b/docs/index.rst index 07c8b91..ed5eda4 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -7,22 +7,24 @@ Here's a simple example: .. plot:: - from retrieval import transit_depth import astropy.units as u import matplotlib.pyplot as plt import numpy as np + from retrieval import Planet + temperatures = np.arange(1000, 3000, 500) * u.K + planet = Planet(1 * u.M_jup, 1 * u.R_jup, 1e-3 * u.bar, 2.2 * u.u) + for temperature in temperatures: - sp = transit_depth(temperature) + sp = planet.transit_depth(temperature) ax = sp.plot(label=temperature) ax.set_xlabel('Wavelength [$\mu$m]') ax.set_ylabel('Transit depth') ax.legend() - plt.tight_layout() plt.show() diff --git a/example/fit_spectrum.py b/example/fit_spectrum.py new file mode 100644 index 0000000..f37fa9d --- /dev/null +++ b/example/fit_spectrum.py @@ -0,0 +1,26 @@ +import sys +sys.path.insert(0, '../') + +from retrieval import Planet + +import numpy as np +from scipy.optimize import fmin_l_bfgs_b +import astropy.units as u + +example_spectrum = np.load('../retrieval/data/example_spectrum.npy') + +planet = Planet(1 * u.M_jup, 1 * u.R_jup, 1e-3 * u.bar, 2.2 * u.u) + + +def minimize(p): + temperature = p[0] * u.K + return np.sum((example_spectrum[:, 1] - + planet.transit_depth(temperature).flux)**2 / + example_spectrum[:, 2]**2) + +initp = [1700] # K + +bestp = fmin_l_bfgs_b(minimize, initp, approx_grad=True, + bounds=[[500, 5000]])[0][0] * u.K + +print(bestp) \ No newline at end of file diff --git a/example/generate_spectrum.py b/example/generate_spectrum.py new file mode 100644 index 0000000..6fa8e98 --- /dev/null +++ b/example/generate_spectrum.py @@ -0,0 +1,26 @@ +import sys +sys.path.insert(0, '../') + +from retrieval import Planet + +import astropy.units as u +import matplotlib.pyplot as plt +import numpy as np + +temperature = 1500 * u.K + +planet = Planet(1 * u.M_jup, 1 * u.R_jup, 1e-3 * u.bar, 2.2 * u.u) +sp = planet.transit_depth(temperature) + +output_path = '../retrieval/data/example_spectrum.npy' + +np.save(output_path, np.vstack([sp.wavelength.value, sp.flux.value, + sp.flux.mean().value / 100 * + np.ones(len(sp.flux))]).T) + +ax = sp.plot(label=temperature) + +ax.set_xlabel('Wavelength [$\mu$m]') +ax.set_ylabel('Transit depth') +ax.legend() +plt.show() \ No newline at end of file diff --git a/example/mcmc.py b/example/mcmc.py new file mode 100644 index 0000000..5c4b321 --- /dev/null +++ b/example/mcmc.py @@ -0,0 +1,42 @@ +import sys +sys.path.insert(0, '../') + +from retrieval import Planet + +import numpy as np +import astropy.units as u +from emcee import EnsembleSampler +from multiprocessing import Pool +import matplotlib.pyplot as plt + +example_spectrum = np.load('../retrieval/data/example_spectrum.npy') + +planet = Planet(1 * u.M_jup, 1 * u.R_jup, 1e-3 * u.bar, 2.2 * u.u) + + +def lnprior(theta): + temperature = theta[0] + + if 500 < temperature < 5000: + return 0 + return -np.inf + + +def lnlikelihood(theta): + temperature = theta[0] * u.K + model = planet.transit_depth(temperature).flux + return -0.5 * np.sum((example_spectrum[:, 1] - model)**2 / + example_spectrum[:, 2]**2) + +nwalkers = 10 +ndim = 1 + +p0 = [[1500 + 10 * np.random.randn()] for i in range(nwalkers)] + +with Pool() as pool: + sampler = EnsembleSampler(nwalkers, ndim, lnlikelihood, pool=pool) + sampler.run_mcmc(p0, 1000) + +plt.hist(sampler.flatchain) +plt.xlabel('Temperature [K]') +plt.show() \ No newline at end of file diff --git a/retrieval/core.py b/retrieval/core.py index d58cc73..632a375 100644 --- a/retrieval/core.py +++ b/retrieval/core.py @@ -1,44 +1,49 @@ import numpy as np import astropy.units as u -from astropy.constants import G, k_B, R_jup, M_jup, R_sun +from astropy.constants import G, k_B from .opacity import water_opacity from .spectrum import Spectrum -__all__ = ['transit_depth'] +__all__ = ['Planet'] gamma = 0.57721 -def transit_depth(temperature): +class Planet(object): """ - Compute the transit depth with wavelength at ``temperature``. - - Parameters - ---------- - temperature : `~astropy.units.Quantity` - - Returns - ------- - sp : `~retrieval.Spectrum` - Transit depth spectrum + Properties of an exoplanet. """ - wavenumber, kappa = water_opacity(temperature) - - g = G * M_jup / R_jup**2 - rstar = 1 * R_sun - - R0 = R_jup - P0 = 1e-3 * u.bar - - mu = 2 * u.u - - scale_height = k_B * temperature / mu / g - tau = P0 * kappa / g * np.sqrt(2.0 * np.pi * R0 / scale_height) - r = R0 + scale_height * (gamma + np.log(tau)) - - depth = (r / rstar) ** 2 - wavelength = wavenumber.to(u.um, u.spectral()) - - return Spectrum(wavelength, depth) + def __init__(self, mass, radius, pressure, mu): + self.mass = mass + self.radius = radius + self.pressure = pressure + self.mu = mu + + def transit_depth(self, temperature, rstar=1 * u.R_sun): + """ + Compute the transit depth with wavelength at ``temperature``. + + Parameters + ---------- + temperature : `~astropy.units.Quantity` + + Returns + ------- + sp : `~retrieval.Spectrum` + Transit depth spectrum + """ + wavenumber, kappa = water_opacity(temperature) + + g = G * self.mass / self.radius**2 + P0 = self.pressure + + scale_height = k_B * temperature / self.mu / g + tau = P0 * kappa / g * np.sqrt(2 * np.pi * self.radius / scale_height) + r = self.radius + scale_height * (gamma + np.log(tau)) + + depth = (r / rstar).decompose() ** 2 + wavelength = wavenumber.to(u.um, u.spectral()) + + return Spectrum(wavelength, depth) diff --git a/retrieval/data/example_spectrum.npy b/retrieval/data/example_spectrum.npy new file mode 100644 index 0000000000000000000000000000000000000000..2b848b102c0d5bc9ffc666d6fef0c1f5cec30d84 GIT binary patch literal 52928 zcmbSSiCaxy7rxE&oS}>fmAOnwteq(`WF8~4424LEDBLJz%sdt)G8Hl=5s{)qCCU)$ z&hx#^ee2iWf5CS;&vTx8+Gp>z*IIj@v)22rO$Yu!sGHq3q|j1VaUI@5i35_ z!qUQbC%zWdLiUH@8>g4l!o?9yuk*gtLiY`xXQ8|n#+GyuKCGw(0sr;WL6x=ed(=Xs z8$VF~Rh%!5ziZ)zNmbX2e`>+5rX&8t?XThGjS}-gtXyr84!l<+L5ZR3LYdd&gyhfiyYV#m? z)L@R59uMwuZjn6(cpY_COV^kOZ;v>PjW*@MyX3WcMix9UcRw>R&WZ;i)@7!aws^Jh znVHa@2X<`-4{Tl5FT2=@2hS(2>D#j_4;+sU3BA#s2S$Ax#<@E3Kt1h8!lPb1@EmWl zhB)(}aaHouw|#if-|ow?8U1;XUYaAU7|4UEy^p_J>BfTxmo~XbhVbAbKgu+GI1e(P zi-Yw%c<`*c%I~Nr4|ckS=XU`fd_P?K^oAD?YSlAb-N*5uD5$E>n+e#zMP};kNjSca zR$BFBY!_MMx@8&Wm67bxy^$+?S{L1r}Loq7Ryb7dpvOepfmVLCJ#(I zct#CT*3;7A6gJLm8qGNbQ+qF1Pou49LXz2U*$_@sjF@6gXS zgoxk0=fQ^&HRrc{;6blw$MP=4JkV9nx%{>SukJTgTS|FQ+;K{Ix34^K-MD6A!8e>& zZ|4cyt9WqA<4JMvUp%mSHc9Zg8tqb`eSGgf9vr`D_0EmQgD2aQ;(rUU-^^b>kBQM< z&l_KkmGNMtRhJcq06WrgU1p+yzKQFyzzV(nRo9^>_k3{P!8WVZ89J3 zHl{wYp2i3Bze5wsrt{&@do!I2v-sds< z{J&3kMdUKnN9MV|lUDNKv&rosn>Bn$)!%=paxEX!A#=y4Zs5c0QB#g@*~o{>vrmUX zFdy{Zu2Y$W@Im{L!su%#A9gQTk(m;X&#QP1TO;^T$bVQiHj)p|ZZuD^*~^FMXM%;* zTt1{c&Ct4kfDfMY_j+*;@}bkbL-z9yVf&-H@h(UCkgRThrs+5z@(!p+7o9{sE&aUf z>KXif{nz~M=lQVuns4E>OMIxlnJCxNcFuxFe8@_Z-G1|w4;|$*J|#Wp!{)64 z$9CuNL3Ut`#4i{7ITywA$m7Gme!NJ#d_JUqE_zT`z=!3q{y_OhJ_HUF)IIvd2h)hX z24}zUVQ%|1_ji<|y<6Y5F8;=c+cUoNNEIIz%}B(Q-cgofv*4NQLrt_}ph+9gb&HJ(D!@p})_! zPJLSVaONdYH$4=-gib5VLu;)gS>J+|rPwKUe@C4lUf*G*ke?@8=mS z+yxN$Idm9*qyWxqE*tzFEr4C4w(qL&!uyrO&XtT4z}nuHOA5UOFmU(CygVNPc&wI` zzMLw6ip|U8pHCNn=k)jdC$j{wEpMsh;au$J!Gd!S76>5QQY60Thx4kL!O!p)K+clv z_>AQOSd_WpYsM-8M2|aKcrO6u2wS=C!FmBaEbnsnVUPeG=&3F~*(`u54;_X--zosx z=ikC!h6-Si&Z3oh;rRT4vigDu)Zdzz5haoM|1YVI6?+A+zf;hoU(o`PU$Uy>#R%Zj zkN3ajaRNBB`Ab;SVF7S;*F4cbE&%&wkFQys6u_$hJ+Cfj1kj=VzG#>8IG>$CVZ$y7 zz#_l7{kS9nXpdXJd}gu$7WhaPE=v_a=echgbZvI`oeXt5?)D3V`F??wWn805pg9UmT?^gg%e!CoR(xLaA@+b*`Zh);4|2 zx?v)O=smZhi_C@4Jtsm}VdXqg)2*0P7Y>;~k;q;2~f?kt_ z5IAM+zXj8TAnf@1@u3+)a8va0&Ydj;ogNj38s-V%tD*H#*F{2*J6TU#yF>`fJKMiW zTqcB`m)F);ti*9+Nv&;w5H@|NI_A3`pWC&sJ~jyF@t?Ui3t^DzU!zWl5Hy>5p7sh8 zLhOrmvfVpyKKmOCvUj2U(YD`PqJ+RbUmpluA#8cg&x$-Cg#Mje3SJ&Wxxc7R>mA1N zCLP052>2@77=JRqPI@7Pf@S7gPQOC?etT+O`$h;?3lam}^Mz1!N9RaPp%7Mf>9y@!kr2ifPU=43 zGxl@x!J()!Aw1rAGv`Z%5H1A;Kkic{g!s3Uw?+OE!h%mb>q`F!VH?MX45&qY)o7+h z3xrTL;=RKUu@D{<7(R2C3!%$#Rp;X>eExFjIB}g2K03<=`!or`Z+9#|r4{wju{Ja zs=Wv@2lq}e?J9yDOYisG>>z?6({p#1_7uUH%Ja#d&LY^T%RPCeuLx3py_#tuhln?&Fw>c;)BMFe|h9vy48lxb0G8*Y~Ul#x`Cny?H?dHWiZt1}2JNVVJ?kjH@CbH+!25 zPZfb+qpI_>G!dNpm90DG4(fI6rRgE1a7A}ugV{bpf&e_-tsIFjGFhi zoA8AQHhX`!3eClNX4HJsek+2>-)8SR{9XhPoHo~Y{wRWc$JnW9pF}VvQWWk{DuPZn zyr}%IA{bYezkE?8&P(*Oi|D5amYm&wZ}%S&^!2HBw5b(=NBoA3*90P{v|D-7Qz8O~ zf{tfC$wlxZ`D^G}wFt^SsNA#~M39s4EdNY1>UDht4Av4uzHtB1LR~S~eC?YPU?_&t z_K!d4nu@{n-tH%f7GfCk>*GGK5kp{OegCTVV%WcDT4sd37`k=+?-0YC^N&|$JBZ<) zVsremUU-GmCv;uJurxPfS86{o+!k(`ICYR1YJ3V5^1))rJ6#@mX1ExRKIBWtNHI9w zo$bXNEry==0z!_B5rd-R_BhW8V)(nPJi2C*7^V*Uy7Kr`F`Qogp*NW!27iHZp=h=k zV*dV`cYeMY{4S>CP4>h8@4wl($zKddMtaY?wL%PIBl1K2*N9=Od+aXj^Ry*<0*z|8|J%QiV(wzX+h`bMTx=ptYw`gR}8=RC3wC+ zAciG@zU!jn#BlKY$*3_$P+!~^ds>f+A%1L+psZ73IN#NJ+|F}iF!BDU>v>5G|2nAd zHYJJSQ2P@TpQni7hrDw^F!ZzneBB z%EiF-8k(|AEe3(GXUUWXG3d=dt8i=)LryoTL6epQ49ue}zvxMz#4gqBj*$fBuX)iL zV;|% zr$p1(yh{R0`YoNT-XnoAT3yabqVc({hXuS?34F`yIO^{qly988|2yH_^?AX%uOzTAenX>K9*+Bd*};zQB~bkAT#ud~ zCD6*>l{2VB0!wW-1&l6}fN!7A9j1Ph!08#5MSeddaQ%18(ZFg6q^(Wy+g&4pODXSN zj|n7jFy7xPMIwQ?32z&pC?sGsV1ukkBY~XH>GGOJrhkqU%0E}dnQBYHi{oZKp=P&`?3{XrQop7a*Uya z6l_l0PwV1{Y|=P3pA%_m%z zj>i6$9xU2Pr0`#KeRhl#8arLfeK<}E6Tf*_mrjtvIql9Ok@>wksJkV3m>=l7hNC52&^E7Bj&k-~Lb!N%|Nq#*63pVhiR3N_6u z&N%v^eyct8CN7b}0nc}wz@<{~iS>#*wph`GwS=-7bw{x1u_5JpXLxLB(IO@i27rAO5yJi|Fkvx zrQkU3-i}i-QZSx6@@L*bl+)GiYhAn)UdTVK?sEj?OrIa)cT5T!Z?=%b2~tpu zQVN^Zr#GojOTl-!F0aoyDWr_wrSZQYg}J-NoIi0%3SErj%ibqSF+M+gUHhsO?xZJM zd!$G~XRXjWI2Gk=i+`mcx9w6q!3wn?z~pM6b9#-x{fW7!ZPlMrQ8ov2+Z!*?`@G3c8!7KCMDROv*D!g z7b$Ev`!7L(_Pxa0{Pn98UYb4j>GoX;b4rfASXCv34RM(zH-DlYH7oWhe@nquQ}N03 zuM{fd7iDv6r10ixVJJQ&9S!( zgneI3+~th@&5b$uwT}$a9J*{B+Fu51#M#=%24cVCGKWiCWiYgOoXylBGB~{E?uoQv zGKlAl&Nd!__n*2QSvOJ!>EBBX-;a{P%A(aS{Rqy(*;O7tMh0c)!d6JfVL!%uBW8Ka zAZYtZ@^F$2UfR~*>@?-SdC6WzPLsjev9I_1n=XT%@~6Tnvt%$^r?b}mIWoAly5EaV z^JUOs@Q^-YZU%C(~gOlEz zPv+q=n7P8@UgS=chcj!eVwVgWAKUI-9VG*!xfg>g_R8Sr%MsST(K3kZ=(;@j0P6Gf zhvgmzaopQ}w(0TsThc4k`G^cE>_ZM-K8F3-otR;FLIzt~iryVNC4*{%txD6gGDyAf zu6X}>8O&-67;zrePm*^h%Ha2K@%g4JsOJ?r6T?%m{gnyV>QZI!X_@ZR&@|MiYs_BF zEg5Y5+}1-iCh`{7Y%mKlYw>eTkEUuWKe%#=!Sj;GVnR$x9hy5QsJaU|*EKTe)?Dii4KkRt{8ZeRCTzd5@$iOL8Avy-nW?KS2Oo^?&g;sdh9tEc zW*`SaLHiCxMsipssVZM*Du*zgX^ZqN`t%AxX5Cwo5!ImFyQdaKz{4(>6@x|e&)Ves^u-Nw4eVfL=Z zmVbTaFz`uj(!l|E|NA==S64X{i-}x4oF#{+{1;OX&6Pt!ue)nufgCj1Cevk$NDeB$-aRjDmP1veo%`Ia za#(V|)y^_h4lByt^7FUL;kU$mQuIza)IE}%9J5;v-JMTwneO-+c`AAM|hPm=mZU`=JWG({i|W(`9}B zSve$;9+`(P$YFTjh|hB_qh0Gyzv+BM4udv*j{B1=hmY}xoNlG!^+4}a+tTFlrl|27 zxh;o9=0R$cbkvVy<-f9fa#+{*`mL*(a_F(cZF2BqIZR0WnoXX`A!v20iDkANCKl=r z`k5n#TwUS`InuN`tchp>AkB^Z8UJWwEQvsoihmP)Op#aOpBhGWI6%cC^ zS;@0iz@n?0wLf=Iz~$g=CfW81nCqaayU|4fdyf0PKh<3UYlIeixsD3Be>lAR*4_$u z?4W;jg^L2?wa(_V`zau}#>s#DKm}Y(J$7oin*zRX$$QvmC_X>f^+i{A1*F@I#1qHc zGyR~^Xay`8R5^{~rGTFnQ=2qn6|i`&DTY1@;P2Sqo9a0$v+6`j;(Jz^dc_MVtz- z8GNnr+ZqL!P7S^LW1RvDWBbqf6R3bG-$w9xn-t)(#dMQ&ivolN?P@gJ6fpey8(&Vi z0_x0<+%Sqz!1c-18kl>#}u&NAvtH;2?gAmG4I}i(+W6$?!Q9;7WV(b?_5;C-9c$HauabLasJ&at}0-J zV_1#!ngRw6U3bkeO##)D)0R5iR)Fn)YD-r@es}qcIrkNiuu56%hGgl0$Hj z0yqP<-#%HQ0NHfok)xE*t!ZCt4UJL$V$Za)>!JuG>XU$p|VOebjV z|EGWlRlO53c<5h&uUviz6_6Kw>ZPev0X-*n^B<{DKxLP*)^%#MC&v{AU}-F0b<;-)GYcCA*G*MIyYVw@$Inng?8F@g$7d@+Prudd z_dF$d>8x;eTcm{RT?WkCyF>|RVr$NRS*C=0W$l&DtCV0b7`q}OKnX8(G_{}BE8+g* zM_iYUN?6=tFmfueI%`}l;kVn@C=ix zRUeda-tV(sWw8?G{C9l8`3~w)_@P`09manQob+7@O;csVU;j{oETBujG1W>)=)m)N z@elh4egg27VA%Je@e7d>G7tQ_=_OM_*anN8xk~g0-*;gX>(CD@Z(b>ER02`?HqUBB zKS-S#{Y_g1zekRSmHI01bb0D5Hc~uRRe|b4e8Yzc*zO)= z2KlHUam}FH7SmLqJfKybK0^g_+Z3b%uk9c78s@7Yz|vB6$xj8QDOrx*OI5IP=#FEO z71-WIM#5&m%J3`|cubG)$H`GaaPkhF{9F|r^4Qex#9P$s zzfP&k3REzxd0^CC(PISOi4V#wD96O|!8ZLRSnC;-AhT`Q``4;`v;Ff5g+&V}N{LoSVq(jtj zxZ3?twYwV51r%H;^HhV>=25#sq6Y6B&(yERso_r6S+8dk)u6M_R`Xzr8UnS19nz<( zp+Zr8_10`PY`>a)F>SsYhSrYJPxC`LySsf)TdIaN9y1+ouT(=>dO&G+q(YRu#8je_A-1sk24JJM#e3bjt(7o9B zoKB1ymJFy|ZX2(LV>@>w^**YG*=y4mdYn+hzM+%je9x$%uUIWcVyPbTWVMoZC)-)$NQ56o6R!S(DwQJL=8>po=BxPzAKd52Hf@NFJ z|4_rzEv8Qjs@32(<&XL}|e!Ee*u{9vB?0r-At^7F^3U(tw+CV4$wK26}EjnmgWF1C9&*ZXal`fh#2& z`xe@3VEDB(Pn&KU=sHgNd!eHSYO}j`yXvGto>ND;xUUAbh+g=O9HfB}77hL}L-2Xo zwWe}+4NRBa9^7Y?28uc;^&-49Frf{|(ZHUmMZ2ARG|)SDh|BJ28W`&uyt#a)2I@4M z77U)JfzB@WzYj0cz?Z%AOvU~h7{ATt>(rGRaOyf@(wzVegigA@*g8-H9YQzbsex@} zi$g0zG~haQ&93p`8psXY+5gTi4fNR)$F<+9fdl*hJ2dch;F2NggBm#fTPuF`5e)=n ztvptqpn=y3+s98iqk$XB!Pjywpqvf|9=RuJAf?0l@bnZ7tUUBd)Hh87bJh24QtxPB z%nl<_ulpJ}r295J`H=>sLwQ-fvoz3S)9Q!UbFly75dj0#0Re8sL_OPLI*5gXNK5NCg>#BW3$jIR5+*J5-X%)Qs7pF5@w?#|X;67F3GHy&xnx=pDA z2Y_$jet7f>D*Yc zF|ZDv-Wa=NSa2O|+J9nr>$W9X-Htk#Db0&L8HxR*RqC(e*1`9tcB_ZR*1@T3 z@z)IwV|(A_pS~ThgHK@@ztT_F!Pl4zBN^^9+aKVVTp7HAH8K176@#yLqf3BYK=IR+=4)3?|KPxd zp7G%78UL-G@!skg->shU-0B&>4coWz+K@BZ#%HT%JhpnqU#n-lwR*-^t7km5dd5$y zXS}p}#z(7XJhXboKdWcFGwi30Z`QzgW(|yA*1&jW4UA9Lz<6X0j6c@Ecw-HWFV?_# zVhxNR*1&jS4U7-gz<6K{jQ`cZcwY^S?}fkHcwP;R-_^i)T@8%S)xdaM`2RNkRs-X0 zH88$b1LJ8mFn(48<7G84K2`(cVKuP)a|7dDH88$a1LIjWFn(17<5e{k@1BZ8BYl1ZQ}h*Pb1^|G%}tK z&cBV{)5v%|jf~IJ$ap-BjK9;!csq@ZuhYbMI!%n9)5Lf=O^lDz#CSN^u8n`w#CSJN zjBnG#cs5OpU(>{RHBF39)5Lf*O^iR&#CS7Jj4#u~crs0lAA{qx@nV`7AEt@%VDPy% z{!0_%y)-etOB3U{G%gU zX=Z$oX2t_)X8ez4#`|bye2-?v^Jr%Lj%LQ|Xl8tlX2#=aX8et2#@lFSe2r$t(`aV= zjAq8mXl8tjX2!#4X8en0#=B@{e2ZqrvuI}gie|>Ez~|cd6fKNL(Zcu>EsQtO!uS#` zKzR}^j33d$co8j(57EMS5G{=V(872REsXEb!gvlXjNj10cnvL#&(Okn3@wbm(872N zEsU?w!gvZTjGxfLcnK|xkI=$+2rZ0%(872JEsSr_!gvNPj9<{gcm*wtPtd}61TBm| z(8_oNt&A_w%6I~;j33a-cmb`9575eZ0IdxFZ)LcDE5rL+8P4Cz@cmYX>$fsIAOGKm z_mErEK3~z5`ID0F@*IOB`-pcUw zR)(W9`8gayadQr#csYkqoSZ`_KF%Q&7v~U)hjR$U!8wHD-}qb`?v3}`@NN#FI5&q- ze49fkuFWA7&*l({V{-__uQ`O`)*M3dY7U_|HHT1qnnNfq%^?(z<`9ZQa|p$sIfUZQ z976GC4xuz0I4_4#e3wHguFD}5&*cz`<8lbaZ#jhGwj4t7S`MK&Er(Ei zmP05m%OMnxQ;Sgi68mYxEwu>6E42v4DYXd2C$$L0CAA2}Bee*{AyLjY{1LBhxFfb}!yB~- z#Tm5-#TT^*#TD^>8=i=Aw&93cgyM%QbZ)p>XYiSdT zXK53PV`&qLUuhGHTWJ%DS7{T9Q{nG6d`g>8TuPfzJW88T97>x|{7IWo+)0~Iyh)o- zoC&XO_>wlExRN%Zc#;mGIFb&b_>m5wxRDN_c##gFIFSya_>c~vxR4H^c#saEIFJsZ z_>T^uxQ`B@c#jUDIFAmY_zw2hhU@4Mis$GMisR@Iir?rEireTAir456iqq&2iqGf} zip%H_ipS^>io@s-iofU(io56#inr(xinHhtim&JpimT`lil^uhilgWdil68Ziks*V zikIjRij(LNijSba+HesaLh%qCLU9m$uIk*p&)+03Utw(6STaVB@w;rMSZ9PKs+IocM zv-JqgW9t!`zec&+^49u<=BxDy%~R_WnxEDuG%u}BXg*q>&^)w0q4{ThLi5i0gyx&| z3C%O>6PjPvCp53DPiQ__pU^zAKB4(zeM0lb`h@0-^$E=r>l2zE)+aPCjB>Q)gY^l` z1M3r-|J5fn@2gK}zE_{nJg+``Kj;&h*VQL9pQ}%39v9_k%irn~nzuC|G+%2#Xr9)9 z(EO|cp?O&YLi4fs|F%4=0ipR<148qz288BY4G7J%8W5UaH6S#vYCvc{)qv1EssW+- zQv*WtrUr!OOAQFklNu13A2lE}FKR$&J`~%x+D(U8#mBHnMyD;g4-Pc$U7uGx^#`ej2x>y!-%tw%N_wC>oD(E4I*-vJH6pYQ*ND)1TO&g2YK;i3pEV-1PS%Lf zdRQYu>t2lrt#36Vw2sw?(0WxPLhDkE2(3RgBDBubh|qdcBSQ0RMugUf8WCCtYD8$g zrxBrbokoP#ZyFJrH!~)*9@Ci6x=Uk1>nn{3t)nz1v|iGf(7H%tLhB!m39WN9CbXW> zn9#aK{N0w{GA6VR(U{PBLt{eg3XR$KvoWD{g2sf_0~!-r_h(FKeV;L*b$rHz*6SG) znujtbwEoV7&^kL4LhI>F2(6nlA+$cugwQ%T6GH3VObD%OGaO{FESyt z?#zVH`Z5zj>&UR*wtSEYp><&dQ1te-!Ub$PRESUdK@!C>u$^lt*=(GeYZD%m}SdF(b4N#f;E;6Ei~VO3VnYA2B1ePQ;ARdJr>0^AzTU z)_0f_TE}5dXuXCxp>-MNgw|h}6Iy3sPG~)aIiYnE=7iQqm=jtDVNPhhgE^sf4d#T_ zFPIZrr(jNKJ%Tx*bqD5z))!b1T1Q|(XuW_1p>+Wkbo~2()1)=L^97B5?}IN6fL$CfNEwq)_JC5wYC+4{9)ajzwd zcP&|*YsunUOBUB!vUt{##j%zwezjzAt0jw9tyr9D#o|*d7MEJFc+`r;p;j#Z#CC0Q zrxlAgtyr9C#o|jV7FSxac+!f+kyb2zv|@3i6^j?GSe$6Z;zKJI7h17+(2B)@R_y(b z?c3r$D;Doru{h6)#dp>$uCr$GoHdK%tXceK&Ehs|7Oz>eIL(^HXVxq(vu5#_HH*Wn zS^QQ^ybFzvE8arFW&bRVVToWFz%=Hn`>QPP#m>}3mb3z8_|0)7bH{czDvEhz|TK9C;TWEMtv9en|qE6*EDZFTwTkB z#{GwVPcFjoA`;g}59PvyF#+GU26ExsDAV(Qc5}f_d&PqnIBxrQ?}Hzm;)3VOj4lST zTv%Q=vmnin3yI2t($1k=xcBs*vAri3fDVoa z*iL>U&EF)I3*Jt}U#73%g3n>2b{RN+p~a6)clL5&dV9BZ_W0bK{Z{RkZ$SByR*j39 z#D)8ldH({gae=q&!qjimxZoEh>t}`Y%QZaA-;H`a+SSHS2me2GVcw}NJ5Ud!B7YBE z#f2Mj&$G2q?(Bt{Qx7h4VeR`Wi7bZ;LR0JZ!-BY=v)$eP#C9$mbMEV(hWaoINE&<$ zpBq^{TYP>5wvQWdwD%bDUwPcE;fTr9;@CbF21Z+xc90G23O|yLc|7Sgc&AK>ePLG7IjDb~@s}a<|rb zE>t}ex+b3BLX=U&lvsZ*ylwH{)_V>Y!tNI5{NZw;wQyDC9Mntx$@3kG(2ljG|0;ju zI9tabE14L9_ET>cYoWfp@5kPWUWRtgfZ=m-9)cS0yRT7SQJYE^be+xx#mU<3n*ChN z>rcqc--UjtI1v9HpX*v}I_n4eT`}i-!xZ$RxA(?;y^Hw@pxa4r;$ zU-Bjee-AsgB=r=^@jUxe>Kl|Zb$plbyU;gLdvduIj{V)X&YWR>KdJaiJyM*dlo!jw4%<@aH5KMub=E?M+5~y@(p* zIFAbhevGoqzmM__h@I4ox4NPqZd?3*|HEMPo5)kOrCgkEnm}cT`naDye_Fy3 zF6P4>TOOjkqfdy3aXzD6y0-+?9zuIAzH@$05*KEzOHN8Vi*k*y=`^913q!(37HMI7 zy%5a^mu+0&n10i_9LI%2%I6)s;Z>p2LFL)Zh1tVw{70i-1>RCF`-|g52ljXPjq)~K zjvg};`_E{SrWwEH!tdAbcXYdl>(6WA6JwNn#`|({!WJ&{obmg22V4jD&3CQ%`h*L0 z^3ywZ7ja?19G&|QUvS~+$BTt;%TXUEKV7`?iVG2CK31OSwTO<&`-jK-8zPHy*a#XM~4$QkE89^7hdH; z{NGhaOt`q;Oj}#&_Y&uSWlq0^h1g%gPN`QG?kDpcQZJw!$7)9eHsEt9W>xu7IM1(# zpUe$Ic@7Wuzdopv3wy>s8t#StoDmgdWvswlY=i;{i#5 zcb=htxNf`FffEf&rrK;XyTS#_1}p7J*#7JBitCXmM?=OUoqag}V}@Jq-Nf-v8ZG_p zb(ag@T|WQ1`Uuy_Sz*PMJT7cs=;0n7!iD-RoX4MWpE2J$$+t#38p_Sr?mL2VUYf-< z_6_C2yP}KdeR2E_MM--4H@VPD`r@8tG#6eyRwVa9eXN_j?wEHBuE!j|7hGH?R}-fs z_DB02>Dc$nAoTkKkw-rcM!C8-mIis8<3f#kjCTs^J=CL5yb;>fiPM-6itQ!Z9gTnD zJlwlyjZ@y>!rh-=rq<&AG9T%~C7p57{s&q#?svG+J*^GYwSoiwA=jgYe2M?Q~#{J-}>A!#^>?iDV{hoE$ z|Cfn5Z_ROi>KA+(zZvIkSui_^`xX7ewKg{CAztkZZ#AO-&6w$XJJu>1cq{iGdyeyI z?5q2A9InrsVdv9lT)_86nP*ufu1n9Dp@;mVaDCJ~YrKa3(EX}eV=~$~eS2m^|4%5l z?fAb+^!K%qwtKQaa^de$y?f(6ee=(B{ebnjzMn;1YeN6hd$_;s;TtaO z8|eKd?;!5S*DotW@VCjxlvsDPpK~9-+2c*4f!|T@+NdzJm*bFxZK$sgmsgifLw^}< zcl5IW?KkdVQt9FAXs4Y0kCM^e6`%ZE)}wyx_k~9E!g2C@RQCRM4E^hblWQvalY5`~ zbLY{XHxDh=8uAqFCLe3|bu$-EEJ}DZ1NSZCg-5;J{`2pklu|p?n`855juXDml_K-g z)3KkNh5eF^@On(=;HA+hXS)Rr?qR5>;lD1O*Ft?-+isP|mT{q^+-+Se%Gs^ym*JKw zE_7~M)a5}E?oUBm$!@f>eNAbnUO1mJZE;6q+*dYR`RY!^ec8t*dEpV%e~HbwXnWLO z%9M%otWk~$Nz)6rW1K0VJ8W`LO-+QGo;rV9hH{}*# zrRpFqNWBL7JV8B+2QO$h6Wbr@wD_D4>dy{7j~{{aj(=HN*%`-o`E%rg8uy(^aW$bH z=wJ6&R7Di8!~II;C#%5mhJ_?HbhwDu^GhymjzPbf%*~jB{#-n3zTZB)aue6&$DzIp zU*!%8UXJ^i;TvvyTxT~8a%XJ{#{SGM{ydBFJ$hC(WWxr0&!vxl7>f4wOLi<=hwsZp zcdy(}#r5p&FtEQ5eov?y*Q{;Eb^px&&)$zHkMGUUG~O-z9=NxDNjDs)4d06?vG!wY@%`&NT~GK0 z_vL=QKIH`BdwY-M#+kodnDS@3#~Jh|y_@SovnA*U`|@6U>)?CX&{6mU_43^C z*UM4MUUhhj>(tp{>mRf)Z_^03Ed^Y-ed?K!?lSDZeqWU{u8+k&!=nHDzPS5kj@drk zAH(LD5udfVE)E$5#P3HxO}k$)>mC=hzixM3i0_@IpYPw6;(FJT?0UBmui{yqj-Nxj zbp5fiWNH6sNF8&(V;#;btNZ8|6L8$EmJ94B;W(dzmWRyx#D$O6@x#WWKCbr^?!NOC z_m`aLes1X3&FyA6N%4E)+slEXySU%}$kb{1xdY`49Y3@y?sv}%bdOK^Pd^53>k@zB zJ}&KkIu-4*X?$wUMx1xhPOt8-&~C*Ce0Llr0!qf-<5~&iASFY4eJsO zB`v{8d$9kI8xGx^t8jlBu;}y`^s|{o18iR5eiJ^{w4lf+8V+tPUEB@*F)q-)@bh-m z<0UFfgZypU}!eTz(?;H&gcI_%NaQzE16VhkY;rqj5)QUTH(a6i%T9Q|c^Ra*Xb5ck& z?8r}w`hxQ|a9(b0V-pPv{zYwHVi^rSUln;;3fzyfKg4+_qJ3w_?>=9F{=!|H;fVX_ zN5?B?#(u^)Cd6!-_Iup_4V2T3ouVP8UEc?xBGiN9w{`w+aG$S{u6Urtb?rHR@n5w2 z%$-eN*LIDD{9QZpUiFNIz~bNUPorNSIAS!YpcU8erM^S0@H-);$9=nLDtsR_+}S^3 zVKm4Fto}0`?YM8xpoltr4=%098_*f|@$zP=^-le0=u^`E$PbkJgNN+KaaxxqWB+!1MvB@l~Ji-_R(;C9k(K? z9^Xr$f;U%CzO@fuc!cwD-Br0R9bJmw4WA2qTf{?ab0{`6g)Op zi0_lVJ+5(4zk)MK1MTqre>~rDPbbvtp-p?=Cgc3POE>Jld_5Yf`#cU#!~NjvjjKjG z&`$%JFZE8($M?;T!(foVl*1|N&-f73|E?>BJ&oE&!+PC)(G57i+V+D&(;IN#y0mL;elr(-Zm1fU zivH1OZ-?0-Xg`k`+m1~}`THr1%Y|9^eI(k}PV^Mp?cM6Ubb2%ddyRL_a*u}ikm?6p z(eBl+qYmuG=cY8-nmD1JUH+EGTu1rDuP!?gw6h)GsM}{siNlwa4r}s&U?d{!5`_G~5Z)T|W}-)4RC&=Ja0Cu<)(p+c7x* zli_EVZ+eCA`+*<2Y{2h1_xtbPh8=yKiYFTxaUgL-}gVm=RL-8 zfBon%>zxU%i;3NjI^w>#BFfawAHUCknblYKsm69XXPi6Y^G9uKHXEaSp;3RVXW+aK z9-lv8@ekbhbe@$g#rt-}*Di4NqQQRRKWLBbhX`Dky5l;ml8rV#j(+VGA-t`|bvr%Q zE9YvzXs|ip__Gr2qMl%B(t`U!{2Ybx16+Rv@eW4)Fz$%t9^QS_5BIUSqscyK=j=$q znZd6yKIqtaPj&&m2fEZkf0SQ6a+_lwzV~{(+i!mRI*xZ>#G4k}{}=oY?Ku7$>$ zNG|%%&@n>0$+`GGSY3ZPvIXClPwvcfMSbNsR8<&ae=Q%5zL|x;HCHTtG~@e@P|1KDMW@@ldQa z+FQM#=SA#y{>JQs-LT)A>wPRLQC}BNc&;w&5e+wg?#!t}|8uy@_kN1{j2*tu`u9K7 z&&=$-_tDOAIsGgaVn34}Kl9#<{v5P7wVRO?-!t;eB^`19UvqPgs27gEapk;ynfN}b zIrA-}FdOeTkGal4Jv>}!)?%iOYDA}OT+#L<+siZ#rq$IPwv(U?a|pF zNPZvvO;4S_*%tj&5%hd!IQmIq%eL75w$V@~DIJ@Qzk7t7^3TD2XVvxMtl{{+%+P&3 zV>;^b-?F#guquR4VNs}5}4yAH3`zC(So%5i_$n=)lM z+RLCr&Ov|c`k4)hSYl?@6U%=jeF5vE9^Zr2go5)Yk#g6l+|FQ)u$u@kq_k zC)A@~veQ@o-9p?u%1WG(j}yCs0xE#-Vb=N2ez>X})7fj#5=nH`GgoxM9^-&5VH-Gn zQ2u(pig(M{LMi=NHi*Stic|)ThWzEIhZSrng)D+4E^SrwY@k5 zd7iWMG#7UJ8U`&`e2xC$l0x&oT<~aCe^m%WmV4fHz^=Y=cP{y%~1 zXWz#n>wvRsx~GdIU%`%}mETyAzjOULua96q9=7|`+u)aDN=Y6h#J}v|BN0NLMBJ7M zdCLd?ggiYEZ;E>3<=(b626j~aePnGm2fa6E^?^wdI_T=BU%lYP4nF*IBKYyrj*P~d zRm?xKRQfcdOn={Myc#6{X{pdEJn96QG{gSr_U zcbu;a`wG!BjJ9c$h&7#meAmH~++HTCL(rLNIb8i_;N#-R?Gt<1p=bY@tPdkz{^2JQ z7qOqHk4<-C7x1F5z~OY5MEqFtwyc3(*)&jl(=~@g==I1x|13%(NaGd8_u-f7Q#BXk zq1()?>n=(nzi&I<3l49Bo*Z6lwgAcUoX?~esfCek<= zVSh<+(T4%##fn-`YgQPE$T)w!`8(p0zI^w=3i2eY;mzq{=wbPfXV*s%mr;`{ii+>x zISD^q_BYVgsge$?$n)~N_IvMPry6g+O9gPA(m@k>8hVF$Y2r^`Gx*j;_`)D?QZFhk zydAt|*8S}El@Y`_tbv;c=lJ3^%VZ9}&7Q8ju82DPV{42}9{g~PS+6(?cu`GJdRE;A zT^Ie+Js!Uks;uQS!a0jsM>;}1NrbbnUSk&Qjof@QY5{Rk*{ROy2tTN%dy94g*R2CD zT^vz27*%DDR9Zs6J*@q0jlLzE<7l(02KFybDUbe(epWL0itIP&1g(C-fe)C69X!S| zin!~wb1-X)lZcyKo}TiE%eed#>wz`sz6b1PywIT&N4FbuNWh|6a9<-E;!;d-{edUe=ZNGFTDePrNbIj#D2rR zuuCgR@Y4wY**q`g+uyk2)pPLY)h`+j>(Gyrcb5mA?LeQ<-YOx7?`g9Nj)-8LvucKB z>)>ncg>)B-10>??uOKc#=&5u{kG4ekuRcdR*+&5VDz)K_0PItxBhh^vc{RJ4Bkv*n zGj|}dE-jiwBw9IU4nbe!WL`S|VmXpH{JO!(1O40KTd&N9fv1Z8D)_bc1c@+Zix!ZUCK08>*7rD0 zlL(^=?%$Z;PemetdBe|^2d;MSz<#~=vY3SMx2DS3abEauSKD?qN#rd_@33nIaQRHb zVKN6i5G43skp?(CcJ<$@WdeTa-cnqLbM3u#x6}y!il(&`>f|R8C#^VGlaZ(M92Z{P z#yL2LIqjt3$LUiulaa8yJ%l~}8}{?eNwp;6NyLL;Hq94b;LlI-spaU8ggP&VF)X+yz1TQ@F^H}@uaMk816^D)SLR)f^WPFzl7K_lL&1l1*s6!KScK7BK@An@+8m(hIzdc1pN?~Nz8wkM?QfjzhQKkj%3yy>)F z49P?tIq;X1pN{>4_LTi$!Z}!OPf6rqU9m4opNF9{wqBR+e0dkV)**SRy*iT6@^${6 zhP)WbiCxiyzcU%G$oru$t2iF2{0jR8B#R|Hz&Qp){N3}n0&m{=w7Y@F?AeF;rl`jr zg*DOb_`Q?QZn|WgdshJcejCK==K@FeJmO>-dVbjt_U%0FI{yQDX=u8~Mi_dSzwLhN z6!Kk7QDJETbuEvTmV+I2awnu6WBN^4_S6j=*xo2H&%2*830OO2;oqss~?I zRUO~P4<39m>2#?R{1W<2mD6Q|XL*Px4X}4eLt3dB@o6~Sk--9gy)bpHYZKMNr`iv ze^2EU_Zo8xuk%SyA3*oAGq3+P=&s|^r1-*6o zx)ji)%k{WZ8bCziS13h~&o? zPaosCm!6o81LCfC>Qqn&;-M)+5th~sKP^;bAIG}W$p-&jLtdCzQ|dIHMBLuTh1^0M z3jebky|97D?sD&fN9gZ#jhB7^Z%%$Mq8s23F;1Gk8^Ft#%T&9aVE5s|)|5Wr`NQ{` zn-35t=BtPJk0Bn?Vgt(@z;Vm7-z%ob|1!@&nWNyd9k;*#o&XQlO0X83!+mYP+k{); z&rGgeHATp?13R})&>=1nH_CkN!Qa0f=Ou-KAI&3Oz3spci)q|vZ#VSoFBUqsV!hwy za=kM+m)m|%Y0>-0i#=Z}mw;c*Q(IRkfdje!2!So&)d5S10DcNKX`eF|IssHS#8I^0k=yB~P%a17^p4Lg4~I%%$A&P8i^TKI7fiKx@or8Bqz zz4cAVjUBo{Wq@N}H~i|Z-c?!_hdJs_wZ4a_L%)<#lntThjHAX5ZoP*2sq<-O2i_>c zwVl7>whrpTxibt;z`q|HrgHzzfe(*%SoNUKryMleejtEEwCbMW;*-Ssf!X&D!=D#M zq;d91fn(EL5{nQ!5d;L-`5akM;e zW8Su~e&5g%b31kV=q^9-jrPC(i6G1yKYh$?R)qb`6cKTIe393oi(+e!;6J^*C_~^# zM0(Kr7}hH;C7pW<-J}~I&E%|!`spbVRu8>X6_genj5sE4_s>&>F8P<3-IN^_R8#AhK*sb+gz>7u>b@1f3u3jG!alL4n+Z_5yDYf051b@d@?oHza z9I3Bx@L?U+HV?Fj9 zdgx#3tNXX%$C$;zsW-sKhlr#~FEz}kstTf29pUGDPv_ib&<|AYZ0o|g&zo?&e?k3B z++`%bj{9%7oU^|lP9h%GDtPJPoX%?|H=gmwn1wC;0bK zw%gGR;5Yo>gJg?+GcG za2Z0k6<&+z1in~a^Tyc04!0@Q9%e@g!rB40IU% zzu-fN=avl8!EgHDqoDRI!E?x)_9RbN^idNnQeT{IBOcA=4<5je2lUx%^o+1h#lJJ% z&|ylyBUn6fZi&>p33?Y$r@qO{Glzk{RQC@5#(salGb~Hhq92obY4#jC`kGosoO&yC z&RXMFe&EiOUvF_gKJPg_GE@j2=`230Iw*yF2_qF<@`*yd9WdwMDVvGD*;`~Qd( zeMS6aj-FT*h(q7ITE^3a@4x#EjoTpJg7z=%6_Lly6I1Ge&^uP8Yo7ZUpa-q!)6O8D zS!6S$bXbu;KKEW~?1evB`M2|ykq9c*KG$|PJio9zxAToF^mXnlRRaDAt;nFDMf_ya zZ{765=c#WC6GnoOpAnOYIO-2n$6Kz5Lzrzt?Q!t> zWjlLCWeN1R%jN5$Snowk>8=#?6DkHH`NP05-9r79NAr08(ensm2=AJ3%S7;o4PT*&fPqytC64}_f#L+{AGRi$8~A`yOX_*pf8yEHk0 z4hpO@D8BxEz>q{7+_FSfu@&{g;=ru_IG#@|bn%bi`(H(!sRMNAmmelHU&MNHCa*!V!se`X}cS_|ml7+uaDLEs>`TJgzqTEyS8fyEX6sn#v~8-hB$e)Z+@c7F8p z|6O;=!~KSP407wF8qjYWlCLOMNMLvEvm`MctLt{}_ zYUFXf_|bpRb1SA}$6MeJ!R--Rmlg^=p{r9{Qs6Ge1kX}LE@tO^B;$f>+Kzy9t3T~$2{!sR7et2bq5I(Z0 zq;3r76v}B0fc}jz9n4>|V%*Y`X`0Hcz z5w0cR#5=jOR{1S-lvl;fCFn7y>kN91T+ko+NpDIJSKXQWjWJkPTGyam5O#bpk~y9O z|Fp9mmY76b0+Z%qBheoO3smO@LB~?`2xN;K>b7@-nEY8#9 zXLkQK`jk;tnS+_QN_(tX{_4b>EKe|y?g{#gHNSjc;Q7Q8JHblq!*#9Ie&hMS5WQz) z0X}CbGYPl^e0^g+6<`HFSo50gp3lYeT$`J}T6IZyPI9^QBjPiA_H9of;`?87#h>^4 zpdTKcuIz&ywieAI8o-?^EkoHD@Jo70{C6KX^5Om(?Q#nG@^EaG7UGk0<%WMA{FANZ z*Vd^4UPw-`+<5M&b84pfE%IpPet@7a&XE)A8d}T=e(}kj)WSN(c*{%zxEJUt?z`g% zzPVgHWYd7>n13P*vV+7kvPi_L+(P0{;G0t7xBJda;Hg7D(--}~U~salA$Vi6 zsl2B&@b`TzlRZWWb;Yo`G+qVIYZ_+;vkDMLZ|jintH4<|SJ1~_sCUQSRl7N(zANuN zbq~C!+`-IN3!c$_q)SRiK54HV8VCf=J65;4Zs|uK-(Sw%55HQ5U2fcaj6@vM|L(kn ziA4AgdJX>tuRXQjUa$i^bkR!wg<$~fFdjN_e4{Tg*?i6mysZ6NU{@pfDj=|;%5V(x zhfCYt7=fcZlhcgJF9illwm$MlHu8F2_2wVveJ)wLzOm} z!h(KFuk!Me74(Htn%Jeq1n|c7K+Dj*m>1J8(kp=1-PWru6mjlNPj3E*gB>Kj0v99j z?T#Dn9eFRoD_^K}E?+{O+!P^b^8$K$$;2eZ7rZ;gIL>AU|H#+gV*@@$)N*gEqrTm> zKSj5w27d$$=_F%chOk;|_A}tyo{KpydQn8HUHY3M#7&srC0y+^>@f97tzJYu@`_&G zbpZPJR;>BE^U!ZwT4n_ep)a^TQ5y(DU*w-+4%rR8`>bbR`~l``myV0Ix{`>L-s^S$ zP;ZL2Uan@1C1H+!@%m%@f2YH&sr`tBwFi{3owI*&O%Rec@o@5`B+2A!Ms?tEBHBy{N^ z%1r?Sz!ATx^Z~>rgFC2o4SBOyrf2*o;%|G@U2zNcmt#HG(gL2*^t`j8fVk}RT$#HF zel{{XpVt69eKz^$E{XnEX8s9<_Bj$^j29QAUZEbJxkx$>|1)~9|BON&ht4!eZ#)m_ zKJ|3!yj$g1KZ&8<)eT$`!r$5DBg`VS(8=E?bq^Ghh~;(C zHO?)#PSbSST0w`N49U)f4&bblm@0!FKcK}SVU6?Jr*3iHhq%P1kDZRou&F^@TUZG!bq$}FfaT5Mr^Bw$w&l!}I<#i44yUh*GjM2ai{RHV5@LbY% zLqfs_d?jJeaJ2=`8$33-@qo8wlRHDg;76PA=l6QvlJNgmYG>y`T+WK&^>Z~ULC zP;ab%M1PlOFx~?wDzO@{BN5f}XXa@9F@H1TJiGCp11ry|Q-6O04|y!#)7&xtHSgS2 z(hFXVzB68r{514%Nh^k*m|F-B66OF?A-z&M;iz*3+k56vN49fJFD%s}|4-Z-JLm%3 zO*4-2f}d==7VKkzx1$$>ZnfdcKc>j^3Gub=7eDYAJT$}TaQ_4F-mTADBzqQp{ogx3 z%pB1dy#K2BP8|5UnUQ<+DV|pyYSCG^3%@5Fn&yIC3`HmYaUqVEK2uU#-^V(81cE!@ zkJ6Cf{T<*37jxt1c0!mRxn81T*CG*bX${QHv!Mqjrw81?<3Ij(8KnWwVe>I6=VAZn z_qLCjZ15gTQ$%D9cp@ZaVY1_r$Ea${tS8`s}g0$(XhL0<)# zaY?}sOC63M9D#Eau^zf&;N^T2x4Jfd|M6?!3sdlyc;7D}9~J1P7`?F^@T879C%*&m z-AI=l^$zE?NRCf`2HZFnJ#EUNAQ9U})_oU&gU#k|xja!1hKFudbHGo1%Z=d{sOR-H zG1F@B=VYypdNC{TmJ*jP1O4`CZ6CFi5&8`mAI_j|)VWWOXP#i6$9apMv#_(>Z&R2B z_`F`AS-1%Nsk4-(T2Bc+Tb4T|27f!BD%IS}1|8|E(e8)q)~eqZ{BVEPu9%!y;Fyyw zOP;d{b&rQXvwJ`0-Dk}Tx58d$k(E2AabKffwf8OfM?d&e%N6jSyF(4{68OPCrO$N* zy!7-%@=hzPD|wRp8Xf98?d)gXMt0P(J^w}b;r_n2$xm#N&s_9WN)BDn0Yx*f0^t{R z{_c2r#NBXShMEEO(jd=ecP{p8sxbdd1^=xZ)Rr^?Pt#?e)TcLtr|o$~YOTRHwTAC| zfD_ZZ%Qm~gGhMRAJU(`)b2DvgUE6@m!&-Bnpu0|8P?;rxzZv$s2Ij+$^?O2nCE}1@ zSBKK_p+oh@{#1%LL+7g5e`6hl{|J^R6X*jZuf96kj=KE9L-FNP)LFlG@-(y1In{hV z%`}Kl$LHpjwV;b+Z4TEz z2%T~{cc9l0e90E_jx_;&(!W#eDp-G=a_n2=R>Y+zUXT&{oGf#7)W`42QY(CJLXRof zUdm>ogI>FQ!I2+2!cR0;p#r#1e3V}<0Df0#Yl!0qPT1oOJ3YZCbiU$tldvoB%+-;P z1ma`WSm=bf5KXE}whqYiGjbyN=E(cw>kpm-C-JO)2YQhQr}y^KtI~iMMu|BytV5}v znj|3$9&MFnTmugD7B}&3%>7rA9lq}hKpk(@o*0Wk|C2Bp^6U+s4=4-LCn=(?Y2I8G z2S2#k8xHB>bA=1;otf}=mb{i_GxSFCZRgTP8O(8V|K-;r|1X;?zGTGt-+2~HAHsbC zQ)eVXV5d!i{tYGYhlSuuMT{!y!=Lj1-edl2O>He%B8cZkYUkqy;J=9>c0NY%)uBc8 zg@)It_gj`_tFi8B;^qOA?`G|338t>d?u0JI*`&QK(0Ot_E0NrTa02#^91Ws)$`wUOowawKz9Zycg-J_PqEj{B9Rlwgvgu#&F-amnn*f&N-x|JAvmXl(n%-(dbi; zB$?I0KW~I}_tMx#5#|I5e4UsSB*QJ3ayQ6YZe?~X;kMx1Ba@{XnM z!Ca~_^Fj>bd&PLWvX(aHttvk4Gq+$r-K6A64dnU5&=ZYr(4+B!te?OqDmf1xx0`It z$;VAZk*9GkOId8icuwCQK3b_3MX1T2Ta^GmbeDhiIK2;bC4=wV@Hy1KEA)rb!E^QA zG71CU`23Ql@@p&b-maX}YZstD?g{z+#=aK%bnGb^SbyIm+24rMlfR-j(=VbgRkiw; z2Odlx{I1(_D2iy0ySjPf{l)(xGkE8`@SI)tW1Q}Woz$6^-91rf6*5kGOGXi2%Vt(i zFCfnUsp`@CVXiBh_943f@kuD^n-Ga2oNYJLj-Wo7=S1E`9bV|QqpIV-5Gtwg|jxX?x?@}T)y4p!nyx66g4G6KcxE88@r+&Zq7bKLLcHC z@4IW&7I@9`aQy3odGvrbTSgCbm+nK$qIZbj!$n`OWW-O39CzYtW7V>LN;EU*~XMAs6ba z=+cq#o7k7LIPB1t>)?rHPVSDI`*?DL-%EG&3=9wbAvugzQ2fPo2K>f9P)ULvFLGwHS*FpnfvV` z^1S*@+Q$3gQCALC#sPqS$iWDaG8@!u8%|9d$Jh7EG@es(jVH!=A zCg{exEN^$j)wge9TF?=Cq^eGPH+1`4y0rLB-S2;bb$=XS{p zx`@Z$rW$8c=z&b}*$0T1M5?ENC-N|TII`^B4)l9173)#>zG|dKo(AXVZeu*y0=vs} zO2@~LU(!dbb(WET9^bEYhKfT+-hDQDejV?PvWndPhF)MhL6v|Qu<@l)Wh+`Vj6SihJl$MzKV6;3O9^ZpljAuhP52kZ7+Powq6 zx>Lq=QW4m1K0B@Sh6i-}vqJ~fao&LM#)F!>;1B2953WEz@o=ZqY}Dhlu-Hr*P59$Z zIURP^!}`wyNMDIC5@E;r+uwx~d>@@fH*AFYQ_yQ& z9q6Ip8wreUh)cRkwx}8ANejGTq-)T@k5hx$H{J(hYW_f*3x8x580!9n9+obN=X?Zv z!V(q*R&GHD{v9D$?X8`YB$ZF)m{7{JMFW(8+&rQV=UxfG`E95Zt z!ujbZ$JOHSJKZC(JZ#Y0$rkJ!5BH;ve&N<00$wa;PN}(LPH0qj(bwAm^R0oeEr!6` z8_jg4?K7CidRxBRhHJr4#wQ`1>jk^F?;vz~NSGiyBi8#xlUK6BguG;C2@QZBP0FSW zihn^roKxD*&xZ32R4g`byhlLrWm024YV}&x_?jGWnErNEpM^wl=zifigmW67mj{j-E`=Uj2vXQ`r=Hz>qqsZ&uriEKmThQMz_3zEp!2Cyyu|gVgROoWi+NjIJ z2Bqcp0^kLP<0{91i{dY_tp31XVZ&A1Xw2t@)5TwV!9LC{tp@vmE8A$1ZN*ktf17ng zH}cO)t1-bJc69N7a8ySg_Ym;GGrlXv z67#(Af@l-Gx8T%L6KU`e^P!188f{+ivp|L0%ze!92X>1$p8*f8+|D<|=R4~>fAt~H zwX04Bo2$a^Yzp_V36!6Xz6#`M{Yevl8UNV?~3)_=CtlS#Jhu=>2q` z>ih`k@e1)A5}>L z-XAB#E~(i=Pc!qhPeN}bJ6W(l-H20%4u==)IYFyOosPJ4#SBQElZPLyXU)~(Ro|A-5*C~_glsE9D-TMyz zD|^D>yrAH_QRryH=EJ?f$>z75v$}MskLrhnwU9^KZpBbV2*E$PggpiHcHNPCgKxnT zYHgcyw%tJ8_gQn7M!s+^tp5rxgX217sc`^8lF#$UH#i0WuGed4S9VWF8>%0GS8KJV53F zG7pe>fXoAA9w74onFq)`K;{8550H6)%mZW|AoBp32gp1?<^eJfka>X217sc`^8lF# i$UH#i0WuGed4S9VWF8>%0GS8KJV53FG7tRUJn%n0OV+yp literal 0 HcmV?d00001 diff --git a/retrieval/spectrum.py b/retrieval/spectrum.py index 41a738f..af38409 100644 --- a/retrieval/spectrum.py +++ b/retrieval/spectrum.py @@ -1,4 +1,5 @@ import matplotlib.pyplot as plt +import numpy as np __all__ = ['Spectrum'] @@ -12,8 +13,9 @@ class Spectrum(object): flux, transit depth, etc. """ def __init__(self, wavelength, flux): - self.wavelength = wavelength - self.flux = flux + sort = np.argsort(wavelength) + self.wavelength = wavelength[sort] + self.flux = flux[sort] def plot(self, ax=None, **kwargs): """ diff --git a/retrieval/tests/test_core.py b/retrieval/tests/test_core.py new file mode 100644 index 0000000..7b6a947 --- /dev/null +++ b/retrieval/tests/test_core.py @@ -0,0 +1,11 @@ +import astropy.units as u + +from ..core import Planet + + +def test_radius(): + temperature = 1500 * u.K + planet = Planet(1 * u.M_jup, 1 * u.R_jup, 1e-3 * u.bar, 2.2 * u.u) + sp = planet.transit_depth(temperature) + + assert abs(sp.flux.mean() - 0.01075) < 1e-5