#!/usr/bin/python3

from math import *
import numpy as np
from scipy.fft import fft, ifft, fftfreq, next_fast_len
from scipy.optimize import brentq as findRoot
import subprocess
from time import sleep

#########################################################################################################
class gnuPlotter:
        def __init__(self):
                self.gpp={}

        def open_gpp(self, gppnamestr=None):
                if (gppnamestr is None) or (gppnamestr not in self.gpp.keys()):
                        gppnew=subprocess.Popen(['gnuplot', '-p'], shell=False, universal_newlines=True, stdin=subprocess.PIPE)
                        #gppnew=subprocess.Popen(['gnuplot', '-p'], shell=False, universal_newlines=True, stdin=subprocess.PIPE, stderr=subprocess.PIPE)
                        gppnew.stdin.write("set term qt size 1900,950; \n")
                        #gppnew.stdin.write("set term x11 size 1900,950; \n")
                if gppnamestr is None:
                        return gppnew
                if gppnamestr not in self.gpp.keys():
                        self.gpp[gppnamestr]=gppnew
                return self.gpp[gppnamestr]

        #def gpp_Raw(self, ldat, rawstr, gppobj=None):
        def gpp_Raw(self, ldat, rawstr, gppstr=None):
                gppCurr=self.open_gpp(gppstr)
                #write ldat to datablock
                gppCurr.stdin.write("$idat << EOD \n")
                for irow in ldat:
                        gppCurr.stdin.write(("{} "*len(irow)+" \n").format(*irow))
                gppCurr.stdin.write("EOD \n")
                #send raw command; if datafile/block is passed to script, use '$idat'
                gppCurr.stdin.write(rawstr)

        #def gpp_Plot(self, ldat, strdat, lstrlbl, gppobj=None, lineWidth=2, toFile=''):
        def gpp_Plot(self, ldat, strdat, lstrlbl, gppstr=None, lineWidth=2, toFile='', config=''):
                gppCurr=self.open_gpp(gppstr)
                gppCurr.stdin.write(config)
                lcol=[r.split("[")[0].strip() for r in strdat.split(",")]
                def get_plotStr(strlbl):
                        llbl=[l.strip().split(":") for l in strlbl.split(",")]
                        nlbl=len(llbl)
                        lidx=[ [ [lcol.index(s) for s in lcol if s.endswith(lbl[ii])] for lbl in llbl] for ii in [0,1] ]
                        lidxpairs=[(lidx[0][i][0],lidx[1][i][0]) for i in range(0,nlbl)]
                        listr=["{}:{}".format(i[0]+1,i[1]+1) for i in lidxpairs]
                        plotStr="plot "
                        #litemStr=["$idat u {} w l lw 2 title '{}'".format(listr[i], llbl[i][1]) for i in range(0,len(listr))]
                        litemStr=["$idat u {} w l lw {} title '{}'".format(listr[i], lineWidth, llbl[i][1]) for i in range(0,len(listr))]
                        return plotStr+", ".join(litemStr)+" \n"
                lplotStr=[get_plotStr(lbl) for lbl in lstrlbl]

                if toFile.endswith(".pdf"):
                        gppCurr.stdin.write("set term pdfcairo enhanced font ',6' \n")
                        #print("set output '{}' \n".format(toFile))
                        gppCurr.stdin.write("set output '{}' \n".format(toFile))

                #write ldat to datablock
                gppCurr.stdin.write("$idat << EOD \n")
                for irow in ldat:
                        gppCurr.stdin.write(("{} "*len(irow)+" \n").format(*irow))
                gppCurr.stdin.write("EOD \n")

                #plot or multiplot
                if len(lstrlbl)==1:
                        gppCurr.stdin.write(lplotStr[0])
                elif len(lstrlbl)==2:
                        gppCurr.stdin.write("set multiplot \n")
                        gppCurr.stdin.write("set origin 0,0 \n  set size 1,0.5 \n")
                        gppCurr.stdin.write(lplotStr[0])
                        gppCurr.stdin.write("set origin 0,0.5 \n  set size 1,0.5 \n")
                        gppCurr.stdin.write(lplotStr[1])
                        gppCurr.stdin.write("unset multiplot \n")
                        #gppCurr.stdin.write(lplotStr[1]+" pause -1 \n")
                elif len(lstrlbl)==3:
                        gppCurr.stdin.write("set multiplot \n")
                        gppCurr.stdin.write("set origin 0,0 \n  set size 1,0.35 \n")
                        gppCurr.stdin.write(lplotStr[0])
                        gppCurr.stdin.write("set origin 0,0.333 \n  set size 1,0.35 \n")
                        gppCurr.stdin.write(lplotStr[1])
                        gppCurr.stdin.write("set origin 0,0.666 \n  set size 1,0.35 \n")
                        gppCurr.stdin.write(lplotStr[2])
                        gppCurr.stdin.write("unset multiplot \n")
                        #gppCurr.stdin.write(lplotStr[1]+" pause -1 \n")
                #if toFile is not None:
                if toFile != '':
                        gppCurr.stdin.write("set output \n")
                return gppCurr


#########################################################################################################
class fftIntegrator:
        def __init__(self, nptspp, fperiod, npad, nprev, ztype2fn, ztype2fn_dc, resultfn, resultbounds):
                self.nPtsPP=nptspp
                self.fPeriod=fperiod
                self.nPtsPS=round(self.nPtsPP*self.fPeriod)
                self.tauLen=1.0
                self.tauScale=1.0/self.fPeriod
                self.nPad=npad
                self.nPrev=nprev
                self.nPerTot=self.nPad+self.nPrev+self.nPad
                self.nPtsTot=self.nPtsPP*self.nPerTot
                self.zType2Fn=ztype2fn
                self.zType2Fn_DC=ztype2fn_dc
                self.resultFn=resultfn
                self.resultBounds=resultbounds
                self.tArr=None
                self.freqArr=None
                self.prevprevArr=None
                self.prevArr=None
                self.padArr=None
                self.gpp={}

        def fftIntegrate(self, arrcurr, flagSavePrev, plotit=False):
                if not "fftInt" in self.gpp.keys():
                        self.gpp["fftInt"]=subprocess.Popen(['gnuplot', '-p'], shell=False, universal_newlines=True, stdin=subprocess.PIPE, stderr=subprocess.PIPE)
                        self.gpp["fftInt"].stdin.write("set term x11 size 1900,950; \n")
                if self.freqArr is None:
                        self.freqArr=fftfreq(self.nPtsTot, d=self.tauScale/(self.nPtsPP))
                if self.tArr is None:
                        tpad_pre = np.linspace(-(self.nPrev+self.nPad)*self.tauLen, -self.nPrev*self.tauLen, round(self.tauScale*self.nPad*self.nPtsPS), endpoint=False)
                        tprev = np.linspace(-self.nPrev*self.tauLen, 0.0, round(self.tauScale*self.nPrev*self.nPtsPS), endpoint=False)
                        tpad_post = np.linspace(0.0, self.nPad*self.tauLen, round(self.tauScale*self.nPad*self.nPtsPS), endpoint=False)     #first period starting t=0 overwritten by current period to integrate
                        self.tArr = np.hstack((tpad_pre, tprev, tpad_post))
                if self.prevArr is None:
                        self.prevArr=np.zeros(round(self.tauScale*self.nPrev*self.nPtsPS))
                        self.padArr=np.zeros(round(self.tauScale*self.nPad*self.nPtsPS))
                idx0=-round(self.tauScale*self.nPad*self.nPtsPS)
                idxf=-round(self.tauScale*(self.nPad-1)*self.nPtsPS)
                idxshift_beg=-round(self.tauScale*(self.nPad+self.nPrev-1)*self.nPtsPS)                 #one period past start of previous history
                idxshift_end=-round(self.tauScale*(self.nPad-1)*self.nPtsPS)                            #to end of current solution, which=idxf
                if self.tArr[idx0]!=0.0 or self.tArr[idxf]!=1.0:
                        print("fftIntegrate: ERROR: check tcurr boundaries, should be 0.0, 1.0 by definition:", idx0, self.tArr[idx0], idxf, self.tArr[idxf])
                        quit()
                ft=np.hstack((self.padArr, self.prevArr, arrcurr, self.padArr[idxf:]))
                result = self.fftIntegrate_Type2(ft, idx0)
                rr=result[idx0:idxf]
                vaoclip=np.less(self.resultFn()+result.real, self.resultBounds[1])*np.greater(self.resultFn()+result.real, self.resultBounds[0])
                #ft= ft*vaoclip  #not the most accurate method, but seems best
                ft= ft-0.9*ft*np.logical_not(vaoclip)
                result=np.clip(result.real, a_min=self.resultBounds[0]-self.resultFn(), a_max=self.resultBounds[1]-self.resultFn())

                if plotit>1:
                        self.gpp["fftInt"].stdin.write("$idat << EOD \n")
                        for i in range(0, len(ft)):
                                ll=[self.tArr[i], ft[i], result[i].real, int(vaoclip[i])]
                                self.gpp["fftInt"].stdin.write(("{} "*len(ll)+" \n").format(*ll))
                        self.gpp["fftInt"].stdin.write("EOD \n")
                        self.gpp["fftInt"].stdin.write("set multiplot \n")
                        self.gpp["fftInt"].stdin.write("set origin 0,0 \n  set size 1,0.5 \n")
                        #self.gpp["fftInt"].stdin.write("plot $idat u 1:3 w l lw 2 title 'result' \n")
                        self.gpp["fftInt"].stdin.write("plot $idat u 1:3 w l lw 2 title 'result', $idat u 1:4 w l lw 2 title 'vaoclip' \n")
                        self.gpp["fftInt"].stdin.write("set origin 0,0.5 \n  set size 1,0.5 \n")
                        self.gpp["fftInt"].stdin.write("plot $idat u 1:2 w l lw 2 title 'ft' \n")
                        self.gpp["fftInt"].stdin.write("unset multiplot \n")
        
                if flagSavePrev:
                        self.prevprevArr=self.prevArr
                        self.prevArr=ft[idxshift_beg:idxshift_end]

                return rr

        def fftIntegrate_Type2(self, ft, idx0):
                fw = fft(ft, norm="forward")
                #print("freq[1]={}  freq[{}]={}  self.tauScale={}    t0=t[{}]={}  sum(fw)={:1.6e}".format(self.freqArr[1], (self.nPtsTot/2)-1, self.freqArr[int(self.nPtsTot/2)-1:int(self.nPtsTot/2)+1], self.tauScale, idx0, self.tArr[idx0], sum(fw)))
                fwlp = np.hstack((0, self.zType2Fn(fw[1:], self.freqArr[1:])))
                ftlp = ifft(fwlp, norm="forward")
                #print("DC components: fw[0]={:1.6e},  ftlp[0]={:1.6e},  self.tauScale={:1.6e}, self.zType2Fn_DC().real={:1.6e}".format(fw[0], ftlp[0], self.tauScale, self.zType2Fn_DC().real))
                ftlp += (self.zType2Fn_DC().real*fw[0]*self.tauScale*(self.tArr-self.tArr[idx0])) - ftlp[idx0]
                return ftlp


