In [1]:
%matplotlib inline
import numpy as np
from astropy import wcs
from astropy import units as u
from astropy.io import fits
from astropy.utils.data import download_file
from astropy.coordinates import SkyCoord
from astrodendro import Dendrogram, pp_catalog
from astrodendro.analysis import PPStatistic
import aplpy
import matplotlib.pyplot as plt
from matplotlib import rc
#rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
#rc('text', usetex=False)
rc('font',**{'family':'Helvetica','sans-serif':['Helvetica']})
rc('text', usetex=True)
Load combined GBT and VLA data. This data set is created using GBT data as image model for the VLA imaging.
Here I also setup some of the variables needed throughout the notebook. The condensations id are identified after exploring the dendrogram analysis off-line (interactively).
In [2]:
data_file='data/B5_VLA_GBT_model_11_mom0.fits'
Tk_file='data/B5_VLA_GBT_model_Tk_QA.fits'
dv_file='data/B5_VLA_GBT_model_dv_QA.fits'
vc_file='data/B5_VLA_GBT_model_vc_QA.fits'
file_scuba2_450='data/B5_450um_ext_v2_regrid_beam.fits'
file_scuba2_850='data/B5_850um_ext_v2_regrid_beam.fits'
# Load NH3 integrated intensity map
hdu =fits.open(data_file)
header=hdu[0].header
data =hdu[0].data
hdu.close()
# rms in integrated intensity map
nh3_rms=2.8e-3 # Jy/beam km/s
scuba2_450_rms=0.23e-3 # Jy/pixel
scuba2_850_rms=2.6e-5 # Jy/pixel
# if keyword si true then the figures will also be saved as PDF and PNG files
save_figure=True
figure_dir='figures/'
#For VLA-only data: 13, 12, and 7
#For VLA+GBT feathered data: 35, 29, and 23
id_Cond1=72
id_Cond2=68
id_Cond3=58
# Extra parameters
distance=250. # pc
beamsize=6.0 # arcsec
pixel_size=0.5 # arcsec
Dendrogram for combined VLA + GBT data. The dendrogram viewer does not work in the iPython Notebook, it works in the regular iPython mode. The commands commented out can be used in a regular iPython session to explore the dendrogram output.
In [3]:
dendro = Dendrogram.compute(data, min_value=4*nh3_rms, min_delta=2*nh3_rms, min_npix=250)#, verbose=True)
#dendro.save_to(data_file.replace('mom0.fits','dendro.fits'))
#v_dendro = dendro.viewer()
#v_dendro.show()
To calculate the catalogue properties I need to define the metadata structure
In [4]:
metadata = {}
metadata['data_unit'] = u.Jy / u.beam ## * u.kilometer / u.second
metadata['spatial_scale'] = pixel_size * u.arcsec
metadata['beam_major'] = beamsize * u.arcsec
metadata['beam_minor'] = beamsize * u.arcsec
dendro_catalogue = pp_catalog(dendro, metadata)
The conversion factor between Jy/beam into Jy/pixel is calculated as $$\textrm{conv_beam2pix} = \frac{8\,\log(2)}{2\pi B_{maj}\times B_{min}} (\textrm{pixel_size})^2$$
In [5]:
conv_beam2pix= (1./(2*np.pi/(8*np.log(2))*metadata['beam_major']*metadata['beam_minor']) * (metadata['spatial_scale'])**2).value
We sefine the mask for each condensation, and we also calculate the cumulative flux function (to calculate mass), and the condensation radius for each flux contour.
The function created calculates the flux over the background.
In [6]:
def B5_Dendrogram_Mask_Flux_Radius( data, dendro_struct, leaf_id):
""" Function that takes the data, dendrogram structure and leaf ID to
calculate the mask, flux, radius and (x,y) peak position for a single leaf id.
Flux is in units of data.units*pixel, and it contains a correction for background emission.
Radius is in pixels units.
i_peak in pixel units.
j_peak in pixel units.
"""
leaf_mask = np.zeros(data.shape, dtype=bool)
leaf = dendro_struct.leaves[leaf_id]
leaf_mask = leaf_mask | leaf.get_mask()
raw_flux = data[np.where( leaf_mask == 1 )]
leaf_flux = np.cumsum(np.sort(raw_flux)[::-1]-np.min(raw_flux))
leaf_radius = np.sqrt(np.arange(1,len(leaf_flux)+1)/np.pi)
j_peak, i_peak = np.unravel_index( np.nanargmax(data*leaf_mask), data.shape)
return leaf_mask, leaf_flux, leaf_radius, i_peak, j_peak
mask_C1, w_C1, r_C1, i_C1, j_C1 = B5_Dendrogram_Mask_Flux_Radius( data*conv_beam2pix, dendro, id_Cond1)
mask_C2, w_C2, r_C2, i_C2, j_C2 = B5_Dendrogram_Mask_Flux_Radius( data*conv_beam2pix, dendro, id_Cond2)
mask_C3, w_C3, r_C3, i_C3, j_C3 = B5_Dendrogram_Mask_Flux_Radius( data*conv_beam2pix, dendro, id_Cond3)
Now we calculate the uncertainty associated with the NH$_3$(1,1) integrated intensity flux.
$$\sigma_W = rms_{NH_3}\times \textrm{conv_beam2pix} \sqrt{N_{pix}} = rms_{NH_3}\times \textrm{conv_beam2pix} \sqrt{\pi R^2}$$
In [7]:
ew_C1=conv_beam2pix*nh3_rms*np.pi * r_C1
ew_C2=conv_beam2pix*nh3_rms*np.pi * r_C2
ew_C3=conv_beam2pix*nh3_rms*np.pi * r_C3
print np.max(w_C1), np.max(w_C2), np.max(w_C3)
print np.max(ew_C1), np.max(ew_C2), np.max(ew_C3)
The condensations' peak position are stored in the variables i_CN and j_CN. We will use astropy to determine the RA and Dec positions for the condensations.
The position of B5-IRS1 is used to calculate the separations.
In [8]:
ra_yso =np.array([ 15.*(3+(47.+41.548/60)/60.)])
dec_yso=np.array([ 32+(51+43.57/60.)/60.])
wcs = wcs.WCS(data_file)
xxyy = np.array( np.transpose([[i_C1,i_C2,i_C3],[j_C1,j_C2,j_C3]] ) )
xxyy_yso= wcs.wcs_world2pix( np.transpose( np.array( [ ra_yso, dec_yso] )), 0)
radec_C = wcs.wcs_pix2world( xxyy, 0)
ra_blob = radec_C[:,0]
dec_blob = radec_C[:,1]
In [9]:
coord_YSO= SkyCoord(ra=ra_yso[0]*u.degree, dec=dec_yso[0]*u.degree, frame='icrs')
coord_C1 = SkyCoord(ra=ra_blob[0]*u.degree, dec=dec_blob[0]*u.degree, frame='icrs')
coord_C2 = SkyCoord(ra=ra_blob[1]*u.degree, dec=dec_blob[1]*u.degree, frame='icrs')
coord_C3 = SkyCoord(ra=ra_blob[2]*u.degree, dec=dec_blob[2]*u.degree, frame='icrs')
print coord_C1.ra.to_string(unit=u.hour,precision=3)+' '+coord_C1.dec.to_string(unit=u.degree, alwayssign=True, precision=2)
print coord_C2.ra.to_string(unit=u.hour,precision=3)+' '+coord_C2.dec.to_string(unit=u.degree, alwayssign=True, precision=2)
print coord_C3.ra.to_string(unit=u.hour,precision=3)+' '+coord_C3.dec.to_string(unit=u.degree, alwayssign=True, precision=2)
In [10]:
sep_C1_YSO=np.around((coord_YSO.separation(coord_C1)).arcsec * distance, decimals=1)
sep_C2_YSO=np.around((coord_YSO.separation(coord_C2)).arcsec * distance, decimals=1)
sep_C3_YSO=np.around((coord_YSO.separation(coord_C3)).arcsec * distance, decimals=1)
sep_C1_C2 =np.around((coord_C1.separation(coord_C2)).arcsec * distance, decimals=1)
sep_C1_C3 =np.around((coord_C1.separation(coord_C3)).arcsec * distance, decimals=1)
sep_C2_C3 =np.around((coord_C2.separation(coord_C3)).arcsec * distance, decimals=1)
print 'B5-IRS1 and B5-Cond1:', sep_C1_YSO, ' au'
print 'B5-IRS1 and B5-Cond2:', sep_C2_YSO, ' au'
print 'B5-IRS1 and B5-Cond3:', sep_C3_YSO, ' au'
print 'B5-Cond1 and B5-Cond2:', sep_C1_C2, ' au'
print 'B5-Cond1 and B5-Cond3:', sep_C1_C3, ' au'
print 'B5-Cond2 and B5-Cond3:', sep_C2_C3, ' au'
Now, we calculate the position of the condensations with respect to B5-IRS1. This will be used later to calculate how bound is the multiple system.
In [11]:
print 'Source Delta x Delta y'
print ' (au) (au)'
print 'B5-Cond1:',np.around(-(i_C1 - xxyy_yso[0,0])*pixel_size*distance, decimals=1), np.around((j_C1 - xxyy_yso[0,1])*pixel_size*distance, decimals=1)
print 'B5-Cond2:',np.around(-(i_C2 - xxyy_yso[0,0])*pixel_size*distance, decimals=1), np.around((j_C2 - xxyy_yso[0,1])*pixel_size*distance, decimals=1)
print 'B5-Cond3:',np.around(-(i_C3 - xxyy_yso[0,0])*pixel_size*distance, decimals=1), np.around((j_C3 - xxyy_yso[0,1])*pixel_size*distance, decimals=1)
We need to calculate the Tk for each condensation
In [12]:
Tk_11 = fits.getdata(Tk_file)
Tk_C1 = np.around( np.nanmean(Tk_11[np.where( mask_C1 == 1 )] ), decimals=1)
Tk_C2 = np.around( np.nanmean(Tk_11[np.where( mask_C2 == 1 )] ), decimals=1)
Tk_C3 = np.around( np.nanmean(Tk_11[np.where( mask_C3 == 1 )] ), decimals=1)
eTk_C1 = np.around( np.nanstd(Tk_11[np.where( mask_C1 == 1 )] ), decimals=1)
eTk_C2 = np.around( np.nanstd(Tk_11[np.where( mask_C2 == 1 )] ), decimals=1)
eTk_C3 = np.around( np.nanstd(Tk_11[np.where( mask_C3 == 1 )] ), decimals=1)
print Tk_C1, eTk_C1
print Tk_C2, eTk_C2
print Tk_C3, eTk_C3
The average kinetic temperature in B5 is:
In [13]:
Tk_mean=int(np.round(np.nanmean(Tk_11)))
Tk_stddev=int(np.round(np.nanstd(Tk_11)))
print Tk_mean, '$\pm$', Tk_stddev
The conversion factor between NH3 Flux and mass is calculated using the SCUBA2 450um data. The SCUBA2 data has previously been regridred to match the NH3 data grid.
Assuming optically thin emission,
$$M_{dust}=
\frac{d^2\,F_\nu}{\kappa_\nu B_\nu(T)} =
0.12M_\odot \left( e^{1.439(\textrm{mm}/\lambda)(10\textrm{K}/T)}-1 \right)
\left( \frac{\kappa_\nu}{0.01 \textrm{cm}^2 \textrm{g}^{-1}} \right)^{-1}
\left( \frac{S_\nu}{ \textrm{Jy}} \right)
\left( \frac{D}{100 \textrm{pc}} \right)^2
\left( \frac{\lambda}{\textrm{mm}} \right)^3
=0.36 M_\odot \left( \frac{S_\nu}{ \textrm{Jy}} \right)~,$$
where $d$ is the distance, $F_\nu$ is the total flux, $\kappa_{\nu}$ is the dust opacity
per unit mass at frequency $\nu$, and
$B_\nu(T)$ is the Planck function at temperature $T$.
We use a dust to gas ratio of 0.01, and a dust opacity per unit mass of
$\kappa_\nu = 0.1\,(\nu/1\,\textrm{THz})^2$,
which at 450$\mu$m gives
$\kappa_{450\mu\textrm{m}}=0.044\textrm{cm}^2\textrm{g}^{-1}$.
We assume a temperature of T=10\,K, which is the typical temperature of starless cores in Perseus, and distance of 250\,pc.
In the case of B5-Cond1, this gives a mass estimate of 0.39\,M$_\odot$.
The 450um SCUBA2 flux over B5-Cond1 is calculated removing the background.
The mass determination uncertainty is due to the uncertainty in the temperature: $$\sigma_M = M \frac{1}{ \left(1-e^{-1.439(\textrm{mm}/\lambda)(10\textrm{K}/T)} \right) }\frac{1.439 \textrm{mm}}{\lambda} \left(\frac{10\textrm{K}}{T}\right) \frac{\sigma_T}{T}$$ for T=9K and $\sigma_T$=1K it gives $\sigma_M/M=0.325$, while for T=10K and $\sigma_T$=1K it gives $\sigma_M/M=0.266$.
In [14]:
scuba2=fits.getdata(file_scuba2_450)
scuba2=np.squeeze(scuba2)
kappa_450=0.044 # cm^2 g^-1
scuba2_wavelength=0.450 # mm
scuba2_raw_flux_C1 = scuba2[np.where( mask_C1 == 1 )]
scuba2_leaf_flux_C1 = np.cumsum(np.sort(scuba2_raw_flux_C1)[::-1]-np.min(scuba2_raw_flux_C1))
Temperature=10.0
print 'Total SCUBA2 450um flux over leaf of B5-Cond1:'
print 'Flux :', np.max(scuba2_leaf_flux_C1),'+-', scuba2_450_rms*np.sqrt(len(scuba2_leaf_flux_C1)), 'Jy'
Jy2Msun=0.12*(np.exp(1.439 *10./Temperature /scuba2_wavelength)-1)*(0.01/kappa_450)*(distance/100.)**2 *(scuba2_wavelength)**3
print 'Conversion factor between SCUBA2 450um flux (Jy) and mass (Msun)',Jy2Msun
Using the conversion factors determined above I calculate the mass in the leaf, and also I convert the radii from pixel units to astronomical units.
In [15]:
conv_ang = pixel_size * distance
conv_mass=np.max(scuba2_leaf_flux_C1) * Jy2Msun /np.max(w_C1)
M_C1=w_C1*conv_mass
M_C2=w_C2*conv_mass
M_C3=w_C3*conv_mass
eM_C1=np.sqrt( (nh3_rms*np.sqrt(np.arange(1,len(M_C1)+1))*conv_mass)**2 + (M_C1*0.266)**2 )
eM_C2=np.sqrt( (nh3_rms*np.sqrt(np.arange(1,len(M_C2)+1))*conv_mass)**2 + (M_C2*0.266)**2 )
eM_C3=np.sqrt( (nh3_rms*np.sqrt(np.arange(1,len(M_C3)+1))*conv_mass)**2 + (M_C3*0.266)**2 )
rp_C1 = r_C1*conv_ang
rp_C2 = r_C2*conv_ang
rp_C3 = r_C3*conv_ang
print np.max(M_C1), conv_mass
To calculate the virial parameter, we need to determine the average linewidth within the condensation. I load the velocity dispersion obtained from the line fitting, and the calculate the average velocity dispersion for each condensation.
For the combined data (GBT+VLA) I use the results of fitting the VLA+GBT data reduction that are of good quality.
In [16]:
def c_sound( A, Tk=10):
return 0.0908537*np.sqrt(Tk/float(A))
cs_h2 =c_sound(2.33)
cs_nh3=c_sound(17.0)
dv_11 = fits.getdata(dv_file)
# Calculate mean and stddev of sigma_v for each condensation
dv_C1 = np.mean(dv_11[np.where( mask_C1 == 1 )] )
dv_C2 = np.mean(dv_11[np.where( mask_C2 == 1 )] )
dv_C3 = np.mean(dv_11[np.where( mask_C3 == 1 )] )
edv_C1 = np.std(dv_11[np.where( mask_C1 == 1 )] )
edv_C2 = np.std(dv_11[np.where( mask_C2 == 1 )] )
edv_C3 = np.std(dv_11[np.where( mask_C3 == 1 )] )
# calculate the total velocity dispersion (dv_nt) and uncertainty of the mean particle for each condensation
dv_tot_C1 = np.sqrt( dv_C1**2 - cs_nh3**2 + cs_h2**2)
dv_tot_C2 = np.sqrt( dv_C2**2 - cs_nh3**2 + cs_h2**2)
dv_tot_C3 = np.sqrt( dv_C3**2 - cs_nh3**2 + cs_h2**2)
edv_tot_C1 = (dv_C1/dv_tot_C1)*edv_C1
edv_tot_C2 = (dv_C2/dv_tot_C2)*edv_C2
edv_tot_C3 = (dv_C3/dv_tot_C3)*edv_C3
# Uncertainty in mass calculation, sigma_M/M = 0.266
dM_M=0.266
# Calculate virial parameter
alpha_C1=928* rp_C1*(u.au).to(u.pc)/M_C1 * dv_tot_C1**2
alpha_C2=928* rp_C2*(u.au).to(u.pc)/M_C2 * dv_tot_C2**2
alpha_C3=928* rp_C3*(u.au).to(u.pc)/M_C3 * dv_tot_C3**2
ealpha_C1 = alpha_C1*np.sqrt( (eM_C1/M_C1)**2 + (2*edv_tot_C1/dv_tot_C1)**2 )
ealpha_C2 = alpha_C2*np.sqrt( (eM_C2/M_C2)**2 + (2*edv_tot_C2/dv_tot_C2)**2 )
ealpha_C3 = alpha_C3*np.sqrt( (eM_C3/M_C3)**2 + (2*edv_tot_C3/dv_tot_C3)**2 )
print dv_tot_C1, edv_tot_C1, dv_tot_C2, edv_tot_C2, dv_tot_C3, edv_tot_C3
print np.round(np.sqrt( dv_C1**2 - cs_nh3**2)/cs_h2, decimals=2), np.round(np.sqrt( dv_C2**2 - cs_nh3**2)/cs_h2, decimals=2), np.round(np.sqrt( dv_C3**2 - cs_nh3**2)/cs_h2, decimals=2)
print np.median(rp_C1), np.median(M_C1), np.median(dv_tot_C1)
print np.median(alpha_C1), np.median(alpha_C2), np.median(alpha_C3)
# minimum radius where alpha<2
factor_round=100
print np.round(rp_C1[np.argmin(np.abs(alpha_C1-2))]/factor_round, decimals=0)*factor_round
print np.round(rp_C2[np.argmin(np.abs(alpha_C2-2))]/factor_round, decimals=0)*factor_round
print np.round(rp_C3[np.argmin(np.abs(alpha_C3-2))]/factor_round, decimals=0)*factor_round
In [17]:
# Figure width for Nature in two column format is 183 mm ~7.2in
fig, ax = plt.subplots(1, 2, figsize=(7.2,3))
color_C1='#E41A1C'
color_C2='#377EB8'
color_C3='#4DAF4A'
draw_error=False
ax[0].plot(rp_C1, M_C1, label='B5-Cond1', color=color_C1)
ax[0].plot(rp_C2, M_C2, label='B5-Cond2', color=color_C2)
ax[0].plot(rp_C3, M_C3, label='B5-Cond3', color=color_C3)
# Add half mass radius
id_mass_C1=np.argmin( np.abs(M_C1 - M_C1[-1]*0.5))
id_mass_C2=np.argmin( np.abs(M_C2 - M_C2[-1]*0.5))
id_mass_C3=np.argmin( np.abs(M_C3 - M_C3[-1]*0.5))
ax[0].scatter(rp_C1[id_mass_C1], M_C1[id_mass_C1], c=color_C1, zorder=11, alpha=0.5)
ax[0].scatter(rp_C2[id_mass_C2], M_C2[id_mass_C2], c=color_C2, zorder=12, alpha=0.5)
ax[0].scatter(rp_C3[id_mass_C3], M_C3[id_mass_C3], c=color_C3, zorder=13, alpha=0.5)
print np.round(rp_C1[id_mass_C1]/100)*100, np.round(rp_C2[id_mass_C2]/100)*100,np.round(rp_C3[id_mass_C3]/100)*100
if draw_error:
ax[0].fill_between(rp_C1, M_C1-eM_C1, M_C1+eM_C1, color=color_C1, alpha=0.1)
ax[0].fill_between(rp_C2, M_C2-eM_C2, M_C2+eM_C2, color=color_C2, alpha=0.1)
ax[0].fill_between(rp_C3, M_C3-eM_C3, M_C3+eM_C3, color=color_C3, alpha=0.1)
# Add legend
ax[0].legend(loc=0, frameon=False, scatterpoints=1, numpoints=1, prop={'size':11}, handlelength=1)
# Add grey area to show beam size
ax[0].axvspan(0,3*conv_ang, color='gray', alpha=0.5, zorder=1)
ax[0].text( 3*conv_ang, 0.2, 'Beam', rotation=270., verticalalignment='center')
ax[0].text( 3*conv_ang*0.5, 0.2, 'size', rotation=270., verticalalignment='center')
# Calculate and plor M(r) for two density radial profiles
## r_MR=np.arange(0.004,0.009,0.0005) # in pc
r_MR=np.arange(800,1900,100.) # in au
M_15=0.07*np.power(r_MR/800., 3-1.5)
M_2 =0.07*np.power(r_MR/800., 3-2)
ax[0].plot( r_MR, M_15, ':', color='black')
ax[0].plot( r_MR, M_2, '--', color='black')
ax[0].text( np.max(r_MR)+20, 0.9*np.max(M_15), "$\\rho\\propto r^{-1.5}$",fontsize=14)
ax[0].text( np.max(r_MR)+20, np.max(M_2), "$\\rho\\propto r^{-2}$",fontsize=14)
#
# Virial parameter
#
ax[1].plot(rp_C1, alpha_C1, label='B5-Cond1', color=color_C1)
ax[1].plot(rp_C2, alpha_C2, label='B5-Cond2', color=color_C2)
ax[1].plot(rp_C3, alpha_C3, label='B5-Cond3', color=color_C3)
ax[1].scatter(rp_C1[id_mass_C1], alpha_C1[id_mass_C1], c=color_C1, zorder=11, alpha=0.5)
ax[1].scatter(rp_C2[id_mass_C2], alpha_C2[id_mass_C2], c=color_C2, zorder=12, alpha=0.5)
ax[1].scatter(rp_C3[id_mass_C3], alpha_C3[id_mass_C3], c=color_C3, zorder=13, alpha=0.5)
print np.round(alpha_C1[id_mass_C1], decimals=2), np.round(alpha_C2[id_mass_C2], decimals=2), np.round(alpha_C3[id_mass_C3], decimals=2)
if draw_error:
ax[1].fill_between(rp_C1, alpha_C1-ealpha_C1, alpha_C1+ealpha_C1, color=color_C1, alpha=0.1)
ax[1].fill_between(rp_C2, alpha_C2-ealpha_C2, alpha_C2+ealpha_C2, color=color_C2, alpha=0.1)
ax[1].fill_between(rp_C3, alpha_C3-ealpha_C3, alpha_C3+ealpha_C3, color=color_C3, alpha=0.1)
ax[1].set_ylim([1,10])
ax[1].legend( frameon=False, scatterpoints=1, numpoints=1, prop={'size':11}, handlelength=1)
ax[1].axhline(y=2, color='black')
ax[1].text( 2000, 2.1, 'Bound', horizontalalignment='center')
ax[1].arrow(2000, 2.0, 0.0, -0.25, fc='k', ec='k', head_width=125, head_length=0.1, width=0.05)
ax[1].axvspan(0, 3*conv_ang, color='gray', alpha=0.5, zorder=1)
ax[1].text( 3*conv_ang, 1.5, 'Beam', rotation=270., verticalalignment='center')
ax[1].text( 3*conv_ang*0.5, 1.5, 'size', rotation=270., verticalalignment='center')
ax[1].set_yscale("log")
plt.setp( ax[0], xlabel="Condensation effective radius, $R_{eff}$ (au)")
plt.setp( ax[0], ylabel="Condensation Mass ($M_{Sun}$)")
plt.setp( ax[1], xlabel="Condensation effective radius, $R_{eff}$ (au)")
ax[1].set_ylabel("Condensation virial parameter, $\\alpha$", labelpad=-6)
ax[0].locator_params(axis = 'both', nbins = 4)
ax[1].locator_params(axis = 'x', nbins = 4)
ax[0].tick_params(which='major', length=6)
ax[1].tick_params(which='minor', length=6)
ax[1].tick_params(which='major', axis='x', length=6)
# Add a) and b) labels
ax[0].text(0.15,0.90, 'a)', color='black', transform=ax[0].transAxes)
ax[1].text(0.15,0.90, 'b)', color='black', transform=ax[1].transAxes)
if draw_error:
ax[0].set_ylim([0.0,0.42])
ax[1].set_ylim([0.9,10])
else:
ax[0].set_xlim([0,3000])
ax[1].set_xlim([0,3000])
ax[0].set_ylim([0.0,0.4])
ax[1].set_ylim([1.0,10])
if save_figure:
plt.savefig(figure_dir+'B5_Mass_Alpha_Radius.pdf', bbox_inches='tight')
plt.savefig(figure_dir+'B5_Mass_Alpha_Radius.eps', bbox_inches='tight')
plt.savefig(figure_dir+'B5_Mass_Alpha_Radius.png', bbox_inches='tight', dpi=200)
In [18]:
print np.max(M_C1), np.max(eM_C1), np.min(alpha_C1), alpha_C1[-1]
print np.max(M_C2), np.max(eM_C3), np.min(alpha_C2), alpha_C2[-1]
print np.max(M_C3), np.max(eM_C3), np.min(alpha_C3), alpha_C3[-1]
The free fall time is calculated as: $$ t_{ff}=\sqrt{\frac{3\pi}{32 G \rho}} = \sqrt{\frac{\pi^2R^3}{8GM}} = 1.6\times 10^6\left(\frac{R}{0.1\textrm{pc}}\right)^{3/2} \left(\frac{M}{0.1M_\odot}\right)^{-1/2}\,\textrm{yr}$$
In [19]:
print 't_ff B5-Cond1=', int(1.6e6 * np.sqrt(np.max(rp_C1*(u.au).to(u.pc))/0.1)**3 * np.sqrt(0.1/np.max(M_C1))), 'yr'
print 't_ff B5-Cond2=', int(1.6e6 * np.sqrt(np.max(rp_C2*(u.au).to(u.pc))/0.1)**3 * np.sqrt(0.1/np.max(M_C2))), 'yr'
print 't_ff B5-Cond3=', int(1.6e6 * np.sqrt(np.max(rp_C3*(u.au).to(u.pc))/0.1)**3 * np.sqrt(0.1/np.max(M_C3))), 'yr'
Therefore, we can use a typical free-fall timescale of 40,000 yr
Here we create the two panel figure that shows the NH$_3$ integrated emission annotated in the first panel, while the second panel is the schematic view of the region.
For the second panel, we need to create the mask with all identified condensations. Also, we need to create an image that shows the typical distance travelled over a free-fall time, for this we use the custom function defined below makeGaussian.
Here we define the makeGaussian function and calculate the distance traveled at the sound speed over a free fall time.
In [20]:
def makeGaussian(size_x, size_y, fwhm = 3, center=None):
""" Make a square gaussian kernel.
size is the length of a side of the square
fwhm is full-width-half-maximum, which
can be thought of as an effective radius.
"""
x = np.arange(0, size_x, 1, float)
xy= np.arange(0, size_y, 1, float)
y = xy[:,np.newaxis]
if center is None:
x0 = size_x // 2
y0 = size_y // 2
else:
x0 = center[0]
y0 = center[1]
return np.exp(-4*np.log(2) * ((x-x0)**2 + (y-y0)**2) / fwhm**2)
v_sound=0.2*u.km / u.s
t_ff = 4e4*u.yr
radius_move_pc=( v_sound * t_ff).to(u.pc)
pix_size=(pixel_size*u.arcsec).to(u.radian).value * (distance * u.pc)
radius_move_pix=(radius_move_pc/pix_size).value
Here we define some of the variables needed to annotate the image
In [21]:
ra_c = ( 3+ (47+38./60.)/60.)*15.
dec_c = (32+(52.+50./60.)/60.)
text_label=['B5-IRS1', 'B5-Cond1', 'B5-Cond2', 'B5-Cond3']
align=['left', 'left', 'center', 'center']
v_align=['center', 'center', 'bottom', 'top']
ra_label =np.array([ (3+(47+38./60.)/60.), (3+(47+36./60.)/60.), (3+(47+44./60.)/60.), (3+(47+44./60.)/60.)]) *15.
dec_label=np.array([32+(51+30/60.)/60., 32+(52+30/60.)/60.,32+(52+15/60.)/60.,32+(51+00/60.)/60.])
dx_arrow=np.array([40,35,-30,-18])/3600.
dy_arrow=np.array([10,-15,-15,22])/3600.
sz=data.shape
G1=np.array([[makeGaussian( sz[1], sz[0], fwhm=2*radius_move_pix, center=xxyy[0,:] )]])
G2=np.array([[makeGaussian( sz[1], sz[0], fwhm=2*radius_move_pix, center=xxyy[1,:] )]])
G3=np.array([[makeGaussian( sz[1], sz[0], fwhm=2*radius_move_pix, center=xxyy[2,:] )]])
G4=np.array([[makeGaussian( sz[1], sz[0], fwhm=2*radius_move_pix, center=xxyy_yso[0,:] )]])
# Create masks with gaussians showing how far they move
mask_hdu_1 = fits.PrimaryHDU(G1, header)
mask_hdu_2 = fits.PrimaryHDU(G2, header)
mask_hdu_3 = fits.PrimaryHDU(G3, header)
mask_hdu_4 = fits.PrimaryHDU(G4, header)
mask_hdu_all= fits.PrimaryHDU(G1+G2+G3+G4, header)
# create masks for condensation boundaries
mask_hdu_C1 = fits.PrimaryHDU(mask_C1.astype('short'), header)
mask_hdu_C2 = fits.PrimaryHDU(mask_C2.astype('short'), header)
mask_hdu_C3 = fits.PrimaryHDU(mask_C3.astype('short'), header)
Now we create the figures
In [22]:
fig=plt.figure(figsize=( 7.2, 5.5))
color_table='Blues_r'
star_size=80
fig0=aplpy.FITSFigure(data_file, hdu=0, figure=fig, subplot=[0.10,0.1,0.45,0.89])
fig1=aplpy.FITSFigure(mask_hdu_all, hdu=0, figure=fig, subplot=[0.55,0.1,0.45,0.89])
fig0.show_colorscale( cmap=color_table, vmin=-0.035, vmax=0.27)
fig1.show_grayscale( vmin=0.5, vmax=1.1, invert=True)
fig0.recenter(ra_c, dec_c, width=0.07, height=0.1) # degrees
fig1.recenter(ra_c, dec_c, width=0.07, height=0.1) # degrees
# circle of motion
fig1.show_contour(mask_hdu_1, colors='green', linewidths=1.0, levels=[0.5], zorder=30,linestyles='dotted',slices=[0, 0])
fig1.show_contour(mask_hdu_2, colors='green', linewidths=1.0, levels=[0.5], zorder=31,linestyles='dotted',slices=[0, 0])
fig1.show_contour(mask_hdu_3, colors='green', linewidths=1.0, levels=[0.5], zorder=32,linestyles='dotted',slices=[0, 0])
fig1.show_contour(mask_hdu_4, colors='green', linewidths=1.0, levels=[0.5], zorder=33,linestyles='dotted',slices=[0, 0])
# filament outline
fig1.show_contour(data_file, colors='black', linewidths=0.5, levels=[0.1], zorder=34)
# Beams
fig0.add_beam()
fig0.beam.set_color('#E31A1C')
fig0.ticks.set_minor_frequency(4)
fig1.add_beam()
fig1.beam.set_color('#E31A1C')
fig1.ticks.set_minor_frequency(4)
# YSO
fig0.show_markers(ra_yso, dec_yso, marker='*', alpha=0.7, layer='lay_yso', c='orange', s=star_size)
fig0.show_markers(ra_blob, dec_blob, marker='o', alpha=0.7, layer='lay_blob', c='orange')
for i in np.arange(len(text_label)):
fig0.add_label(ra_label[i], dec_label[i], text_label[i], horizontalalignment=align[i], color='yellow',zorder=61+i, verticalalignment=v_align[i])
fig0.show_arrows(ra_label[i], dec_label[i], dx_arrow[i], dy_arrow[i], alpha=0.5, color='red')
# Add YSO
fig1.show_markers(ra_yso, dec_yso, marker='*', alpha=0.7, layer='lay_yso', c='orange', s=star_size)
ang_sep = 1e4/distance / 3600. # from arcsec -> deg
fig0.add_scalebar(ang_sep)
fig0.scalebar.set(color='white')
fig0.scalebar.set_label('10,000 au')
fig1.add_scalebar(ang_sep)
fig1.scalebar.set(color='black')
fig1.scalebar.set_label('10,000 au')
# Add contour with condensation boundaries
fig0.show_contour(mask_hdu_C1, colors='red', linewidths=0.5)
fig0.show_contour(mask_hdu_C2, colors='red', linewidths=0.5)
fig0.show_contour(mask_hdu_C3, colors='red', linewidths=0.5)
fig1.show_contour(mask_hdu_C1, colors='red', linewidths=0.5)
fig1.show_contour(mask_hdu_C2, colors='red', linewidths=0.5)
fig1.show_contour(mask_hdu_C3, colors='red', linewidths=0.5)
#
fig0.axis_labels.set_xtext('Right Ascension (J2000)')
fig0.axis_labels.set_ytext('Declination (J2000)')
fig1.axis_labels.set_xtext('')
fig1.axis_labels.set_ytext('')
fig0.tick_labels.set_xformat('hh:mm:ss')
fig0.tick_labels.set_yformat('dd:mm')
fig1.tick_labels.set_xformat('hh:mm:ss')
fig1.tick_labels.set_yformat('dd:mm')
fig1.tick_labels.hide()
fig0.ticks.set_color('black')
fig0.set_nan_color('gray')
fig1.ticks.set_color('black')
fig1.set_nan_color('gray')
fig0.add_label(0.85,0.95, 'Barnard 5\nJVLA NH$_3$(1,1)', relative=True, color='white')
fig0.add_label(0.05,0.95, 'a)', relative=True, color='white')
fig1.add_label(0.05,0.95, 'b)', relative=True, color='black')
if save_figure:
plt.savefig(figure_dir+'B5_sketch_dendro.pdf', bbox_inches='tight')
plt.savefig(figure_dir+'B5_sketch_dendro.png', bbox_inches='tight')
Make SCUBA2 450um annotated figure
In [23]:
fig=plt.figure(figsize=( 3.5, 5.5))
fig1=aplpy.FITSFigure(file_scuba2_450, hdu=0, figure=fig)
fig1.recenter(ra_c, dec_c, width=0.07, height=0.1)
cont_lev_450=np.arange(3*scuba2_450_rms,45e-3,3*scuba2_450_rms) # in Jy/pixel
fig1.show_colorscale( cmap='Reds_r', vmin=-0.0001, vmax=0.0014)
fig1.show_contour( file_scuba2_450, levels=cont_lev_450, colors='black', kernel='box', smooth=5)
# Add beams
fig1.add_beam()
fig1.beam.set_color('white')
fig1.add_label(0.78,0.93, 'Barnard 5\nSCUBA2 450$\\mu$m', relative=True, color='white')
fig1.ticks.set_minor_frequency(4)
# YSO
fig1.show_markers(ra_yso, dec_yso, marker='*', alpha=0.9, layer='lay_yso', c='gold', zorder=50)
fig1.show_markers(ra_blob, dec_blob, marker='o', alpha=0.9, layer='lay_blob', c='gold', zorder=51)
for i in np.arange(len(text_label)):
fig1.add_label(ra_label[i], dec_label[i], text_label[i], horizontalalignment=align[i], color='yellow',zorder=51+i, verticalalignment=v_align[i])
fig1.show_arrows(ra_label[i], dec_label[i], dx_arrow[i], dy_arrow[i], alpha=0.5, color='blue',zorder=51+i)
# Scalebar 10,000 au
ang_sep = 1e4/distance / 3600. # from arcsec -> deg
fig1.add_scalebar(ang_sep)
fig1.scalebar.set(color='white')
fig1.scalebar.set_label('10,000 au')
fig1.axis_labels.set_xtext('Right Ascension (J2000)')
fig1.axis_labels.set_ytext('Declination (J2000)')
fig1.tick_labels.set_xformat('hh:mm:ss')
fig1.tick_labels.set_yformat('dd:mm')
fig1.ticks.set_color('white')
fig1.set_nan_color('gray')
if save_figure:
plt.savefig(figure_dir+'B5-450-anotated.pdf', bbox_inches='tight')
plt.savefig(figure_dir+'B5-450-anotated.png', bbox_inches='tight')
Optional annotated figure with SCUBA2 450um and 850um maps
In [24]:
fig=plt.figure(figsize=( 7.2, 5.5))
fig1=aplpy.FITSFigure(file_scuba2_450, hdu=0, figure=fig, subplot=[0.10,0.1,0.45,0.89], slices=[0])
fig2=aplpy.FITSFigure(file_scuba2_850, hdu=0, figure=fig, subplot=[0.55,0.1,0.45,0.89], slices=[0])
fig1.recenter(ra_c, dec_c, width=0.07, height=0.1)
fig2.recenter(ra_c, dec_c, width=0.07, height=0.1)
cont_lev_450=np.arange(3*scuba2_450_rms,45e-3,3*scuba2_450_rms) # in Jy/pixel
fig1.show_colorscale( cmap='Reds_r', vmin=-0.0001, vmax=0.0014)
fig1.show_contour( file_scuba2_450, levels=cont_lev_450, colors='black', kernel='box', smooth=5, slices=[0])
cont_lev_850=np.arange(3*scuba2_850_rms,45e-3,3*scuba2_850_rms) # in Jy/pixel
fig2.show_colorscale( cmap='Greens_r', vmin=-0.08e-3, vmax=0.3e-3)
fig2.show_contour( file_scuba2_850, levels=cont_lev_850, colors='black', kernel='box', smooth=5, slices=[0])
fig1.add_beam()
fig1.beam.set_color('white')
fig2.add_beam()
fig2.beam.set_color('white')
fig1.add_label(0.825,0.925, 'Barnard 5\nSCUBA2 450$\\mu$m', relative=True, color='white')
fig2.add_label(0.825,0.925, 'Barnard 5\nSCUBA2 850$\\mu$m', relative=True, color='white')
fig1.add_label(0.075,0.925, 'a)', relative=True, color='white')
fig2.add_label(0.075,0.925, 'b)', relative=True, color='white')
fig1.ticks.set_minor_frequency(4)
fig2.ticks.set_minor_frequency(4)
# YSO
blue_color='#67a9cf'
orange_color=(0.96875,0.7890625,0.08203125)
fig1.show_markers(ra_yso, dec_yso, marker='*', alpha=0.95, layer='lay_yso', edgecolor='black', facecolor=blue_color, zorder=50, s=90)
fig1.show_markers(ra_blob, dec_blob, marker='o', alpha=0.95, layer='lay_blob',edgecolor='black', facecolor=blue_color, zorder=51, s=30)
fig2.show_markers(ra_yso, dec_yso, marker='*', alpha=0.95, layer='lay_yso', edgecolor='black', facecolor=orange_color,zorder=52, s=90)
fig2.show_markers(ra_blob, dec_blob, marker='o', alpha=0.95, layer='lay_blob',edgecolor='black', facecolor=orange_color,zorder=53, s=30)
for i in np.arange(len(text_label)):
fig1.add_label(ra_label[i], dec_label[i], text_label[i], horizontalalignment=align[i], color='yellow',zorder=51+i, verticalalignment=v_align[i])
fig1.show_arrows(ra_label[i], dec_label[i], dx_arrow[i], dy_arrow[i], alpha=0.5, color='blue',zorder=51+i)
# Scalebar 10,000 au
ang_sep = 1e4/distance / 3600. # from arcsec -> deg
fig1.add_scalebar(ang_sep)
fig1.scalebar.set(color='white')
fig1.scalebar.set_label('10,000 au')
#labels and tick labels
fig1.axis_labels.set_xtext('Right Ascension (J2000)')
fig1.axis_labels.set_ytext('Declination (J2000)')
fig2.axis_labels.set_xtext('')
fig2.axis_labels.set_ytext('')
fig1.tick_labels.set_xformat('hh:mm:ss')
fig1.tick_labels.set_yformat('dd:mm')
fig2.tick_labels.set_xformat('hh:mm:ss')
fig2.tick_labels.set_yformat('dd:mm')
fig2.tick_labels.hide()
fig1.ticks.set_color('white')
if save_figure:
plt.savefig(figure_dir+'B5-450-850-anotated.pdf', bbox_inches='tight')
plt.savefig(figure_dir+'B5-450-850-anotated.eps', bbox_inches='tight')
plt.savefig(figure_dir+'B5-450-850-anotated.png', bbox_inches='tight')
The centroid velocity (v$_c$) map is loaded into the vc_map variable.
In [25]:
vc_map=fits.getdata(vc_file)
The centroid velocity to be used towards the condensations and the protostar is calculated as the average value in a box of size equals to the beamsize. The uncertainty is estimated by the standard deviation of the values within the box. The mean and standard deviation are rounded to the same precision.
We compare these values to the value at the center of the box.
In [26]:
Vc_YSO=vc_map[ np.round(xxyy_yso[0,1]), np.round(xxyy_yso[0,0])]
Vc_C1=vc_map[ np.round(xxyy[0,1]), np.round(xxyy[0,0])]
Vc_C2=vc_map[ np.round(xxyy[1,1]), np.round(xxyy[1,0])]
Vc_C3=vc_map[ np.round(xxyy[2,1]), np.round(xxyy[2,0])]
In [27]:
mVc_YSO=np.round(np.mean(vc_map[ np.round(xxyy_yso[0,1])-6:np.round(xxyy_yso[0,1])+6, np.round(xxyy_yso[0,0])-6:np.round(xxyy_yso[0,0])+6]), decimals=2)
mVc_C1 =np.round(np.mean(vc_map[ np.round(xxyy[0,1])-6:np.round(xxyy[0,1])+6, np.round(xxyy[0,0])-6:np.round(xxyy[0,0])+6]), decimals=2)
mVc_C2 =np.round(np.mean(vc_map[ np.round(xxyy[1,1])-6:np.round(xxyy[1,1])+6, np.round(xxyy[1,0])-6:np.round(xxyy[1,0])+6]), decimals=2)
mVc_C3 =np.round(np.mean(vc_map[ np.round(xxyy[2,1])-6:np.round(xxyy[2,1])+6, np.round(xxyy[2,0])-6:np.round(xxyy[2,0])+6]), decimals=2)
dVc_YSO=np.round(np.std(vc_map[ np.round(xxyy_yso[0,1])-6:np.round(xxyy_yso[0,1])+6, np.round(xxyy_yso[0,0])-6:np.round(xxyy_yso[0,0])+6]), decimals=2)
dVc_C1 =np.round(np.std(vc_map[ np.round(xxyy[0,1])-6:np.round(xxyy[0,1])+6, np.round(xxyy[0,0])-6:np.round(xxyy[0,0])+6]), decimals=2)
dVc_C2 =np.round(np.std(vc_map[ np.round(xxyy[1,1])-6:np.round(xxyy[1,1])+6, np.round(xxyy[1,0])-6:np.round(xxyy[1,0])+6]), decimals=2)
dVc_C3 =np.round(np.std(vc_map[ np.round(xxyy[2,1])-6:np.round(xxyy[2,1])+6, np.round(xxyy[2,0])-6:np.round(xxyy[2,0])+6]), decimals=2)
print 'center mean stddev'
print Vc_YSO, mVc_YSO, dVc_YSO
print Vc_C1, mVc_C1, dVc_C1
print Vc_C2, mVc_C2, dVc_C2
print Vc_C3, mVc_C3, dVc_C3
For the pair B5-IRS1 and B5-Cond2, we calculate the ratio between kinetic and gravitational energy as a function of final stellar mass for B5-Cond2. Also, plot this.
In [28]:
def Check_Bound_Binary(m1, m2, v1, v2, r):
# m1 in Msun
# m2 in Msun
# v1 V_LSR for m1
# v2 V_LSR for m2
# r in AU
#
v_cm=(m1*v1 + v2*m2)/(m1+m2)
dv1=v1-v_cm
dv2=v2-v_cm
U = -0.887*(m1*m2)/(r*1e-3)
K = m1*dv1*dv1*0.5 + m2*dv2*dv2*0.5
return K+U, -K/U
def Check_Bound_Binary_v3d(m1, m2, v1, v2, r):
# m1 in Msun
# m2 in Msun
# v1 V_LSR for m1
# v2 V_LSR for m2
# r in AU
#
v_cm=(m1*v1 + v2*m2)/(m1+m2)
dv1=v1-v_cm
dv2=v2-v_cm
U = -0.887*(m1*m2)/(r*1e-3)
K = m1*dv1*dv1*0.5 + m2*dv2*dv2*0.5
return 3*K + U, -3*K / U
# Mass range
Ma = 0.1 # Msun
Mbmin = 0.01 # Msun
Mbmax = 0.2 # Msun
dMb = 0.01 # Msun
# Velocities
Va = mVc_YSO # km/s
Vb = mVc_C2 # km/s
separation= sep_C2_YSO # 3570.85 # AU
M_fin_C2=np.max(M_C2)
Mb_range= np.arange(Mbmin, Mbmax+dMb, dMb)
Energy_1 = np.zeros(len(Mb_range))
Energy_2 = np.zeros(len(Mb_range))
Energy_3 = np.zeros(len(Mb_range))
Energy_4 = np.zeros(len(Mb_range))
KU_1 = np.zeros(len(Mb_range))
KU_2 = np.zeros(len(Mb_range))
KU_3 = np.zeros(len(Mb_range))
KU_4 = np.zeros(len(Mb_range))
for i in np.arange(len(Mb_range)):
Energy_1[i], KU_1[i] = Check_Bound_Binary(Ma, Mb_range[i], Va, Vb, separation)
Energy_2[i], KU_2[i] = Check_Bound_Binary_v3d(Ma, Mb_range[i], Va, Vb, separation)
Energy_3[i], KU_3[i] = Check_Bound_Binary(Ma, Mb_range[i], Va, Vb, separation*np.sqrt(2))
Energy_4[i], KU_4[i] = Check_Bound_Binary_v3d(Ma, Mb_range[i], Va, Vb, separation*np.sqrt(2))
color=['#D7191C', '#FDAE61', '#ABD9E9', '#2C7BB6']
labels=['a) 1D results', 'b) 3D velocities', 'c) 3D separation', 'd) 3D separation and velocities']
####
#### Energy ratio
####
fig=plt.figure( figsize=( 5.0, 4.5) )
ax= plt.subplot()
plt.plot( Mb_range, KU_1, label=labels[0], color=color[0])#, s=10**2)
plt.plot( Mb_range, KU_2, label=labels[1], color=color[1])#, s=10**2)
plt.plot( Mb_range, KU_3, label=labels[2], color=color[2])#, s=10**2)
plt.plot( Mb_range, KU_4, label=labels[3], color=color[3])#, s=10**2)
plt.setp( ax, xlabel="Final stellar mass from condensation B5-Cond2 (M$_{Sun}$)")
plt.setp( ax, ylabel="Kinetic Energy/$|$Gravitational Energy$|$")
ax.axvline(x=M_fin_C2*0.5, ymin=0.03, ymax=0.2, linestyle='-', color='gray')
ax.axvline(x=M_fin_C2*0.3, ymin=0.03, ymax=0.2, linestyle='-', color='gray')
ax.axvline(x=M_fin_C2*0.15, ymin=0.03, ymax=0.2, linestyle='-', color='gray')
ax.text( M_fin_C2*0.495, 0.0015, '$\\epsilon_{cond}$=0.5', horizontalalignment='right')
ax.text( M_fin_C2*0.295, 0.0015, '$\\epsilon_{cond}$=0.3', horizontalalignment='right')
ax.text( M_fin_C2*0.145, 0.0015, '$\\epsilon_{cond}$=0.15', horizontalalignment='right')
ax.axvline(x=0.076, ymin=0.24, ymax=0.58, linestyle='-', color='black', linewidth=2.0)
ax.arrow(0.073, 0.0095, -0.025, 0.0, fc='k', ec='k', head_width=0.00075, head_length=0.0075, width=0.00005)
ax.arrow(0.079, 0.0095, +0.025, 0.0, fc='k', ec='k', head_width=0.00075, head_length=0.0075, width=0.00005)
ax.text( 0.072, 0.0105, 'Brown Dwarfs', zorder=30, horizontalalignment='right')
ax.text( 0.080, 0.0105, 'Stars', zorder=31, horizontalalignment='left')
ax.legend( loc=1, frameon=False, scatterpoints=1, numpoints=1)
plt.xlim(0.0, 0.2)
plt.ylim(0.0, 0.033)
ax.locator_params(axis = 'both', nbins = 5)
#ax.locator_params(axis = 'x', nbins = 4)
if save_figure:
plt.savefig(figure_dir+'B5_KU_binary_line.pdf', bbox_inches='tight')
plt.savefig(figure_dir+'B5_KU_binary_line.png', bbox_inches='tight')
In [29]:
print M_fin_C2
For the pair B5-IRS1 and B5-Cond2, we calculate the ratio between kinetic and gravitational energy as a function of final stellar mass for B5-Cond2. Also, plot this.
In [30]:
def plot_virial_ratio(Vel_Cond, Vel_YSO, Mass_Cond, Mass_YSO, linear_separation, ax_i):
"""
Function to calculate the virial ratio of one Condensation and the YSO.
It also receives an axis object and plots the virial ratios.
"""
color=['#D7191C', '#FDAE61', '#ABD9E9', '#2C7BB6']
labels=['a) 1D results', 'b) 3D velocities', 'c) 3D separation', 'd) 3D separation and velocities']
# Mass range
Mbmin = 0.01 # Msun
Mbmax = 0.2 # Msun
dMb = 0.01 # Msun
Mb_range= np.arange(Mbmin, Mbmax+dMb, dMb)
Energy_1 = np.zeros(len(Mb_range))
Energy_2 = np.zeros(len(Mb_range))
Energy_3 = np.zeros(len(Mb_range))
Energy_4 = np.zeros(len(Mb_range))
#
KU_1 = np.zeros(len(Mb_range))
KU_2 = np.zeros(len(Mb_range))
KU_3 = np.zeros(len(Mb_range))
KU_4 = np.zeros(len(Mb_range))
#
for i in np.arange(len(Mb_range)):
Energy_1[i], KU_1[i] = Check_Bound_Binary( Mass_YSO, Mb_range[i], Vel_YSO, Vel_Cond, linear_separation)
Energy_2[i], KU_2[i] = Check_Bound_Binary_v3d(Mass_YSO, Mb_range[i], Vel_YSO, Vel_Cond, linear_separation)
Energy_3[i], KU_3[i] = Check_Bound_Binary( Mass_YSO, Mb_range[i], Vel_YSO, Vel_Cond, linear_separation*np.sqrt(2))
Energy_4[i], KU_4[i] = Check_Bound_Binary_v3d(Mass_YSO, Mb_range[i], Vel_YSO, Vel_Cond, linear_separation*np.sqrt(2))
ax_i.plot( Mb_range, KU_1, label=labels[0], color=color[0])
ax_i.plot( Mb_range, KU_2, label=labels[1], color=color[1])
ax_i.plot( Mb_range, KU_3, label=labels[2], color=color[2])
ax_i.plot( Mb_range, KU_4, label=labels[3], color=color[3])
ax_i.axvline(x=Mass_Cond*0.5, ymin=0.03, ymax=0.2, linestyle='-', color='gray')
ax_i.axvline(x=Mass_Cond*0.3, ymin=0.03, ymax=0.2, linestyle='-', color='gray')
ax_i.axvline(x=Mass_Cond*0.15, ymin=0.03, ymax=0.2, linestyle='-', color='gray')
ax_i.text( Mass_Cond*0.495/Mbmax, 0.05, '$\\epsilon_{cond}$=0.5', horizontalalignment='right', transform=ax_i.transAxes)
ax_i.text( Mass_Cond*0.295/Mbmax, 0.05, '$\\epsilon_{cond}$=0.3', horizontalalignment='right', transform=ax_i.transAxes)
ax_i.text( Mass_Cond*0.145/Mbmax, 0.05, '$\\epsilon_{cond}$=0.15', horizontalalignment='right', transform=ax_i.transAxes)
ax_i.axvline(x=0.076, ymin=0.24, ymax=0.58, linestyle='-', color='black', linewidth=2.0)
fig_mp, ax_mp = plt.subplots(3, 1, figsize=(5.3,12), sharex=True)
M_fin_C1=np.max(M_C1)
M_fin_C2=np.max(M_C2)
M_fin_C3=np.max(M_C3)
plot_virial_ratio(mVc_C1, mVc_YSO, M_fin_C1, 0.1, sep_C1_YSO, ax_mp[0])
plot_virial_ratio(mVc_C2, mVc_YSO, M_fin_C2, 0.1, sep_C2_YSO, ax_mp[1])
plot_virial_ratio(mVc_C3, mVc_YSO, M_fin_C3, 0.1, sep_C3_YSO, ax_mp[2])
#ax_mp[1].text( 0.078, 0.0165, 'Brown Dwarf limit', rotation=270., zorder=30)#, transform=ax_mp[1].transAxes)
#
ax_mp[1].arrow(0.073, 0.0095, -0.025, 0.0, fc='k', ec='k', head_width=0.00075, head_length=0.0075, width=0.00005)
ax_mp[1].arrow(0.079, 0.0095, +0.025, 0.0, fc='k', ec='k', head_width=0.00075, head_length=0.0075, width=0.00005)
ax_mp[1].text( 0.072, 0.0105, 'Brown Dwarfs', zorder=30, horizontalalignment='right')
ax_mp[1].text( 0.080, 0.0105, 'Stars', zorder=31, horizontalalignment='left')
#
ax_mp[0].text( 0.75, 0.5, 'B5-Cond1', transform=ax_mp[0].transAxes, size='x-large')
ax_mp[1].text( 0.75, 0.5, 'B5-Cond2', transform=ax_mp[1].transAxes, size='x-large')
ax_mp[2].text( 0.75, 0.5, 'B5-Cond3', transform=ax_mp[2].transAxes, size='x-large')
# Fine-tune figure; make subplots close to each other and hide x ticks for
# all but bottom plot.
fig_mp.subplots_adjust(hspace=0.05)
plt.setp([a.get_xticklabels() for a in fig_mp.axes[:-1]], visible=False)
plt.setp( ax_mp[2], xlabel="Final stellar mass from condensation (M$_{Sun}$)")
plt.setp( ax_mp[2], ylabel="Kinetic Energy/$|$Gravitational Energy$|$")
ax_mp[0].legend( loc=1, frameon=False, scatterpoints=1, numpoints=1)
plt.xlim(0.0, 0.2)
#plt.ylim(0.0, 0.033)
ax_mp[0].locator_params(axis = 'both', nbins = 5)
ax_mp[1].locator_params(axis = 'both', nbins = 5)
ax_mp[2].locator_params(axis = 'both', nbins = 5)
Print V$_{LSR}$, R$_{eq}$, Flux NH$_3$, Flux 450$\mu$m SCUBA2, and Mass estimates for B5-IRS1, B5-Cond1, B5-Cond2, and B5-Cond3. The numbers are rounded off to be used in Table.
In [31]:
print 'Source Vlsr Req Flux NH3 Flux SCUBA2 Mass'
print 'B5-IRS1 :', mVc_YSO,'$\pm$',dVc_YSO, 'km\,s$^{-1}$'
print 'B5-Cond1:', mVc_C1,'$\pm$',dVc_C1, 'km\,s$^{-1}$', np.round(np.max(rp_C1), decimals=4), ' pc' , np.round(np.max(w_C1), decimals=3), '$\pm$', np.round(np.max(ew_C1), decimals=3), np.round(np.max(scuba2_leaf_flux_C1), decimals=3),'$\pm$', np.round(scuba2_450_rms *np.sqrt(len(scuba2_leaf_flux_C1)), decimals=3), np.round(Jy2Msun*np.max(scuba2_leaf_flux_C1), decimals=3),'$\pm$',np.round(dM_M*Jy2Msun*np.max(scuba2_leaf_flux_C1), decimals=3)
print 'B5-Cond2:', mVc_C2,'$\pm$',dVc_C2, 'km\,s$^{-1}$', np.round(np.max(rp_C2), decimals=4), ' pc' , np.round(np.max(w_C2), decimals=3), '$\pm$', np.round(np.max(ew_C2), decimals=3), ' ---------- ', np.round(np.max(M_C2), decimals=2), '$\pm$', np.round(np.max(eM_C2), decimals=2)
print 'B5-Cond3:', mVc_C3,'$\pm$',dVc_C3, 'km\,s$^{-1}$', np.round(np.max(rp_C3), decimals=4), ' pc' , np.round(np.max(w_C3), decimals=3), '$\pm$', np.round(np.max(ew_C3), decimals=3), ' ---------- ', np.round(np.max(M_C3), decimals=2), '$\pm$', np.round(np.max(eM_C3), decimals=2)
In [32]:
print np.round(np.min(alpha_C1), decimals=2), np.round(np.min(alpha_C2), decimals=2), np.round(np.min(alpha_C3), decimals=2)
print np.round(np.min(alpha_C1), decimals=2), np.round(np.min(alpha_C2), decimals=2), np.round(np.min(alpha_C3), decimals=2)
In [32]: