In [1]:

    

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 [2]:

    

solmass=1.99e30


    

In [3]:

    

G


    

    
Out[3]:





6.67408e-11

In [4]:

    

parsec


    

    
Out[4]:





3.0856775813057292e+16

In [5]:

    

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


    

In [6]:

    

print uconv


    

    




3.29680980573e-35

In [7]:

    

4.65e-20/np.sqrt(solmass)


    

    
Out[7]:





3.2962976032887622e-35

In [8]:

    

4.0/3.0*pi


    

    
Out[8]:





4.1887902047863905

In [9]:

    

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


    

In [ ]:

    




    

In [ ]:

    




    

In [29]:

    

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


    

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 [33]:

    

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


    

In [34]:

    

popt


    

    
Out[34]:





array([ 1.,  1.,  1.])

In [35]:

    

pcov


    

    
Out[35]:





array([[ inf,  inf,  inf],
       [ inf,  inf,  inf],
       [ inf,  inf,  inf]])

In [36]:

    

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, linexpy2(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_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 [37]:

    

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


    

In [38]:

    

popt2
pcov2


    

    
Out[38]:





array([[ inf,  inf,  inf],
       [ inf,  inf,  inf],
       [ inf,  inf,  inf]])

In [41]:

    

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, linexpy2(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_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 [42]:

    

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


    

In [44]:

    

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, linexpy2(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_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 [5]:

    

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


    

In [8]:

    

print popt3
print pcov3
print perr


    

    




[   365.94513368    638.77723908  14721.29122626]
[[  2.15500329e+19  -3.75079532e+19   9.53865951e+18]
 [ -3.75079532e+19   6.52828032e+19  -1.66020910e+19]
 [  9.53865951e+18  -1.66020910e+19   5.16418170e+19]]
[  4.64220129e+09   8.07977742e+09   7.18622411e+09]

In [11]:

    

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, linexpy2(xfine, popt3[0], popt3[1],popt3[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_y2.pdf')
plt.savefig('gasdwarf_fit_y2.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 [45]:

    

plt.hist([1,2,3,4])


    

    
Out[45]:





(array([ 1.,  0.,  0.,  1.,  0.,  0.,  1.,  0.,  0.,  1.]),
 array([ 1. ,  1.3,  1.6,  1.9,  2.2,  2.5,  2.8,  3.1,  3.4,  3.7,  4. ]),
 <a list of 10 Patch objects>)






    

In [ ]: