Figure showing DM velocity curve

rotation.py

"""
Rotation velocity of DM
=======================
"""

import matplotlib
import numpy as np
import matplotlib.pyplot as plt

from matplotlib.ticker import MaxNLocator
from scipy.interpolate import spline
from numpy.random import normal


matplotlib.rc('text', usetex=False)
plt.xkcd()
matplotlib.rc('font', **{'size': 25})


fig, ax = plt.subplots(figsize=(8, 8))

r = np.linspace(0, 50, 10)

def kepler(r, match=5., scale=1.):
    """
    """
    x = r / match
    
    if x >= 1:
        return scale * x**-0.5
    else:
        return scale * x**0.5
 
def DM(r, match=5., scale=1):
    """
    """
    x = r / match
    
    if x >= 1:
        return scale
    else:
        return scale * x**0.5
        
dense = np.linspace(0, 50, 1000)
K = spline(r, map(kepler, r), dense)
D = spline(r, map(DM, r), dense)

r_data = np.linspace(0, 50, 30) 
mu_data = spline(r, map(DM, r), r_data) 
data = [normal(mu, 0.05 * abs(mu)) for mu in mu_data]

plt.plot(dense, K, color="Crimson", lw=4, label="Keplerian")
plt.plot(dense, D, color="RoyalBlue", lw=4, label="DM halo")
plt.scatter(r_data, data, color="Olive", label="Data", zorder=100)

plt.xlim([0, 50])
plt.ylim([0, 1.2])
plt.setp(ax, yticklabels=[], xticklabels=[])

plt.xlabel("Distance from center")
plt.ylabel("Rotation speed")
ax.legend(loc='lower right')
        
plt.savefig("rotation.png")