####################--
# File: eqcirc3.py
# Equivalent Circuit Parameter Estimator for Piezoelectric Structures
# Author: <NAME>
# Associate Professor of Mechanical Engineering
# 282 Toomey Hall
# 400 W. 13th St.
# Rolla, MO 65409-0050
# Email: <EMAIL>
# Original release: eqcirc1.py Version 0.1.0 3-29-2015
# Modified and renamed to eqcirc2.py to input impedance 4-17-2015:DSS
# Now also automatically determines the number of points to process.
#
# Modified for partial Python 3 compatibility, moved results printing
# into functions, refined plot formatting, and added
# flag to set the saved figure format. 4-26-2016: DSS
# Added time stamp on figures. 7-7-2016 <NAME>
# Added phase model. 7-8-2016: DSS
# Moved time stamp plot plot title and appended it to the plotfile, and
# added LaTeX formatting to the figure annotations 7-13-2016: DSS
# Added error bar functionality, and changed rmserr to calculate the
# more conservative sqrt(variance/(n-4)) estimate --
# reflecting the four degrees of freedom in the model.
# The error bars are ±rmserr centered on the model
# evaluated at the data point frequency. 1-2-2017: DSS
####################--

"""
    
    This program calculates the equivalent
    circuit parameters from frequency-impedance magnitude
    magnitude data stored in the standard
    HP4294A Impedance Analyser output data
    format.
    
    An image of the equivalent circuit may be
    found here:
    http://web.mst.edu/~stutts/piezoequivcircuit0.png
    
    The program first calculates the approximate
    equivalent circuit parameters for a single
    resonance-antiresonance frequency pair.
    It then uses the Levenberg-Marquardt (LM)
    nonlinear least squares algorithm to
    optimize the model in the least squares
    sense.  The LM algorithm is invoked
    via a call to leastsq(rez, z0, args=(yy, xx),
    full_output=1) from the scipy.optimize library.
    
    See: http://docs.scipy.org/doc/scipy-0.14.0/
    reference/generated/scipy.optimize.leastsq.html
    for more information.
    
    eqcirc3.py calculates the following outputs stdout:
    
    (1) fr (the resonance frequency)
    (2) fa (the anti-resonance frequency)
    (3) C0 (the parallel capacitance)
    (4) R1 (the motional resistance)
    (5) L1 (the motional inductance)
    (6) C1 (the motional capacitance)
    (7) Q (the series R1L1C1 resonance
    quality factor = 1/2zeta)
    (8) RMS Deviation
    
    A graph of the data and model is also produced.
    
    Example call: python eqcirc3.py inputdatafile.txt
    
    The graph may be saved in PNG format, and the text
    may be redirected from stdout to a file like so:
    
    python eqcirc3.py inputdatafile.txt > outdata.txt
    
    # This code is copyrighted by the author, but released under the MIT
    # license:
    
    Copyright (c) 2015 -- eqcirc3.py
    
    S&T and the University of Missouri Board of Curators
    license to you the right to use, modify, copy, and distribute this
    code subject to the MIT license:
    
    Permission is hereby granted, free of charge, to any person obtaining a copy
    of this software and associated documentation files (the "Software"), to deal
    in the Software without restriction, including without limitation the rights
    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
    copies of the Software, and to permit persons to whom the Software is
    furnished to do so, subject to the following conditions:
    
    The above copyright notice and this permission notice shall be included
    in all copies or substantial portions of the Software.
    
    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
    THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
    DEALINGS IN THE SOFTWARE.
    
    The author kindly requests that any publications benefiting from the use
    of this software include the following citation:
    
    @Misc{eqcirc3_2015,
    author =   {<NAME>.},
    title = {{eqcirc3.py}: {Equivalent Circuit Parameter Estimator
    for Piezoelectric Structures.}},
    howpublished = {\\url{https://github.com/MSTESG/EQCIRC3.git}},
    year = {2016}}
    
    """


#from pylab import *
import sys
from scipy.optimize import leastsq
import numpy as np
from numpy import array
from numpy import sqrt
import matplotlib.pyplot as plt
import time #to allow time stamp on output
# Test for Python version:
cur_version = sys.version_info
# Initialize some lists:
ydat = []
x = []
xx = []
yy = []
zdat = []
f = 0.0
# Define functions:
def y(f, z):  # Admittance model
    return 0.2e1 * np.pi * f * np.sqrt(0.4e1 * z[0] ** 2*z[3]**2*
    2 * np.pi ** 2 * f ** 2 + (-0.4e1 *z[0]*z[3]*z[2] * np.pi**2*f**2
    + z[0]+z[3])**2)*((-0.4e1 * z[3]*z[2]*np.pi ** 2*f**2+0.1e1)**2
    + 0.4e1*z[1]**2*z[3]**2*np.pi**2*f**2)**(-0.1e1/0.2e1)

# Real admittance:
Y_R = lambda f,cg,cg5,cg3,cg1: cg5 * cg1 ** 2 / (16 * cg ** 2 *
cg1 ** 2 * cg3 ** 2 * np.pi ** 4 * f ** 4 + 4 * cg ** 2 * cg1 ** 2 *
np.pi ** 2 * cg5 ** 2 * f ** 2 - 8 * cg ** 2 * cg1 * cg3 * np.pi ** 2
* f ** 2 - 8 * cg * cg1 ** 2 * cg3 * np.pi ** 2 * f ** 2 + cg ** 2 +
2 * cg * cg1 + cg1 ** 2)

# Imaginary admittance:
Y_I = lambda f,cg,cg5,cg3,cg1: -0.1e1 / np.pi / f * (16 * cg * cg1 ** 2
* cg3 ** 2 * np.pi ** 4 * f ** 4 + 4 * cg * cg1 ** 2 * np.pi ** 2
* cg5 ** 2 * f ** 2 - 8 * cg * cg1 * cg3 * np.pi ** 2 * f ** 2 -
4 * cg1 ** 2 * cg3 * np.pi ** 2 * f ** 2 + cg + cg1) / (16 * cg ** 2 *
cg1 ** 2 * cg3 ** 2 * np.pi ** 4 * f ** 4 + 4 * cg ** 2 * cg1 ** 2 *
np.pi ** 2 * cg5 ** 2 * f ** 2 - 8 * cg ** 2 * cg1 * cg3 *
np.pi ** 2 * f ** 2 - 8 * cg * cg1 ** 2 * cg3 * np.pi ** 2 * f ** 2
+ cg ** 2 + 2 * cg * cg1 + cg1 ** 2) / 2
phi = lambda f, cg,cg5,cg3,cg1:180*np.arctan2(Y_I(f,cg,cg5,cg3,cg1),
Y_R(f,cg,cg5,cg3,cg1))/np.pi

def C0_i(Ymin, Ymax, fr, fa):  # Parallel capacitance estimate
    return np.sqrt(0.2e1*(fa ** 2 - fr**2)*Ymin**2/np.pi**2/fa**4
    + 0.2e1*np.sqrt((fa ** 2 - f ** 2)**2/np.pi**4/fa**8*Ymin**4
    + 0.4e1*Ymin**2*Ymax**2/np.pi**4/fa**4))/0.4e1


def R1_i(Ymin, Ymax, fr, fa, C0):  # Motional resistance estimate
    return (-0.4e1*np.pi**2*fr**2*C0**2+Ymax**2)**(-0.1e1/0.2e1)


def L1_i(fr, fa, C0):  # Motional inductance estimate
    return 0.1e1 / np.pi ** 2 / (fa ** 2 - fr ** 2) / C0 / 0.4e1


def C1_i(fr, fa, C0):  # Motional capacitance estimate
    return (fa ** 2 / fr ** 2 - 1) * C0


def rez(z, ydat, f):  # Residual function
    return ydat - y(f, z)
#
# Comment out py3print() if you are using Python 2.x.x
#
def py3print():
    '''
        Author: <NAME>
        4-22-2016
        This function uses the print() function
        according to the Python 3.x.x requirements.
        '''
    print( "Ymax = ", Ymax, " at fr = ", fr, "\n")
    print( "Ymin = ", Ymin, " at fa = ", fa, "\n")
    print( "fr = ", fr, "\n")
    print( "fa = ", fa, "\n")
    # Initial estimates:
    print("C0i = ", C0i,"\n")
    print("R1i = ", R1i,"\n")
    print("L1i = ", L1i,"\n")
    print("C1i = ", C1i,"\n")
    print("Qi = ", Qi,"\n")
    # Optimal estimates:
    print( "C0 = ", C0, "\n")
    print( "R1 = ", R1, "\n")
    print( "L1 = ", L1, "\n")
    print( "C1 = ", C1, "\n")
    print( "Q = ", Q, "\n")
    print("RMS Diviation = ", rmserr,"\n")

