%pylab inline
import SB_model
from pymc import MCMC, Matplot

M = MCMC(SB_model)

plt.plot(M.r.value, M.SB.value, c='gray', label='model')
plt.scatter(M.r.value, M.mags.value, c='r', label='observations')

plt.gca().invert_yaxis()
plt.xlabel('radius (")')
plt.ylabel('Surface Brightness (mag)')
plt.legend(loc='best')


# reasonable values for this model
# this will take a few seconds to run, so be patient!
M.sample(iter=5e4, burn=1000, thin=100)

# plot function takes the model (or a single parameter) as an argument:
Matplot.plot(M)

print 'R_effective (bulge) / R_effective (disk) =', \
       M.stats()['r_e_B']['mean'] / M.stats()['r_e_D']['mean']
print 'Effective surface brightness of the bulge: \n', \
       '    Best-fit value:', M.stats()['M_e_B']['mean'], \
       '\n    95% Confidence interval:', M.stats()['M_e_B']['quantiles'][2.5], \
        'to', M.stats()['M_e_B']['quantiles'][97.5]

for i in range(50):
    plt.plot(M.r.value, M.trace('SB')[i], c='gray', alpha=.25)

plt.scatter(M.r.value, M.mags.value, c='r') 
plt.gca().invert_yaxis()
plt.xlabel('radius (")')
plt.ylabel('Surface Brightness (mag)')

# our best-fit parameters
r_e_B = M.stats()['r_e_B']['mean']
r_e_D = M.stats()['r_e_D']['mean']
M_e_B = M.stats()['M_e_B']['mean']
M_e_D = M.stats()['M_e_D']['mean']
n     = M.stats()['n']['mean']
r     = M.r.value

from SB_model import sersic
# the sersic profile defined in SB_model.py expects brightnesses
#  in flux units, not magnitudes, so we have to convert M_e_B and
#  M_e_D before feeding them in.
plt.plot(r, M.stats()['SB']['mean'], c='k')
plt.plot(r, -2.5*np.log10(sersic( r, r_e_B, n, 10**(-.4*M_e_B) )), 'g:', label='bulge')
plt.plot(r, -2.5*np.log10(sersic( r, r_e_D, 1., 10**(-.4*M_e_D) )), 'g--', label='disk')
plt.scatter(M.r.value, M.mags.value, c='r', label='obs')

plt.gca().invert_yaxis()
plt.axis( [0, 6000, 26, 14] )
plt.ylabel('Surface Brightness (mag)')
plt.xlabel('Radius (")')
plt.legend(loc='best')