Deriving the stellar mass density around the solar system from Gaia data

In [1]:
import astropy.units as u
from astropy.coordinates.sky_coordinate import SkyCoord
from astropy.units import Quantity
from astroquery.gaia import Gaia

Created TAP+ (v1.0.1) - Connection:
	Host: gea.esac.esa.int
	Use HTTPS: False
	Port: 80
	SSL Port: 443
In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

# Suppress warnings. Comment this out if you wish to see the warning messages
import warnings
warnings.filterwarnings('ignore')
In [9]:
# Get the data from Gaia archive
job = Gaia.launch_job_async("SELECT TOP 3000000 source_id,ra,ra_error,dec,dec_error,parallax,parallax_error,phot_g_mean_mag,bp_rp,radial_velocity,radial_velocity_error,phot_variable_flag,teff_val,a_g_val FROM gaiadr2.gaia_source  WHERE (parallax>=2e1)" \
, dump_to_file=True)
print (job)
r = job.get_results()

Jobid: 1525673124602O
Phase: COMPLETED
Owner: None
Output file: async_20180507150524.vot
Results: None
In [18]:
# Plot the cumulative stellar number counts vs. the distance from the earth
plt.hist(1000.0/r['parallax'],color='r', bins=np.logspace(-1,2, 50), alpha=0.3, cumulative=True)
x=np.logspace(-1.0, 2.0, 30)

def y (x,numdens):
    return (4.0/3.0)*np.pi*(x**3)*numdens

plt.plot(x,y(x,0.1),label = "n = 0.1 per pc^3")
plt.plot(x,y(x,1.0),label = "n = 1 per pc^3")
plt.plot(x,y(x,10.0),label = "n = 10 per pc^3")

plt.legend(loc="top left")
plt.xscale('log')
plt.yscale('log')
plt.xlim(1,50)
plt.ylim(1e0,1e7)
plt.show()

If the stellar number density is n=0.1 per pc^3, then

In [21]:
#the stellar number density is
print 0.1*2e33/(4.0/3.0*np.pi*(3e18)**3) , '[g]'

1.76838825658e-24 [g]