notebook.community

Edit and run



In [7]:

    
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from astropy.table import Table
from matplotlib import rc
from matplotlib.ticker import MultipleLocator
from statsmodels.stats.weightstats import DescrStatsW

rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})
rc('text', usetex=True)

#%matplotlib inline
%matplotlib qt5



In [14]:

    
Lup_Ms=Table.read('Lup_Ms.fit')
Lup_Fmm=Table.read('Lup_Fmm.fit')
Lup_Fmm.show_in_browser()

Lup_Ms.remove_rows([13,14,35,53])      #removing non observed sources

name=Lup_Fmm['Name']
Ms=Lup_Ms['Mass']
err_Ms=Lup_Ms['e_Mass']
Md=Lup_Fmm['MDust']
err_Md=Lup_Fmm['e_MDust']
Fmm=Lup_Fmm['F890']
rms=Lup_Fmm['rms']
#Lup_Fmm.info



In [9]:

    
nnan=np.isnan(Ms)        #20 obscured sources
notnan=np.logical_not(nnan)

Reproduce Distribution of Stellar Masses in Fig 9 of Mortier+11 and assign randomly Mstar values to the nan values of Ms based on that distribution.



In [10]:

    
x=np.array([-1.2,-0.9,-0.6,-0.3,0,0.3,0.6])
y=np.array([0,5,12,9,1,1,0])
p=y/28.     #41 is the  number of x elements used in Mortier+11?
fig, ax = plt.subplots(1, 2, figsize=(5, 5))
ax[0].step(x,y)
ax[0].set_xlim(-1.9,0.5)
ax[0].set_ylim(0,15)
ax[1].step(x,p)
ax[1].set_xlim(-1.9,0.5)
ax[1].set_ylim(0,1)

plt.show()

logM=np.random.choice(x, 20, p=p)
M=10**logM

weighted_stats = DescrStatsW(x, weights=p, ddof=0)
std=weighted_stats.std
e_M=M*std
print M
print e_M
print std









    



[ 0.25118864  1.99526231  0.25118864  0.50118723  0.25118864  0.12589254
  0.25118864  0.25118864  0.25118864  0.50118723  0.50118723  0.25118864
  0.25118864  0.50118723  0.50118723  0.50118723  0.50118723  0.25118864
  0.25118864  0.12589254]
[ 0.0699222   0.55541181  0.0699222   0.13951314  0.0699222   0.03504412
  0.0699222   0.0699222   0.0699222   0.13951314  0.13951314  0.0699222
  0.0699222   0.13951314  0.13951314  0.13951314  0.13951314  0.0699222
  0.0699222   0.03504412]
0.278365308359



In [12]:

    
Ms[nnan]=10**logM
err_Ms[nnan]=e_M

notdelta= Fmm < 3*rms          #upperlimits
delta=np.logical_not(notdelta)        #observed sources
print len(Fmm),len(Fmm[delta]), len(Ms[nnan])
upp=np.empty(89,dtype=str)
upp[notdelta]='<'



In [15]:

    
T=Table([name,Ms,err_Ms,upp,Md,err_Md], 
        names=('Name','M*','err_M*','l_Md','Md','err_Md'))       #with the 20 obscured sources
T2=Table([name[notnan],Ms[notnan],err_Ms[notnan],upp[notnan],Md[notnan],err_Md[notnan]], 
        names=('Name','M*','err_M*','l_Md','Md','err_Md'))      #without the 20 obscured sources
#T.show_in_browser()
#T2.show_in_browser()

T2.write('Lupus_Tab.fit', format='fits', overwrite='True')



In [201]:



In [ ]: