from matplotlib.pyplot import *
from numpy import *

f, Vpeak, cycles = 50., 1., 1.5
fs, Tlim = 2.*f, cycles/f
Vavg, Vrms = Vpeak*2./pi, Vpeak/sqrt(2.)
t = linspace(0,Tlim,fs*cycles)
w = 2.*pi*f # 50Hz AC
signal = lambda x,p: sin(w*x+p*2.*pi/3.)
def halfWave(time):
        s1, s2, s3 = signal(time,0.), signal(time,1.), signal(time,2.)
        if s1 > s2 and s1 > s3:
                if s2 > s3:
                        return s1, s2
                else:
                        return s1, s3
        elif s2 > s1 and s2 > s3:
                if s1 > s3:
                        return s2, s1
                else:
                        return s2, s3
        else:
                if s1 > s2:
                        return s3, s1
                else:
                        return s3, s2
 
def fullWave(time):
        s1, s2, s3 = abs(signal(time,0.)), abs(signal(time,1.)), \
                abs(signal(time,2.))
        if s1 > s2 and s1 > s3:
                if s2 > s3:
                        return s1, s2
                else:
                        return s1, s3
        elif s2 > s1 and s2 > s3:
                if s1 > s3:
                        return s2, s1
                else:
                        return s2, s3
        else:
                if s1 > s2:
                        return s3, s1
                else:
                        return s3, s2

xTickPts = []
for time in t:
        s1, s2, s3 = abs(signal(time,0.)), abs(signal(time,1.)), \
                abs(signal(time,2.))
        if s1 == s2:
                xTickPts = append(xTickPts, time)
                print time
        elif s2 == s3:
                xTickPts = append(xTickPts, time)
                print time
        elif s3 == s1:
                xTickPts = append(xTickPts, time)
                print time
 
def myAxes(this):
        this.grid(True)
        this.set_xlim(0,Tlim)
        this.set_xticks(arange(0,cycles+.25,.25)/f)
        this.set_xticklabels([])
        this.set_ylabel(r"Voltage (V)",fontsize=12)
	this.set_ylim(-2.*Vpeak-.1,2.*Vpeak+.1)
	this.set_yticks([-1.73*Vpeak,-Vpeak,0,Vpeak,1.73*Vpeak])
	this.set_yticklabels([r"$-\sqrt{3}V_{\mathrm{peak} }$",r"$-V_{\mathrm{peak} }$",\
        	r"0",r"$V_{\mathrm{peak} }$",r"$\sqrt{3}V_{\mathrm{peak} }$"])

fig = figure(figsize=(7,12))
ax = fig.add_subplot(311)
ax.plot(t,signal(t,0),'b',linewidth=2,label=r"$\phi=0^\circ$")
ax.plot(t,signal(t,1),'r',linewidth=2,label=r"$\phi=120^\circ$")
ax.plot(t,signal(t,2),'g',linewidth=2,label=r"$\phi=240^\circ$")
myAxes(ax)
ax.set_title(r'3-Phase signals',fontsize=12)
ax.legend(loc=1, \
        bbox_to_anchor=(.8,.35),\
        frameon=False,handletextpad=.05)

ax = fig.add_subplot(312)
S, H = [], []
for time in t:
        s, h = halfWave(time)
        S = append(S,s)
        H = append(H,h)
ax.plot(t,S,'k',linewidth=2.)
ax.plot(t,signal(t,0),'b--',linewidth=1.)
ax.plot(t,signal(t,1),'r--',linewidth=1.)
ax.plot(t,signal(t,2),'g--',linewidth=1.)
myAxes(ax)
ax.set_title(r"Half-wave rectification", fontsize=12) 
 
ax = fig.add_subplot(313)
S, H = [], []
for time in t:
        s, h = fullWave(time)
        S = append(S,s)
        H = append(H,h)
ax.plot(t,S+H,'k',linewidth=2.)
ax.plot(t,(signal(t,0)),'b--',linewidth=1.)
ax.plot(t,(signal(t,1)),'r--',linewidth=1.)
ax.plot(t,(signal(t,2)),'g--',linewidth=1.)
myAxes(ax)
ax.set_title(r"Full-wave rectification", fontsize=12) 
 
myLabel = []
for i in arange(0,cycles+.25,.25):
   myLabel = append(myLabel,r"%.2fT"%i)
#    myLabel = append(myLabel,r"${}^{%.1fT}_{\pi/%.1f}$"%(i,(i*2)))
 
ax.set_xticklabels(myLabel,fontsize=10)
ax.set_xlabel(r"Time",fontsize=14)
 
#fig.suptitle("3-phase AC rectification",fontsize=16)
 
fig.savefig("3_phase_rectification_2.svg",bbox_inches="tight",\
        pad_inches=.15)