#########################################################################################################
# UCC28070
#########################################################################################################
class UCC28070:
        def __init__(self, plotpdf=False):
                #######################
                # DESIGN PARAMETERS
                #######################
                #       Voltage Source AC V RMS
                self.Vac_rms=240
                self.Fac=60
                #       Voltage Output DC V
                self.Vout=600
                self.Vovp=self.Vout*1.06  #by definition
                self.Vdrop=20  #constraints on Cout; max Vdrop in time thalf; start with Vripple instead?
                #       Power Total W RMS
                self.Ptot=1600
                #       Frequency Switching PWM
                self.Fpwm=60*3334   #use int multiple of Fac
                self.Fsync=2*self.Fpwm
                self.Dmax=0.98
                #######################
                # INTERNAL PARAMETERS
                #######################
                #       Regulation Voltages
                self.Vreg_vsense=3.0            #not adjustable
                self.Vreg_vinac=3.0             #Disable Synthesis --> adjustable: optimal ~2.7774 = sqrt(Kvff_max)
                #self.Vreg_vinac=2.7774             #Disable Synthesis --> adjustable: optimal ~2.7774 = sqrt(Kvff_max)
                self.DVvao=3.6                  #Vvao range, 5-1=4.0 max, leave room for power loss, vac fluctuations
                #self.Vvao_bounds=[0.0, 5.0]
                self.Vreg_currsense=3.2         #max cm volts input 3.6V, leave room for peaks
                self.Vreg_currpeak=4.0          #must have room for ILxPeak, up to 4V
                self.DVramp=4.0                 #not adjustable
                self.Vref=6.0
                #       Voltage Amplifier Transconductance VAO
                self.Gmv=70.0e-6   #uS
                #       Current Amplifier Transconductance CAOx
                self.Gmc=100.0e-6  #uS
                #       Current Multiplier IMO
                self.Kmult=17.0e-6 #uA
                #######################
                # COMPONENTS
                #######################
                #       Inductors
                self.Lx=540.0e-6
                #       VSENSE divider
                self.RAvs=3.0e6;  self.RBvs=self.get_RB(self.RAvs, self.Vreg_vsense, self.Vout)   #0.2mA --> only 120mW
                self.Cvs=1.0e-9
                #       VINAC divider (div ratio should match vsense ratio, but with NO synthesis and Kvff error, we can in fact rescale vinac a bit to cancel error, best match actual power output
                # ??? NOT full 3V range, most useful for 120VAC applications!! more accurately detect zero crossing and better signal/noise!! hooray for no current synthesis)
                # BUT can I use full scale and simply rescale Rimo to compensate for rescaled kvff? I think yes...
                self.RAac=3.0e6;  self.RBac=self.get_RB(self.RAac, self.Vreg_vinac, self.Vout)
                #self.RAac=3.0e6;  self.RBac=self.get_RB(self.RAac, self.Vreg_vinac, self.Vac_rms*sqrt(2.0))
                self.Cac=1.0e-9
                #       Dither / Sync
                self.Rrdm=0  #GND:no sync no dither; >5V:sync to pulse on RDM;
                self.Ccdr=0  #dither off: CDR to VREF;
                #       IMO current mult output; Rimo=((1/2)*Iin_pk*Rsequiv)/Imo_max; Imo_max=17uA*(Vac_rms_lopk)*(Vvao-1)/kvff; kvff=0.398, Vvao=5V #self.Rimo=17500 #using 3V/171uA
                self.Rimo=self.get_Rimo()
                self.Rsequiv=self.get_Rsequiv()
                #       Current Synthesis downslope; VRsynth = 3.0 pm0.09 V; TI info USELESS. propto L*(RB/(RA+RB))/tdrop, tdrop=1/(Rsequiv*0.1nF), Rsequiv=3V/Ipkavg
                self.Enable_Synthesis=0    #disable synthesis if full-time current sensing used
                self.Rsyn=self.get_Rsyn()
                #       Peak Current limit; Vpklimit = 3.3 pm0.03 V at CSx rising
                self.RApk=6.8e3;  self.RBpk=self.get_RB(self.RApk, self.Vreg_currpeak, self.Vref)  #choose to draw <<1mA from Vref=6V #self.RBpk=6.6e3 #3.5 when div app to Vref=6V
                #       Oscillator Rate; Rrt(kOhm) = 15000 / fsync(kHz); seems internal C=66.6666pF #Duty Cycle Max; Rdmx = (Rrt)*(2*Dmax-1-Dsync); Dsync=tsync/(1/fsync); tsync>=200nS; internal osc -->Dsync=0
                self.PWsyncExt=250.0e-9         #external sync pulse-width (not period); minimum from datasheet -> 200.0e-9
                self.Rrt, self.Rdmx = self.get_Rrt_Rdmx()
                #       Vout filter
                self.Cout=self.get_Cout(self.Vdrop)
                #       VAO compensation type 2
                self.fracDVvao=0.02  #fraction of DVvao equivalent to Vripple
                self.Cpv=self.get_Cpv()
                self.Fpv=self.get_Fpv()
                self.Fzv=(1/6)*self.get_Fpv()  #set zero freq to fraction of pole freq
                self.Rzv, self.Czv = self.get_Rzv_Czv()
                #       CAOx compensation type 2
                self.Fzc=(1/12)*self.Fpwm
                self.Gpsc=abs(self.get_Gpsc(self.Fzc))  #Gain: power stage current (magnitude)
                self.Rzc, self.Czc = self.get_Rzc_Czc()
                self.Fpc=(1/3)*self.Fpwm
                self.Cpc=self.get_Cpc()
                #       Vref,Vcc filters
                self.Cvref=1.0e-6       #minimum 22nF
                self.Cvcc=10.0e-6       #minumum 0.1uF
                #       Soft Start; recommend Css >= Czv to limit overshoot; try optimal finite time instead
                self.Css=self.get_Css(0.5) #300ms too fast?
                #######################
                # STATE VARIABLES
                #######################
                self.Clock=0
                self.Rload=225.0        #ideal pure resistive load
                self.EnablePh1=True     #False-> disable second phase (Ph1)
                self.SkipHalfMod=1      #1->no skip  n->skip n
                self.SkipHalf=False     #dynamic, computed in VLoop

                self.gppPlotter=gnuPlotter()
                self.PlotPDF=plotpdf
                self.PlotPrefix="pdfplots/"

                self.VL0 = self.VoltageLoop(self)
                self.CL0 = self.CurrentLoop(self)
                self.PM0 = self.PowerMeter()

        #########################################################################################################
        # OSCILLATOR : Rate, Duty Max
        def get_DsyncExt(self, tsync, tpropdly=100.0e-9):
                return self.PWsyncExt*self.Fsync        #no prop dly
                #return (self.PWsyncExt-tpropdly)*self.Fsync        #with prop dly

        def get_Rrt_Rdmx(self):
                rrt = (15.0e9 / self.Fsync)
                rdmx = rrt * (2.0*self.Dmax - 1.0 - self.get_DsyncExt(self.PWsyncExt))
                return rrt, rdmx

        def get_RB(self, ra, numer, denom):
                kfrac=numer/denom
                return (kfrac*ra)/(1.0-kfrac)

        def get_Css(self, tss):
                return tss*(10.0e-6/2.25) #10uA charging cap to 2.25V in ~tss sec

        #########################################################################################################
        # VSENSE, VAO : Voltage Sense, Voltage Amplifier
        def get_KVsense(self):
                return self.RBvs/(self.RAvs+self.RBvs)

        def get_Iimo(self, vinac, vvao):
                return (self.Kmult)*vinac*(vvao-1.0)/self.get_Kvff()

        def get_Vimo(self, iimo):
                return iimo*self.Rimo*(1.0+float(self.EnablePh1==False))*(1.0+float(self.SkipHalf==True))

        def get_Cout(self, vdrop):
                return float( ((1.0/(2.0*self.Fac))*self.Ptot)/((1/2)*((self.Vout)**2 - (self.Vout-vdrop)**2)) )

        def get_Vripple(self):
                return self.get_ICout(sqrt(2)*self.Vac_rms, 2*self.Ptot)*(1/(2*pi*2*self.Fac*self.Cout))

        # VAO : Voltage Amplifier Compensation
        def get_Cpv(self):
                ztmp=(self.fracDVvao*self.DVvao)/(self.get_Vripple()*self.get_KVsense()*self.Gmv)
                #print("ztmp",ztmp)
                return 1.0/(2*pi*(2*self.Fac)*ztmp)

        def get_Fpv(self):
                return (1/(2*pi))*sqrt( (self.get_KVsense()*self.Gmv*self.Ptot)/(self.DVvao*self.Vout*self.Cout*self.get_Cpv()) )

        def get_Rzv_Czv(self):
                rzv=1/(2*pi*self.get_Cpv()*self.get_Fpv())
                czv=1/(2*pi*rzv*self.Fzv)
                return rzv, czv

        # Gcv : Voltage Compensation Network Gain
        def get_Gcv(self, f):
                return self.get_KVsense() * self.Gmv * self.get_Zcomp_type2(1.0, f, self.Rzv, self.Czv, self.Cpv)

        # Gpsv : Voltage Loop Power Stage Gain
        def get_Gpsv(self, f):
                return ((self.Ptot/self.Vout)/(1j*2*pi*f*self.Cout))/self.DVvao

        # GTOTv : TOTAL VOLTAGE LOOP GAIN
        def get_GTOTv(self, f):
                return self.get_Gcv(f) * self.get_Gpsv(f)


        #########################################################################################################
        # IMO : Current Multiplier Output
        def get_IMO_max(self):
                #if Vac increases or decreases, what is possible change in current? within kvff level?
                return (self.Kmult)*(self.get_Vvinac(sqrt(2)*self.Vac_rms))*(self.DVvao)/self.get_Kvff()

        def get_Rimo(self):
                return float(self.Vreg_currsense/self.get_IMO_max())

        def get_KVvinac(self):
                return self.RBac/(self.RAac+self.RBac)

        def get_Vvinac(self, vinac):
                return vinac*self.get_KVvinac()

        def get_Vinac_Correction_Factor(self):
                vpk=self.get_Vvinac(sqrt(2)*self.Vac_rms)
                kvpk=self.get_Kvff()
                result=[vpk, kvpk, (vpk**2)/2.0, ((vpk**2)/2.0)/kvpk, sqrt(2.0*kvpk), sqrt(2.0*kvpk)/vpk]
                if self.PlotPDF:
                        #print("result:",result)
                        def fk(k):
                                return 3.0*k/1000.0
                        self.kvffDat = [ [ fk(k), self._get_Kvff(fk(k)), fk(k)**2/2.0, (fk(k)**2/2.0)/self._get_Kvff(fk(k)), vpk, kvpk, (vpk**2/2.0)/(kvpk), 1.0] for k in range(1,1000) ]
                        self.kvff_str= "Vvinac, Kvff, a*Pvinac, Pratio, VPeak, KVPeak, KVPeakRatio, Perfect"
                        self.kvff_lblstr=["Vvinac:Kvff, Vvinac:a*Pvinac, VPeak:Kvff","Vvinac:Pratio, Vvinac:KVPeakRatio, VPeak:Pratio, Vvinac:Perfect"]
                        self.gppPlotter.gpp_Plot(self.kvffDat, self.kvff_str, self.kvff_lblstr, toFile=self.get_FileName("kvff_ratios.pdf"))
                        #quit()
                return result  #Kvff discrete levels cause error in IMO current command, fix this

        def get_Kvff(self):
                # NOTE this will be irrelevant when current synthesis is disabled thanks to full-time inductor current sensing via hall effect sensors
                #  which has huge advantages for low VAC applications, we can scale Vinac divider for full accuracy zero detect and best signal/noise
                return self._get_Kvff(self.get_Vvinac(sqrt(2)*self.Vac_rms))

        def _get_Kvff(self, v, dvin_sign=1):
                #klst=[[0,0], [1.0, 0.398], [1.2, 0.600 ], [1.4, 0.839], [1.65, 1.156], [1.95, 1.604], [2.25, 2.199], [2.6, 2.922], [3.0*33, 3.857]]
                return np.piecewise(v, \
                        [v<1.0, (v>=1.0)&(v<1.2), (v>=1.2)&(v<1.4), (v>=1.4)&(v<1.65), (v>=1.65)&(v<1.95), (v>=1.95)&(v<2.25), (v>=2.25)&(v<2.6), v>=2.6], \
                        [0.398, 0.600, 0.839, 1.156, 1.604, 2.199, 2.922, 3.857] )

        # CAOx : Current Synthesizer, Amplifiers
        def get_Rsequiv(self):
                return self.Vreg_currsense/self.get_ILxTop(self.Vac_rms, self.Ptot/2)

        def get_Rsyn(self):
                return float(self.Enable_Synthesis)*self.Lx*self.get_KVsense()/((0.1e-9)*self.Rsequiv) #appears 100pF cap inside

        def get_Rzc_Czc(self):
                rzc=1.0/(self.Gmc*abs(self.get_Gpsc(self.Fzc)))
                czc=1.0/(2*pi*self.Fzc*rzc)
                return rzc, czc

        def get_Cpc(self):
                return 1/(2*pi*self.Fpc*self.Rzc)

        # Gcc : Current Compensation Network Gain
        def get_Gcc(self, f):
                return self.Gmc * self.get_Zcomp_type2(1.0, f, self.Rzc, self.Czc, self.Cpc)

        # CAOx : Current Loop Power Stage Gain
        def get_Gpsc(self, f):
                return (self.Vout*self.Rsequiv)/(1j*2*pi*f*self.Lx*self.DVramp)

        # GTOTc : TOTAL CURRENT LOOP GAIN
        def get_GTOTc(self, f):
                return self.get_Gcc(f) * self.get_Gpsc(f)

        # AC voltage, current
        def get_Vinac_Sin(self, theta):         # Pi Periodic, PosDef for real theta
                return sqrt(2)*self.Vac_rms*np.fabs(np.sin(theta))

        def get_Iin_Sin(self, theta, prms):     # Pi Periodic, PosDef for real theta
                return sqrt(2)*(prms/self.Vac_rms)*np.fabs(np.sin(theta))

        # Idealized approximate quantities
        def get_Duty(self, vin):
                return (self.Vout-vin)/self.Vout

        def get_Iin(self, vin, pout):
                return pout/vin

        def get_ILx(self, vin, poutLx):
                return self.get_Iin(vin, poutLx)

        def get_ILxTop(self, vin, poutLx):
                return self.get_ILx(sqrt(2)*vin, 2*poutLx)

        def get_ILxPeak(self, fracsign=1):
                return (fracsign*self.get_dILxTop(self.Vac_rms)/2.0) + self.get_ILxTop(self.Vac_rms, self.Ptot/2)

        def get_dILx(self, vin, freq):
                return vin*(self.get_Duty(vin)/freq)/(self.Lx)

        def get_dILxTop(self, vin):
                return self.get_dILx(sqrt(2)*vin, self.Fpwm)

        def get_LxMin_dILx(self, vin, poutLx, freq, fracilx):
                return sqrt(2)*poutLx / ( (self.get_ILxTop(vin, poutLx)**2)*fracilx*freq )

        def get_LxMin_CCM(self, vin, poutLx, freq):
                return (vin**2)/(2*poutLx*freq)

        def get_PoutMin_CCM(self, vin, freq):
                return 2*(vin**2)/(2*self.Lx*freq)

        def get_ICout(self, vin, ptot):
                return 2*self.get_ILx(vin, ptot/2)*(1.0-self.get_Duty(vin))

        def _get_ILx_Vals(self, n, pout):
                nmax=(self.Fpwm/self.Fac)/2.0
                ph_on, ph_mid = pi*((n)/nmax), pi*((n+0.5)/nmax)
                v, i = self.get_Vinac_Sin(ph_mid), self.get_Iin_Sin(ph_mid, pout)
                d = self.get_Duty(v)
                ph_off, ph_end = pi*((n+d)/nmax), pi*((n+1)/nmax)
                dt_on, dt_off = d/self.Fpwm, (1-d)/self.Fpwm
                iL, diL = i, (v*dt_on)/(2*self.Lx)
                return ph_on, ph_off, ph_end, v, i, d, dt_on, dt_off, iL-diL, iL+diL

        def get_FileName(self, fstr):
                if self.PlotPDF:
                        return self.PlotPrefix+fstr
                else:
                        return ''

        def plotBode_VoltageLoop(self):
                ftmp=[p/100.0 for p in range(0, 300, 1)]; flo=1
                gptmp=[ self.get_Z_lst(flo*pow(10,f), self.get_GTOTv(flo*pow(10,f))) for f in ftmp]
                self.gppPlotter.gpp_Raw(gptmp, "call 'gpBode.gp' '$idat' '" + self.get_FileName("bodePlot_VoltageLoop.pdf") + "' \n")

        def plotBode_CurrentLoop(self):
                ftmp=[p/100.0 for p in range(301, 600, 1)]
                gptmp=[ self.get_Z_lst(pow(10,f), self.get_GTOTc(pow(10,f))) for f in ftmp]
                self.gppPlotter.gpp_Raw(gptmp, "call 'gpBode.gp' '$idat' '" + self.get_FileName("bodePlot_CurrentLoop.pdf") + "' \n")


        #########################################################################################################
        class PowerMeter:
                def __init__(self):
                        self.dictE={}
                        self.dictT={}

                def initLbl(self, lbl):
                        self.dictE[lbl]=0.0
                        self.dictT[lbl]=0.0

                def add_E_T(self, lbl, e, t):
                        if lbl not in self.dictE.keys():
                                self.initLbl(lbl)
                        self.dictE[lbl] += e
                        self.dictT[lbl] += t

                def Report(self):
                        print("PowerMeter: ")
                        for k in self.dictE.keys():
                                print(k, self.dictE[k]/self.dictT[k], self.dictE[k], self.dictT[k])


        #########################################################################################################
        class VoltageLoop:
                def __init__(self, parent):
                        self.parent=parent
                        self.Vvao_bounds=[0.0, 5.0]
                        self.nptsIvao=20
                        self.fcycIvao=2*self.parent.Fac
                        self.thalf=(1.0/self.fcycIvao)
                        self.dt=self.thalf*(1.0/self.nptsIvao)
                        self.Tau = np.linspace(0, np.pi, self.nptsIvao, endpoint=False)
                        self.Vinac=self.parent.get_Vinac_Sin(self.Tau)
                        self.gppPlotter=gnuPlotter()
                        self.fftIntegrator_Ivao = fftIntegrator(self.nptsIvao, self.fcycIvao, 20, 120, self.get_ZType2_VAO, self.get_ZType2_VAO_DC, self.get_Vvao, self.Vvao_bounds)

                        self.Vvao=0
                        self.QCss=0
                        self.cntThalf=0
                        self.Clock=0

                def addClockThalf(self):
                        self.cntThalf+=1
                        self.Clock+=1.0/(self.parent.Fac*2.0)

                def get_Vss(self):
                        return self.QCss/self.parent.Css

                def get_Vsense(self, vout):
                        return vout*self.parent.get_KVsense()

                def get_VsenseRef(self):
                        return min(3.0, self.get_Vss())

                def get_Ivao(self, v):  #v==Vsense; note -sign since interest in current coming out of vao
                        isat = 30.0e-6
                        iscp = 100.0e-6*float(self.get_Vss()>4.0)
                        vlo, vmid, vhi = -0.42857, -0.2, 0.42857
                        dv = v-self.get_VsenseRef()
                        rnew = -np.piecewise(dv, \
                                [(dv<vlo),      (dv>=vlo)&(dv<vmid),            (dv>=vmid)&(dv<vhi),            (dv>=vhi)], \
                                [-isat-iscp,    lambda dum: dum*self.parent.Gmv-iscp,  lambda dum: dum*self.parent.Gmv,       isat] \
                        )
                        return rnew

                def get_Vvao(self):
                        return self.Vvao

                def get_ZType2_VAO(self, bk, f):
                        return self.parent.get_Zcomp_type2(bk, f, self.parent.Rzv, self.parent.Czv, self.parent.Cpv)
        
                def get_ZType2_VAO_DC(self): #hack that returns approx integration constant for DC component: 1/(cp+cz)
                        return self.parent.get_Zcomp_type2(1, 1/(2*pi*1j), 0, self.parent.Czv, self.parent.Cpv)

                #-----------------------------------------------------------------------------------------------------------------------------------
                def stepTHALF(self, n, useCurrentLoop, plotIterDat=True, plotThalfDat=False, nskip=0):
                        if(self.Vvao<=1.0):
                                self.Vvao=1.0001
                        lvIterDat=[]
                        self.lvIterDat_str= "self.Clock, self.parent.CL0.QCOUT.get_VCout(), 100*self.get_Vss(), deltaVao, IvaoAvg, 100*self.Vvao, self.parent.Rload, 0.5*ploadavg, PinacAvg, PcoutAvg, IcoutAvg"
                        self.lvIterDat_lblstr=["Clock:VCout(), Clock:ploadavg, Clock:PinacAvg, Clock:PcoutAvg, Clock:Vvao, Clock:Vss()"]
                        self.lvThalfDat_str= "Tphase, nphon*Iinac, Vsense, nphon*Pinac, nphon*Pcout, Duty, Vimo, Ivao, IIvao.real"
                        self.lvThalfDat_lblstr=["Tphase:Ivao", "Tphase:IIvao.real", "Tphase:Duty"]
                        self.gppPlotter.open_gpp("lvThalfDat")

                        if not useCurrentLoop:
                                Duty=self.get_Duty(self.Vinac)

                        for icnt in range(-nskip, n):
                                self.QCss+=(10.0e-6*self.thalf*(self.get_Vss()<5.0))
                                if not useCurrentLoop:
                                        # assumption-> load current constant
                                        iloadavg=self.parent.CL0.QCOUT.get_VCout()/self.parent.Rload
                                        ploadavg=self.parent.CL0.QCOUT.get_VCout()*iloadavg
                                        Vimo=self.parent.get_Vimo(self.parent.get_Iimo(self.parent.get_Vvinac(self.Vinac), self.Vvao))    #assuming Vvao fluctuations small compared to Vvinac over thalf
                                        Iinac=((self.Vvao>1.0)*Vimo/self.parent.Vreg_currsense)*self.parent.get_Iin(sqrt(2)*self.parent.Vac_rms, 2*self.parent.Ptot/2)  # assumption-> vinac power commanded by Vimo instantly, and duty known
                                        Icout=(1.0-Duty)*(Iinac/Duty)*(float(not self.parent.SkipHalf))           #current from one inductor into cout averaged over dt
                                        Icout_net=(1.0+float(self.parent.EnablePh1))*Icout-iloadavg              #n inductors minus load current
                                        DVCout=np.cumsum(Icout_net*dt)/self.parent.Cout        #variation in VCout vs dt
                                        VCout=self.parent.CL0.QCOUT.get_VCout()+DVCout
                                        self.parent.CL0.QCOUT.QCout += (1.0+float(self.parent.EnablePh1))*self.thalf*np.mean(Icout) - self.thalf*iloadavg  #some redundancy, more than one way to calc this
                                        Pinac=self.Vinac*Iinac
                                        Pcout=VCout*Icout

                                if useCurrentLoop:
                                        self.parent.CL0.StatsCL.__init__()
                                        self.parent.CL0.stepPWM(0, 1667, plotDatPks=(icnt==n-1), debug=False)
                                        phTest, Duty, Vimo, Iinac, Icout, Pinac, Pcout, VCout = self.parent.CL0.StatsCL.getInterpAll(self.Tau)       #Iinac,Icout changed def here to instantaneous
                                        iloadavg=np.mean(VCout/self.parent.Rload)
                                        ploadavg=np.mean(VCout*iloadavg)

                                Ivao=self.get_Ivao(self.get_Vsense(VCout))      
                                IIvao=self.fftIntegrator_Ivao.fftIntegrate(Ivao, flagSavePrev=True)  #True:save integrand to history
                                self.Vvao=max(self.Vvao_bounds[0], min(self.Vvao + IIvao[-1].real, self.Vvao_bounds[1]))

                                if plotThalfDat:  #realtime display, never to file
                                        nphon=(1.0+float(self.parent.EnablePh1))
                                        Tphase, Vsense = self.thalf*self.Tau/np.pi, self.get_Vsense(VCout)
                                        self.lvThalfDat=np.vstack(eval("["+self.lvThalfDat_str+"]"))
                                        #self.gppPlotter.gpp_Plot(self.lvThalfDat.transpose(), self.lvThalfDat_str, self.lvThalfDat_lblstr, self.gppPlotter.gpp["lvThalfDat"], lineWidth=2.0)
                                        self.gppPlotter.gpp_Plot(self.lvThalfDat.transpose(), self.lvThalfDat_str, self.lvThalfDat_lblstr, "lvThalfDat", lineWidth=2.0)

                                if icnt>=0:
                                        deltaVao, IvaoAvg, PinacAvg, PcoutAvg, IcoutAvg  = IIvao[-1], np.mean(Ivao), np.mean(Pinac), np.mean(Pcout), np.mean(Icout)
                                        lvIterDat.append(eval("["+self.lvIterDat_str+"]"))
                                        self.addClockThalf()

                                self.parent.SkipHalf=(self.parent.SkipHalfMod>1) and not ((self.cntThalf%self.parent.SkipHalfMod)==0)

                        if plotIterDat:
                                #self.gppPlotter.gpp_Plot(lvIterDat, self.lvIterDat_str, self.lvIterDat_lblstr, lineWidth=1.0, toFile=self.parent.get_FileName("lvIterDat.pdf"))
                                gppconfig="set grid ytics back ; \n set grid lw 0.3 lt 1 lc rgb 0x80888888 ; \n"
                                self.gppPlotter.gpp_Plot(lvIterDat, self.lvIterDat_str, self.lvIterDat_lblstr, lineWidth=1.0, toFile=self.parent.get_FileName("lvIterDat.pdf"), config=gppconfig)


        #########################################################################################################
        class CurrentLoop:
                Vcaox_bounds=[0.0, 6.0]
                nptsIcaox=200           #puts 1st stage root finding error < 0.5%
                Vramp_min=0.7           #not adjustable

                def __init__(self, parent):
                        self.parent=parent
                        self.fcycIvao=2*self.parent.Fac
                        self.thalf=(1.0/(2*self.parent.Fac))
                        self.fcycIcaox=self.parent.Fpwm
                        self.tpwm=(1.0/self.parent.Fpwm)
                        self.npwmMax=round(self.parent.Fpwm/(2*self.parent.Fac))
                        self.ntonPh=np.pi/self.npwmMax
                        self.Tau = np.linspace(0.0, 1.0, self.nptsIcaox, endpoint=False)
                        self.lcDatPks=[]
                        self.lcDatRT=[]

                        self.gppPlotter=gnuPlotter()
                        self.Vimo=None
                        self.StatsCL=self.statsCL()

                        self.QCOUT = self.QCout_Tau_Obj(self, 0.0)
                        self.IL0 = self.ILx_Tau_Obj(self, self.QCOUT)
                        self.CA0 = self.CAx_CSx_Obj(self, self.IL0)

                        self.fftIntegrator_Icaox = fftIntegrator(self.nptsIcaox, self.fcycIcaox, 20, 20, self.get_ZType2_CAOx, self.get_ZType2_CAOx_DC, self.CA0.get_Vcaox, self.Vcaox_bounds)
                        self.Vramp = self.get_Vramp(self.Tau, self.Vramp_min, self.parent.DVramp)
                        self.Clock=self.parent.Clock
                        self.cntTPWM=0


                def addClockTPWM(self, n=1):
                        self.cntTPWM+=n
                        self.Clock+=n*1.0/(self.parent.Fpwm)

                def get_Vramp(self, tau, vmin, dv):
                        return vmin + tau*dv    #tau in [0,1]


                #########################################################################################################
                class statsCL:
                        def __init__(self):
                                self.nPh0=[]
                                self.TauOn=[]
                                self.Vimo=[]
                                self.ILx_Vinac=[]
                                self.ILx_VCout=[]
                                self.Pinac=[]
                                self.Pcout=[]
                                self.VCoutAvg=[]

                        def addStats(self, nph0, tauon, vimo, ilxvinac, ilxvcout, pinac, pcout, vcoutavg):
                                self.nPh0.append(nph0)
                                self.TauOn.append(tauon)
                                self.Vimo.append(vimo)
                                self.ILx_Vinac.append(ilxvinac)
                                self.ILx_VCout.append(ilxvcout)
                                self.Pinac.append(pinac)
                                self.Pcout.append(pcout)
                                self.VCoutAvg.append(vcoutavg)

                        def getPCoutAvg(self):
                                return np.mean( np.array(self.VCoutAvg)*np.array(self.ILx_VCout)*(1.0-np.array(self.TauOn)) )

                        def getInterpAll(self, taunew):
                                r=[ np.interp(taunew, self.nPh0, arrcurr) for arrcurr in [self.nPh0, self.TauOn, self.Vimo, self.ILx_Vinac, self.ILx_VCout, self.Pinac, self.Pcout, self.VCoutAvg] ]
                                return r[0],r[1],r[2],r[3],r[4],r[5],r[6],r[7]


                #########################################################################################################
                class CAx_CSx_Obj:
                        def __init__(self, parentcl, ilxobj, vcaox=0.0):
                                self.parentCL=parentcl
                                self.ILxobj=ilxobj
                                self.Vcsx=None
                                self.Icaox=None
                                self.IIcaox=None
                                self.PWMx=None
                                self.solnFound=False
                                self.TauCurr=None
                                self.TauRoot=None
                                self.TauRootFin=None
                                self.ErrAvg_PWMx=0.0
                                self.Vcaox=vcaox
                                self.Vcaox_prev=0.0
                                self.tStamp=-1.0

                        def reset(self):
                                self.solnFound, self.TauRoot, self.TauRootFin, self.Vcaox_Next = False, None, None, None

                        def get_Vcsx(self, ilx):
                                return ilx*self.parentCL.parent.Rsequiv

                        def get_Icaox(self, vimo, vcsx):
                                return self.parentCL.parent.Gmc*(vimo-vcsx)

                        def get_Vcaox(self):
                                return self.Vcaox

                        def get_PWMx(self):
                                return self.parentCL.Vramp-(self.get_Vcaox() + self.IIcaox.real)

                        def configure(self, nph, tauroot, flagSavePrev):
                                self.ILxobj.configure(nph, tauroot)
                                self.Vcsx=self.get_Vcsx(self.ILxobj.get_ILxCurrArr(self.parentCL.Tau))
                                self.Icaox=self.get_Icaox(self.parentCL.Vimo, self.Vcsx)
                                self.IIcaox=self.parentCL.fftIntegrator_Icaox.fftIntegrate(self.Icaox, flagSavePrev)
                                self.PWMx=self.parentCL.Vramp-(self.get_Vcaox()+self.IIcaox.real)
                                self.tStamp=self.parentCL.Clock

                        def peakLimitTrigger(self):
                                return (self.Vcsx > self.parentCL.parent.Vreg_currpeak).any()

                        def has_Root(self, arr):
                                return ( (arr > 0.0).any() and (arr < 0.0).any() )

                        def has_Root_PWMx(self):
                                return self.has_Root(self.PWMx)

                        def has_Root_PeakLimit(self):
                                return self.has_Root(self.Vcsx[:round(self.parentCL.nptsIcaox*self.TauRoot)] - self.parentCL.parent.Vreg_currpeak)

                        def get_FnRoot_PWMxArr(self):
                                return self.PWMx[round(self.parentCL.nptsIcaox*self.TauRoot)]

                        def get_TauRoot_PeakLimit(self, tlo=-0.1, thi=1.1):
                                result=None
                                try:
                                        result=findRoot( lambda t: np.interp(t, self.parentCL.Tau, self.Vcsx - self.parentCL.parent.Vreg_currpeak), tlo, thi )
                                except:
                                        pass
                                return result

                        def get_TauRoot_PWMx(self, tlo=-0.1, thi=1.1):
                                result=None
                                try:
                                        result = findRoot( lambda t: np.interp(t, self.parentCL.Tau, self.PWMx), tlo, thi )
                                except:
                                        pass
                                return result

                        def get_TauRoot_Any(self, tlo=-0.1, thi=1.1):
                                rootpwm, rootpklim = None, None
                                if self.has_Root_PWMx():
                                        rootpwm = self.get_TauRoot_PWMx(tlo, thi)
                                if self.has_Root_PeakLimit():
                                        rootpklim = self.get_TauRoot_PeakLimit(tlo, thi)
                                elif (rootpwm is None) and self.peakLimitTrigger():
                                        rootpklim = 0.0
                                if rootpklim is None:                   #easy, no peaklimit, simple root
                                        result = rootpwm, False
                                elif rootpwm is not None:               #both roots found, use earliest
                                        result = np.fmin(rootpwm, rootpklim), (rootpklim<=rootpwm)
                                        #print("Both Roots Found. rootpwm, rootpklim, result == root, is_peak_limited", rootpwm, rootpklim, result)
                                else:                                   #no roots found, returning None or zero
                                        result = rootpklim, True
                                        print("No Roots Found. rootpwm, rootpklim, result == root, is_peak_limited", rootpwm, rootpklim, result)
                                return result

                        def get_FnRoot_PWMx(self, tauroot):
                                return np.interp(tauroot, self.parentCL.Tau, self.PWMx)

                        def get_IIcaox_Last(self):
                                return self.IIcaox[-1].real

                        def save_Vcaox_Next(self):
                                self.Vcaox_Next=max(self.parentCL.Vcaox_bounds[0], min(self.Vcaox + self.get_IIcaox_Last(), self.parentCL.Vcaox_bounds[1]))

                        def apply_Vcaox_Next(self):
                                self.Vcaox_prev=self.get_Vcaox()
                                self.Vcaox=self.Vcaox_Next


                #########################################################################################################
                class ILx_Tau_Obj:
                        def __init__(self, parentcl, qcobj, ilx=0.0):
                                self.parentCL=parentcl
                                self.Lx=self.parentCL.parent.Lx         #convenience copy
                                self.QCobj=qcobj
                                self.ILx=ilx
                                self.Aac=sqrt(2)*self.parentCL.parent.Vac_rms
                                self.Wac=2*np.pi*self.parentCL.parent.Fac
                                self.Kac=self.Aac/(self.Wac*self.Lx)
                                self.ILxCurrArr=None

                        def configure(self, nph0, tauOn):
                                self.nPh0 = nph0
                                self.TauOn=tauOn
                                self.QCobj.configure(self.TauOn)
                                self.ILxCurrArr=None    #must regenerate after config
                                self.set_ILxCurrArr(self.parentCL.Tau)
                                self.ILx_Next=None

                        def get_DILx_Vinac(self, tauI, tauF, vcout=0): #using positive current convention
                                #       i(t) ~= -I0 - ( Vinac - VCout )*t/L
                                #               Vinac ~= A*sin(ph0 + dph/2);  dph=tpwm*TauRoot  --> mid interval voltage approx, matters little at small dt
                                #               VCout ~= QCout/Cout
                                dphI, dphF = self.parentCL.ntonPh*tauI, self.parentCL.ntonPh*tauF
                                dt=self.parentCL.tpwm*(tauF-tauI)
                                vinac=self.Aac*np.sin(self.nPh0+dphI+dphF/2.0)
                                return (vinac-vcout)*dt/self.Lx
                        def zget_DILx_Vinac(self, tauI, tauF, vcout=0): #using positive current convention
                                #return ( (sqrt(2)*self.parentCL.parent.Vac_rms/(2*np.pi*self.parentCL.parent.Fac))*(np.cos(theta)-np.cos(theta+self.parentCL.ntonPh*tau)) - vcout*self.parentCL.tpwm*tau )/self.Lx
                                #                       i(t) =  -C0 + K*cos(ph0+w*t) + VCout*t/L        assuming VCout=0 when t<0; since cos decreasing, i getting more negative as desired
                                #                       -C0 = -I0 - K*cos(ph0)                          sets initial current to I0 < 0
                                dphI=self.parentCL.ntonPh*tauI
                                dphF=self.parentCL.ntonPh*tauF
                                dt=self.parentCL.tpwm*(tauF-tauI)
                                dCos=-(np.cos(self.nPh0+dphF) - np.cos(self.nPh0+dphI))
                                return self.Kac*dCos - vcout*dt/self.Lx

                        def _get_ILx_On(self, tauF, flagICout=False):
                                if not flagICout:
                                        return self.ILx + self.get_DILx_Vinac(0.0, tauF)
                                else:
                                        return 0.0

                        def _get_ILx_Off(self, tauF):
                                result = self._get_ILx_On(self.TauOn) + self.get_DILx_Vinac(self.TauOn, tauF, self.parentCL.QCOUT.get_VCout_Tau(self.TauOn))
                                return np.fmax(0.0, result)

                        def get_ILx_Tau(self, tau, flagICout=False):     #prevent negative currents for root finding, but retain ability to explicitly check for DCM operation
                                result = np.piecewise(tau, \
                                        [tau<0.0,  (tau>=0.0)&(tau<=self.TauOn),        (tau>self.TauOn)&(tau<=1.0),            (tau>1.0)],\
                                        [0,         lambda tau: self._get_ILx_On(tau, flagICout),   lambda tau: self._get_ILx_Off(tau),     0] )
                                return np.fmax(result, 0.0)

                        def get_ILxAvgOn(self):
                                return (self.get_ILx_Tau(0.0)+self.get_ILx_Tau(self.TauOn))/2.0
                        def get_ILxAvgOff(self):
                                return (self.get_ILx_Tau(1.0)+self.get_ILx_Tau(self.TauOn))/2.0
                        def get_ILx_TauOn(self):
                                return self.get_ILx_Tau(self.TauOn)

                        def get_DE_Vinac(self, tauI, tauF, vcout=0.0):
                                #       E(t)-E(0) ~= Vinac*( I0*t + (1/2)*(Vinac - VCout)*t**2/L )
                                dphI=self.parentCL.ntonPh*tauI
                                dphF=self.parentCL.ntonPh*tauF
                                dt=self.parentCL.tpwm*(tauF-tauI)
                                vinac=self.Aac*np.sin(self.nPh0+dphI+dphF/2.0)
                                i0 = self.get_ILx_Tau(tauI)
                                return vinac*dt*(i0 + 0.5*(vinac-vcout)*dt/self.Lx)
                        def zget_DE_Vinac(self, tauI, tauF, vcout=0.0):
                                # E(t)-E(0) = (A/w)*C(t)*I0 + (1/2)*K*(A/w)*(C(t))**2 - VCout*K*t*C(t) - VCout*(K/w)*S(t)       # Don't forget there are 2 inductors!
                                #                       C(t) = -(cos(ph0+w*t)-cos(ph0)) > 0
                                #                       S(t) = (sin(ph0+w*t)-sin(ph0))
                                #                       K  = (A/(w*L))
                                dphI=self.parentCL.ntonPh*tauI
                                dphF=self.parentCL.ntonPh*tauF
                                dt=self.parentCL.tpwm*(tauF-tauI)
                                dCos=-(np.cos(self.nPh0+dphF) - np.cos(self.nPh0+dphI))
                                dSin=(np.sin(self.nPh0+dphF) - np.sin(self.nPh0+dphI))
                                i0 = self.get_ILx_Tau(tauI)
                                return (self.Aac/self.Wac)*dCos*(i0 + 0.5*self.Kac*dCos) - vcout*self.Kac*(dt*dCos + dSin/self.Wac)

                        def get_DE_VCout(self, tauI, tauF, vcout):
                                #       E(t)-E(0) ~= -VCout*( I0*t + (1/2)*(Vinac - VCout)*t**2/L )
                                dphI=self.parentCL.ntonPh*tauI
                                dphF=self.parentCL.ntonPh*tauF
                                dt=self.parentCL.tpwm*(tauF-tauI)
                                vinac=self.Aac*np.sin(self.nPh0+dphI+dphF/2.0)
                                i0 = self.get_ILx_Tau(tauI)
                                #return -vcout*dt*(i0 + 0.5*(vinac-vcout)*dt/self.Lx)
                                return np.fmin(-vcout*dt*(i0 + 0.5*(vinac-vcout)*dt/self.Lx), 0.0)      #DIODE --> ALWAYS NEGATIVE
                        def zget_DE_VCout(self, tauI, tauF, vcout):
                                # E(t)-E(0) = -VCout*(I0 + K*cos(ph0))*t + VCout*(K/w)*S(t) + (1/2)*VCout**2*(t**2)/L
                                dphI=self.parentCL.ntonPh*tauI
                                dphF=self.parentCL.ntonPh*tauF
                                dt=self.parentCL.tpwm*(tauF-tauI)
                                dSin=(np.sin(self.nPh0+dphF) - np.sin(self.nPh0+dphI))
                                i0 = self.get_ILx_Tau(tauI)
                                return vcout*( (-i0 - self.Kac*np.cos(self.nPh0+dphI) + 0.5*vcout*dt/self.Lx)*dt + self.Kac*dSin/self.Wac )

                        def set_ILx(self, ilx):
                                self.ILx=max(0.0, ilx)

                        def save_ILx_Next(self):
                                self.ILx_Next=self.get_ILx_Tau(1.0)

                        def apply_ILx_Next(self):
                                tOn=self.parentCL.tpwm*(self.TauOn)
                                deVinacOn=self.get_DE_Vinac(0.0, self.TauOn)
                                tOff=self.parentCL.tpwm*(1.0-self.TauOn)
                                vcoutRoot=self.QCobj.get_VCout_Tau(self.TauOn)
                                deVinacOff=self.get_DE_Vinac(self.TauOn, 1.0, vcoutRoot)
                                self.parentCL.parent.PM0.add_E_T("DE_Vinac_On", deVinacOn, tOn)
                                self.parentCL.parent.PM0.add_E_T("DE_Vinac_Off", deVinacOff, tOff)
                                self.parentCL.parent.PM0.add_E_T("DE_Vinac_Tot", deVinacOn, tOn)
                                self.parentCL.parent.PM0.add_E_T("DE_Vinac_Tot", deVinacOff, tOff)
                                de_ilx_vcout = self.get_DE_VCout(self.TauOn, 1.0, vcoutRoot)
                                self.parentCL.parent.PM0.add_E_T("DE_ILx_Net", deVinacOn, tOn)
                                self.parentCL.parent.PM0.add_E_T("DE_ILx_Net", deVinacOff, tOff)
                                self.parentCL.parent.PM0.add_E_T("DE_ILx_Net", de_ilx_vcout, self.parentCL.tpwm*(1.0-self.TauOn)*0)
                                self.set_ILx(self.ILx_Next)

                        def set_ILxCurrArr(self, tauarr):
                                self.ILxCurrArr = self.get_ILx_Tau(tauarr)

                        def get_ILxCurrArr(self, tauarr):
                                #self.set_ILxCurrArr(tauarr)
                                return self.ILxCurrArr

                        def get_ILxCout_Arr(self):
                                return self.get_ILx_Tau(self.parentCL.Tau, flagICout=True)


                #########################################################################################################
                class QCout_Tau_Obj:
                        def __init__(self, parentcl, qcout=0.0):
                                self.parentCL=parentcl
                                self.QCout=qcout
                                self.Cout=self.parentCL.parent.Cout  #convenience copy

                        def configure(self, tauOn):
                                self.TauOn=tauOn
                                self.QCout_Next=None

                        def _get_ICout_On(self, tau):
                                tRC=(self.parentCL.parent.Rload*self.Cout)
                                return -(self.QCout/tRC)*np.exp(-self.parentCL.tpwm*tau/(tRC))

                        def _get_ICout_Off(self, tau):
                                return self._get_ICout_On(0*tau)

                        def get_ICout_Tau(self, tau, flagPosOk=False):
                                result = np.piecewise(tau, \
                                        [tau<0.0,  (tau>=0.0)&(tau<=self.TauOn),        (tau>self.TauOn)&(tau<=1.0),            (tau>1.0)],\
                                        [0,         lambda tau: self._get_ICout_On(tau),   lambda tau: self._get_ICout_Off(tau),     0] )
                                if flagPosOk:
                                        return result
                                else:
                                        return np.fmin(result, 0.0)

                        def get_DQCout_Rload_dT(self, q0, t):
                                # Q(t)-Q(0) = Q0*( exp(-t/(Rload*Cout)) - 1 )
                                return q0*(np.exp(-t/(self.parentCL.parent.Rload*self.Cout))-1.0)

                        def get_DE_Rload_dQ(self, qI, qF):
                                # E(t)-E(0) = (1/2)*( Q(t)**2 - Q(0)**2 )/C
                                return 0.5*(qF**2-qI**2)/self.Cout

                        def get_DE_Rload_dT(self, q0, t):
                                # DECout(t) = ELx - t*V(0)**2/Rload ; assumption: constant current
                                return -t*((q0/self.Cout)**2)/self.parentCL.parent.Rload

                        def _get_QCout_On(self, tau):
                                return self.QCout + self.get_DQCout_Rload_dT(self.QCout, self.parentCL.tpwm*tau)

                        def _get_QCout_Off(self, tau):
                                # DECout(t) = ELx - t*V(0)**2/Rload
                                # QCout_new = sqrt( QCout**2 + 2*DECout*Cout )
                                de_ilx_vcout = -(1.0+float(self.parentCL.parent.EnablePh1))*self.parentCL.IL0.get_DE_VCout(self.TauOn, tau, self._get_QCout_On(self.TauOn)/self.Cout)
                                de_cout_rload = self.get_DE_Rload_dT(self.get_QCout_TauOn(), self.parentCL.tpwm*(tau-self.TauOn))
                                de_cout_net = de_ilx_vcout + de_cout_rload      #notice this does change voltage a bit
                                qnew2=np.fmax(0.0, (self._get_QCout_On(self.TauOn)**2 + 2.0*de_cout_net*self.Cout))
                                return np.sqrt(qnew2)

                        def get_QCout_Tau(self, tau):     #prevent negative charge, ever useful in cap context?
                                result = np.piecewise(tau, \
                                        [tau<0.0,  (tau>=0.0)&(tau<=self.TauOn),          (tau>self.TauOn)&(tau<=1.0),              (tau>1.0)],\
                                        [0,         lambda tau: self._get_QCout_On(tau),   lambda tau: self._get_QCout_Off(tau),     0] )
                                return np.fmax(result, 0.0)
                        def get_VCout_Tau(self, tau):
                                return self.get_QCout_Tau(tau)/self.Cout
                        def get_QCout_TauOn(self):
                                return self.get_QCout_Tau(self.TauOn)
                        def get_VCout_TauOn(self):
                                return self.get_QCout_Tau(self.TauOn)/self.Cout

                        def get_PCout_Lx(self):
                                return -self.parentCL.IL0.get_DE_VCout(self.TauOn, 1.0, self.get_VCout_TauOn())/self.parentCL.tpwm
                                #return -self.parentCL.IL0.get_DE_VCout(self.TauOn, 1.0, self._get_QCout_On(self.TauOn)/self.Cout)/self.parentCL.tpwm

                        def get_VCoutAvg(self):
                                return (self.get_VCout_Tau(0.0)+self.get_VCout_Tau(self.TauOn)+self.get_VCout_Tau(1.0))/3.0

                        def add_QCout_dE(self, dE):
                                # QCout_new = sqrt( QCout**2 + 2*DECout*Cout )
                                qnew2=self.QCout**2 + 2*dE*self.Cout
                                qnew2=max(0.0, qnew2)
                                self.set_QCout(self.QCout + sqrt(qnew2))

                        def set_QCout(self, qcout):
                                self.QCout=max(0.0, qcout)

                        def get_VCout(self):
                                return self.QCout/self.Cout

                        def save_QCout_Next(self):
                                self.QCout_Next=self.get_QCout_Tau(1.0)

                        def apply_QCout_Next(self):
                                de_ilx_vcout = -(1.0+float(self.parentCL.parent.EnablePh1))*self.parentCL.IL0.get_DE_VCout(self.TauOn, 1.0, self._get_QCout_On(self.TauOn)/self.Cout)
                                de_cout_rload = self.get_DE_Rload_dT(self.get_QCout_TauOn(), self.parentCL.tpwm*(1.0-self.TauOn))
                                de_cout_net = de_ilx_vcout + de_cout_rload
                                self.parentCL.parent.PM0.add_E_T("de_ilx_vcout", de_ilx_vcout, self.parentCL.tpwm)
                                self.parentCL.parent.PM0.add_E_T("de_cout_rload", de_cout_rload, self.parentCL.tpwm*(1.0-self.TauOn))
                                self.parentCL.parent.PM0.add_E_T("de_cout_net", de_cout_net, self.parentCL.tpwm)
                                dq = self.get_DQCout_Rload_dT(self.QCout, self.parentCL.tpwm*self.TauOn)
                                deon = self.get_DE_Rload_dQ(self.QCout, self.QCout+dq)
                                deoff = self.get_DE_Rload_dT(self.get_QCout_TauOn(), self.parentCL.tpwm*(1.0-self.TauOn))
                                self.parentCL.parent.PM0.add_E_T("DE_QCout_Rload_On", deon, self.parentCL.tpwm*(self.TauOn))
                                self.parentCL.parent.PM0.add_E_T("DE_QCout_Rload_Off", deoff, self.parentCL.tpwm*(1.0-self.TauOn))
                                self.parentCL.parent.PM0.add_E_T("DE_QCout_Rload_Tot", deon+deoff, self.parentCL.tpwm)
                                self.set_QCout(self.QCout_Next)


                #-----------------------------------------------------------------------------------------------------------------------------------
                def stepPWM(self, n0=None, n=1, plotDatPks=False, plotDatTau=False, plotDatRT=False, debug=False, dtau=1.0/(2.0*nptsIcaox), eps=1.0e-6):
                        self.gppPlotter.open_gpp("lcDatTau")
                        self.lcDatTau=np.array([]).reshape(7,0)
                        self.ILxCurr1, self.ILxCout1, self.ILxCurrTot, self.ILxCoutTot = self.Tau, self.Tau, self.Tau, self.Tau,
                        self.lcDatTau_str= "self.Clock+self.tpwm*self.Tau, self.IL0.ILxCurrArr, self.IL0.get_ILxCout_Arr(), self.ILxCurr1, self.ILxCout1, self.ILxCurrTot, self.ILxCoutTot"
                        self.lcDatTau_lblstr=["Tau:ILxCurrArr, Tau:ILxCurr1, Tau:ILxCurrTot","Tau:get_ILxCout_Arr(), Tau:ILxCout1, Tau:ILxCoutTot"]
                        self.gppPlotter.open_gpp("lcDatRT")
                        self.lcDatRT_str= "self.Tau, self.IL0.ILxCurrArr, self.CA0.Vcsx, self.Vimo, self.CA0.Icaox, self.CA0.IIcaox.real, self.CA0.PWMx"
                        self.lcDatRT_lblstr=["Tau:ILxCurrArr","Tau:PWMx"]
                        self.gppPlotter.open_gpp("lcDatPks")
                        self.lcDatPks_str= "self.Clock,  self.IL0.get_ILx_Tau(0.0), self.IL0.get_ILx_TauOn(), self.IL0.get_ILx_Tau(1.0),"+ \
                                "self.QCOUT.get_VCout_Tau(0.0), self.QCOUT.get_VCout_TauOn(), self.QCOUT.get_VCout_Tau(1.0),"+ \
                                "self.CA0.Vcaox_prev, self.CA0.get_Vcaox(), self.CA0.get_IIcaox_Last(), self.CA0.Icaox[-1],"+ \
                                "ncnt, self.CA0.get_FnRoot_PWMxArr(), self.CA0.get_FnRoot_PWMx(self.CA0.TauRootFin), self.CA0.TauRoot, 2*self.CA0.TauRootFin, self.CA0.ErrAvg_PWMx, "#+ \
                        self.lcDatPks_lblstr=["Clock:VCout_TauOn(), Clock:VCout_Tau(1.0)", "Clock:ILx_Tau(0.0), Clock:ILx_TauOn(), Clock:ILx_Tau(1.0), Clock:TauRootFin"]

                        if n0 is not None:
                                self.cntTPWM=n0
                        if n0==0:
                                self.parent.PM0.__init__()
                        for ncnt in range(self.cntTPWM, self.cntTPWM+n):
                                self.CA0.reset()
                                nPh0 = (ncnt%self.npwmMax)*self.ntonPh  #this + delta must be <= Pi; recall pain around integral of fabs(sin) at Pi
                                Vinac=self.parent.get_Vinac_Sin(nPh0 + self.ntonPh*self.Tau)
                                self.Vimo=self.parent.get_Vimo(self.parent.get_Iimo(self.parent.get_Vvinac(Vinac), self.parent.VL0.Vvao))    #assuming Vvao fluctuations small compared to Vvinac over tpwm
                                if (self.CA0.get_Vcaox()>self.Vramp_min) and (not self.parent.SkipHalf):         #begins with state On
                                        if debug:               print("  ON, trying Dmax for root find; ")
                                        self.CA0.TauRoot=self.parent.Dmax
                                else:                                           #off entire cycle
                                        if debug:               print("  OFF entire cycle. Accepting 0 solution as final... ILx:", self.IL0.ILx, self.QCOUT.QCout)
                                        self.CA0.TauRoot = 0.0  #WOW USING INTEGER 0 HERE BROKE IT IN SNEAKIEST SHITTY WAY, not impressed, python FAHHHHKKK
                                        self.CA0.solnFound = True
                                self.CA0.configure(nPh0, self.CA0.TauRoot, flagSavePrev=False)  #False:not committed to this integrand yet, don't save in history
                                if not self.CA0.solnFound:
                                        if self.CA0.has_Root_PWMx():
                                                self.CA0.TauRoot, self.CA0.PeakLimit = self.CA0.get_TauRoot_Any()
                                                if self.CA0.PeakLimit:
                                                        #print("PEAK LIMIT:", ncnt, self.CA0.TauRoot)
                                                        self.CA0.solnFound=True
                                                        self.CA0.TauRootFin=self.CA0.TauRoot
                                                elif (self.CA0.TauRoot < self.parent.Dmax) and (self.CA0.TauRoot > 0.0):        #Good Root Found, refine and calc actual solution
                                                        self.CA0.configure(nPh0, self.CA0.TauRoot, flagSavePrev=False)
                                                        self.CA0.TauRootFin = self.CA0.get_TauRoot_PWMx(self.CA0.TauRoot-5*dtau, self.CA0.TauRoot+5*dtau)
                                                        if ( fabs(self.CA0.get_FnRoot_PWMx(self.CA0.TauRootFin)) < min(fabs(self.CA0.get_FnRoot_PWMxArr()),eps) ) and ( fabs(fabs(self.CA0.TauRootFin) - fabs(self.CA0.TauRoot)) < dtau ):     #accept new tau_root
                                                                pwmx_root=fabs(self.CA0.get_FnRoot_PWMx(self.CA0.TauRootFin))
                                                        else:
                                                                print("ERROR?? Should Never Happen. keeping old tau_root: pwmx_root_ini={:1.6e}, pwmx_root_fin={:1.6e}, tau_root={:1.6e}, tau_root_fin={:1.6e}".format(self.CA0.get_FnRoot_PWMxArr(), self.CA0.get_FnRoot_PWMx(self.CA0.TauRootFin), self.CA0.TauRoot, self.CA0.TauRootFin))
                                                                self.CA0.TauRootFin=None
                                                                pwmx_root=fabs(self.CA0.get_FnRoot_PWMxArr())
                                                                quit()
                                                        self.CA0.ErrAvg_PWMx=(self.CA0.ErrAvg_PWMx + pwmx_root)/2.0
                                                        self.CA0.solnFound = True
                                        if not self.CA0.solnFound:
                                                if self.CA0.peakLimitTrigger():
                                                        self.CA0.TauRoot=0.0
                                                else:
                                                        self.CA0.TauRoot=self.parent.Dmax
                                                if debug or 1:
                                                        if self.CA0.TauRoot > self.parent.Dmax:         print("  OOB: Root exceeds bounds: Dmax limiting; tau_root=", ncnt, self.CA0.TauRoot, end=" \n")
                                                        else:                                           pass #print("  No Root found in PWMx: peak limited-->zero; tau_root=", ncnt, self.CA0.TauRoot, end=" \n")
                                                self.CA0.solnFound = True
                                if self.CA0.solnFound:
                                        if self.CA0.TauRootFin is None:
                                                self.CA0.TauRootFin=self.CA0.TauRoot
                                        self.CA0.configure(nPh0, self.CA0.TauRootFin, flagSavePrev=True)
                                        #self.CA0.update_Vcaox_with_IIcaox()
                                        self.StatsCL.addStats(nPh0, self.CA0.TauRootFin, np.mean(self.Vimo), self.IL0.get_ILxAvgOn(), self.IL0.get_ILxAvgOff(), np.mean(Vinac)*self.IL0.get_ILxAvgOn(), self.QCOUT.get_PCout_Lx(), self.QCOUT.get_VCoutAvg())
                                        #WHY HERE? some variables lose scope outside loop
                                        self.lcDatPks.append(eval("["+self.lcDatPks_str+"]"))
                                        if plotDatTau:
                                                self.lcDatTau=np.hstack([self.lcDatTau, np.vstack(eval("["+self.lcDatTau_str+"]"))])
                                else:
                                        print("ERROR: no root found or assumed. stopping")
                                        quit()
                                #since some are dependent, must save all new values first
                                self.CA0.save_Vcaox_Next()
                                self.IL0.save_ILx_Next()
                                self.QCOUT.save_QCout_Next()
                                #now ok to apply changes
                                self.CA0.apply_Vcaox_Next()
                                self.IL0.apply_ILx_Next()               #also updates power meter
                                self.QCOUT.apply_QCout_Next()           #also updates power meter

                                self.addClockTPWM()

                                if plotDatRT:      #realtime display, never to file
                                        #self.lcDatRT=list( eval("zip(*["+self.lcDatRT_str+"])") )
                                        self.lcDatRT=np.vstack(eval("["+self.lcDatRT_str+"]"))
                                        #self.gppPlotter.gpp_Plot(self.lcDatRT.transpose(), self.lcDatRT_str, self.lcDatRT_lblstr, self.gppPlotter.gpp["lcDatRT"])
                                        self.gppPlotter.gpp_Plot(self.lcDatRT.transpose(), self.lcDatRT_str, self.lcDatRT_lblstr, "lcDatRT")

                        if plotDatTau: #contortions needed to dup currents pi out of phase and add them
                                self.lcDatTau[3]=np.interp(self.lcDatTau[0]-self.tpwm/2.0, self.lcDatTau[0], self.lcDatTau[1])
                                self.lcDatTau[4]=np.interp(self.lcDatTau[0]-self.tpwm/2.0, self.lcDatTau[0], self.lcDatTau[2])
                                self.lcDatTau[5]=self.lcDatTau[1]+self.lcDatTau[3]
                                self.lcDatTau[6]=self.lcDatTau[2]+self.lcDatTau[4]
                                tmpname=self.parent.get_FileName("ilxTau_{}_{}.pdf".format(n, self.cntTPWM))
                                #self.gppPlotter.gpp_Plot(self.lcDatTau.transpose(), self.lcDatTau_str, self.lcDatTau_lblstr, self.gppPlotter.gpp["lcDatTau"], lineWidth=0.5, toFile=tmpname)
                                self.gppPlotter.gpp_Plot(self.lcDatTau.transpose(), self.lcDatTau_str, self.lcDatTau_lblstr, "lcDatTau", lineWidth=0.5, toFile=tmpname)

                        if plotDatPks:
                                #self.gppPlotter.gpp_Plot(self.lcDatPks, self.lcDatPks_str, self.lcDatPks_lblstr, self.gppPlotter.gpp["lcDatPks"], lineWidth=1.0, toFile=self.parent.get_FileName("lcDatPks.pdf"))
                                self.gppPlotter.gpp_Plot(self.lcDatPks, self.lcDatPks_str, self.lcDatPks_lblstr, "lcDatPks", lineWidth=1.0, toFile=self.parent.get_FileName("lcDatPks.pdf"))

                        self.parent.PM0.Report()
                        return self.lcDatPks[-1]


                def get_ZType2_CAOx(self, bk, f):
                        return self.parent.get_Zcomp_type2(bk, f, self.parent.Rzc, self.parent.Czc, self.parent.Cpc)

                def get_ZType2_CAOx_DC(self):
                        # should be integration coefficient, simply 1/(cz+cp), to calc DC component integral approx over dt: dv = dq/Ctot ~= dt*DCcomponentOfCurrent/(cp+cz)
                        return self.parent.get_Zcomp_type2(1, 1/(2*pi*1j), 0, self.parent.Czc, self.parent.Cpc)


        #########################################################################################################
        # FILTERS : Complex Impedance, Bode Plots
        def get_Zcomp_type2(self, bk, f, rz, cz, cp):
                omega = 2*pi*f
                zz = rz + 1.0/(omega*cz*1j)
                return bk / ( (omega*cp*1j) + (1.0/zz) )

        def get_Z_lst(self, f, res):
                return [f, res, abs(res), atan(res.imag/res.real)]

        def get_impedance_div2R2C(self, omega, ra, rb, ca, cb):
                tmpa = 1.0 / ( (omega*ca*1j) + (1.0/ra) )
                tmpb = 1.0 / ( (omega*cb*1j) + (1.0/rb) )
                res = tmpb / (tmpa + tmpb)
                return [omega, res, abs(res), atan2(res.imag, res.real)]


