import numpy as np
from numpy import sqrt, pi, exp, linspace, loadtxt, arctan
from lmfit import  Model

import matplotlib.pyplot as plt

#import data
data = loadtxt('lab3_data1.txt')
frequency = data[:, 0]
Vin = data[:,1]
Vout = data[:,4]
phase = data[:, 5]
VoutVin = data[:,6]

 #plotting Vout vs. frequency (High Pass) experimental data
 plt.figure(1)
 plt.plot(frequency*(2*np.pi) ,VoutVin, 'ro')
 plt.xlabel('Angular Frequency')
 plt.ylabel('Vout/Vin')
 plt.show()

 #plotting phase vs. frequency (High Pass) experimental data
 plt.figure(2)
 plt.plot(frequency*(2*np.pi),phase, 'rx')
 plt.xlabel('Angular Frequency (Hz)')
 plt.ylabel('Phase (rad)')
 plt.show()

#expected data
F = np.arange(50,20000,50)
RC = 470*0.00000022
P = np.arctan(1/(RC*F*2*np.pi))

plt.figure(3)
plt.plot(F,P,'bx',frequency, phase, 'rx')
plt.xlabel('Frequency')
plt.ylabel('Phase Expected')
plt.show()

import numpy as np
from scipy.optimize import curve_fit
from matplotlib import pylab



 #Now fitting measured data for amplitude ratio to model, allowing RC
  to vary

 ratio1 = (F*RC)/((1+(F*RC)**2)**.5)

 def ratio(frequency, RC):
      return ((frequency*RC)/(1+(frequency*RC)**2)**.5)
 model = Model(ratio)


result =model.fit(ratio, frequency=frequency)

plt.figure(4)
plt.plot(frequency,ratio,'g.')
plt.xlabel('angular frequency')
plt.ylabel('voltage ratio')
plt.show()
plt.plot(frequency, result.init_fit, 'r--', linewidth=1)
plt.show()
plt.plot(frequency, result.best_fit, 'b-', linewidth=2)
plt.show()
plt.legend()
plt.show()

print result.fit_report()
fig=result.plot()
fig.show()
#---------------------------------------------------------------
#Now fitting measured data for phase to model, allowing RC to vary
ratio1 = (F*RC)/((1+(F*RC)**2)**.5)

def phasefit(frequency, RC):
      return (np.arctan(1/(frequency*RC)))
model = Model(phasefit)

result =model.fit(phase, frequency=frequency)

plt.figure(5)
plt.plot(frequency,phase,'g.')
plt.xlabel('angular frequency')
plt.ylabel('phase')
plt.show()
plt.plot(frequency, result.init_fit, 'r--', linewidth=1)
plt.show()
plt.plot(frequency, result.best_fit, 'b-', linewidth=2)
plt.show()
plt.legend()
plt.show()

print result.fit_report()
fig=result.plot()