flux_mass_conversion.py

import numpy as np
import scipy.constants as sc

M_Sun = 1.98844e30           # [M_Sun] = kg
M_Earth = 5.9723e24          # [M_Earth] = kg
M_Jup = 1.8986e27            # [M_Jup] = kg
M_Sat = 5.6846e26            # [M_Sat] = kg


def Bnu(nu, T, RJ=False):
    if RJ:
        # Planck function Rayleigh-Jeans [W m-2 Hz-1 sr-1].
        Bnu_RJ = 2 * np.power(nu, 2.) * sc.k * T / sc.c**2
        return Bnu_RJ

    # Full Planck function [W m-2 Hz-1 sr-1].
    Bnu = 2 * sc.h * np.power(nu, 3) / np.power(sc.c, 2)
    Bnu /= (np.exp(sc.h * nu / (sc.k * T)) - 1.)
    return Bnu


def flux_mass(nu, Fnu, T, distance=100., RJ=False, verbose=False):
    """Mass of solids assuming optically thin continuum emission Fnu.
    Inputs: nu [Hz], Fnu [Jy], T [K], distance [pc].  
    Output: mass [kg]. """

    Fnu *= 1e-26                # [Jy]->[W m-2 Hz-1]

    # Dust opacity (per unit of solids mass)
    # kappa = 0.02                # [cm2/gr] (Beckwith)
    kappa = 0.1*(nu/1e12)       # [cm2 gr-1] (total mass, gas + dust)

    kappa *= (1000/100**2.)     # --> [m2 kg-1]
    kappa *= 100.               # gas-to-dust mass ratio

    M_dust_RJ = Fnu * np.power(distance * sc.parsec, 2)
    M_dust_RJ /= (kappa * Bnu(nu, T, RJ=True))

    M_dust = Fnu * np.power(distance * sc.parsec, 2.)
    M_dust /= (kappa * Bnu(nu, T))

    # Difference between full Planck and Rayleigh-Jeans.
    # print("Bnu_RJ/Bnu = ", (Bnu(nu,T,RJ=True)/Bnu(nu,T)), "at nu =", (nu), "Hz")

    if verbose:
        print("M_dust (Planck)\t=", (M_dust/M_Earth), "M_Earth")
        print("\t\t=", (M_dust/M_Jup), "M_Jup")
        print("\t\t=", (M_dust/M_Sun), "M_Sun")
        print("M_dust (RJ)\t=", (M_dust_RJ/M_Earth), "M_Earth")
        print("\t\t=", (M_dust_RJ/M_Jup), "M_Jup")
        print("\t\t=", (M_dust_RJ/M_Sun), "M_Sun")
    if RJ:
        return M_dust_RJ
    else:
        return M_dust


# Continuum FU/EX ors from Antonio's sample
Tdust = 20                      # [k]
nu = 225.5e9                    # [Hz]

print("Mass of solids in")

print(" V582 Aur =", '{:.3f}'.format(
    flux_mass(nu, 5.3e-3, Tdust, distance=2575)/M_Earth), "M_Earth")
print(" V900 Mon =", '{:.3f}'.format(
    flux_mass(nu, 9.8e-3, Tdust, distance=1500)/M_Earth), "M_Earth")
print(" UZ Tau E =", '{:.3f}'.format(
    flux_mass(nu, 134.e-3, Tdust, distance=131)/M_Earth), "M_Earth")
print(" GM Cha   =", '{:.3f}'.format(
    flux_mass(nu, 10.4e-3, Tdust, distance=160)/M_Earth), "M_Earth")

print("Using Rayleigh-Jeans approximation:")

print(" V582 Aur =", '{:.3f}'.format(
    flux_mass(nu, 5.3e-3, Tdust, distance=2575, RJ=True)/M_Earth), "M_Earth")
print(" V900 Mon =", '{:.3f}'.format(
    flux_mass(nu, 9.8e-3, Tdust, distance=1500, RJ=True)/M_Earth), "M_Earth")
print(" UZ Tau E =", '{:.3f}'.format(
    flux_mass(nu, 134.e-3, Tdust, distance=131, RJ=True)/M_Earth), "M_Earth")
print(" GM Cha   =", '{:.3f}'.format(
    flux_mass(nu, 10.4e-3, Tdust, distance=160, RJ=True)/M_Earth), "M_Earth")

Mass of solids in
V582 Aur = 1054.959 M_Earth
V900 Mon = 661.933 M_Earth
UZ Tau E = 69.032 M_Earth
GM Cha = 7.992 M_Earth
Using Rayleigh-Jeans approximation:
V582 Aur = 795.149 M_Earth
V900 Mon = 498.915 M_Earth
UZ Tau E = 52.031 M_Earth
GM Cha = 6.024 M_Earth