# acpower.py
#
# Python module for doing alternating current phasor computations
#
# (c) Samuel May 7/11/09
#

__doc__ = """
A python module containing classes and functions for quickly making power
engineering calculations at the Python prompt.

Put in your current directory and load with 'import acpower' after starting
python. You may then create phasors with something like:

     $ python
     Python 2.6.4 (r264:75706, Oct 27 2009, 06:25:13)
     [GCC 4.4.1] on linux2
     Type \"help\", \"copyright\", \"credits\" or \"license\" for more information.
     >>> import acpower
     >>> v = acpower.Phasor(1+3j)
     >>> v
     phasor with:
     polar coods        =  3.162 /_ 71.565
     rectangular coords =  1.000+3.000j

You can call the helper functions with 'acpower.function(args)'
e.g:

     >>> z = acpower.Phasor([1.414,44.977])
     >>> z
     phasor with:
     polar coods        =  1.414 /_ 44.977
     rectangular coords =  1.000+0.999j
     >>> total = acpower.parallel(v,z)
     >>> total
     phasor with:
     polar coords       =  1.000 /_ 53.113
     rectangular coords =  0.600+0.800j

The python operators are overloaded, so you can use +,-,*,and / like you
would with an int or a float:

     >>> v/z
     phasor with:
     polar coods        =  2.236 /_ 0.464
     rectangular coords =  2.000+1.001j

By default all angles are in degrees. You can change this to radians with:

      >>> acpower.angle_unit = \"rad\"

eg:

      >>> z
      phasor with:
      polar coords       =  1.414 /_ 44.977
      rectangular coords =  1.000+0.999j
      >>> acpower.angle_unit = \"rad\"
      >>> z
      phasor with:
      polar coords       =  1.414 /_ 0.785
      rectangular coords =  1.000+0.999j

and back with:

      >>> acpower.angle_unit = \"deg\"
"""

import math

# Angle representation extension. Internally, everything is in radians, so we
# only worry when reading in an angle (should we convert it to radians?), and
# when we print in out (should we convert it to degrees?)
angle_unit = "deg"

def deg2rad(theta):
    """Convert the given angle from degrees to radians."""
    return (theta/180.0)*math.pi

def rad2deg(theta):
    """Convert the given angle from radians to degrees."""
    return (theta/math.pi)*180.0

def is_phasor(p):
    """Returns true or false depending on whether the given object is a phasor."""
    return hasattr(p,"magnitude") and hasattr(p,"phase") and hasattr(p,"rect")

def is_scalar(p):
    """Returns true or false depending on whether the given object is a phasor."""
    return not is_phasor(p)

class Phasor:
    """An object with magnitude and phase. Contains two internal representaions,
    a complex number in rectangular form, and a magnitude and phase. Possible
    ways to create a phasor:
    
    No arguments:
    
        >>> p = acpower.Phasor()   - sets everything to zero
             
    With an ordinary complex number:
        
        >>> p = acpower.Phasor(1+3j)
    
    With a list or tuple of magnitude and phase:
        
        >>> p = acpower.Phasor([1.414,0.785])"""

    def __init__(self,z=0+0j):
        
        # see if it's a list or tuple, otherwise, treat it as a complex number
        if hasattr(z,"__getitem__"):
            self.magnitude = z[0]
            # compensate if they have given us the phase in degrees
            if (angle_unit == "deg"):
                z[1] = deg2rad(z[1])
            self.phase = z[1]
            self.__update_rect()
        # no? it must be a complex number
        else:
            self.rect = z
            self.__update_polar()
        # if not, exception is thrown and you deal with the error


    def __update_polar(self):
        """Sets this phasor's magnitude and phase from its rectangular complex
        number representation. Private function for messing about with the
        internal representation - DO NOT USE"""
        # save a little work by using some slightly more sophisticated functions
        # from the python math library. 'hypot(x,y) calculates sqrt(x^2 + y^2),
        # i.e. the norm or modulus.
        self.magnitude = math.hypot(self.rect.real,self.rect.imag)
        # 'atan2(x,y)' calculates inverse tan of x/y, i.e. argument.
        self.phase = math.atan2(self.rect.imag,self.rect.real)

    def __update_rect(self):
        """Sets this phasor's rectangular complex number representation from its
        magnitude and phase. Private funtion for messing about with the internal
        representation - DO NOT USE"""
        # go backwards from the modulus and argument, using cos and sin math
        # functions.
        self.rect = complex(self.magnitude*math.cos(self.phase),
                            self.magnitude*math.sin(self.phase))

    ## overloaded arithmetic methods ##

    # this is the good thing about having the two internal reps.  we can use the
    # x+iy form for addition and subtraction, and the polar form for
    # multiplication and division.

    # the reason we have the __blah__ function names is because we're
    # overwriting the python methods for '+', '-', '*', and '/', allowing us to
    # treat phasors just likes ints or doubles

    # also trying a tweak to add multiplication/division by scalars. Not sure if
    # I want to extend this to addition/subtraction.
    
    def __add__(self,other):
        """Returns the phasor that is the sum of two phasor arguments."""
        p = Phasor()
        p.rect = self.rect + other.rect
        p.__update_polar()
        return p
    
    def __sub__(self,other):
        """Returns the phasor that is the second argument subracted
        from the first."""
        p = Phasor()
        p.rect = self.rect - other.rect
        p.__update_polar()
        return p
    
    def __mul__(self,other):
        """Returns the product of two phasors."""
        p = Phasor()
        if is_scalar(other):
            other = Phasor([other,0])
        
        p.magnitude = self.magnitude * other.magnitude
        p.phase = self.phase + other.phase
            
        # keep theta element of [-pi,pi]
        if p.phase > math.pi:
            p.phase -= 2*math.pi
        elif p.phase < -math.pi:
            p.phase += 2*math.pi
        
        p.__update_rect()
        return p
    
    def __div__(self,other):
        """Returns the result of dividing the first phasor by the second."""
        p = Phasor()
        if is_scalar(other):
            other = Phasor([other,0])
        
        p.magnitude = self.magnitude / other.magnitude
        p.phase = self.phase - other.phase
        # keep theta element of [-pi,pi]
        if p.phase > math.pi:
            p.phase -= 2*math.pi
        elif p.phase < -math.pi:
            p.phase += 2*math.pi

        # update x+iy form
        p.__update_rect()
        return p

    # this is the string representation you see. (syntax is like printf)
    def __repr__(self):
        # print in degrees if that is the external angle representation
        if (angle_unit == "deg"):
            angle_rep = rad2deg(self.phase)
        else:
            angle_rep = self.phase
        
        string = """phasor with:
polar coords       = % .3f /_ %.3f
rectangular coords = % .3f%+.3fj""" % (self.magnitude,angle_rep,self.rect.real,self.rect.imag)
        return string

## end of Phasor class ##

# useful formulae for working with impedances. (you really should know these)

# parallel impedances
#
# 1     1  1  1  1
# -  =  -  -  -  -  ...
# z     z  z  z  z
#  p     2  2  2  2
#
def parallel(*impedances):
    """Return the equivalent total impedance of a list of impedances in parallel."""
    if not impedances:
        return Phasor()
    
    one = Phasor(1)
    total = Phasor()
    for z in impedances:
        total = total + one/z
        
    return one/total

# series impedances
# 
# z  = z + z + z + z  ...
#  s    1   2   3   4
#
def series(*impedances):
    """Return the equivalent total impedance of a list of impedances in series."""
    i = Phasor()
    for z in impedances:
        i = i + z
    return i

# ac power measurements

def apparent_power(V,I):
    """Return the apparent power of the phasors voltage V and current I."""
    return (V.magnitude*I.magnitude)/2

def power_factor(V,I):
    """Return the power factor of the phasors voltage V and current I."""
    return math.cos(V.phase - I.phase)

def average_power(V,I):
    """Return the average power of the phasors voltage V and current I."""
    # magnitude = Vm*Im/2
    a = apparent_power(V,I)
    # power factor
    pf = power_factor(V,I)
    return a*pf

def reactive_power(V,I):
    """Return the reactive power of the phasors voltage V and current I."""
    a = apparent_power(V,I)
    # sine of difference of phases
    s = math.sin(V.phase - I.phase)
    return a*s

def complex_power(V,I):
    """Return the complex power of the phasors voltage V and current I."""
    return complex(average_power(V,I),reactive_power(V,I))