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)  

# define the fit function (NOTE THE INDENT!)
def fit_function(R, V0):
    return  (V0/R)

# set the array of initial value(s)
initial_values = (5.0)

# call the minimization procedure (NOTE THE ARGUMENTS)
pars, covm = curve_fit(fit_function, R, I, initial_values, numpy.sqrt(dI**2+((V0/(R+initial_values)**2)*dR)**2))

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

# calculate the resulting chisquare and the number of dof
chisq = ((I - fit_function(R, pars))**2/(dI**2+((V0/(R+initial_values)**2)*dR)**2)).sum()
ndof = len(R) - len(pars)

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

# bellurie 
pylab.rc('font',size=16)
pylab.xlabel('R  [kohm]')
pylab.ylabel('I  [mA]')
pylab.xscale('log'); pylab.yscale('log')
pylab.title('Fit of a subset without internal resistance')

# data plot
pylab.errorbar(R,I,dI,dR,linestyle = '', color = 'black', marker = '.')

func_grid = numpy.logspace(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/fig1.pdf')

# show the plot on screen
pylab.show()