import matplotlib.pyplot as plt
import numpy as np

read in a dataset of stars, get their g-band magnitude

from astropy.table import Table
gaia_data=Table.read('http://burro.case.edu/Academics/Astr323/HW/HW2/Gaia_bLT10_pGT5.fits')
g_mag=gaia_data['phot_g_mean_mag']

Plot a histogram of magnitude from 2-20 in steps of half a mag

magbins=np.arange(2,20.5,0.5)
plt.hist(g_mag,bins=magbins)
plt.xlabel('Magnitude')
plt.ylabel('N')

Text(0, 0.5, 'N')

now lets do it in a way that lets us work with the counts in each bin

Fig,(linPlot,logPlot) = plt.subplots(1,2,figsize=(8,4))

# use np.histogram to actually give us the counts and bins
# in variables that we can work with
N,bin_edges = np.histogram(g_mag,bins=magbins)

# calculate log N and bin centers
logN = np.log10(N)
bin_cent = 0.5*(bin_edges[1:]+bin_edges[:-1])

# make a plot of N(m)
linPlot.scatter(bin_cent,N)
linPlot.set_xlabel('Magnitude')
linPlot.set_ylabel('N')

# make a plot of logN(m)
logPlot.scatter(bin_cent,logN)
logPlot.set_xlabel('Magnitude')
logPlot.set_ylabel('log(N)')

Fig.tight_layout()

Now try fitting a slope to logN over the magnitude range 3-9

fit_range = np.logical_and(bin_cent>3,bin_cent<9)

coeff,cov=np.polyfit(bin_cent[fit_range],logN[fit_range],1,cov=True)
coeff_err = np.sqrt(np.diag(cov))
print('    slope = {:.3f} +/- {:.3f}'.format(coeff[0],coeff_err[0]))
print('intercept = {:.3f} +/- {:.3f}'.format(coeff[1],coeff_err[1]))
polynomial=np.poly1d(coeff)


# plot the data
plt.scatter(bin_cent,logN)
plt.xlabel('Magnitude')
plt.ylabel('log(N)')

# overplot your fit
xfit=np.arange(3,9)
plt.plot(xfit,polynomial(xfit),color='red',lw=3)

    slope = 0.395 +/- 0.018
intercept = 0.151 +/- 0.111

[<matplotlib.lines.Line2D at 0x108e247d0>]