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 [3]:

    

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 [7]:

    

UppSco=Table.read('UpperScorpius_Tab.fit')
#UppSco.show_in_browser()

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


    

In [63]:

    

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

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


    

In [82]:

    

d=145.        # distance of Chamaleon in pc
e_d=20.       # error on distance d in pc
r=1.          #dust to gass mass ratio
k=2.3*(338./225.4)**(0.4)        #cm^2/g at 887um
lam=0.88     #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!
e_Md=eMd(ecal_Fmm,Fmm,Md)  #error on Md which include the error on flux
add_Md=e_Md+2.*e_d*Md/d    #error on Md which include the error on flux and on the distance.

#Md[notdelta]=3.*e_Md[notdelta]

logMd=np.log10(Md)
print len(Md)
prim=UppSco['Type']!='Debris/Ev. Trans.'     # only primordial disks not debris/evolved


    

    




106

In [83]:

    

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

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


    

In [ ]: