In [ ]:

    

import splat
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import matplotlib as mpl
from astropy.visualization import ZScaleInterval
import matplotlib
import wisps
%matplotlib inline


    

In [ ]:

    

df_ids=wisps.UCD_SPECTRA.grism_id


    

In [ ]:

    

df=pd.DataFrame()


    

In [ ]:

    

df['grism_id']=df_ids.apply(lambda x: x.replace('g141', 'G141')).values


    

In [ ]:

    

from dask.distributed import Client, progress
client = Client(processes=False, threads_per_worker=4,
                n_workers=1, memory_limit='2GB',  silence_logs='error')
client


    

In [ ]:

    

res=wisps.get_multiple_sources(df_ids, source_type='source')


    

In [ ]:

    

df['spectra']=res


    

In [ ]:

    

#remove extra problematic spectra
problems=['cosmos-02-g141_04153', 'cosmos-05-g141_07469']


    

In [ ]:

    

df=(df[ ~df.grism_id.str.lower().isin(problems)]).reset_index(drop=True)


    

In [ ]:

    

df['spt']=df['spectra'].apply(lambda x: wisps.make_spt_number(x.spectral_type[0]))


    

In [ ]:

    

df['spt_unc']=df['spectra'].apply(lambda x: wisps.make_spt_number(x.spectral_type[1]))


    

In [ ]:

    

df=df.sort_values(by='spt').reset_index(drop=True)


    

In [ ]:

    

s=df.spectra.values[0]


    

In [ ]:

    

df['distance']=df.spectra.apply(lambda x: x.distance['val'].value)


    

In [ ]:

    

df['RA']=df.spectra.apply(lambda x: x.coords.ra.value)
df['DEC']=df.spectra.apply(lambda x: x.coords.dec.value)


    

In [ ]:

    

df['snr1']=df.spectra.apply(lambda x: x.snr['snr1'])


    

In [ ]:

    

fdf=df[(df.spt >=17.) ]


    

In [ ]:

    

len(fdf)


    

In [ ]:

    

#distance distributions
from astropy.coordinates import SkyCoord
import astropy.coordinates as astrocoord
import astropy.units as u
coords=SkyCoord(ra=fdf.RA.values*u.deg, dec=fdf.DEC.values*u.deg, distance=fdf.distance.values*u.pc)


    

In [ ]:

    

galoc=coords.transform_to(astrocoord.Galactocentric(galcen_distance=8.3*u.kpc))


    

In [ ]:

    

#Rsun=8300.
#Zsun=27.
#r=np.sqrt( (coords.distance.value * np.cos( coords.galactic.b.radian ) )**2 + Rsun * (Rsun - 2 * coords.distance.value * np.cos( coords.galactic.b.radian ) * np.cos( coords.galactic.l.radian ) ) )
#z=Zsun+ coords.distance.value * np.sin( coords.galactic.b.radian - np.arctan( Zsun / Rsun) )
x, y, z=galoc.cartesian.xyz
r=(x**2+y**2)**0.5


    

In [ ]:

    

fig, (ax, ax1)=plt.subplots(figsize=(11, 4), ncols=2)
sc=ax.scatter(r, z, c=fdf.spt.values, s=50, cmap='viridis',  marker='*')
br=plt.colorbar(sc)
#ax.ylabel('z (pc)', fontsize=20)
#ax.xlabel('r (pc)', fontsize=20)
#plt.grid()

ax.minorticks_on()

sc=ax1.scatter(coords.cartesian.x.value, coords.cartesian.y.value, c=fdf.spt.values, s=50, cmap='viridis', marker='*')


#plt.grid()
ax.minorticks_on()

ax.set_xlim([6000, 11000])
ax.set_ylim([-5000, 6000])

#ax1.set_xlim([-11000, -5000])
#ax1.set_ylim([-3000, 3000])

ax.set_ylabel('Z (pc)', fontsize=20)
ax.set_xlabel('R (pc)', fontsize=20)
ax1.set_ylabel('X (pc)', fontsize=20)
ax1.set_xlabel('Y (pc)', fontsize=20)

ax.minorticks_on()
ax1.minorticks_on()


br.set_ticks([20, 25, 30, 35, 40])
br.set_ticklabels(['L0', 'L5', 'T0', 'T5', 'Y0'])
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/candidate_distances.pdf', bbox_inches='tight')


    

In [ ]:

    

mags=wisps.Annotator.reformat_table(pd.DataFrame.from_records(fdf.spectra.apply(lambda x: x.mags).values))


    

In [ ]:

    

#MAKE THIS COLOR ABSOLUTE MAG from your mag relations


    

In [ ]:

    

abs_mag_relation_140=wisps.POLYNOMIAL_RELATIONS['sp_F140W']
abs_mag_relation_110=wisps.POLYNOMIAL_RELATIONS['sp_F110W']
abs_mag_relation_160=wisps.POLYNOMIAL_RELATIONS['sp_F160W']


    

In [ ]:

    

dhgjk=fdf.sort_values(by='distance')[['grism_id', 'distance', 'spt']]


    

In [ ]:

    

latc_df=pd.DataFrame()


    

In [ ]:

    

latc_df['Shortname']=fdf.spectra.apply(lambda x: x.shortname.upper())
latc_df['designation']=fdf.spectra.apply(lambda x: x.designation.upper())
latc_df['grism id']=fdf.spectra.apply(lambda x: x.name.upper())
latc_df['snrj']=fdf.spectra.apply(lambda x: int(round(x.snr['snr1'])))
latc_df['spt']=fdf.spectra.apply(lambda x: splat.typeToNum(x.spectral_type[0]) +r'$\pm$'+ str(x.spectral_type[1]))
latc_df['ra']=fdf.spectra.apply(lambda x: x.ra.value)
latc_df['dec']=fdf.spectra.apply(lambda x: x.dec.value)
latc_df['f110']=fdf.spectra.apply(lambda x:  str(round(x.mags['F110W'][0],3)) +r'$\pm$'+ str(round(x.mags['F110W'][1],3)))
#latc_df['f110_er']=fdf.spectra.apply(lambda x: round(x.mags['F110W'][1], 1))
latc_df['f140']=fdf.spectra.apply(lambda x:  str(round(x.mags['F140W'][0],3)) +r'$\pm$'+  str(round(x.mags['F140W'][1],3)))
#latc_df['f140_er']=fdf.spectra.apply(lambda x: round(x.mags['F140W'][1], 1))
latc_df['f160']=fdf.spectra.apply(lambda x:   str(round(x.mags['F160W'][0],3)) + r'$\pm$'+ str(round(x.mags['F160W'][1],3)))
#latc_df['f160_er']=fdf.spectra.apply(lambda x: round(x.mags['F160W'][1], 1))
latc_df['distance']=fdf.spectra.apply(lambda x: str(int(round(x.distance['val'].value, 0)))  + r'$\pm$'+ 
                                                 str((int(round(x.distance['er'].value, 0)))))


    

In [ ]:

    

latc_df=latc_df.replace('nan$\\pm$nan', 'nodata')


    

In [ ]:

    

#get the right sequence of coplumns

col_list=['designation','grism id', 'ra', 'dec', 'f110', 'f140', 'f160', 'snrj', 'spt', 'distance']


    

In [ ]:

    

latc_df[col_list]


    

In [ ]:

    

latc_df[col_list].to_latex(wisps.LIBRARIES+'/candidates.tex',
             header=True, index=False, na_rep=' ')


    

In [ ]:

    

df[df.grism_id=='goodsn-33-G141_09283'].spt.apply(wisps.make_spt_number)>=17


    

In [ ]:

    

'goodsn-33-G141_09283'.lower() in latc_df['grism id'].apply(lambda x: x.lower()).values


    

In [ ]:

    

fghvj


    

In [ ]:

    

#%%capture
fold='/Users/caganze/research/wisps/figures/ltwarfs/'

ids=0
for idx, row in fdf.iterrows():
    try:
        s=row.spectra
        print (s)
        filename=fold+'spectrum'+str(ids)+'.jpeg'
        s.plot(save=True, filename=filename)
        ids=ids+1
    except:
        s=wisps.Source(filename=row.grism_id.replace('g141', 'G141'),is_ucd=False )
        print (s)
        filename=fold+'spectrum'+str(ids)+'.jpeg'
        s.plot(save=True, filename=filename)
        ids=ids+1


    

In [ ]:

    

from astropy.coordinates import SkyCoord


    

In [ ]:

    

#indices_df=pd.DataFrame.from_records(df.spectra.apply(lambda x: x.indices).values)


    

In [ ]:

    

#for k in indices_df.columns: df[k]= indices_df[k]


    

In [ ]:

    

#df.to_hdf(wisps.OUTPUT_FILES+'/true_spectra_cands.hdf', key='with_indices')


    

In [ ]:

    

import astropy.units as u


    

In [ ]:

    

sk=SkyCoord(ra=fdf.RA.values*u.deg, dec=fdf.DEC.values*u.deg, distance=fdf.distance.values*u.pc)


    

In [ ]:

    

norm = matplotlib.colors.Normalize(vmin=20., vmax=37.0)


    

In [ ]:

    

from matplotlib import cm


    

In [ ]:

    

import wisps.simulations as wispsim


    

In [ ]:

    

df['pointing']=df.grism_id.apply(wispsim.get_proper_pointing)


    

In [ ]:

    

pntcoords=SkyCoord([x.coord for x in wisps.OBSERVED_POINTINGS if not x.name in df.pointing.values])


    

In [ ]:

    

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='mollweide')

ax.scatter(pntcoords.galactic.l.wrap_at(180*u.deg).radian, pntcoords.galactic.b.wrap_at(180*u.deg).radian, c='grey', marker='x', alpha=0.5)
c=ax.scatter(sk.galactic.l.wrap_at(180*u.deg).radian, sk.galactic.b.wrap_at(180*u.deg).radian,
           c=np.vstack(fdf.spt.values)[:,0], cmap='viridis')


plt.xlabel("l", fontsize=18)
plt.ylabel("b", fontsize=18)
plt.grid()
plt.minorticks_on()
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/candidate_skymap.pdf', bbox_inches='tight')


    

In [ ]:

    

def plot_vals(values, ax, title):
    values=wisps.drop_nan(np.sort(values))
    kernel = stats.gaussian_kde(values)
    height = kernel.pdf(values)
    mode_value = values[np.argmax(height)]
    h=ax.hist(values, bins='auto', normed=True, histtype='step')
    ax.plot(values, height)
    ax.axvline(mode_value, c='k')
    ax.set_xlabel(title, fontsize=18)


    

In [ ]:

    

#fig, ax=plt.subplots(ncols=3,figsize=(12, 4))
#plot_vals((wisps.OBSERVED_POINTINGS[np.random.randint(533)]).mags['F110'], ax[0], 'F110W')
#plot_vals((wisps.OBSERVED_POINTINGS[np.random.randint(533)]).mags['F140'], ax[1], 'F140W')
#plot_vals((wisps.OBSERVED_POINTINGS[np.random.randint(533)]).mags['F160'], ax[2], 'F160W')
#plt.tight_layout()
#plt.savefig(wisps.OUTPUT_FIGURES+'/mag_limits_illustration.pdf', bbox_inches='tight')


    

In [ ]:

    

df['designation']=df.spectra.apply(lambda x: x.designation)


    

In [ ]:

    

df.shape


    

In [ ]:

    

df.to_pickle(wisps.OUTPUT_FILES+'/true_spectra_cands.pkl')


    

In [ ]:

    

dfn=df[df.spt>=17]


    

In [ ]:

    

import seaborn as sns


    

In [ ]:

    

fig, ax=plt.subplots()

sns.distplot(fdf.spt,  ax=ax, hist_kws={'log': True},  norm_hist=False, kde=False)
#h=ax.hist(dfn.spt, log=True, bins=38-17)
#dfn['spt'].plot(kind='hist')

ax.set_xticks([17, 20, 25, 30, 35, 40])
ax.set_xticklabels(['M7', 'L0', 'L5', 'T0', 'T5', 'Y0'])

ax.set_xlabel('Spectral Type', fontsize=18)
ax.set_ylabel('Number', fontsize=18)

plt.minorticks_on()
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/sptdistribution.pdf')


    

In [ ]:

    

fig, ax=plt.subplots()

sns.distplot(fdf.distance.apply(np.log10),  ax=ax, hist_kws={'log': True}, norm_hist=False, kde=False)
#h=ax.hist(dfn.spt, log=True, bins=38-17)
#dfn['spt'].plot(kind='hist')

#ax.set_xticks([17, 20, 25, 30, 35, 40])
#ax.set_xticklabels(['M7', 'L0', 'L5', 'T0', 'T5', 'Y0'])

ax.set_xlabel('Log Distance (pc)', fontsize=18)
ax.set_ylabel('Number', fontsize=18)

plt.minorticks_on()
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/distdistribution.pdf')


    

In [ ]:

    

s=fdf.spectra.iloc[-1]


    

In [ ]:

    

%%capture 
res=splat.classifyByTemplate(s.splat_spectrum, range=[[1.15, 1.65]], spt=[34,38], force=True )


    

In [ ]:

    

res


    

In [ ]:

    

fig, ax=plt.subplots()
ax.plot((res['spectra'][0]).wave, (res['spectra'][0]).flux*0.78)
ax.plot(s.splat_spectrum.wave, s.splat_spectrum.flux)
ax.set_xlim([1.1, 1.7])


    

In [ ]:

    

evol_models=pd.read_csv(wispsim.EVOL_MODELS_FOLDER+'//'+'phillips2020'.lower()+'.csv')


    

In [ ]:

    




    

In [ ]:

    

evol_models.columns


    

In [ ]:

    

import splat


    

In [ ]:

    

pols=wisps.POLYNOMIAL_RELATIONS


    

In [ ]:

    

color_pol=pols['colors']


    

In [ ]:

    

fdf['f110']=fdf.spectra.apply(lambda x: x.mags['F110W'][0])
fdf['f110_er']=fdf.spectra.apply(lambda x: x.mags['F110W'][1])

fdf['f140']=fdf.spectra.apply(lambda x: x.mags['F140W'][0])
fdf['f140_er']=fdf.spectra.apply(lambda x: x.mags['F140W'][1])

fdf['f160']=fdf.spectra.apply(lambda x: x.mags['F160W'][0])
fdf['f160_er']=fdf.spectra.apply(lambda x: x.mags['F160W'][1])


    

In [ ]:

    

spdf=splat.searchLibrary()


    

In [ ]:

    

spdf['SHORTNAME']=spdf.DESIGNATION.apply(lambda x: splat.designationToShortName)


    

In [ ]:

    

spdf=spdf[~ (spdf.PARALLAX.isna() | spdf.SHORTNAME.isin(wisps.FORBIDDEN_LIST))]


    

In [ ]:

    

spdf['absF110W']=pols['sp_F110W'](spdf['SPEX_TYPE'].apply(wisps.make_spt_number))
spdf['absF140W']=pols['sp_F140W'](spdf['SPEX_TYPE'].apply(wisps.make_spt_number))
spdf['absF160W']=pols['sp_F160W'](spdf['SPEX_TYPE'].apply(wisps.make_spt_number))


    

In [ ]:

    

spdf['relF110W']=spdf['J_2MASS']-color_pol['j_f110'][0](spdf['SPEX_TYPE'].apply(wisps.make_spt_number))
spdf['relF140W']=spdf['J_2MASS']-color_pol['j_f140'][0](spdf['SPEX_TYPE'].apply(wisps.make_spt_number))
spdf['relF160W']=spdf['H_2MASS']-color_pol['h_f160'][0](spdf['SPEX_TYPE'].apply(wisps.make_spt_number))


    

In [ ]:

    

absf110= pols['sp_F110W'](np.vstack(fdf.spt.values)[:,0])
absf140= pols['sp_F140W'](np.vstack(fdf.spt.values)[:,0])
absf160= pols['sp_F160W'](np.vstack(fdf.spt.values)[:,0])


    

In [ ]:

    

spts0ssss=np.vstack(fdf.spt.values)[:,0]


    

In [ ]:

    

plt.plot(spts0ssss, color_pol['j_f110'][0](spts0ssss))
plt.plot(spts0ssss, color_pol['j_f140'][0](spts0ssss))
plt.plot(spts0ssss, color_pol['h_f160'][0](spts0ssss))
#plt.plot(np.vstack(latc_df.spt.values)[:,0], absf110)
#plt.plot(np.vstack(latc_df.spt.values)[:,0], absf160)
#plt.gca().invert_xaxis()


    

In [ ]:

    

fig, ax=plt.subplots(ncols=3, figsize=(12, 4))

ax[0].errorbar(spdf['relF160W']-spdf['relF110W'], spdf['absF160W'],  yerr=spdf['H_2MASS_E'],
               xerr= (color_pol['j_f160'][1]**2+spdf['J_2MASS_E']**2)**0.5, fmt='+', ms=5)

ax[1].errorbar(spdf['relF160W']-spdf['relF140W'], spdf['absF160W'], yerr=spdf['H_2MASS_E'],
                xerr= (color_pol['j_f160'][1]**2+spdf['J_2MASS_E']**2)**0.5, fmt='+',  ms=5)

ax[2].errorbar(spdf['relF140W']-spdf['relF110W'], spdf['absF140W'],yerr=spdf['J_2MASS_E'],
                xerr= (color_pol['j_f140'][1]**2+spdf['J_2MASS_E']**2)**0.5, fmt='+',  ms=5)



    

ax[0].errorbar(fdf['f160']-fdf['f110'], absf160, yerr=(pols['sigma_sp_F160W']**2+fdf['f160_er']**2)**0.5, 
               xerr=(fdf['f160_er']**2+fdf['f110_er']**2)**0.5, fmt='o')

ax[1].errorbar(fdf['f160']-fdf['f140'], absf160, yerr=(pols['sigma_sp_F160W']**2+fdf['f160_er']**2)**0.5, 
               xerr=(fdf['f160_er']**2+fdf['f140_er']**2)**0.5, fmt='o')

ax[2].errorbar(fdf['f140']-fdf['f110'], absf140, yerr=(pols['sigma_sp_F140W']**2+fdf['f140_er']**2)**0.5,  
               xerr=(fdf['f140_er']**2+fdf['f110_er']**2)**0.5,fmt='o')


ax[0].set_xlabel('F160W-F110W')
ax[0].set_xlabel('abs F160W')

ax[0].axhline(pols['sp_F110W'](20), c='k')
ax[1].axhline(pols['sp_F140W'](20), c='k')
ax[2].axhline(pols['sp_F160W'](20), c='k')

ax[0].set_xlabel('F160W-F110W', fontsize=18)
ax[0].set_ylabel('abs F160W', fontsize=18)

ax[1].set_xlabel('F160W-F140W', fontsize=18)
ax[1].set_ylabel('abs F160W', fontsize=18)

ax[2].set_xlabel('F140W-F110W', fontsize=18)
ax[2].set_ylabel('abs F140W', fontsize=18)

for a in ax:
    a.invert_yaxis()
    #a.set_xlim([-2., 1.0])

plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/color_color_plots.pdf')


    

In [ ]:

    

pols['sp_F110W'](20)


    

In [ ]: