In [1]:
import wisps
import wisps.simulations as wispsim
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import bisect, numba
from tqdm import tqdm
%matplotlib inline
In [2]:
OBSERVED_ =pd.read_pickle(wisps.OUTPUT_FILES+'/bayesian_observed_pointings.pkl')
In [3]:
LF=wisps.LUMINOSITY_FUCTION
LFDES=wisps.DES_LUMINOSITY_FUCTION
MAGLIMITS=wisps.MAG_LIMITS
PNTS=OBSERVED_
CANDS=wisps.datasets['candidates']
DIST_DATAFRAME=pd.DataFrame.from_records([x.samples for x in OBSERVED_])
SIMULATED_DIST=wispsim.simulate_spts()
In [4]:
def drop_nan(x):
x=np.array(x)
return x[(~np.isnan(x)) & (~np.isinf(x)) ]
In [5]:
fig, ax=plt.subplots()
for idx in tqdm(np.arange(13)):
spts=drop_nan(SIMULATED_DIST['spts'][idx][:,0])
#ages=drop_nan(SIMULATED_DIST['ages'][idx][:,0])
h=plt.hist(spts, histtype='step', bins='auto')
In [6]:
fig, ax=plt.subplots()
for idx in tqdm(np.arange(13)):
spts=drop_nan(SIMULATED_DIST['ages'][idx])
#ages=drop_nan(SIMULATED_DIST['ages'][idx][:,0])
h=plt.hist(spts, histtype='step', bins='auto')
In [7]:
fig, ax=plt.subplots()
spt_dists={}
import matplotlib
from matplotlib import cm
norm = matplotlib.colors.Normalize(vmin=20., vmax=37.0)
for idx in tqdm(np.arange(20., 37.)):
spt_dists[idx]=np.log10(np.concatenate(DIST_DATAFRAME[idx].values))
h=ax.hist(np.log10(np.concatenate(DIST_DATAFRAME[idx].values)),color=cm.viridis(norm(idx)),\
histtype='step', bins='auto', alpha=0.5,
normed=True)
sm = plt.cm.ScalarMappable(cmap=cm.viridis, norm=norm)
sm.set_array([])
br=plt.colorbar(sm, ticks=np.arange(20., 40., 5 ))
br.set_ticklabels(['L0', 'L5', 'T0', 'T5', 'Y0'])
ax.set_xlabel('Log distance (pc)')
ax.set_ylabel('Density')
plt.minorticks_on()
plt.tight_layout()
ax.yaxis.set_ticks_position('both')
ax.xaxis.set_ticks_position('both')
plt.savefig(wisps.OUTPUT_FIGURES+'/distance_distribution.pdf')
In [8]:
from astropy.coordinates import SkyCoord
import astropy.units as u
In [9]:
pnts=wisps.OBSERVED_POINTINGS
In [10]:
#big_skycoord= []
In [11]:
#for p in pnts:
# big_skycoord.append(SkyCoord(ra=p.coord.ra, dec=p.coord.dec, distance=np.concatenate([p.samples[k] for k in p.samples.keys()])*u.pc))
In [12]:
#skcoord=SkyCoord(big_skycoord)
In [13]:
#x, y, z=skcoord.galactic.cartesian.xyz
In [14]:
#r=(x**2+y**2)**0.5
In [15]:
import seaborn as sns
In [16]:
import wisps
import numpy as np
from wisps.simulations import custom_volume_correction, Pointing
from wisps import get_distance
from astropy.coordinates import SkyCoord
import astropy.units as u
import matplotlib.pyplot as plt
%matplotlib inline
In [17]:
wisps.MAG_LIMITS['wisps']
Out[17]:
In [18]:
coord=SkyCoord(ra=50*u.deg, dec=50*u.deg)
In [61]:
p=wispsim.BayesianPointing(coord=PNTS[0].coord)
In [62]:
p.survey='wisps'
p.mag_limits=wisps.MAG_LIMITS['wisps']
p.compute_distance_limits()
p.computer_volume()
In [63]:
h=plt.hist(p.random_draw(10000, 200, 10000 ), bins='auto', normed=True)
plt.xlabel('distance (pc)')
Out[63]:
In [64]:
vals, cdf=p.cdf(100, 2000)
In [65]:
#plt.plot(vals, cdf)
plt.plot(vals,cdf)
plt.xlabel('distance (pc)')
Out[65]:
In [66]:
#DIST_DATAFRAME.applymap(np.nanmin)
In [67]:
lms=pd.DataFrame.from_records([x.dist_limits for x in OBSERVED_])
In [68]:
big_skycoord= []
for p in pnts:
big_skycoord.append(SkyCoord(ra=p.coord.ra, dec=p.coord.dec, distance=np.concatenate([p.samples[20] for k in p.samples.keys()])*u.pc))
small_skycoord= []
for p in pnts:
small_skycoord.append(SkyCoord(ra=p.coord.ra, dec=p.coord.dec, distance=np.concatenate([p.samples[37] for k in p.samples.keys()])*u.pc))
In [69]:
cood=SkyCoord(big_skycoord)
cood2=SkyCoord(small_skycoord)
In [70]:
x, y, z=cood.galactic.cartesian.xyz
x2, y2, z2=cood2.galactic.cartesian.xyz
In [71]:
idx=np.arange(len(cood))
In [72]:
ints=np.random.choice(idx, size=10000)
In [73]:
xr=x[ints]
yr=y[ints]
zr=z[ints]
xr2=x2[ints]
yr2=y2[ints]
zr2=z2[ints]
In [74]:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
# plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
#ax.scatter(xr, yr, zr, s=1., alpha=0.1, color='r' )
ax.scatter(xr2, yr2, zr2, s=1., alpha=0.5, color='b')
#ax.
ax.view_init(20, 185)
plt.show()
In [75]:
@numba.vectorize
def proper_motion(vt, d):
"""
mu in arcsec/y
d in pc
vt in
"""
return vt/(4.74*d)
@numba.vectorize("float64(float64)")
def probability_of_detection_smearing(d):
#draw a uniform
vt=np.random.uniform(1, 100)
mu= proper_motion(vt, d)
return np.log10(mu)
In [76]:
r2=(xr2**2+yr2**2)**0.5
d2=(xr2**2+yr2**2+zr2**2)**0.5
#sl=probability_of_detection_smearing(d2.value)
f=plt.scatter(r2, zr2, s=1., alpha=0.5)
plt.ylabel('z')
plt.xlabel('r')
Out[76]:
In [77]:
r=(xr**2+yr**2)**0.5
d=(xr**2+yr**2+zr**2)**0.5
sl=probability_of_detection_smearing(d.value)
f=plt.scatter(r, zr, s=1., alpha=0.5)
plt.ylabel('z')
plt.xlabel('r')
Out[77]:
In [78]:
scipy.f_test.cdf?
In [ ]: