#!/usr/bin/env python
# coding: utf-8

# ## Importing Necessary Packages for carrying out the computation
# 
# ### I have followed [this](https://mmas.github.io/least-squares-fitting-numpy-scipy) tutorial by Modesto Mas. My most of the python codes are inspired from this one.

# In[19]:


import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sc
import scipy.optimize as sp


# ### Plotting $\sin(x)$ and $\exp(x)$

# In[3]:


x0 = np.linspace(0,10, 100)
y01 = np.sin(x0)
y02 = np.exp(x0)
plt.figure(1,(8,8))
plt.plot(x0,y01,'ro')


# In[4]:


plt.figure(1,(8,8))
plt.plot(x0,y02,'g')


# In[5]:


plt.figure(1,(8,8))
plt.plot(x0,y01,x0,y02/20000.)


# ## Fitting a sample Curve
# 
# Let us create a set of sata which seems to be linear.
# 

# In[6]:


fun0 = np.poly1d([6,3]) #creates a polynomial function 6x+3
print(np.poly1d(fun0))
print(fun0(5)) #evaluates function at x=5
roots=fun0.r #finds the root
print(roots)


# In[7]:


x1 = np.linspace(0,20,50)
y1 = fun0(x1) + 10*np.random.normal(size=len(x1))
plt.figure(1,(8,8))
plt.plot(x1,y1,'r.')


# ## fitting can be made using different algorithms. Common method is least square fit

# In[8]:


#method 1 - polyfit from numpy
m1, c1 = np.polyfit(x1,y1,1)
m1,c1


# In[9]:


#method 2 - liear regression from scipy
m2,c2, r,p,er = sc.linregress(x1,y1)
m2,c2


# In[10]:


plt.plot(x1,y1,'r.')
x1_new = np.linspace(0,20,200)
y1_new = np.polyval([m1,c1],x1_new) #polyval creates new function using the fitting parameters we got
plt.plot(x1_new,y1_new)


# ## Fitting Polynomials
# 
# Using ```poly1d``` function, a polynomial function $5x^2-x+3$ is created

# In[11]:


fun1 = np.poly1d([5,-1,3]) #creates a polynomial function of 5x^2-x+3
print(np.poly1d(fun1))


# In[12]:


x2 = np.linspace(0,10,20)
y2 = fun1(x2)+25*np.random.normal(size=len(x2))
plt.plot(x2,y2,'.')


# In[104]:


m21,m22, m23 = np.polyfit(x2,y2,2)
x2_new = np.linspace(0,10,200)
y2_new = np.polyval([m21,m22,m23],x2_new)
plt.plot(x2,y2,'.')
plt.plot(x2_new,y2_new)
m21, m22, m23


# In[144]:


def fun2(x3,a,b):
    return a*np.sin(b*np.pi*x3)
p0 = [0.5,2]
x3 =  np.linspace(0,1,30)
y3 = fun2(x3,*p0) + 0.5*np.random.normal(size=len(x3))


# In[145]:


x3_new = np.linspace(0,1,200)
opt,conv = sp.curve_fit(fun2,x3,y3,[-10,4])


# In[146]:


opt


# In[147]:


plt.plot(x3,y3,'or')
plt.plot(x3_new, fun2(x3_new,*popt))


# The above one isn't probably a best fit. It can be corrected by using appropriate initial values.

# ## Fit and find Hubble parameter
# 
# ### A text file with data, distance(d) and velocities of galaxies (v) can be downloaded from [here](data_sdss.txt). 
# 
# We have to **unpack** it to read as separate columns and store in z_val and v_val. Now we plot it and fit using Linear regression method. According to Hubble's Law,
# 
# $v= Hd$
# 
# Thus slope of *v* versus *d* plot should give Hubble parameter.
# 

# In[119]:


d, v = np.loadtxt("data_sdss.txt",unpack=True) #Import data file "data_sdss.txt" 
#and save the columns in two variables


# In[153]:


plt.figure(1,(10,6))
plt.xlabel('d (Mpc)') #X label
plt.ylabel('v (km/s)') #Y label
plt.scatter(d,v)
slope, intercept, r_value, p_value, std_err = sc.linregress(d,v)
plt.plot([0, 2000], [intercept, slope * 2000 + intercept], '-k', lw=2) #here we plot x axis from 0 to 2000, 
#y axis as intercept + slope * x value, which gives mx+c, equation of line.
print("Actual value of Hubble Parameter = 69.8")
print('Obtained Value of Hubble Parameter = %1.2f'%slope)


# ## Thus we got the slope, Hubble parameter which is near to the currently estimated value

# # Resources
# 
# ## Reference : [Least Square Fit Tutorial by Modesto Mas](https://mmas.github.io/least-squares-fitting-numpy-scipy)
# 
# ## [SDSS modified text file with distance and velocity of galaxies data](data_sdss.txt)
# 
# ## [Download the whole code as python file](fitcode.py)
# 
# ## [Download the ipython Notebook for Jupyter](fitcode.ipyb)

# In[ ]:





# In[ ]: