# Buck converter, simulation, exact, phase plane
#
# Vin = 10 V
# D = 0.5, fixed
# L = 10 uH
# C = 10 uF
# Iout = 2 A or Iout = 5 A to simulate transients
# Ts = 10 us, fs = 100 kHz
# ideal switch, ideal diode



from pylab import *



def phaseplane(periods, points_per_period, ic_mc, ic_jl, Jout, w0Ts):

    n = periods * points_per_period
    mc = empty(n + 1)
    jl = empty(n + 1)
    per = linspace(0, periods, n + 1)
    delta_phi = w0Ts / points_per_period
    perpart = modf(per)[0]
    s = perpart < 0.5
    mc[0] = ic_mc
    jl[0] = ic_jl
    r = sqrt((ic_mc - 1)**2 + (ic_jl - Jout)**2)
    p = arctan2(ic_jl - Jout, ic_mc - 1)
    for i in range(1, n + 1):
        if s[i - 1] == True and s[i] == True:
            p -= delta_phi
            mc[i] = 1 + r * cos(p)
            jl[i] = Jout + r * sin(p)
        elif s[i - 1] == False and s[i] == False:
            p -= delta_phi
            mc[i] = 0 + r * cos(p)
            jl[i] = Jout + r * sin(p)
        elif s[i - 1] == False and s[i] == True:
            r = sqrt((mc[i - 1] - 1)**2 + (jl[i - 1] - Jout)**2)
            p = arctan2(jl[i - 1] - Jout, mc[i - 1] - 1)
            p -= delta_phi
            mc[i] = 1 + r * cos(p)
            jl[i] = Jout + r * sin(p)
        elif s[i - 1] == True and s[i] == False:
            r = sqrt((mc[i - 1] - 0)**2 + (jl[i - 1] - Jout)**2)
            p = arctan2(jl[i - 1] - Jout, mc[i - 1] - 0)
            p -= delta_phi
            mc[i] = 0 + r * cos(p)
            jl[i] = Jout + r * sin(p)
        else:
            print('this is not possible!')
    return (mc, jl, per, s)



def averagedphaseplane(periods, points_per_period, ic_mc, ic_jl, Jout, w0Ts):

    n = periods * points_per_period
    mc = empty(n + 1)
    jl = empty(n + 1)
    per = linspace(0, periods, n + 1)
    delta_phi = w0Ts / points_per_period
    mc[0] = ic_mc
    jl[0] = ic_jl
    r = sqrt((ic_mc - 0.5)**2 + (ic_jl - Jout)**2)
    p = arctan2(ic_jl - Jout, ic_mc - 0.5)
    for i in range(1, n + 1):
        p -= delta_phi
        mc[i] = 0.5 + r * cos(p)
        jl[i] = Jout + r * sin(p)
    return (mc, jl, per)
    

   
close('all')
 
# 100 points per period, LRA initial condition
ic_mc = 0.5
delta_jl = 0.125 # by linear ripple approximation
Jout = 0.2
ic_jl = Jout - delta_jl
w0Ts = 1
mc, jl, per, s = phaseplane(10, 100, ic_mc, ic_jl, Jout, w0Ts)

figure(1, figsize = (6,6))
plot(mc, jl, label = 'LRA')
plot(1, Jout, 'r.')
plot(0, Jout, 'b.')
xlim(-0.1, 1.1)
ylim(-0.1, 1.1)
xlabel('mc')
ylabel('jl')
legend()
show()

input('press Enter to continue . . . ')

# 100 points per period, exact initial condition
ic_mc = 0.5
delta_jl = 0.5 * tan(0.25) # by the exact solution
Jout = 0.2
ic_jl = Jout - delta_jl
w0Ts = 1
mc, jl, per, s = phaseplane(10, 100, ic_mc, ic_jl, Jout, w0Ts)
plot(mc, jl, label = 'exact')
xlim(-0.1, 1.1)
ylim(-0.1, 1.1)
xlabel('mc')
ylabel('jl')
legend()
show()

input('press Enter to continue . . . ')



close('all')
 
# 10000 points per period, LRA initial condition
ic_mc = 0.5
delta_jl = 0.125 # by linear ripple approximation
Jout = 0.2
ic_jl = Jout - delta_jl
w0Ts = 1
mc, jl, per, s = phaseplane(10, 10000, ic_mc, ic_jl, Jout, w0Ts)

figure(1, figsize = (6,6))
plot(mc, jl, label = 'LRA')
plot(1, Jout, 'r.')
plot(0, Jout, 'b.')
xlim(-0.1, 1.1)
ylim(-0.1, 1.1)
xlabel('mc')
ylabel('jl')
legend(loc = 'upper right')
show()

input('press Enter to continue . . . ')

# 10000 points per period, exact initial condition
ic_mc = 0.5
delta_jl = 0.5 * tan(0.25) # by the exact solution
Jout = 0.2
ic_jl = Jout - delta_jl
w0Ts = 1
mc, jl, per, s = phaseplane(10, 10000, ic_mc, ic_jl, Jout, w0Ts)
plot(mc, jl, label = 'exact')
xlim(-0.1, 1.1)
ylim(-0.1, 1.1)
legend(loc = 'upper right')
show()

input('press Enter to continue . . . ')

# averaged model
ic_mc = 0.5
Jout = 0.2
ic_jl = Jout
w0Ts = 1
mcav, jlav, perav = averagedphaseplane(10, 10000, ic_mc, ic_jl, Jout, w0Ts)
plot(mcav, jlav, '.', label = 'avg')
xlim(-0.1, 1.1)
ylim(-0.1, 1.1)
xlabel('mc')
ylabel('jl')
legend(loc = 'upper right')
show()

input('press Enter to continue . . . ')

close('all')
figure(1)
plot(per, jl, label = 'jl')
plot(per, mc, label = 'mc')
legend(loc = 'lower left')
xlabel('t/ Ts')
ylabel('mc, jl')
xlim(0, 10)
ylim(-0.1, 0.6)
show()

input('press Enter to continue . . . ')

plot(perav, jlav, label = 'jlav')
plot(perav, mcav, label = 'mcav')
legend(loc = 'lower left')
xlabel('t/ Ts')
ylabel('mc, jl')
xlim(0, 10)
ylim(-0.1, 0.6)
show()

input('press Enter to continue . . . ')



close('all')
 
# 100 points per period, LRA initial condition
ic_mc = 0.5
delta_jl = 0.125 # by linear ripple approximation
Jout = 0.5
ic_jl = 0.2 - delta_jl
w0Ts = 1
mc, jl, per, s = phaseplane(100, 10000, ic_mc, ic_jl, Jout, w0Ts)

figure(1, figsize = (6,6))
plot(mc, jl, label = 'swm')
plot(1, Jout, 'r.')
plot(0, Jout, 'b.')
xlim(-0.1, 1.1)
ylim(-0.1, 1.1)
xlabel('mc')
ylabel('jl')
legend(loc = 'upper right')
show()

input('press Enter to continue . . . ')

# averaged model
ic_mc = 0.5
Jout = 0.5
ic_jl = 0.2
w0Ts = 1
mcav, jlav, perav = averagedphaseplane(100, 10000, ic_mc, ic_jl, Jout, w0Ts)
plot(mcav, jlav, label = 'avg')
xlim(-0.1, 1.1)
ylim(-0.1, 1.1)
xlabel('mc')
ylabel('jl')
legend(loc = 'upper right')
show()

input('press Enter to continue . . . ')

close('all')
figure(1)
plot(per, jl, label = 'jl_sw')
plot(per, mc, label = 'mc_sw')
legend(loc = 'lower left')
xlabel('t/ Ts')
ylabel('mc, jl')
xlim(0, 100)
ylim(-0.1, 1.1)
show()

input('press Enter to continue . . . ')

figure(2)
plot(perav, jlav, label = 'jl_av')
plot(perav, mcav, label = 'mc_av')
legend(loc = 'lower left')
xlabel('t/ Ts')
ylabel('mc, jl')
xlim(0, 100)
ylim(-0.1, 1.1)
show()



input('press Enter to continue . . . ')

close('all')