import matplotlib.pyplot as plt
import numpy as np
import time

# cs is a class for computing the response of passive circuits with R, L, and C

# Starting and ending frequency and number of frequencies
# Making an instance creates a frequency array with settable properties
# This frequency array is used in the various methods without the user
# needing to explicitly include it.
class cs():
    def __init__(self, fs = 1., fe = 20000., fnum = 1000., fext = 'none'):
        if type(fext) == str:
            self.f = np.linspace(fs, fe, num = fnum)
        else:
            self.f = fext
            if self.f[0] == 0.:
                self.f[0] = self.f[1]

# Method to return impedance of a capacitor
    def fZC(self, C):
        return -1.j/(2.*np.pi*self.f*C)

# Method to return impedance of an inductor
    def fZL(self, L):
        return 1.j*2*np.pi*self.f*L

# Method to return impedance of a resistor (as complex array)
    def fZR(self, R):
       ZR = np.zeros(self.f.size, dtype = np.complex128)
       ZR.real = np.linspace(R, R, self.f.size)
       return ZR

# Method to parallel 2 or more impedances
    def par(self, Zf, *Zs):
        admit = 1./Zf
        for Z in Zs:
            admit += 1./Z
        return 1./admit

# Method to find the response of a voltage divider
    def vdiv(self, Zser, Zshu):
        return Zshu/(Zser + Zshu)

# Method to make impedance of parallel R, L, C
    def fRLC(self, R, L, C):
        return self.par(self.fZL(L), self.par(self.fZR(R)), self.fZC(C))

# Method to make impedance of a network to simulate speaker impedance 
# Rb, Lb and Cb are in parallel and are the base resonance
# Lb and Cb have resistors in series
# Lh0 and Lh1 are two inductors in series with resistors across them that
# make the impedance rise with frequency with approximate allowance
# for eddy current losses.
# LH0 and LH1 also have resistances in series with them
    def fZs(self, 
        Rb = 100.,        # Resistance across LB (Set large for no function) 
        Lb = 12.e-3,     # Big choke for bass resonance
        Cb = 380.e-6,      # Big C for bass resonance
        RLb = .5,          # Resistance in series with Lb
        Rcb = .05,         # Resistance in series with Cb 
        Rs = 5.8,          # Main dissipator
        Lh0 = .45e-3,      # Choke for high frequencies, number 0
        Rh0 = 6.5,         # Resistance across Lh0 
        Rs0 = .5,          # Resistance in series with Lh0
        Lh1 = .45e-3,      # Choke for high frequencies, number 1
        Rh1 = 50.,         # Resistance across Lh1
        Rs1 = .5           # Resistance in series with Lh1
        ):        
        ZLser = self.fZL(Lb) + RLb
        ZRb = self.fZR(Rb)
        ZCser =self.fZC(Cb) + Rcb
        Zb = self.par(ZRb, ZLser, ZCser) + Rs
        return Zb + self.par(self.fZL(Lh0) + Rs0, Rh0) + \
        self.par(self.fZL(Lh1) + Rs1, Rh1)

# Method to make a plot of the results of vdviv
    def plres(self, res, str):
        plt.semilogx(self.f, 20.*np.log10(np.absolute(res)))
        plt.grid()
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Response (0 = unity gain)')
        plt.title('{0:s}  {1:s}'.format(str, time.ctime()))

# Method to plot an impedance
    def plIm(o, Zs, str):
        pstr = 'Bl: real, Gr: imag  '
        plt.semilogx(o.f, Zs.real)
        plt.semilogx(o.f, Zs.imag)
        plt.grid()
        plt.grid(which = 'minor')
        plt.xlim(o.f[0], o.f[-1])
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Impedance')
        plt.title(pstr +'{0:s}  {1:s}'.format(str, time.ctime()))

# Method to plot the magnitude of an impedance
    def plMag(o, Zs, str):
        pstr = 'Bl: real, Gr: imag  '
        plt.loglog(o.f, np.absolute(Zs))
        yl = plt.ylim()
        plt.ylim(np.amin(np.absolute(Zs)) - 1., np.amax(np.absolute(Zs)) + 1.)
        plt.grid()
        plt.grid(which = 'minor')
        plt.xlim(o.f[0], o.f[-1])
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Msgnitude of Impedance')
        plt.title(pstr +'{0:s}  {1:s}'.format(str, time.ctime()))

# for computing impedance versus frequency from  pickup parameters,
# including the the eddy current terms a derived in "Mutual...  "
# Default values are provided for everything.  Actaul values can
# be supplied when making the object, or modified later.
class mimp():
    def __init__(self,
        Lc = 2.15,
        Rc = 7000.,
        k = .5,
        Rse = 1.e5, 
        C = 1.e-10, 
        fs = 1.,
        fe = 20000.,
        fnum = 1000.,
        fext = 'none'):
        self.Lc = Lc
        self.Rc = Rc
        self.k = k
        self.Rse = Rse
        self.C = C
        if type(fext) == str:
            self.f = np.linspace(fs, fe, num = fnum)
        else:
            self.f = fext
            if self.f[0] == 0.:
                self.f[0] = self.f[1]

#  This is the series combination of Lc and Rc.
    def mZLR(self):
        return 1j*2.*np.pi*self.f*self.Lc + self.Rc

# This is the eddy current part (from eq. 7 in Mutual....).
    def mZec(self):
        return (2.*np.pi*self.f*self.k*self.Lc)**2/ \
            (1j*2.*np.pi*self.f*self.Lc + self.Rse)

# This is the pickup impedance, no C.
    def mZnc(self):
        return self.mZLR() + self.mZec()

# This is the pickup impedance with C.
    def mZpu(self):
        return self.parC(self.mZLR() + self.mZec())

# This puts C in parallel.
    def parC(self, Z):
        self.Zc = -1.j/(2.*np.pi*self.f*self.C)
        return Z*self.Zc/(Z + self.Zc)

# Routine to compute equation 2, page 484, Radiotron... 4th edition
def bypass(f, Ck, Rk, Rprime, mu, rp):
    omega = 2*np.pi*f
    oTerm = omega*Ck*Rk
    muTerm = (mu + 1)*Rk/(Rprime + rp)
    return np.sqrt((1 + oTerm**2)/(1 + muTerm**2 + oTerm**2))

def phalf(f, r):
    half = (r[0] + r[-1])/2.
    halfloc = np.argmin(np.fabs(r - half))
    plt.semilogx(f[halfloc], r[halfloc], 'ko')

def pby(f,r):
    plt.semilogx(f, r)
    phalf(f, r)