notebook.community

Edit and run



In [183]:

    
from scipy.optimize import curve_fit
import astropy.io.ascii as ascii
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
from scipy.constants import *



In [198]:

    
solmass=1.99e30



In [184]:

    
G









    Out[184]:





6.67408e-11



In [197]:

    
parsec









    Out[197]:





3.0856775813057292e+16



In [213]:

    
uconv=np.sqrt(G*1e-12/(parsec*solmass))



In [214]:

    
print uconv









    



3.29680980573e-35



In [215]:

    
4.65e-20/np.sqrt(solmass)









    Out[215]:





3.2962976032887622e-35



In [217]:

    
4.0/3.0*pi









    Out[217]:





4.1887902047863905



In [2]:

    
tab1=ascii.read("Table1.csv")
tab2=ascii.read("Table2.csv")
tab3=ascii.read("Table3.csv")



In [9]:

    
tab1









    Out[9]:




<Table length=21>

radius rvel Theoryval Ydiffval rvelerr
float64 float64 float64 float64 float64
0.215827338129 12.9 9.0 3.9 1.8
0.44964028777 17.1 10.8 6.3 1.8
0.647482014388 23.7 11.85 11.85 1.8
0.89928057554 26.4 12.0 14.4 1.8
1.18705035971 27.6 12.0 15.6 1.8
1.36690647482 27.0 11.55 15.45 1.8
1.61870503597 27.6 11.55 16.05 1.8
1.8345323741 31.2 11.55 19.65 1.8
2.05035971223 32.7 11.4 21.3 1.8
2.26618705036 33.3 11.55 21.75 1.8
2.55395683453 37.2 12.15 25.05 1.8
2.76978417266 39.0 12.6 26.4 1.8
3.23741007194 41.7 13.2 28.5 1.8
3.72302158273 43.2 13.5 29.7 1.8
4.20863309353 45.3 14.1 31.2 1.8
4.64028776978 47.1 14.4 32.7 1.8
5.10791366906 49.2 14.7 34.5 1.8
5.59352517986 50.4 14.7 35.7 1.8
6.04316546763 50.7 14.7 36.0 1.8
6.52877697842 51.0 14.4 36.6 1.95
7.01438848921 52.35 14.1 38.25 2.7



In [4]:

    
plt.figure()
plt.errorbar(tab1['radius'],tab1['rvel'],yerr=tab1['rvelerr'],fmt='o',capsize=3)
plt.plot(tab1['radius'],tab1['Theoryval'],'-')
plt.plot(tab1['radius'],tab1['Ydiffval'],'--')
plt.title('Gas Dominated Dwarf (NGC 3741)')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('gasdwarf.pdf')
plt.savefig('gasdwarf.jpeg')



In [49]:

    
def linexpy(x, M, a):
    return (M*np.exp(-x/10.0)/np.sqrt(x))*(1.0+a*(1+x/6.0)*np.exp(-x/6.0))



In [50]:

    
popt, pcov = curve_fit(linexpy, tab1['radius'], tab1['rvel'])



In [51]:

    
popt









    Out[51]:





array([ 868.58810348,   -0.98697009])



In [52]:

    
pcov









    Out[52]:





array([[  7.34537177e+02,  -2.77073525e-02],
       [ -2.77073525e-02,   3.24922269e-06]])



In [53]:

    
plt.figure()
plt.errorbar(tab1['radius'],tab1['rvel'],yerr=tab1['rvelerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 7.5, 100)  # define values to plot the function for
plt.plot(xfine, linexpy(xfine, popt[0], popt[1]), 'r-')
plt.title('Gas Dominated Dwarf best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('gasdwarf_fit_y.pdf')
plt.savefig('gasdwarf_fit_y.jpeg')









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in divide
  from ipykernel import kernelapp as app



In [11]:

    
tab2









    Out[11]:




<Table length=31>

Radius Rvelocity Theory Veldiff Velerr
float64 float64 float64 float64 float64
0.616438356164 77.4725274725 51.6483516484 25.8241758242 10.4395604396
1.43835616438 100.0 89.010989011 10.989010989 3.2967032967
2.19178082192 109.89010989 100.549450549 9.34065934066 3.2967032967
2.87671232877 118.131868132 96.1538461538 21.978021978 3.2967032967
3.69863013699 121.428571429 86.8131868132 34.6153846154 3.2967032967
4.45205479452 121.428571429 81.3186813187 40.1098901099 3.2967032967
5.20547945205 117.582417582 74.1758241758 43.4065934066 3.2967032967
5.61643835616 116.483516484 72.5274725275 43.956043956 3.2967032967
6.09589041096 115.934065934 70.3296703297 45.6043956044 3.2967032967
6.84931506849 115.934065934 66.4835164835 49.4505494505 3.2967032967
... ... ... ... ...
15.0684931507 117.582417582 51.6483516484 65.9340659341 3.2967032967
15.8904109589 115.934065934 50.5494505495 65.3846153846 3.2967032967
16.5753424658 114.835164835 50.0 64.8351648352 3.2967032967
17.397260274 114.835164835 48.9010989011 65.9340659341 3.2967032967
18.1506849315 114.835164835 47.8021978022 67.032967033 3.2967032967
19.5205479452 114.835164835 46.1538461538 68.6813186813 3.2967032967
20.9589041096 115.934065934 45.0549450549 70.8791208791 3.2967032967
21.8493150685 114.835164835 43.956043956 70.8791208791 4.3956043956
22.6712328767 118.131868132 43.4065934066 74.7252747253 4.3956043956
23.4246575342 116.483516484 42.3076923077 74.1758241758 10.4395604396



In [12]:

    
plt.figure()
plt.errorbar(tab2['Radius'],tab2['Rvelocity'],yerr=tab2['Velerr'],fmt='o',capsize=3)
plt.plot(tab2['Radius'],tab2['Theory'],'-')
plt.plot(tab2['Radius'],tab2['Veldiff'],'--')
plt.title('Disk-Dominated Spiral (NGC6503)')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('diskspiral.pdf')
plt.savefig('diskspiral.jpeg')



In [71]:

    
def linexpy(x, M, a):
    return (M*np.exp(-x/2.0)/np.sqrt(x))*(1.0+a*(1+x/6.0)*np.exp(-x/6.0))



In [72]:

    
popt2, pcov2 = curve_fit(linexpy, tab2['Radius'], tab2['Rvelocity'])



In [73]:

    
popt2
pcov2









    Out[73]:





array([[  1.60045358e+07,  -1.18068680e+01],
       [ -1.18068680e+01,   4.67971564e-05]])



In [74]:

    
plt.figure()
plt.errorbar(tab2['Radius'],tab2['Rvelocity'],yerr=tab2['Velerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 25, 100)  # define values to plot the function for
plt.plot(xfine, linexpy(xfine, popt2[0],popt2[1]), 'r-')
plt.title('Disk-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('diskspiral_fit_y.pdf')
plt.savefig('diskspiral_fit_y.jpeg')









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in divide
  from ipykernel import kernelapp as app



In [16]:

    
print popt2
print pcov2









    



[ 4.1327069]
[[ 0.04702947]]



In [17]:

    
tab3









    Out[17]:




<Table length=18>

radius rvel Rtheory veldiff velerr
float64 float64 float64 float64 float64
1.59482758621 255.833333333 276.666666667 -20.8333333333 14.1666666667
1.76724137931 241.666666667 246.666666667 -5.0 14.1666666667
2.80172413793 224.166666667 220.0 4.16666666667 6.66666666667
3.87931034483 220.833333333 199.166666667 21.6666666667 6.66666666667
5.04310344828 219.166666667 184.166666667 35.0 6.66666666667
6.12068965517 219.166666667 170.0 49.1666666667 6.66666666667
7.19827586207 218.333333333 158.333333333 60.0 6.66666666667
8.23275862069 216.666666667 150.0 66.6666666667 5.0
9.35344827586 213.333333333 142.5 70.8333333333 5.0
10.474137931 211.666666667 135.833333333 75.8333333333 5.0
11.5086206897 207.5 130.0 77.5 5.0
12.4568965517 205.833333333 125.0 80.8333333333 5.0
13.7068965517 206.666666667 119.166666667 87.5 5.0
14.7413793103 205.833333333 115.0 90.8333333333 5.0
15.775862069 205.0 111.666666667 93.3333333333 5.0
16.9396551724 205.0 108.333333333 96.6666666667 5.0
17.974137931 204.166666667 104.166666667 100.0 5.0
19.0517241379 205.0 100.833333333 104.166666667 5.0



In [18]:

    
plt.figure()
plt.errorbar(tab3['radius'],tab3['rvel'],yerr=tab3['velerr'],fmt='o',capsize=3)
plt.plot(tab3['radius'],tab3['Rtheory'],'-')
plt.plot(tab3['radius'],tab3['veldiff'],'--')
plt.title('Bulge-Dominated Spiral (NGC7814)')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('bulgespiral.pdf')
plt.savefig('bulgespiral.jpeg')



In [104]:

    
def linexpy(x, M, a):
    return (M*np.exp(-x/1.0)/np.sqrt(x))*(1.0+a*(1+x/2.0)*np.exp(-x/2.0))



In [105]:

    
popt3, pcov3 = curve_fit(linexpy, tab3['radius'], tab3['rvel'])



In [106]:

    
plt.figure()
plt.errorbar(tab3['radius'],tab3['rvel'],yerr=tab3['velerr'],fmt='o',capsize=3)
xfine = np.linspace(1., 20, 100)  # define values to plot the function for
plt.plot(xfine, linexpy(xfine, popt3[0],popt3[1]), 'r-')
plt.title('Bulge-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('bulgespiral_fit_y.pdf')
plt.savefig('bulgespiral_fit_y.jpeg')



In [21]:

    
print popt3
print pcov3









    



[ 6.15506884]
[[ 0.06964819]]



In [91]:

    
popt, pcov = curve_fit(line, tab1['radius'][2:], tab1['Ydiffval'][2:])



In [90]:

    
tab1['radius'][2:]









    Out[90]:




<Column name='radius' dtype='float64' length=19>

0.647482014388
0.89928057554
1.18705035971
1.36690647482
1.61870503597
1.8345323741
2.05035971223
2.26618705036
2.55395683453
2.76978417266
3.23741007194
3.72302158273
4.20863309353
4.64028776978
5.10791366906
5.59352517986
6.04316546763
6.52877697842
7.01438848921



In [93]:

    
plt.figure()
plt.errorbar(tab1['radius'],tab1['Ydiffval'],yerr=tab1['rvelerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 7.5, 100)  # define values to plot the function for
plt.plot(xfine, line(xfine, popt[0], popt[1]), 'r-')
plt.title('Gas Dominated Dwarf best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Diff between theory & observed (in km/s/Mpc)')
plt.savefig('gasdwarf_diff_fit.pdf')
plt.savefig('gasdwarf_diff_fit.jpeg')



In [95]:

    
popt2, pcov2 = curve_fit(line, tab2['Radius'][4:], tab2['Veldiff'][4:])
plt.figure()
plt.errorbar(tab2['Radius'],tab2['Veldiff'],yerr=tab2['Velerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 25, 100)  # define values to plot the function for
plt.plot(xfine, line(xfine, popt2[0], popt2[1]), 'r-')
plt.title('Disk-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Diff between theory & observed (in km/s/Mpc)')
plt.savefig('diskspiral_diff_fit.pdf')
plt.savefig('diskspiral_diff_fit.jpeg')



In [96]:

    
popt3, pcov3 = curve_fit(line, tab3['radius'][2:], tab3['veldiff'][2:])
plt.figure()
plt.errorbar(tab3['radius'],tab3['veldiff'],yerr=tab3['velerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 20, 100)  # define values to plot the function for
plt.plot(xfine, line(xfine, popt3[0], popt3[1]), 'r-')
plt.title('Bulge-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Diff between theory & observed (in km/s/Mpc)')
plt.savefig('bulgespiral_diff_fit.pdf')
plt.savefig('bulgespiral_diff_fit.jpeg')



In [205]:

    
def linexpy1(x, M, a, l):
    return np.sqrt((M*x**2)*(1.0+a*(1.0+x/l)*np.exp(-x/l)))



In [206]:

    
popt, pcov = curve_fit(linexpy1, tab1['radius'], tab1['rvel'], sigma=tab1['rvelerr'])
perr = np.sqrt(np.diag(pcov))









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in sqrt
  from ipykernel import kernelapp as app



In [207]:

    
print popt
print pcov
print perr









    



[ 65.15482987  10.72570235   0.90868142]
[[  5.41861061e+01  -2.00173220e+00  -5.62382493e-01]
 [ -2.00173220e+00   3.37406238e+00  -9.52051082e-02]
 [ -5.62382493e-01  -9.52051082e-02   1.13174223e-02]]
[ 7.36112125  1.8368621   0.10638337]



In [208]:

    
plt.figure()
plt.errorbar(tab1['radius'],tab1['rvel'],yerr=tab1['rvelerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 7.5, 100)  # define values to plot the function for
plt.plot(xfine, linexpy1(xfine, popt[0], popt[1],popt[2]), 'r-')
plt.title('Gas Dominated Dwarf best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('gasdwarf_fit_omega_expy1.pdf')
plt.savefig('gasdwarf_fit_omega_expy1.jpeg')



In [170]:

    
popt2, pcov2 = curve_fit(linexpy1, tab2['Radius'], tab2['Rvelocity'],sigma=tab2['Velerr'])
perr = np.sqrt(np.diag(pcov2))



In [171]:

    
print popt2
print pcov2
print perr









    



[   298.19379512    240.53110697  13989.21916121]
[[  2.18538137e+19  -1.77738406e+19  -7.99310377e+18]
 [ -1.77738406e+19   1.44555736e+19   6.50084018e+18]
 [ -7.99310372e+18   6.50084014e+18   2.30125830e+19]]
[  4.67480627e+09   3.80204860e+09   4.79714321e+09]



In [172]:

    
plt.figure()
plt.errorbar(tab2['Radius'],tab2['Rvelocity'],yerr=tab2['Velerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 25, 100)  # define values to plot the function for
plt.plot(xfine, linexpy1(xfine, popt2[0], popt2[1],popt2[2]), 'r-')
plt.title('Disk-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('diskspiral_fit_y1.pdf')
plt.savefig('diskspiral_fit_y1.jpeg')









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in divide
  from ipykernel import kernelapp as app



In [176]:

    
popt3, pcov3 = curve_fit(linexpy1, tab3['radius'], tab3['rvel'],sigma=tab3['velerr'])
perr = np.sqrt(np.diag(pcov3))









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in sqrt
  from ipykernel import kernelapp as app



In [177]:

    
plt.figure()
plt.errorbar(tab3['radius'],tab3['rvel'],yerr=tab3['velerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 20, 100)  # define values to plot the function for
plt.plot(xfine, linexpy1(xfine, popt3[0], popt3[1], popt3[2]), 'r-')
plt.title('Bulge-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('bulgespiral_fit_y1.pdf')
plt.savefig('bulgespiral_fit_y1.jpeg')









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in divide
  from ipykernel import kernelapp as app



In [178]:

    
print popt3
print pcov3
print perr









    



[   447.11146951    787.75655669  20877.25087467]
[[ -1.84404813e+21   3.25959608e+21  -5.58621281e+20]
 [  3.25959608e+21  -5.76176208e+21   9.87436120e+20]
 [ -5.58621281e+20   9.87436120e+20   0.00000000e+00]]
[ nan  nan   0.]



In [218]:

    
def linexpy1(x, M, a, l):
    return 4.0/3.0*pi*x*np.sqrt(M*(1.0+a*(1.0+x/l)*np.exp(-x/l)))



In [219]:

    
popt, pcov = curve_fit(linexpy1, tab1['radius'], tab1['rvel'])
perr = np.sqrt(np.diag(pcov))



In [220]:

    
print popt
print pcov
print perr









    



[  3.44300528  10.872621     0.96286125]
[[ 0.16150153 -0.12883633 -0.03299331]
 [-0.12883633  3.27119715 -0.08946629]
 [-0.03299331 -0.08946629  0.01249255]]
[ 0.40187253  1.80864511  0.1117701 ]



In [221]:

    
plt.figure()
plt.errorbar(tab1['radius'],tab1['rvel'],yerr=tab1['rvelerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 7.5, 100)  # define values to plot the function for
plt.plot(xfine, linexpy1(xfine, popt[0], popt[1],popt[2]), 'r-')
plt.title('Gas Dominated Dwarf best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('gasdwarf_fit_y1.pdf')
plt.savefig('gasdwarf_fit_y1.jpeg')



In [222]:

    
popt2, pcov2 = curve_fit(linexpy1, tab2['Radius'], tab2['Rvelocity'])
perr = np.sqrt(np.diag(pcov2))



In [223]:

    
print popt2
print pcov2
print perr









    



[  1.92936299  64.84514769   2.04013059]
[[  4.29240736e-02  -7.35351350e-01  -1.28818061e-02]
 [ -7.35351350e-01   1.37691172e+02  -1.00761148e+00]
 [ -1.28818061e-02  -1.00761148e+00   1.92920203e-02]]
[ 0.40187253  1.80864511  0.1117701 ]



In [225]:

    
plt.figure()
plt.errorbar(tab2['Radius'],tab2['Rvelocity'],yerr=tab2['Velerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 25, 100)  # define values to plot the function for
plt.plot(xfine, linexpy1(xfine, popt2[0], popt2[1], popt2[2]), 'r-')
plt.title('Disk-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('diskspiral_fit_1y.pdf')
plt.savefig('diskspiral_fit_1y.jpeg')



In [227]:

    
popt3, pcov3 = curve_fit(linexpy1, tab3['radius'], tab3['rvel'])
perr = np.sqrt(np.diag(pcov3))



In [228]:

    
print popt3
print pcov3
print perr









    



[   9.39833112  105.57892479    1.47807102]
[[  1.84167176e+00  -1.30802581e+01  -6.33046984e-02]
 [ -1.30802581e+01   7.09713521e+02  -2.12070946e+00]
 [ -6.33046984e-02  -2.12070946e+00   1.68663834e-02]]
[  1.35708208  26.64044897   0.12987064]



In [230]:

    
plt.figure()
plt.errorbar(tab3['radius'],tab3['rvel'],yerr=tab3['velerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 20, 100)  # define values to plot the function for
plt.plot(xfine, linexpy1(xfine, popt3[0], popt3[1],popt3[2]), 'r-')
plt.title('Bulge-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('bulgespiral_fit_1y.pdf')
plt.savefig('bulgespiral_fit_1y.jpeg')



In [ ]:



In [54]:

    
def linexp2(x, a, b):
    return a * np.exp(b/x)



In [55]:

    
popt, pcov = curve_fit(linexp2, tab1['radius'], tab1['Ydiffval'])



In [56]:

    
print popt
print pcov









    



[ 42.74594111  -1.25279963]
[[ 2.8838049  -0.16093615]
 [-0.16093615  0.0126969 ]]



In [57]:

    
plt.figure()
plt.errorbar(tab1['radius'],tab1['Ydiffval'],yerr=tab1['rvelerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 7.5, 100)  # define values to plot the function for
plt.plot(xfine, linexp2(xfine, popt[0], popt[1]), 'r-')
plt.title('Gas Dominated Dwarf best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('gasdwarf_fit_omega_exp2.pdf')
plt.savefig('gasdwarf_fit_omega_exp2.jpeg')









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in divide
  from ipykernel import kernelapp as app



In [58]:

    
popt2, pcov2 = curve_fit(linexp2, tab2['Radius'], tab2['Veldiff'])



In [59]:

    
print popt2
print pcov2









    



[ 82.94523877  -3.59263564]
[[ 7.08181719 -0.70943047]
 [-0.70943047  0.09845183]]



In [60]:

    
plt.figure()
plt.errorbar(tab2['Radius'],tab2['Veldiff'],yerr=tab2['Velerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 25, 100)  # define values to plot the function for
plt.plot(xfine, linexp2(xfine, popt2[0], popt2[1]), 'r-')
plt.title('Disk-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('diskspiral_fit_omega_exp2.pdf')
plt.savefig('diskspiral_fit_omega_exp2.jpeg')









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in divide
  from ipykernel import kernelapp as app



In [61]:

    
popt3, pcov3 = curve_fit(linexp2, tab3['radius'], tab3['veldiff'])



In [62]:

    
plt.figure()
plt.errorbar(tab3['radius'],tab3['veldiff'],yerr=tab3['velerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 20, 100)  # define values to plot the function for
plt.plot(xfine, linexp2(xfine, popt3[0], popt3[1]), 'r-')
plt.title('Bulge-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('bulgespiral_fit_omega_exp2.pdf')
plt.savefig('bulgespiral_fit_omega_exp2.jpeg')









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in divide
  from ipykernel import kernelapp as app



In [ ]:



In [22]:

    
def linexp3(x, M, a):
    return (G*M/np.sqrt(x))*(1.0+a*(1+x/4.0)*np.exp(-x/4.0))



In [23]:

    
popt, pcov = curve_fit(linexp3, tab1['radius'], tab1['rvel'])



In [24]:

    
print popt
print pcov









    



[ 287.43626659   -0.95528935]
[[  1.65489916e+02  -5.75197523e-02]
 [ -5.75197523e-02   5.21181591e-05]]



In [25]:

    
plt.figure()
plt.errorbar(tab1['radius'],tab1['rvel'],yerr=tab1['rvelerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 7.5, 100)  # define values to plot the function for
plt.plot(xfine, linexp3(xfine, popt[0], popt[1]), 'r-')
plt.title('Gas Dominated Dwarf best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('gasdwarf_fit_omega_exp3.pdf')
plt.savefig('gasdwarf_fit_omega_exp3.jpeg')









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in divide
  from ipykernel import kernelapp as app



In [26]:

    
popt2, pcov2 = curve_fit(linexp3, tab2['Radius'], tab2['Rvelocity'])



In [27]:

    
print popt2
print pcov2









    



[ 517.97809376   -0.83464614]
[[  1.45806674e+02  -9.94948656e-02]
 [ -9.94948656e-02   2.40080975e-04]]



In [28]:

    
plt.figure()
plt.errorbar(tab2['Radius'],tab2['Rvelocity'],yerr=tab2['Velerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 25, 100)  # define values to plot the function for
plt.plot(xfine, linexp3(xfine, popt2[0], popt2[1]), 'r-')
plt.title('Disk-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('diskspiral_fit_omega_exp3.pdf')
plt.savefig('diskspiral_fit_omega_exp3.jpeg')









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in divide
  from ipykernel import kernelapp as app



In [29]:

    
popt3, pcov3 = curve_fit(linexp3, tab3['radius'], tab3['rvel'])



In [30]:

    
plt.figure()
plt.errorbar(tab3['radius'],tab3['rvel'],yerr=tab3['velerr'],fmt='o',capsize=3)
xfine = np.linspace(0., 20, 100)  # define values to plot the function for
plt.plot(xfine, linexp3(xfine, popt3[0], popt3[1]), 'r-')
plt.title('Bulge-dominated Spiral best fit')
plt.xlabel('Radius(in kpc)')
plt.ylabel('Rotational velocity (in km/s/Mpc)')
plt.savefig('bulgespiral_fit_omega_exp3.pdf')
plt.savefig('bulgespiral_fit_omega_exp3.jpeg')









    



/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in divide
  from ipykernel import kernelapp as app



In [122]:

    
help(curve_fit)









    



Help on function curve_fit in module scipy.optimize.minpack:

curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=True, bounds=(-inf, inf), method=None, jac=None, **kwargs)
    Use non-linear least squares to fit a function, f, to data.
    
    Assumes ``ydata = f(xdata, *params) + eps``
    
    Parameters
    ----------
    f : callable
        The model function, f(x, ...).  It must take the independent
        variable as the first argument and the parameters to fit as
        separate remaining arguments.
    xdata : An M-length sequence or an (k,M)-shaped array
        for functions with k predictors.
        The independent variable where the data is measured.
    ydata : M-length sequence
        The dependent data --- nominally f(xdata, ...)
    p0 : None, scalar, or N-length sequence, optional
        Initial guess for the parameters.  If None, then the initial
        values will all be 1 (if the number of parameters for the function
        can be determined using introspection, otherwise a ValueError
        is raised).
    sigma : None or M-length sequence, optional
        If not None, the uncertainties in the ydata array. These are used as
        weights in the least-squares problem
        i.e. minimising ``np.sum( ((f(xdata, *popt) - ydata) / sigma)**2 )``
        If None, the uncertainties are assumed to be 1.
    absolute_sigma : bool, optional
        If False, `sigma` denotes relative weights of the data points.
        The returned covariance matrix `pcov` is based on *estimated*
        errors in the data, and is not affected by the overall
        magnitude of the values in `sigma`. Only the relative
        magnitudes of the `sigma` values matter.
    
        If True, `sigma` describes one standard deviation errors of
        the input data points. The estimated covariance in `pcov` is
        based on these values.
    check_finite : bool, optional
        If True, check that the input arrays do not contain nans of infs,
        and raise a ValueError if they do. Setting this parameter to
        False may silently produce nonsensical results if the input arrays
        do contain nans. Default is True.
    bounds : 2-tuple of array_like, optional
        Lower and upper bounds on independent variables. Defaults to no bounds.        
        Each element of the tuple must be either an array with the length equal
        to the number of parameters, or a scalar (in which case the bound is
        taken to be the same for all parameters.) Use ``np.inf`` with an
        appropriate sign to disable bounds on all or some parameters.
    
        .. versionadded:: 0.17
    method : {'lm', 'trf', 'dogbox'}, optional
        Method to use for optimization.  See `least_squares` for more details.
        Default is 'lm' for unconstrained problems and 'trf' if `bounds` are
        provided. The method 'lm' won't work when the number of observations
        is less than the number of variables, use 'trf' or 'dogbox' in this
        case.
    
        .. versionadded:: 0.17
    jac : callable, string or None, optional
        Function with signature ``jac(x, ...)`` which computes the Jacobian
        matrix of the model function with respect to parameters as a dense
        array_like structure. It will be scaled according to provided `sigma`.
        If None (default), the Jacobian will be estimated numerically.
        String keywords for 'trf' and 'dogbox' methods can be used to select
        a finite difference scheme, see `least_squares`.
    
        .. versionadded:: 0.18
    kwargs
        Keyword arguments passed to `leastsq` for ``method='lm'`` or
        `least_squares` otherwise.
    
    Returns
    -------
    popt : array
        Optimal values for the parameters so that the sum of the squared error
        of ``f(xdata, *popt) - ydata`` is minimized
    pcov : 2d array
        The estimated covariance of popt. The diagonals provide the variance
        of the parameter estimate. To compute one standard deviation errors
        on the parameters use ``perr = np.sqrt(np.diag(pcov))``.
    
        How the `sigma` parameter affects the estimated covariance
        depends on `absolute_sigma` argument, as described above.
    
        If the Jacobian matrix at the solution doesn't have a full rank, then
        'lm' method returns a matrix filled with ``np.inf``, on the other hand
        'trf'  and 'dogbox' methods use Moore-Penrose pseudoinverse to compute
        the covariance matrix.
    
    Raises
    ------
    ValueError
        if either `ydata` or `xdata` contain NaNs, or if incompatible options
        are used.
    
    RuntimeError
        if the least-squares minimization fails.
    
    OptimizeWarning
        if covariance of the parameters can not be estimated.
    
    See Also
    --------
    least_squares : Minimize the sum of squares of nonlinear functions.
    stats.linregress : Calculate a linear least squares regression for two sets
                       of measurements.
    
    Notes
    -----
    With ``method='lm'``, the algorithm uses the Levenberg-Marquardt algorithm
    through `leastsq`. Note that this algorithm can only deal with
    unconstrained problems.
    
    Box constraints can be handled by methods 'trf' and 'dogbox'. Refer to
    the docstring of `least_squares` for more information.
    
    Examples
    --------
    >>> import numpy as np
    >>> from scipy.optimize import curve_fit
    >>> def func(x, a, b, c):
    ...     return a * np.exp(-b * x) + c
    
    >>> xdata = np.linspace(0, 4, 50)
    >>> y = func(xdata, 2.5, 1.3, 0.5)
    >>> ydata = y + 0.2 * np.random.normal(size=len(xdata))
    
    >>> popt, pcov = curve_fit(func, xdata, ydata)
    
    Constrain the optimization to the region of ``0 < a < 3``, ``0 < b < 2``
    and ``0 < c < 1``:
    
    >>> popt, pcov = curve_fit(func, xdata, ydata, bounds=(0, [3., 2., 1.]))



In [ ]:

radius	rvel	Theoryval	Ydiffval	rvelerr
float64	float64	float64	float64	float64
0.215827338129	12.9	9.0	3.9	1.8
0.44964028777	17.1	10.8	6.3	1.8
0.647482014388	23.7	11.85	11.85	1.8
0.89928057554	26.4	12.0	14.4	1.8
1.18705035971	27.6	12.0	15.6	1.8
1.36690647482	27.0	11.55	15.45	1.8
1.61870503597	27.6	11.55	16.05	1.8
1.8345323741	31.2	11.55	19.65	1.8
2.05035971223	32.7	11.4	21.3	1.8
2.26618705036	33.3	11.55	21.75	1.8
2.55395683453	37.2	12.15	25.05	1.8
2.76978417266	39.0	12.6	26.4	1.8
3.23741007194	41.7	13.2	28.5	1.8
3.72302158273	43.2	13.5	29.7	1.8
4.20863309353	45.3	14.1	31.2	1.8
4.64028776978	47.1	14.4	32.7	1.8
5.10791366906	49.2	14.7	34.5	1.8
5.59352517986	50.4	14.7	35.7	1.8
6.04316546763	50.7	14.7	36.0	1.8
6.52877697842	51.0	14.4	36.6	1.95
7.01438848921	52.35	14.1	38.25	2.7