In [1]:
import numpy as np
import wisps
import wisps.simulations as wispsim
import wisps.simulations.effective_numbers as eff
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
import astropy.units as u
from tqdm import tqdm
import pandas as pd
import matplotlib as mpl
%matplotlib inline
In [2]:
wispsim.make_pointings()
In [3]:
from wisps import drop_nan
In [4]:
baraffe_data=eff.simulation_outputs()["baraffe2003"]
saumon_data=eff.simulation_outputs()["saumon2008"]
In [5]:
import seaborn as sns
In [6]:
cmap_teff=sns.diverging_palette(100, 300, s=80, l=55, n=19, as_cmap=True)
In [7]:
pnts=wisps.OBSERVED_POINTINGS
In [8]:
volumes=[]
dlimits=[]
for pnt in pnts:
vs=[]
dls=[]
for g in wispsim.SPGRID:
vsx=[]
for h in wispsim.HS:
vsx.append((pnt.volumes[h])[g])
dls.append(pnt.dist_limits[g])
vs.append(vsx)
volumes.append(vs)
dlimits.append(dls)
volumes=np.array(volumes)
dlimits=np.array(dlimits)
In [9]:
galc=SkyCoord([x.coord.galactic for x in pnts])
In [10]:
vls250=volumes[:,0, :,]
In [11]:
fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(111, projection="mollweide")
c=ax.scatter(galc.l.wrap_at(180*u.degree).radian,galc.b.wrap_at(90*u.degree).radian, marker='+', cmap='viridis')
ax.set_xlabel('l (deg)', fontsize=18)
ax.set_ylabel('b (deg)', fontsize=18)
plt.grid()
plt.savefig(wisps.OUTPUT_FIGURES+'/fields_skymap.pdf', bbox_inches='tight')
In [12]:
from matplotlib.colors import Normalize
import seaborn as sns
cnorm=Normalize(wispsim.HS[0]/100, wispsim.HS[-1]/100)
In [13]:
from scipy import integrate
In [14]:
volume_fx=np.vectorize(wispsim.custom_volume)
In [60]:
import seaborn as sns
cmap= sns.color_palette("coolwarm", 8, as_cmap=True)
In [56]:
ds=np.logspace(0, 3.5, 1000)
fig, ax=plt.subplots(figsize=(8, 5))
for idx, h in tqdm(enumerate(wispsim.HS)):
plt.plot(ds, np.log10(volume_fx(0.,np.pi/4, 0, ds,h)), color=cmap(cnorm(h/100)), label=r'h ={}'.format(h))
plt.plot(ds, np.log10(ds**3), c='k', label=r'd$^3$')
plt.ylabel(r'Log Veff (pc$^3$)', fontsize=18)
plt.xlabel('d (pc)', fontsize=18)
plt.legend(fontsize=14)
In [17]:
vsunif=np.nansum((dlimits[:, :,0])**3-(dlimits[:, :,1])**3, axis=0)*4.1*(u.arcmin**2).to(u.radian**2)
In [18]:
VOLUMES=np.nansum(volumes.T, axis=2)*4.1*(u.arcmin**2).to(u.radian**2)
In [54]:
fig, ax=plt.subplots(figsize=(8, 5))
for idx, h in enumerate(wispsim.HS):
plt.plot(wispsim.SPGRID, np.log10(VOLUMES[idx]), color=cmap(cnorm(h/100)),
label=r'h ={}'.format(h), linewidth=3)
#plt.plot(wispsim.SPGRID, np.log10(vsunif), label=r'd$^3$' )
plt.ylabel(r'Log Veff (pc$^3$)', fontsize=18)
plt.xlabel('SpT', fontsize=18)
plt.legend(fontsize=18)
plt.minorticks_on()
plt.tight_layout()
ax.set_xticks([20, 25, 30, 35, 40])
ax.set_xticklabels(['L0', 'L5', 'T0', 'T5', 'Y0'])
plt.savefig(wisps.OUTPUT_FIGURES+'/simulation_volumes.pdf', bbox_inches='tight')
In [20]:
cum_volumes=np.cumsum(volumes, axis=1)
cum_volumes.shape
Out[20]:
In [21]:
np.nansum((cum_volumes[:, :, 0]), axis=1).shape
Out[21]:
In [22]:
steps=np.arange(533)
In [23]:
wispsim.SPGRID
Out[23]:
In [24]:
import splat
In [25]:
fig, ((ax, ax1), (ax2, ax3))=plt.subplots(figsize=(10, 8), ncols=2, nrows=2)
for idx, h in enumerate(wispsim.HS):
ax.step(steps, np.log10(np.cumsum(volumes[:, 0, idx])), color=cmap(cnorm(h/100)))
ax1.step(steps, np.log10(np.cumsum(volumes[:, 9, idx])), color=cmap(cnorm(h/100)))
ax2.step(steps, np.log10(np.cumsum(volumes[:,-6 , idx])), color=cmap(cnorm(h/100)))
ax3.step(steps, np.log10(np.cumsum(volumes[:, -1, idx])), color=cmap(cnorm(h/100)), label=r'h ={}'.format(h))
ax.set_ylabel(r'Log Cumulative Veff (pc$^3$)', fontsize=18)
ax2.set_ylabel(r'Log Cumulative Veff (pc$^3$)', fontsize=18)
for a in [ax, ax1, ax2, ax3]:
a.set_xlabel('Number of Pointings', fontsize=18)
a.minorticks_on()
ax.set_title('{} UCDs'.format(splat.typeToNum(wispsim.SPGRID[0])), fontsize=18)
ax1.set_title('{} UCDs'.format(splat.typeToNum(wispsim.SPGRID[9])), fontsize=18)
ax2.set_title('{} UCDs'.format(splat.typeToNum(wispsim.SPGRID[-6])), fontsize=18)
ax3.set_title('{} UCDs'.format(splat.typeToNum(wispsim.SPGRID[-1])), fontsize=18)
ax3.legend(fontsize=15)
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/cumulative_volumes.pdf', bbox_inches='tight')
In [26]:
import numba
In [27]:
@numba.jit
def get_kde(r, z):
kde=wisps.kernel_density([r, z])
dens=kde.pdf([r, z])
return dens
In [28]:
data=eff.simulation_outputs()['baraffe2003']
In [29]:
r100=data[150]['r'].flatten()
z100=data[150]['z'].flatten()
d100=data[150]['d'].flatten()
In [30]:
r600=data[400]['r'].flatten()
z600=data[400]['z'].flatten()
d600=data[400]['d'].flatten()
In [31]:
data[200].keys()
Out[31]:
In [32]:
dens1=wispsim.density_function(r100, z100)
dens2=wispsim.density_function(r600, z600)
#dens1=data[150]['sl']
#dens2=data[400]['sl']
In [33]:
fig, ax=plt.subplots()
ax.scatter(data[400]['spts'], np.log10(data[400]['snrj']), c=data[400]['sl'], cmap='viridis')
Out[33]:
In [34]:
fig, (ax, ax1)=plt.subplots( figsize=(10, 10), ncols=2, nrows=2)
h=ax[0].hist(np.log10(d100.flatten()), color='#FF4136',
bins='auto', histtype='step', label="h = 150 pc")
h=ax[0].hist(np.log10(d600.flatten()), color='#0074D9', bins='auto',
histtype='step', label="h = 400 pc")
h=ax[1].hist(z100.flatten(), color='#FF4136',
bins='auto', histtype='step', label="h = 150 pc")
h=ax[1].hist(z600.flatten(), color='#0074D9', bins='auto',
histtype='step', label="h = 400 pc")
c1=ax1[0].scatter(r100, z100, s=1., c=np.log10(dens1), cmap='viridis', alpha=0.1)
ax1[0].set_title("h = 150 pc", fontsize=18)
ax1[1].set_title("h = 400 pc", fontsize=18)
c=ax1[1].scatter(r600, z600,s=1., c=np.log10(dens2), cmap='viridis', alpha=0.1)
cbar0=plt.colorbar(c, ax=ax1[1], orientation='horizontal')
cbar1=plt.colorbar(c1, ax=ax1[0], orientation='horizontal')
cbar0.ax.set_xlabel(r'Log $ \rho/ \rho_0$ ', fontsize=18)
cbar1.ax.set_xlabel(r'Log $ \rho/ \rho_0$ ', fontsize=18)
ax[0].legend(fontsize=15, loc='lower left')
ax[0].set_xlabel('Log Distance (pc)', fontsize=18)
ax[0].set_ylabel('N', fontsize=18)
ax[1].set_xlabel('Z (pc)', fontsize=18)
ax[1].set_ylabel('N', fontsize=18)
ax[1].set_yscale('log')
#ax1[0].legend(fontsize=15, loc='lower left')
ax1[0].set_xlabel('R (pc)', fontsize=18)
ax1[0].set_ylabel('Z (pc)', fontsize=18)
#ax1[1].legend(fontsize=15, loc='lower left')
ax1[1].set_xlabel('R (pc)', fontsize=18)
ax1[1].set_ylabel('Z (pc)', fontsize=18)
for a in np.concatenate([ax, ax1]):
a.minorticks_on()
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/simulation_distances.jpeg', bbox_inches='tight', rasterized=True)
In [35]:
maglimits=pd.DataFrame([x.mag_limits for x in pnts])
exposure_times=[np.nanmean(x.exposure_time) for x in pnts]
In [36]:
maglimits['survey']=[x.survey for x in pnts]
maglimits['exp']=np.log10(exposure_times)
In [37]:
#magpolw=wispsim.MAG_LIMITS[survey][][0]
#magpolh=wispsim.MAG_LIMITS[survey][key][0]
#magsctt=MAG_LIMITS[survey][key][1]
#maglt=np.nanmean(np.random.normal(magpol(np.log10(pnt.exposure_time)), magsctt, 100))
In [38]:
import seaborn as sns
In [39]:
wisps_pnts=[x for x in pnts if x.name.startswith('par')]
In [40]:
hst3d_pnts=[x for x in pnts if not x.name.startswith('par')]
In [41]:
pols={'wisps':{}, 'hst3d':{}}
In [42]:
wisps.MAG_LIMITS
Out[42]:
In [43]:
pnt_dicts={'wisps':wisps_pnts, 'hst3d':hst3d_pnts}
In [44]:
p=wisps_pnts[0]
In [45]:
#hst3d_pnts
In [46]:
fig, ax=plt.subplots(ncols=3, figsize=(10, 6))
for mag_key, a in zip(['F110', 'F140', 'F160'],ax):
for survey in pnt_dicts.keys():
pts=pnt_dicts[survey]
nsbools=np.array([len(x.mags[mag_key]) for x in pts])
f110ls=np.array([x.mag_limits[mag_key] for x in pts])
expts=np.array([np.log10(x.exposure_time) for x in pts ])
nless=nsbools<100
nans=np.isnan(f110ls)
p0=np.nan
p0_scatter=np.nan
if len(f110ls[~np.logical_or(nans, nless)]) !=0:
p0=wisps.fit_polynomial( expts[~(nans|nless) ], f110ls[~(nans | nless)], n=1, sigma_clip=False, sigma=False)
p0_scatter=np.nanstd( f110ls[~np.logical_or(nans, nless)])
a.plot(p0(expts[~nans]), expts[~nans])
pols[survey][mag_key]=(p0, p0_scatter)
c=a.scatter(f110ls, expts, c=nless.astype(int).astype(float))
a.set_xlabel(mag_key+'W', fontsize=18)
a.set_ylabel('Log Exposure Time', fontsize=18)
#a.set_xlim([20, 26])
plt.tight_layout()
In [47]:
pols
Out[47]:
In [48]:
#import pickle
#output = open(wisps.OUTPUT_FILES+'/magnitude_cutoff.pkl', 'wb')
#pickle.dump(pols, output)
#output.close()
In [ ]:
In [49]:
fig, ax=plt.subplots(ncols=3, figsize=(12, 4))
sns.scatterplot(x='F110', y='exp', hue='survey', s=15., ax=ax[0], data=maglimits)
sns.scatterplot(x='F140', y='exp', hue='survey', s=15., ax=ax[1], data=maglimits)
sns.scatterplot(x='F160', y='exp', hue='survey', s=15., ax=ax[2], data=maglimits)
(ax[0]).plot(((wispsim.MAG_LIMITS['wisps'])['F110'][0])(maglimits.exp),maglimits.exp)
(ax[1]).plot(((wispsim.MAG_LIMITS['wisps'])['F140'][0])(maglimits.exp),maglimits.exp)
(ax[2]).plot(((wispsim.MAG_LIMITS['wisps'])['F160'][0])(maglimits.exp),maglimits.exp)
#(ax[1]).plot(((wispsim.MAG_LIMITS['wisps'])['F140'][0])(maglimits.exp),maglimits.exp)
#(ax[2]).plot(((wispsim.MAG_LIMITS['wisps'])['F160'][0])(maglimits.exp),maglimits.exp)
#ax[1].scatter(maglimits.F140, exposure_times, c='#FF851B', s=15.)
#ax[2].scatter(maglimits.F160, exposure_times, c='#FF4136', s=15.)
for a,k in zip(ax, maglimits.columns):
#a.set_xlim([20, 25])
#a.set_yscale('log')
a.minorticks_on()
a.set_xlabel('PS Limiting '+ k+'W', fontsize=18)
a.set_ylabel('Log Grism Exposure Time (s)', fontsize=18)
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/mag_limit.pdf', bbox_inches='tight')
In [50]:
asfbfgn
In [ ]:
maglimits[maglimits.F160<15.]
In [ ]:
p=pnts[np.random.randint(533)]
In [ ]:
p.name
In [ ]:
import scipy
In [ ]:
vals=p.mags['F160']
vals=vals[~np.isnan(vals)]
In [ ]:
len(vals)
In [ ]:
kde=scipy.stats.gaussian_kde(vals)
In [ ]:
frq, edges=np.histogram(vals, bins=int(np.ptp(vals)/0.3), normed=True)
In [ ]:
firstd=np.diff(frq, n=1, prepend=frq[0])
sencd=np.diff(firstd, n=1, prepend=firstd[0])
In [ ]:
len(edges[: -1]), len(np.gradient(frq)), len(frq), len(firstd)
In [ ]:
fig, ax=plt.subplots(figsize=(10, 6))
plt.plot(np.sort(vals), kde(np.sort(vals)), c='b')
plt.bar(edges[:-1], frq, width=np.diff(edges), align="edge", label='0.5 hist', color='#7FDBFF',
edgecolor='#7FDBFF', fill=True)
plt.axvline(edges[:-1][firstd==0][0], c='k')
plt.axvline(edges[:-1][firstd==0][-3], c='k')
plt.axvline(edges[:-1][firstd==0][-2], c='k')
plt.axvline(edges[:-1][firstd==0][-1], c='k')
#plt.axvline(edges[:-1][np.argmin(sencd)], c='#2ECC40', label='sec')
plt.axvline(vals[np.argmax(kde(vals))], c='#2ECC40', linestyle='--', label='kde')
plt.legend(fontsize=18)
plt.xlabel('F160W', fontsize=18)
plt.minorticks_on()
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/mag_limit_illustration.pdf')
In [ ]:
x=np.arange(10).astype(int)
In [ ]:
ys=np.random.random(10)
In [ ]:
bs=(x==3)
In [ ]:
x[bs], ys[bs]
In [ ]:
a=np.empty(23)
a[:]=np.nan
In [ ]:
a
In [ ]: