In [1]:

    

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from astropy.table import Table
from astropy import units as u
from astropy.modeling.blackbody import blackbody_nu


    

In [2]:

    

def Mdisk(F,d,r,k,Td,lam):
    #Eq 2 Andrews et al. 2013 
    #F is *mm* Flux in *Jy*
    #d is the distance in *pc*
    #r is the dust to mass ratio *adimensional*
    #k is the dust opacity in *cm^2/g*
    #lam is the wavelenght in *mm*
    c=2.99792458e10      # Light velocity *cm/s*
    h=6.6260755e-27 # Planck constant in erg*s (cgs)
    kb=1.3806488e-16    # Boltzmann constant in erg/K (cgs)
    
    ni=c/(lam*0.1) #in Hz
    B_ni=1e23*2*h*ni**3/c**2*(1./(np.exp(h*ni/kb/Td)-1)) #  Planck function in *Jy* --> 1Jy=1e-23erg/s/cm^2/Hz
#    B_ni=1e23*blackbody_nu(ni, Td)      #  Planck function in *Jy* --> 1Jy=1e-23erg/s/cm^2/Hz
    
#    return np.log10(F)+2*np.log10(d*3.08568e18)-np.log10(r*k*1.988e33)-np.log10(B_ni)  #log(Disk mass/*Msol*) 
    return F*(d*3.086e18)**2/k/1.988e33/r/B_ni       #Disk mass in *Msol*
    
def eMd(errF,F,Md):
    return Md*errF/F


    

In [4]:

    

ChaI=Table.read('ChaI_Fmm_Ms.fit')
ChaI_Md=Table.read('ChaI_Mdust.fit')
#ChaI_Md.show_in_browser()

name=ChaI['_2MASS']
logMs=ChaI['logM_']
up_logMs=ChaI['b_logm__lc']          #upper boundary in logM*
down_logMs=ChaI['b_logM_']        #Lower boundary in logM*
Fmm_tab=ChaI['Fnu']      #Integrated flux density at 887um
e_Fmm=ChaI['e_Fnu']
ecal_Fmm=np.sqrt((e_Fmm)**2+(0.1*Fmm_tab)**2)       #Error on Fmm including the calibration error(10%) in quadratic sum

print name[1],Fmm_tab[1],logMs[1]
print 10**ChaI_Md['logd20']


    

    




J10555973-7724399 34.1 -0.13
 logd20 
 [Mgeo] 
--------
 10.9648
 10.8143
 7.26106
 9.31108
 28.9734
  7.4817
 138.357
 76.3836
 6.65273
 5.15229
     ...
0.334195
0.097051
0.097051
 3.16228
 6.35331
0.776247
0.097051
 1.24451
0.772681
0.161808
Length = 446 rows

In [101]:

    

notdelta= Fmm_tab < 3*e_Fmm        #upperlimits
delta=np.logical_not(notdelta)      #observed sources

upp=np.empty(len(delta),dtype=str)
upp[notdelta]='<'

Fmm=Fmm_tab.copy()
Fmm[notdelta]=3*e_Fmm[notdelta]      # upper limit is about 3 sigma


    

In [121]:

    

d=160        #pc distance of Chamaleon in pc
r=1          #dust to gass mass ratio
k=2.3*(338./230)**(0.4)        #cm^2/g at 887um
lam=0.887    #wavelenght in mm

Td=20        #Temperature of the disk in K

Md=Mdisk(Fmm*1e-3,d,r,k,Td,lam)/3e-6         #Fmm has to be in *Jy*; Md in Earth masses!
logMd=np.log10(Md)

e_Md=eMd(ecal_Fmm,Fmm,Md)  #error on Md

print name[0], logMd[0]


    

    




J10533978-7712338 0.1873052873

In [111]:

    

T=Table([name,logMs,up_logMs,down_logMs,upp, Md,e_Md, Fmm, e_Fmm,ecal_Fmm], 
        names=('Name','logM*','up_logM*','down_logM*','l_Md','Md','e_Md', 'Fmm','e_Fmm','e+cal_Fmm'))    

T.write('ChaI_Tab.fit', format='fits', overwrite='True')


    

In [ ]: