import pylab
import numpy
from scipy.optimize import curve_fit

# data load
f,V,df,dV=pylab.loadtxt('D:/stash/dida13_14/lab2_13_14/bode/hp.txt',unpack=True)

# set the parameter and its uncertainty
V0 = 9.0
dV0 = 0.5

# create the array Aj and DeltaAj
A = V/V0
dA = (dV/V+dV0/V0)*A

# create the corresponding arrays measured in dB
Adb = 20*numpy.log10(A)
dAdb = (20/numpy.log(10))*(dA/A)

# duplicate f df Adb and dAdb for killing reasons
fkill = f
dfkill = df
Adbkill = Adb
dAdbkill = dAdb

# limit the range by killing data points if f > f_lim_high
f_lim_high = 1.5e3
index_kill=[]
for counter in range (0,len(f)):
    if f[counter] > f_lim_high:
        index_kill.append(counter)

fkill = numpy.delete(fkill,index_kill)
dfkill = numpy.delete(dfkill,index_kill)
Adbkill = numpy.delete(Adbkill,index_kill)    
dAdbkill = numpy.delete(dAdbkill,index_kill) 

# define the fit function (NOTE THE INDENT!)
def fit_function(f, fTdb):
    return  (-fTdb+20*numpy.log10(f))

# set the array of initial value(s) (IN THE FORM OF A LIST)
initial_values = [20*numpy.log10(5.e2)]

# call the minimization procedure (NOTE THE ARGUMENTS)
pars, covm = curve_fit(fit_function, fkill, Adbkill, initial_values, dAdbkill)

# calculate the resulting chisquare and the number of dof
chisq = (((Adbkill - fit_function(fkill, pars[0]))/dAdbkill)**2).sum()
ndof = len(fkill) - len(pars)

# go back to the physical parameter and its uncertainty
fT = 10**(pars/20)
dfT = numpy.sqrt(covm.diagonal())*(pars/20)*10**(pars/20-1)

# print the results on the console in an "expert" mode
print('Chisquare/ndof = %f/%d' % (chisq, ndof))
print('Cut off frequency: ', fT)
print('Uncertainty on cut off frequency: ', dfT)

# bellurie 
pylab.rc('font',size=16)
pylab.xlabel('$f$  [Hz]')
pylab.ylabel('$A$  [dB]')
pylab.xscale('log')
pylab.xlim(8,1.1e4)
pylab.ylim(-50,2.)
pylab.minorticks_on()
pylab.title('High-pass: Bode plot')

# data plot (NOTE THE ORDER OF ARGUMENTS)
pylab.errorbar(f,Adb,dAdb,df,linestyle = '', color = 'black', marker = '.')

func_grid = numpy.logspace(1, 4, 100)
pylab.plot(func_grid, fit_function(func_grid, pars[0]), color = 'black')

# plot the dotted line for zero
zeroline=func_grid*0
pylab.plot(func_grid,zeroline, '--')

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

# show the plot on screen
pylab.show()