##################################
# main ###########################
##################################
def main():
#        if(len(sys.argv)>1):
#        l=sys.argv.copy()
#        l.remove(l[0])
        pfc0=UCC28070(plotpdf=True)

        def formatCustom(a):
                if type(a) is float:
                        return "{:20.6}".format(a)
                elif (type(a) is int) or (type(a) is bool):
                        return "{:20}".format(a)
                else:
                        return "{}".format(a)
        with open("logs/log_Vac{}_Ptot{}_Lx{}_Fpwm{}_RBac{}_Cout{}.dat".format(\
                pfc0.Vac_rms, pfc0.Ptot, round(pfc0.Lx*1.0e6), pfc0.Fpwm, round(pfc0.RBac), round(pfc0.Cout*1.0e6) ), 'w') as logfile:
                for k in pfc0.__dict__.keys():
                        #print(k,":",pfc0.__dict__[k], file=logfile)
                        print("{:20} : {}".format(k, formatCustom(pfc0.__dict__[k])), file=logfile)
                print("ICout at Pmax:", pfc0.get_ICout(sqrt(2)*pfc0.Vac_rms, 2*pfc0.Ptot), "Vripple:", pfc0.get_Vripple(), file=logfile)
                print("Vvinac max:", pfc0.get_Vvinac(sqrt(2)*pfc0.Vac_rms), "kvff:", pfc0._get_Kvff(pfc0.get_Vvinac(sqrt(2)*pfc0.Vac_rms)), file=logfile)
                print("get_Vinac_Correction_Factor():", pfc0.get_Vinac_Correction_Factor(), file=logfile)
                print("IMO_max:", pfc0.get_IMO_max(), "Vimo_max:", pfc0.Rimo*pfc0.get_IMO_max(), file=logfile)
                print("Duty at top:", pfc0.get_Duty(sqrt(2)*pfc0.Vac_rms), " at peak voltage:", sqrt(2)*pfc0.Vac_rms, file=logfile)
                print("ILxTop:", pfc0.get_ILxTop(pfc0.Vac_rms,pfc0.Ptot/2), "dILxTop:", pfc0.get_dILxTop(pfc0.Vac_rms), "ILxPeak:", pfc0.get_ILxPeak(), file=logfile)
                print("dILxTop/ILxTop:", pfc0.get_dILxTop(pfc0.Vac_rms)/pfc0.get_ILxTop(pfc0.Vac_rms,pfc0.Ptot/2), file=logfile)
                print("Rsequiv:", pfc0.Rsequiv, " --> need Vreg_currpeak >= ", pfc0.Vreg_currsense*pfc0.get_ILxPeak()/pfc0.get_ILxTop(pfc0.Vac_rms, pfc0.Ptot/2), file=logfile)
                print("LxMin_dILx:", pfc0.get_LxMin_dILx(pfc0.Vac_rms, pfc0.Ptot/2, pfc0.Fpwm, 0.3), "for fracilx=0.3", file=logfile)
                print("LxMin_CCM at Pmax/4:", pfc0.get_LxMin_CCM(pfc0.Vac_rms, pfc0.Ptot/(2*4), pfc0.Fpwm), file=logfile)
                print("LxMin_CCM at Pmax/3:", pfc0.get_LxMin_CCM(pfc0.Vac_rms, pfc0.Ptot/(2*3), pfc0.Fpwm), file=logfile)
                print("PoutMin_CCM:", pfc0.get_PoutMin_CCM(pfc0.Vac_rms, pfc0.Fpwm), "with Lx:", pfc0.Lx, file=logfile)
                #print("_get_ILx_Vals at n=833:", pfc0._get_ILx_Vals(833, pfc0.Ptot/2.0), file=logfile)
        #quit()

        pdiv=1
        pfc0.Rload=225.0*pdiv
        pfc0.EnablePh1=True
        pfc0.SkipHalfMod=1
        pfc0.Css=pfc0.get_Css(0.3)
        pfc0.CL0.QCOUT.QCout=340*pfc0.CL0.QCOUT.Cout
        #pfc0.CL0.QCOUT.QCout=600*pfc0.CL0.QCOUT.Cout
        pfc0.VL0.QCss=1.8*pfc0.Css
        #pfc0.VL0.QCss=4.7*pfc0.Css
        pfc0.VL0.Vvao=1.0 + (float(pfc0.SkipHalfMod))*pfc0.DVvao/pdiv
        pfc0.CL0.CA0.Vcaox=4.6/3
        pfc0.CL0.IL0.ILx=0.0

        iIn, iFin = 0, round(1667/2)
        if 0:
                iStep=33
                for ii in [i for i in range(iIn, int(iFin), iStep)]:
                        print("ii:",ii)
                        pfc0.CL0.stepPWM(n=iStep, plotDatPks=0, plotDatTau=1, debug=False)
                quit()

        if 1:
                pfc0.VL0.stepTHALF(60, useCurrentLoop=True)
                quit()

        if 0:
                pfc0.CL0.stepPWM(iIn, iFin-iIn, plotDatPks=1, plotDatTau=1, plotDatRT=0, debug=False)
                quit()

        pfc0.plotBode_VoltageLoop()
        pfc0.plotBode_CurrentLoop()


##################################
##################################
##################################

if __name__ == '__main__':
    main()






######################################################
#       Voltage Compensation Network Gain: Gcv(f):
#
#series Rz Cz:
#               Zser = Rz + 1/(j*2*pi*f*Cz)
#               Zser = ( Rz*j*2*pi*f*Cz + 1 )/(j*2*pi*f*Cz)
#parallel Zs Cp:
#               Zpar = 1 / ( 1/Zser + (j*2*pi*f*Cp) )
#               Zpar = 1 / ( (j*2*pi*f*Cz) / ( Rz*j*2*pi*f*Cz + 1 ) + (j*2*pi*f*Cp) )
#               Zpar = ( Rz*j*2*pi*f*Cz + 1 ) / ( (j*2*pi*f*Cz)  + (j*2*pi*f*Cp) * ( Rz*j*2*pi*f*Cz + 1 ) )
#               Zpar = ( Rz*j*2*pi*f*Cz + 1 ) / ( (j*2*pi*f*Cz) + (j*2*pi*f*Cp) + (j*2*pi*f*Cp)*(Rz*j*2*pi*f*Cz) )
#               Zpar = ( Rz*j*2*pi*f*Cz + 1 ) / ( (j*2*pi*f)*(Cz+Cp) + (j*2*pi*f)*(Cz+Cp)*(j*2*pi*f*Cz*Cp*Rz)/(Cz+Cp) )
#               Zpar = ( Rz*j*2*pi*f*Cz + 1 ) / ( (j*2*pi*f)*(Cz+Cp) * ( 1 + (j*2*pi*f*Cz*Cp*Rz)/(Cz+Cp) )
#
# Gcv(f) = dVvao/dVout = Kdivvs * Gmv * Zpar
#
######################################################
#       Voltage Loop Power Stage Gain: Gpsv(f):
#
# Gpsv(f) = (Ptot/Vout)*(1/(j*2*pi*f*Cout))/DVvao
#
######################################################
#
#
#
#
#
########################################################################################################
########################################################################################################
# ENERGY BASED APPROACH TO FIND FAST PWM TIMESCALE INDUCTOR CURRENT VALUES
# control strategy: average current must follow
# Iin(t)=A*sin(omega*t) where A=sqrt(2)*(Prms/Vrms), omega*t in [0,pi]
# the average current must track the sin wave,
#
########################################################################################################
########################################################################################################
# when general, v(t) is A*sin(w*t) - VCout, so
# integral di/dt = integral -(A*sin(w*t)-VCout)/L
#       --> i(t) = -C + ((A/w)*cos(w*t)+VCout*t)/L
# but to satisfy initial condition I(0) = -I0,
#       i(0) = -I0 = -C + ((A/w)*cos(ph0)+VCout*(0)/L) ; assuming VCout off until t=0;
#               -I0 = -C + ((A/w)*cos(ph0))/L
#               -C = -I0 - ((A/w)*cos(ph0))/L
# at time t later:
#       i(t) = -C + ((A/w)*cos(ph0+w*t)+VCout*t)/L
#       i(t) = -I0 - ((A/w)*cos(ph0))/L + ((A/w)*cos(ph0+w*t) + VCout*t)/L
######################################################################################################## general form, exact for all t:
#       i(t) = -I0 - (A/(w*L))*(cos(ph0) - cos(ph0+w*t)) + VCout*t/L
####################################################################################################################################### 
#       SIMPLEST: for t<<1/Wac
#       i(t) ~= -I0 - ( Vinac - VCout )*t/L
#               Vinac ~= A*sin(ph0 + dph/2);  dph=tpwm*TauRoot  --> mid interval voltage approx
#               VCout ~= QCout/Cout
####################################################################################################################################### 
#       since cos is decreasing, the current is more negative (into the terminal) from positive ac voltage applied, OK.
#
# interpreting signs: the current from the positive AC terminal is INTO the terminal, hence negative;
#       this current flows through the inductor, out of the VCout terminal, doing POSITIVE work on Cout. OK.
#       the positive VCout on the opposite side changes the current in the opposite direction, OK.
# at left terminal, current flowing in against the positive potential does positive work on the inductor,
# at right terminal, current flowing out with the positive potential does negative work on the inductor (doing positive work on capacitor).
#
#------------------------------------------------------------------------------------------------------------------------------
# NOTE CONVENTION: in software, sign convention is opposite of above, ILx is POSITIVE always.
#       +VCout with positive current inductor-->cap does positive work
#       consistent with positive AC voltage applied to opposite side of inductor, increasing current in positive direction
#------------------------------------------------------------------------------------------------------------------------------
#
# Integrate I*V to get energy - simple right? just watch
# expand current to convenient form,
#       i(t) =  -C0 + K*cos(ph0+w*t) + VCout*t/L  ;
#                                                        K  = (A/(w*L))
#                                                       -C0 = -I0 - K*cos(ph0)
# consider work integrals on left,right separately:
# v_left(t) = A*sin(ph0+w*t)           v_right(t)= VCout
# LEFT SIDE: positive work done on inductor --> E(t) = - integral v_left(t) * i(t)
# E(t) = -integral ( A*sin(ph0+w*t) ) * ( -C0 + K*cos(ph0+w*t) + VCout*t/L  )dt
# E(t) = -integral ( -A*sin(ph0+w*t)*C0 + A*sin(ph0+w*t)*K*cos(ph0+w*t) + A*sin(ph0+w*t)*VCout*t/L )dt
# E(t) = - ( (A/w)*cos(ph0+w*t)*C0 - (1/2)*(A/w)*cos**2(ph0+w*t)*K + ( integral t*sin(ph0+w*t)*A*VCout/L dt ) ) ;       #TRICKY, cos**2 not -sin**2; integration constants gonna getchya
#       by parts: integral t*sin(ph0+w*t) dt = -(1/w)*t*cos(ph0+w*t) - integral -(1/w)*cos(ph0+w*t)*1 dt
#       by parts: integral t*sin(ph0+w*t) dt = -(1/w)*t*cos(ph0+w*t) + (1/w**2)*sin(ph0+w*t)
# E(t) = - ( (A/w)*cos(ph0+w*t)*C0 - (1/2)*(A/w)*cos**2(ph0+w*t)*K + (A*VCout/L)*( -(1/w)*t*cos(ph0+w*t) + (1/w**2)*sin(ph0+w*t) ) )
# E(t) = -(A/w)*cos(ph0+w*t)*C0 + (1/2)*(A/w)*cos**2(ph0+w*t)*K + VCout*(A/(w*L))*t*cos(ph0+w*t) - VCout*(1/w)*(A/(w*L))*sin(ph0+w*t)
# E(t)-E(0) = -(A/w)*cos(ph0+w*t)*C0 + (1/2)*(A/w)*cos**2(ph0+w*t)*K + VCout*K*t*cos(ph0+w*t) - VCout*(1/w)*K*sin(ph0+w*t)                       #MUST EVAL E(t)-E(0);
# E(t)-E(0) = (A/w)*C(t)*C0 + (1/2)*(A/w)*(cos**2(ph0+w*t)-cos**2(ph0))*K - VCout*K*t*C(t) - VCout*(1/w)*K*S(t)                    # C(t) = -(cos(ph0+w*t)-cos(ph0)) > 0,   S(t) = (sin(ph0+w*t)-sin(ph0))
# E(t)-E(0) = (A/w)*C(t)*(I0+K*cos(ph0)) + (1/2)*(A/w)*(cos**2(ph0+w*t)-cos**2(ph0))*K - VCout*K*t*C(t) - VCout*(1/w)*K*S(t)            # -C0 = -I0 - K*cos(ph0)
# E(t)-E(0) = (A/w)*C(t)*I0 - (A/w)*(cos(ph0+w*t)-cos(ph0))*K*cos(ph0) + (1/2)*(A/w)*(cos**2(ph0+w*t)-cos**2(ph0))*K - VCout*K*t*C(t) - VCout*(1/w)*K*S(t)
# E(t)-E(0) = (A/w)*C(t)*I0 - (A/w)*(cos(ph0+w*t)*cos(ph0) - cos**2(ph0))*K + (1/2)*(A/w)*(cos**2(ph0+w*t) - cos**2(ph0))*K - VCout*K*t*C(t) - VCout*(1/w)*K*S(t)
# E(t)-E(0) = (A/w)*C(t)*I0 - (A/w)*(cos(ph0+w*t)*cos(ph0))*K + (1/2)*(A/w)*(cos**2(ph0+w*t) + cos**2(ph0))*K - VCout*K*t*C(t) - VCout*(1/w)*K*S(t)
# E(t)-E(0) = (A/w)*C(t)*I0 + (1/2)*K*(A/w)*( -2*cos(ph0+w*t)*cos(ph0) + cos**2(ph0+w*t) + cos**2(ph0) ) - VCout*K*t*C(t) - VCout*(1/w)*K*S(t)
# E(t)-E(0) = (A/w)*C(t)*I0 + (1/2)*K*(A/w)*(cos(ph0) - cos(ph0+w*t))**2 - VCout*K*t*C(t) - VCout*(1/w)*K*S(t)
########################################################################################################
# WORK ON INDUCTOR FROM VINAC  (these need review - not sure if individual term signs are correct?)
########################################################################################################
# E(t)-E(0) = (A/w)*C(t)*I0 + (1/2)*K*(A/w)*(C(t))**2 - VCout*K*t*C(t) - VCout*(K/w)*S(t)       # Don't forget there are 2 inductors!
#                       C(t) = -(cos(ph0+w*t)-cos(ph0)) > 0
#                       S(t) = (sin(ph0+w*t)-sin(ph0))
#                       K  = (A/(w*L))
#                       -C0 = -I0 - K*cos(ph0)                          sets initial current to I0 < 0
#                       i(t) =  -C0 + K*cos(ph0+w*t) + VCout*t/L        assuming VCout=0 when t<0; since cos decreasing, i getting more negative as desired
#######################################################################################################################################
#       SIMPLEST: for t<<1/Wac
#       E(t)-E(0) ~= Vinac*( I0*t + (1/2)*(Vinac - VCout)*t**2/L )
#               Vinac ~= A*sin(ph0 + dph/2);  dph=tpwm*TauRoot  --> mid interval voltage approx
#               VCout ~= QCout/Cout
####################################################################################################################################### 
# E(t)-E(0) looks sensible, positive definite when VCout=0, tends to 0 as t-->0, OK.
# continue,
# RIGHT SIDE: negative work done by (by definition negative) current flowing out of inductor into +VCout --> E(t) = integral v_right(t) * i(t)
# with v_right(t) = VCout when transistor is off (zero otherwise)
# E(t) = integral ( VCout ) * ( -C0 + K*cos(ph0+w*t) + VCout*t/L  )dt
# E(t) = VCout * integral ( -C0 + K*cos(ph0+w*t) + VCout*t/L  )dt
# E(t) = VCout * ( -C0*t + (1/w)*K*sin(ph0+w*t) + (1/2)*VCout*t**2/L )                 #MUST EVAL E(t)-E(0);
# E(t)-E(0) = VCout * ( -C0*t + (1/w)*K*S(t) + (1/2)*VCout*(t**2)/L )
# E(t)-E(0) = -VCout*C0*t + VCout*(K/w)*S(t) + (1/2)*VCout**2*(t**2)/L
########################################################################################################
# WORK ON INDUCTOR FROM VCOUT
########################################################################################################
# E(t)-E(0) = -VCout*(I0 + K*cos(ph0))*t + VCout*(K/w)*S(t) + (1/2)*VCout**2*(t**2)/L
#######################################################################################################################################
#       SIMPLEST: for t<<1/Wac
#       E(t)-E(0) ~= -VCout*( I0*t + (1/2)*(Vinac - VCout)*t**2/L )
#               Vinac ~= A*sin(ph0 + dph/2);  dph=tpwm*TauRoot  --> mid interval voltage approx
#               VCout ~= QCout/Cout
####################################################################################################################################### 
# Note this minus the energy consumed by the load is the energy left to charge the capacitor.
#       negative in initial current, good.
#       linear decrease in current yields +quadratic term.
#       work done by Vinac also adds a +S(t) term;
#               note this term cancels when two work expressions are added to get total work on inductor
#
# WORK ON CAPACITOR FROM INDUCTOR
#       is simply the negative of the previous expression MINUS the load energy demand during that period
#
# WORK ON LOAD FROM VCOUT
#       while the inductor is grounded with Vinac working on it,
#       the load is simply connected to VCout and during this time can be simply modelled as
#       an exponential decay solution with time constant Rload*Cout, since for Cout
# dQ/dt = - V(t)/Rload ; V(t) =  Q(t)/Cout
# dQ/dt = - Q(t)/(Rload*Cout) -->
# Q(t) = Q0*exp(-t/(Rload*Cout)) ; Q0 = Q(0) = V(0)*Cout
########################################################################################################
# Q(t)-Q(0) = Q0*( exp(-t/(Rload*Cout)) - 1 )
########################################################################################################
#       SIMPLEST: no need to change the following, already good.
########################################################################################################
# and the work done by the cap on the load during this time changes the cap energy as
# E(t)-E(0) = (1/2)*( Q(t)**2 - Q(0)**2 )/C
#
# While the inductor is connected to the cap, quickly dumping energy into it,
# the load continues to demand energy at a rate
# Pload = V(t)**2/Rload
# diverting some of the energy of the inductor away from the capacitor.
##### Since the inductor current switching time is about 1e5 times faster than Cout/Rload timescale,
##### it may be ok to assume constant voltage and do energy jump calculations,
##### rather than develop full complicated power flow diffyQ calculations. so,
# for a known time t, we can calculate the energy demand as
# Eload = Pload * t = t*V(t)**2/Rload
# and knowing how much energy ELx the inductor is transferring from above, we find
# the energy left for the cap as
# DECout = ELx - Eload
# DECout(t) = ELx - t*V(t)**2/Rload
#       which may be positive or negative.  If V(t) is changing only slightly, we can use V(0) = Q(0)/Cout,
######################################################## Don't forget 2 inductors contribute!
# DECout(t) = ELx - t*V(0)**2/Rload
########################################################
#       to update QCout via ECout_new = ECout + DECout,
#       (1/2)*QCout_new**2/Cout = (1/2)*QCout**2/Cout + DECout
#       QCout_new**2 = QCout**2 + 2*DECout*Cout
########################################################
# QCout_new = sqrt( QCout**2 + 2*DECout*Cout )
########################################################
# we can use the new Q to revise the V(t) we used, or just let it be an artificial inefficiency.
# we will be slightly underestimating the voltage and the energy increase during inductor dump.
# approx 1e-5 fractional error in charge, 10s of millivolts or less.
#
########################################################################################################
########################################################################################################








# transconductance weirdness
# like a resistor in series, with current defined by voltage:
# iout = C * vin ; C units 1/r
# r(vin) = vin/iout
# but if iout into non zero va,
# r(vin) = (vin-va)/iout = (vin-va)/(C*vin)
# r(vin) = (1-(va/vin))/C
# which becomes like negative resistance for va>vin

#-----------------------------------------------------------------------------------------------------------------------------------------------


