In [1]:
import numpy as np
import splat
import wisps
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
import wisps.simulations as wispsim
In [3]:
from astropy.coordinates import SkyCoord
In [4]:
J1_NIPS=24 #limiting mag in NISP
J2_NIPS=27 #best case scenario
In [5]:
def get_abs_mag(spt):
#magnitude from spt and relation
j, h=wisps.absolute_magnitude_jh(spt)
return j
In [6]:
spgrid=wispsim.SPGRID
In [7]:
abs_mags=np.array([get_abs_mag(s ) for s in spgrid])
In [8]:
dmaxs1=wisps.get_distance(abs_mags, np.ones_like(abs_mags)*J1_NIPS)
dmaxs2=wisps.get_distance(abs_mags, np.ones_like(abs_mags)*J2_NIPS)
In [9]:
import astropy.units as u
In [10]:
EDFS=SkyCoord(l=24.6*u.deg, b=-82.0*u.deg , frame='galactic').galactic
EDF_CH=SkyCoord("3:32:28.0 -27:48:30" , obstime="J2000", unit=u.deg).galactic
In [11]:
#compute volume correction terms
EDFS_vcs1=np.array([ wispsim.custom_volume(EDFS.l.radian,EDFS.b.radian, 1., dm, 300.) for dm in dmaxs1])
EDF_CH_vcs1=np.array([ wispsim.custom_volume(EDF_CH.l.radian, EDFS.b.radian,1, dm, 300.) for dm in dmaxs1])
EDFS_vcs2=np.array([ wispsim.custom_volume(EDFS.l.radian,EDFS.b.radian, 1., dm, 300.) for dm in dmaxs2])
EDF_CH_vcs2=np.array([ wispsim.custom_volume(EDF_CH.l.radian, EDFS.b.radian,1, dm, 300.) for dm in dmaxs2])
In [12]:
EDFS_CH_vols1=EDF_CH_vcs1*(20*(u.deg**2)).to(u.radian**2)
EDF_vols1=EDFS_vcs1*(10*(u.deg**2)).to(u.radian**2)
EDFS_CH_vols2=EDF_CH_vcs2*(20*(u.deg**2)).to(u.radian**2)
EDF_vols2=EDFS_vcs2*(10*(u.deg**2)).to(u.radian**2)
In [13]:
np.round(spgrid).astype(int)
Out[13]:
In [14]:
#the empirical luminosity function
import wisps.simulations.euclid as eucl
In [ ]:
def get_numbers( model, h=300, field='fornax'):
ds=eucl.expected_numbers(model, field=field, h=h)
simdf=pd.DataFrame.from_records(ds)
cutdf=(simdf[simdf.j.between(0, 27)]).reset_index(drop=True)
NORM = 0.63*(10**-3)/ len(cutdf.teff[np.logical_and(cutdf.teff>=1650, cutdf.teff <=1800)])
NSIM=dict(zip(wispsim.SPGRID,np.zeros(len(wispsim.SPGRID))))
#rounded spectral type
cutdf['spt_r']=cutdf.spt.apply(np.round)
for g in cutdf.groupby('spt_r'):
NSIM[g[0]]=np.nansum((g[1]).prob*NORM)
return NSIM
In [ ]:
%%capture
n_south_brf=get_numbers('baraffe2003', h=300, field='south')
n_south_marl=get_numbers('marley2019', h=300, field='south')
n_south_phil=get_numbers('phillips2020', h=300, field='south')
n_south_saum=get_numbers('saumon2008', h=300, field='south')
In [ ]:
%%capture
n_forn_brf=get_numbers('baraffe2003', h=300, field='fornax')
n_forn_marl=get_numbers('marley2019', h=300, field='fornax')
n_forn_phil=get_numbers('phillips2020', h=300, field='fornax')
n_forn_saum=get_numbers('saumon2008', h=300, field='fornax')
In [ ]:
import seaborn as sns
sns.set_palette('Set2')
In [ ]:
fig, (ax, ax1)=plt.subplots(ncols=2, figsize=(10, 4), sharey=False, sharex=True)
ax.step(spgrid, np.array([n_south_brf[k] for k in wispsim.SPGRID])*EDFS_CH_vols2)
ax.step(spgrid, np.array([n_south_marl[k] for k in wispsim.SPGRID])*EDFS_CH_vols2)
ax.step(spgrid, np.array([n_south_phil[k] for k in wispsim.SPGRID])*EDFS_CH_vols2)
ax.step(spgrid, np.array([n_south_saum[k] for k in wispsim.SPGRID])*EDFS_CH_vols2)
ax1.step(spgrid, np.array([n_forn_brf[k] for k in wispsim.SPGRID])*EDFS_CH_vols2)
ax1.step(spgrid, np.array([n_forn_marl[k] for k in wispsim.SPGRID])*EDFS_CH_vols2)
ax1.step(spgrid, np.array([n_forn_phil[k] for k in wispsim.SPGRID])*EDFS_CH_vols2)
ax1.step(spgrid, np.array([n_forn_saum[k] for k in wispsim.SPGRID])*EDFS_CH_vols2)
ax.set_title(r'Euclid South', fontsize=18)
ax1.set_title(r'Euclid Fornax', fontsize=18)
for a in [ax, ax1]:
a.minorticks_on()
a.set_xlabel('SpT', fontsize=18)
a.legend()
#a.set_yscale('log')
a.set_ylabel(r'N', fontsize=18)
a.set_xticks([15, 20, 25, 30, 35, 40])
a.set_xticklabels(['M5', 'L0', 'L5', 'T0', 'T5', 'Y0'])
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/euclid_predictions.pdf')
In [ ]: