import pylab
import numpy
from scipy.optimize import curve_fit

# data load
R,I,dR,dI=pylab.loadtxt('D:/stash/dida14_15/lab2_14_15/fitte_ohm/data.txt',unpack=True)

# limit the range by killing data points if R value is less than R_lim_low
R_lim_low = 0.3
index_kill=[]
for counter in range (0,len(R)):
    if R[counter] < R_lim_low:
        index_kill.append(counter)

R = numpy.delete(R,index_kill)
I = numpy.delete(I,index_kill)
dR = numpy.delete(dR,index_kill)    
dI = numpy.delete(dI,index_kill) 

# calculate the log10 of the data and the associated uncertainties
x = numpy.log10(R)
y = numpy.log10(I)
dy = dI/(I*numpy.log(10))
dx = dR/(R*numpy.log(10))
 

# define the fit function (NOTE THE INDENT!)
def fit_function(x, a):
    return  (a - x)

# set the array of initial value(s)
initial_values = (numpy.log10(5.))

# call the minimization procedure (NOTE THE ARGUMENTS)
pars, covm = curve_fit(fit_function, x, y, initial_values, numpy.sqrt(dy**2+dx**2))

# extract the estimated standard deviation of the parameter(s) 
delta_a = numpy.sqrt(covm.diagonal())

# calculate the resulting chisquare and the number of dof
chisq = ((y - fit_function(x, pars))**2/(dy**2+dx**2)).sum()
ndof = len(x) - len(pars)

# calculate V_0 and delta V_0 from a and delta a
V0 = 10**pars
deltaV0 = delta_a * (10**pars) *pars/10

# print the results on the console (with bellurie)
print('a = %f +- %f' % (pars, delta_a))
print('V_0 = %f +- %f' % (V0, deltaV0))
print('Chisquare/ndof = %f/%d' % (chisq, ndof))

# bellurie 
pylab.rc('font',size=16)
pylab.xlabel('log(R/R0) , with R0 = 1 kohm')
pylab.ylabel('log(I/I0) , with I0 = 1 mA')

pylab.title('Numerical fit of the linearized subset')

# data plot
pylab.errorbar(x,y,dy,dx,linestyle = '', color = 'black', marker = '.')

func_grid = numpy.linspace(numpy.log10(R_lim_low), 3, 100)
pylab.plot(func_grid, fit_function(func_grid, pars), color = 'black')

# save the plot as a pdf file somewhere (in my own directories!)
pylab.savefig('D:/stash/dida14_15/lab2_14_15/lin_vs_num/fig2.pdf')

# show the plot on screen
pylab.show()