# Set the desired resolution:
res = 300# Use a larger value for PNG
#plottype = ''# Defaults to PNG
plottype = 'EPS'
# Input data file on command line:
infile = sys.argv[1]
data = open(infile, "r")  # get array out of input file
numline = 0
# Count the number of lines in the data file:
for line in data:
    numline +=1

# Calculate the number of magnitude data points:

data.seek(0) # Reset file pointer to the beginning
nummagpts = (numline - 1 - 26)/2
linecount = 0
# read the 21st through total lines from the data file
# and fill x,y lists with floating point numbers:
if cur_version[0]==3:# This is necesary due to the change in the type
    for line in data:# returned by the map function in Python 3.x.x.
        if linecount > 20 and linecount < nummagpts+21:# relative to 2.x.x.
            freqs = list(map(float, (line[0:31]).split()))
            impedances = list(map(float, (line[0:31]).split()))
            x.append(freqs[0])
            zdat.append(impedances[1])
        linecount += 1
else:
    for line in data:
        if linecount > 20 and linecount < nummagpts+21:
            x.append(map(float, (line[0:31]).split())[0])
            zdat.append(map(float, (line[0:31]).split())[1])
        linecount += 1

data.close()# close data file
xx = array(x)
zin = array(zdat)

yy = 1/zin

# Locate Ymax, Ymin, and initial guesses for fr and fa:
zmin = min(zdat)
zminidx = zdat.index(zmin)  # index of min impedance (max admittance)
zmax = max(zdat)
zmaxidx = zdat.index(zmax)  # index of max impedance (min admittance)
Ymax = max(yy)
Ymin = min(yy)
fr = x[zminidx]  # initial guess for resonance frequency
fa = x[zmaxidx]  # initial guess for antiresonance frequency
imax = len(x)

# Estimate initial parameter values:

C0i = C0_i(Ymin, Ymax, fr, fa)
R1i = R1_i(Ymin, Ymax, fr, fa, C0_i(Ymin, Ymax, fr, fa))
L1i = L1_i(fr, fa, C0_i(Ymin, Ymax, fr, fa))
C1i = C1_i(fr, fa, C0_i(Ymin, Ymax, fr, fa))
Qi = 1/(R1i*np.sqrt(C1i/L1i))
# Create initial guess array:
z0 = [C0i, R1i, L1i, C1i]

# Find the best values:
output = leastsq(rez, z0, args=(yy, xx), full_output=1)

C0 = output[0][0]
R1 = output[0][1]
L1 = output[0][2]
C1 = output[0][3]
Q = 1 / (R1 * np.sqrt(C1 / L1))
fr = 1 / np.sqrt(L1 * C1) / 0.2e1 / np.pi
fa = np.sqrt((C0 + C1) / C0 / C1 / L1) / np.pi / 0.2e1

# Put the optimal values in a list:
coeffs = [C0, R1, L1, C1]

# Calculate RMS error: output is a tuple containing a dictionary
var = np.inner(output[2]['fvec'],output[2]['fvec'])
rmserr = sqrt(var/(len(output[2]['fvec'])-4))  # Standard error estimate
# Print the results to std out:
py3print()

# Calculate plot annotation positions:
delx = (fa-fr)/10.0
dely = (Ymax-Ymin)/20.0
noteymax = 0.65*Ymax
# Plot the model and the data:
plt.figure(figsize=(8,7),dpi=res)
plt.subplot(211)
plt.plot(xx, y(xx, coeffs), 'r-', label='model')
plt.plot(xx, yy, 'go', label='data')
#plt.errorbar(xx, yy, xerr=0, yerr=rmserr)# If accounting for x uncertainty
plt.errorbar(xx, y(xx, coeffs), yerr=rmserr, linestyle='None')
plt.annotate(r"$f_r$ = "+'{: 3.3e}'.format(fr),xy=(fa-delx,noteymax))
plt.annotate(r"$f_a$ = "+'{: 3.3e}'.format(fa),xy=(fa-delx,noteymax-1.4*dely))
plt.annotate(r"$C_0$ = "+'{: 3.3e}'.format(C0),xy=(fa-delx,noteymax-2.8*dely))
plt.annotate(r"$R_1$ = "+'{: 3.3e}'.format(R1),xy=(fa-delx,noteymax-4.4*dely))
plt.annotate(r"$L_1$ = "+'{: 3.3e}'.format(L1),xy=(fa-delx,noteymax-5.8*dely))
plt.annotate(r"$C_1$ = "+'{: 3.3e}'.format(C1),xy=(fa-delx,noteymax-7.2*dely))
plt.annotate('Q = '+'{: 3.3e}'.format(Q),xy=(fa-delx,noteymax-8.6*dely))
plt.annotate('RMS Dev. = '+'{: 3.2e}'.format(rmserr),xy=(fa-delx,noteymax-10*dely))

# Add date and time in plot title:
loctime = time.asctime(time.localtime(time.time()))
plt.suptitle('Data File ='+infile+':  '+loctime)
#plt.annotate('Date and time = '+localtime,xy=(fa,noteymax+dely))
print ("Date and Time =", loctime, "\n")


legend = plt.legend(loc='upper right', shadow=True, fontsize='large')
plt.xlabel(r"$f$ (Hz)")
plt.ylabel(r"$\mathscr{Y}$ (A/V)")
plt.grid(True)
# Put a nice background color on the legend:
legend.get_frame().set_facecolor('#00FFCC')
plt.subplot(212)
plt.plot(xx, phi(xx, C0, R1, L1, C1), 'r-', label='model')
plt.xlabel(r"$f$ (Hz)")
plt.grid(True)
plt.ylabel(r"$\phi$ (degrees)")
legend = plt.legend(loc='upper right', shadow=True, fontsize='large')
legend.get_frame().set_facecolor('#00FFCC')
if plottype=='PNG' or plottype=='':# Default to PNG
    # Save plot as PNG:
    plotname = infile.split('.')[0]+"trmodel"+loctime.replace(':','-')+'.PNG'
    plt.savefig(plotname,format='png', dpi=res)
else:# Save plot as EPS:
    plotname = infile.split('.')[0]+"trmodel"+loctime.replace(':','-') + '.eps'
    plt.savefig(plotname,format='eps', dpi=res)

plt.show()


####################--
# File: eqcirc3.py
# Equivalent Circuit Parameter Estimator for Piezoelectric Structures
# Author: <NAME>
# Associate Professor of Mechanical Engineering
# 282 Toomey Hall
# 400 W. 13th St.
# Rolla, MO 65409-0050
# Email: <EMAIL>
# Original release: eqcirc1.py Version 0.1.0 3-29-2015
# Modified and renamed to eqcirc2.py to input impedance 4-17-2015:DSS
# Now also automatically determines the number of points to process.
#
# Modified for partial Python 3 compatibility, moved results printing
# into functions, refined plot formatting, and added
# flag to set the saved figure format. 4-26-2016: DSS
# Added time stamp on figures. 7-7-2016 <NAME>
# Added phase model. 7-8-2016: DSS
# Moved time stamp plot plot title and appended it to the plotfile, and
# added LaTeX formatting to the figure annotations 7-13-2016: DSS
# Added error bar functionality, and changed rmserr to calculate the
# more conservative sqrt(variance/(n-4)) estimate --
# reflecting the four degrees of freedom in the model.
# The error bars are ±rmserr centered on the model
# evaluated at the data point frequency. 1-2-2017: DSS
####################--

"""
    
    This program calculates the equivalent
    circuit parameters from frequency-impedance magnitude
    magnitude data stored in the standard
    HP4294A Impedance Analyser output data
    format.
    
    An image of the equivalent circuit may be
    found here:
    http://web.mst.edu/~stutts/piezoequivcircuit0.png
    
    The program first calculates the approximate
    equivalent circuit parameters for a single
    resonance-antiresonance frequency pair.
    It then uses the Levenberg-Marquardt (LM)
    nonlinear least squares algorithm to
    optimize the model in the least squares
    sense.  The LM algorithm is invoked
    via a call to leastsq(rez, z0, args=(yy, xx),
    full_output=1) from the scipy.optimize library.
    
    See: http://docs.scipy.org/doc/scipy-0.14.0/
    reference/generated/scipy.optimize.leastsq.html
    for more information.
    
    eqcirc3.py calculates the following outputs stdout:
    
    (1) fr (the resonance frequency)
    (2) fa (the anti-resonance frequency)
    (3) C0 (the parallel capacitance)
    (4) R1 (the motional resistance)
    (5) L1 (the motional inductance)
    (6) C1 (the motional capacitance)
    (7) Q (the series R1L1C1 resonance
    quality factor = 1/2zeta)
    (8) RMS Deviation
    
    A graph of the data and model is also produced.
    
    Example call: python eqcirc3.py inputdatafile.txt
    
    The graph may be saved in PNG format, and the text
    may be redirected from stdout to a file like so:
    
    python eqcirc3.py inputdatafile.txt > outdata.txt
    
    # This code is copyrighted by the author, but released under the MIT
    # license:
    
    Copyright (c) 2015 -- eqcirc3.py
    
    S&T and the University of Missouri Board of Curators
    license to you the right to use, modify, copy, and distribute this
    code subject to the MIT license:
    
    Permission is hereby granted, free of charge, to any person obtaining a copy
    of this software and associated documentation files (the "Software"), to deal
    in the Software without restriction, including without limitation the rights
    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
    copies of the Software, and to permit persons to whom the Software is
    furnished to do so, subject to the following conditions:
    
    The above copyright notice and this permission notice shall be included
    in all copies or substantial portions of the Software.
    
    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
    THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
    DEALINGS IN THE SOFTWARE.
    
    The author kindly requests that any publications benefiting from the use
    of this software include the following citation:
    
    @Misc{eqcirc3_2015,
    author =   {<NAME>.},
    title = {{eqcirc3.py}: {Equivalent Circuit Parameter Estimator
    for Piezoelectric Structures.}},
    howpublished = {\\url{https://github.com/MSTESG/EQCIRC3.git}},
    year = {2016}}
    
    """


#from pylab import *
import sys
from scipy.optimize import leastsq
import numpy as np
from numpy import array
from numpy import sqrt
import matplotlib.pyplot as plt
import time #to allow time stamp on output
# Test for Python version:
cur_version = sys.version_info
# Initialize some lists:
ydat = []
x = []
xx = []
yy = []
zdat = []
f = 0.0
# Define functions:
def y(f, z):  # Admittance model
    return 0.2e1 * np.pi * f * np.sqrt(0.4e1 * z[0] ** 2*z[3]**2*
    2 * np.pi ** 2 * f ** 2 + (-0.4e1 *z[0]*z[3]*z[2] * np.pi**2*f**2
    + z[0]+z[3])**2)*((-0.4e1 * z[3]*z[2]*np.pi ** 2*f**2+0.1e1)**2
    + 0.4e1*z[1]**2*z[3]**2*np.pi**2*f**2)**(-0.1e1/0.2e1)

# Real admittance:
Y_R = lambda f,cg,cg5,cg3,cg1: cg5 * cg1 ** 2 / (16 * cg ** 2 *
cg1 ** 2 * cg3 ** 2 * np.pi ** 4 * f ** 4 + 4 * cg ** 2 * cg1 ** 2 *
np.pi ** 2 * cg5 ** 2 * f ** 2 - 8 * cg ** 2 * cg1 * cg3 * np.pi ** 2
* f ** 2 - 8 * cg * cg1 ** 2 * cg3 * np.pi ** 2 * f ** 2 + cg ** 2 +
2 * cg * cg1 + cg1 ** 2)

# Imaginary admittance:
Y_I = lambda f,cg,cg5,cg3,cg1: -0.1e1 / np.pi / f * (16 * cg * cg1 ** 2
* cg3 ** 2 * np.pi ** 4 * f ** 4 + 4 * cg * cg1 ** 2 * np.pi ** 2
* cg5 ** 2 * f ** 2 - 8 * cg * cg1 * cg3 * np.pi ** 2 * f ** 2 -
4 * cg1 ** 2 * cg3 * np.pi ** 2 * f ** 2 + cg + cg1) / (16 * cg ** 2 *
cg1 ** 2 * cg3 ** 2 * np.pi ** 4 * f ** 4 + 4 * cg ** 2 * cg1 ** 2 *
np.pi ** 2 * cg5 ** 2 * f ** 2 - 8 * cg ** 2 * cg1 * cg3 *
np.pi ** 2 * f ** 2 - 8 * cg * cg1 ** 2 * cg3 * np.pi ** 2 * f ** 2
+ cg ** 2 + 2 * cg * cg1 + cg1 ** 2) / 2
phi = lambda f, cg,cg5,cg3,cg1:180*np.arctan2(Y_I(f,cg,cg5,cg3,cg1),
Y_R(f,cg,cg5,cg3,cg1))/np.pi

def C0_i(Ymin, Ymax, fr, fa):  # Parallel capacitance estimate
    return np.sqrt(0.2e1*(fa ** 2 - fr**2)*Ymin**2/np.pi**2/fa**4
    + 0.2e1*np.sqrt((fa ** 2 - f ** 2)**2/np.pi**4/fa**8*Ymin**4
    + 0.4e1*Ymin**2*Ymax**2/np.pi**4/fa**4))/0.4e1


def R1_i(Ymin, Ymax, fr, fa, C0):  # Motional resistance estimate
    return (-0.4e1*np.pi**2*fr**2*C0**2+Ymax**2)**(-0.1e1/0.2e1)


def L1_i(fr, fa, C0):  # Motional inductance estimate
    return 0.1e1 / np.pi ** 2 / (fa ** 2 - fr ** 2) / C0 / 0.4e1


def C1_i(fr, fa, C0):  # Motional capacitance estimate
    return (fa ** 2 / fr ** 2 - 1) * C0


def rez(z, ydat, f):  # Residual function
    return ydat - y(f, z)
#
# Comment out py3print() if you are using Python 2.x.x
#
def py3print():
    '''
        Author: <NAME>
        4-22-2016
        This function uses the print() function
        according to the Python 3.x.x requirements.
        '''
    print( "Ymax = ", Ymax, " at fr = ", fr, "\n")
    print( "Ymin = ", Ymin, " at fa = ", fa, "\n")
    print( "fr = ", fr, "\n")
    print( "fa = ", fa, "\n")
    # Initial estimates:
    print("C0i = ", C0i,"\n")
    print("R1i = ", R1i,"\n")
    print("L1i = ", L1i,"\n")
    print("C1i = ", C1i,"\n")
    print("Qi = ", Qi,"\n")
    # Optimal estimates:
    print( "C0 = ", C0, "\n")
    print( "R1 = ", R1, "\n")
    print( "L1 = ", L1, "\n")
    print( "C1 = ", C1, "\n")
    print( "Q = ", Q, "\n")
    print("RMS Diviation = ", rmserr,"\n")

# Set the desired resolution:
res = 300# Use a larger value for PNG
#plottype = ''# Defaults to PNG
plottype = 'EPS'
# Input data file on command line:
infile = sys.argv[1]
data = open(infile, "r")  # get array out of input file
numline = 0
# Count the number of lines in the data file:
for line in data:
    numline +=1

# Calculate the number of magnitude data points:

data.seek(0) # Reset file pointer to the beginning
nummagpts = (numline - 1 - 26)/2
linecount = 0
# read the 21st through total lines from the data file
# and fill x,y lists with floating point numbers:
if cur_version[0]==3:# This is necesary due to the change in the type
    for line in data:# returned by the map function in Python 3.x.x.
        if linecount > 20 and linecount < nummagpts+21:# relative to 2.x.x.
            freqs = list(map(float, (line[0:31]).split()))
            impedances = list(map(float, (line[0:31]).split()))
            x.append(freqs[0])
            zdat.append(impedances[1])
        linecount += 1
else:
    for line in data:
        if linecount > 20 and linecount < nummagpts+21:
            x.append(map(float, (line[0:31]).split())[0])
            zdat.append(map(float, (line[0:31]).split())[1])
        linecount += 1

data.close()# close data file
xx = array(x)
zin = array(zdat)

yy = 1/zin

# Locate Ymax, Ymin, and initial guesses for fr and fa:
zmin = min(zdat)
zminidx = zdat.index(zmin)  # index of min impedance (max admittance)
zmax = max(zdat)
zmaxidx = zdat.index(zmax)  # index of max impedance (min admittance)
Ymax = max(yy)
Ymin = min(yy)
fr = x[zminidx]  # initial guess for resonance frequency
fa = x[zmaxidx]  # initial guess for antiresonance frequency
imax = len(x)

# Estimate initial parameter values:

C0i = C0_i(Ymin, Ymax, fr, fa)
R1i = R1_i(Ymin, Ymax, fr, fa, C0_i(Ymin, Ymax, fr, fa))
L1i = L1_i(fr, fa, C0_i(Ymin, Ymax, fr, fa))
C1i = C1_i(fr, fa, C0_i(Ymin, Ymax, fr, fa))
Qi = 1/(R1i*np.sqrt(C1i/L1i))
# Create initial guess array:
z0 = [C0i, R1i, L1i, C1i]

# Find the best values:
output = leastsq(rez, z0, args=(yy, xx), full_output=1)

C0 = output[0][0]
R1 = output[0][1]
L1 = output[0][2]
C1 = output[0][3]
Q = 1 / (R1 * np.sqrt(C1 / L1))
fr = 1 / np.sqrt(L1 * C1) / 0.2e1 / np.pi
fa = np.sqrt((C0 + C1) / C0 / C1 / L1) / np.pi / 0.2e1

# Put the optimal values in a list:
coeffs = [C0, R1, L1, C1]

# Calculate RMS error: output is a tuple containing a dictionary
var = np.inner(output[2]['fvec'],output[2]['fvec'])
rmserr = sqrt(var/(len(output[2]['fvec'])-4))  # Standard error estimate
# Print the results to std out:
py3print()

# Calculate plot annotation positions:
delx = (fa-fr)/10.0
dely = (Ymax-Ymin)/20.0
noteymax = 0.65*Ymax
# Plot the model and the data:
plt.figure(figsize=(8,7),dpi=res)
plt.subplot(211)
plt.plot(xx, y(xx, coeffs), 'r-', label='model')
plt.plot(xx, yy, 'go', label='data')
#plt.errorbar(xx, yy, xerr=0, yerr=rmserr)# If accounting for x uncertainty
plt.errorbar(xx, y(xx, coeffs), yerr=rmserr, linestyle='None')
plt.annotate(r"$f_r$ = "+'{: 3.3e}'.format(fr),xy=(fa-delx,noteymax))
plt.annotate(r"$f_a$ = "+'{: 3.3e}'.format(fa),xy=(fa-delx,noteymax-1.4*dely))
plt.annotate(r"$C_0$ = "+'{: 3.3e}'.format(C0),xy=(fa-delx,noteymax-2.8*dely))
plt.annotate(r"$R_1$ = "+'{: 3.3e}'.format(R1),xy=(fa-delx,noteymax-4.4*dely))
plt.annotate(r"$L_1$ = "+'{: 3.3e}'.format(L1),xy=(fa-delx,noteymax-5.8*dely))
plt.annotate(r"$C_1$ = "+'{: 3.3e}'.format(C1),xy=(fa-delx,noteymax-7.2*dely))
plt.annotate('Q = '+'{: 3.3e}'.format(Q),xy=(fa-delx,noteymax-8.6*dely))
plt.annotate('RMS Dev. = '+'{: 3.2e}'.format(rmserr),xy=(fa-delx,noteymax-10*dely))

# Add date and time in plot title:
loctime = time.asctime(time.localtime(time.time()))
plt.suptitle('Data File ='+infile+':  '+loctime)
#plt.annotate('Date and time = '+localtime,xy=(fa,noteymax+dely))
print ("Date and Time =", loctime, "\n")


legend = plt.legend(loc='upper right', shadow=True, fontsize='large')
plt.xlabel(r"$f$ (Hz)")
plt.ylabel(r"$\mathscr{Y}$ (A/V)")
plt.grid(True)
# Put a nice background color on the legend:
legend.get_frame().set_facecolor('#00FFCC')
plt.subplot(212)
plt.plot(xx, phi(xx, C0, R1, L1, C1), 'r-', label='model')
plt.xlabel(r"$f$ (Hz)")
plt.grid(True)
plt.ylabel(r"$\phi$ (degrees)")
legend = plt.legend(loc='upper right', shadow=True, fontsize='large')
legend.get_frame().set_facecolor('#00FFCC')
if plottype=='PNG' or plottype=='':# Default to PNG
    # Save plot as PNG:
    plotname = infile.split('.')[0]+"trmodel"+loctime.replace(':','-')+'.PNG'
    plt.savefig(plotname,format='png', dpi=res)
else:# Save plot as EPS:
    plotname = infile.split('.')[0]+"trmodel"+loctime.replace(':','-') + '.eps'
    plt.savefig(plotname,format='eps', dpi=res)

plt.show()


####################-- # File: eqcirc3.py # Equivalent Circuit Parameter Estimator for Piezoelectric Structures # Author: <NAME> # Associate Professor of Mechanical Engineering # 282 Toomey Hall # 400 W. 13th St. # Rolla, MO 65409-0050 # Email: <EMAIL> # Original release: eqcirc1.py Version 0.1.0 3-29-2015 # Modified and renamed to eqcirc2.py to input impedance 4-17-2015:DSS # Now also automatically determines the number of points to process. # # Modified for partial Python 3 compatibility, moved results printing # into functions, refined plot formatting, and added # flag to set the saved figure format. 4-26-2016: DSS # Added time stamp on figures. 7-7-2016 <NAME> # Added phase model. 7-8-2016: DSS # Moved time stamp plot plot title and appended it to the plotfile, and # added LaTeX formatting to the figure annotations 7-13-2016: DSS # Added error bar functionality, and changed rmserr to calculate the # more conservative sqrt(variance/(n-4)) estimate -- # reflecting the four degrees of freedom in the model. # The error bars are ±rmserr centered on the model # evaluated at the data point frequency. 1-2-2017: DSS ####################-- This program calculates the equivalent     circuit parameters from frequency-impedance magnitude     magnitude data stored in the standard     HP4294A Impedance Analyser output data     format.          An image of the equivalent circuit may be     found here:     http://web.mst.edu/~stutts/piezoequivcircuit0.png          The program first calculates the approximate     equivalent circuit parameters for a single     resonance-antiresonance frequency pair.     It then uses the Levenberg-Marquardt (LM)     nonlinear least squares algorithm to     optimize the model in the least squares     sense.  The LM algorithm is invoked     via a call to leastsq(rez, z0, args=(yy, xx),     full_output=1) from the scipy.optimize library.          See: http://docs.scipy.org/doc/scipy-0.14.0/     reference/generated/scipy.optimize.leastsq.html     for more information.          eqcirc3.py calculates the following outputs stdout:          (1) fr (the resonance frequency)     (2) fa (the anti-resonance frequency)     (3) C0 (the parallel capacitance)     (4) R1 (the motional resistance)     (5) L1 (the motional inductance)     (6) C1 (the motional capacitance)     (7) Q (the series R1L1C1 resonance     quality factor = 1/2zeta)     (8) RMS Deviation          A graph of the data and model is also produced.          Example call: python eqcirc3.py inputdatafile.txt          The graph may be saved in PNG format, and the text     may be redirected from stdout to a file like so:          python eqcirc3.py inputdatafile.txt > outdata.txt          # This code is copyrighted by the author, but released under the MIT     # license:          Copyright (c) 2015 -- eqcirc3.py          S&T and the University of Missouri Board of Curators     license to you the right to use, modify, copy, and distribute this     code subject to the MIT license:          Permission is hereby granted, free of charge, to any person obtaining a copy     of this software and associated documentation files (the "Software"), to deal     in the Software without restriction, including without limitation the rights     to use, copy, modify, merge, publish, distribute, sublicense, and/or sell     copies of the Software, and to permit persons to whom the Software is     furnished to do so, subject to the following conditions:          The above copyright notice and this permission notice shall be included     in all copies or substantial portions of the Software.          THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,     FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL     THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER     LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING     FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER     DEALINGS IN THE SOFTWARE.          The author kindly requests that any publications benefiting from the use     of this software include the following citation:          @Misc{eqcirc3_2015,     author =   {<NAME>.},     title = {{eqcirc3.py}: {Equivalent Circuit Parameter Estimator     for Piezoelectric Structures.}},     howpublished = {\\url{https://github.com/MSTESG/EQCIRC3.git}},     year = {2016}} #from pylab import * #to allow time stamp on output # Test for Python version: # Initialize some lists: # Define functions: # Admittance model # Real admittance: # Imaginary admittance: # Parallel capacitance estimate # Motional resistance estimate # Motional inductance estimate # Motional capacitance estimate # Residual function # # Comment out py3print() if you are using Python 2.x.x # Author: <NAME>         4-22-2016         This function uses the print() function         according to the Python 3.x.x requirements. # Initial estimates: # Optimal estimates: # Set the desired resolution: # Use a larger value for PNG #plottype = ''# Defaults to PNG # Input data file on command line: # get array out of input file # Count the number of lines in the data file: # Calculate the number of magnitude data points: # Reset file pointer to the beginning # read the 21st through total lines from the data file # and fill x,y lists with floating point numbers: # This is necesary due to the change in the type # returned by the map function in Python 3.x.x. # relative to 2.x.x. # close data file # Locate Ymax, Ymin, and initial guesses for fr and fa: # index of min impedance (max admittance) # index of max impedance (min admittance) # initial guess for resonance frequency # initial guess for antiresonance frequency # Estimate initial parameter values: # Create initial guess array: # Find the best values: # Put the optimal values in a list: # Calculate RMS error: output is a tuple containing a dictionary # Standard error estimate # Print the results to std out: # Calculate plot annotation positions: # Plot the model and the data: #plt.errorbar(xx, yy, xerr=0, yerr=rmserr)# If accounting for x uncertainty # Add date and time in plot title: #plt.annotate('Date and time = '+localtime,xy=(fa,noteymax+dely)) # Put a nice background color on the legend: # Default to PNG # Save plot as PNG: # Save plot as EPS: