Project

General

Profile

Feature #1316 » phasor.rb

Phasor.rb: Complex number class using polar coordinates - peterhil (Peter Hillerström), 03/26/2009 04:58 PM

 
# coding: utf-8
#
# phasor.rb -
# $Release Version: 0.1 $
# $Revision: 0.4 $
# $Date: 2009/03/26 07:34:05 $
# by Peter Hillerström
#
# ----
#
# === Author
#
# Written by Peter Hillerström.
#
# === Attributions and license
#
# This class is adapted from the Ruby Complex number class
# written by Keiju ISHITSUKA (SHL Japan Inc.)
#
# === Usage
#
# phasor.rb implements the Phasor class for complex numbers in polar form.
# Some methods in other Numeric classes (including Complex) are redefined or
# added to allow greater interoperability with Complex numbers and Phasors.
#
# Phase vectors can be created in the following manner:
# - <tt>Phasor(radius, theta)</tt>
#
# Additionally, note the following:
# - <tt>Phasor::I</tt> (the mathematical constant <i>i</i>)
# - <tt>Numeric#phase</tt> (e.g. <tt>5.phase -> Phasor(0, 5.0)</tt>)
#
# The following +Math+ module methods are defined to do calculations using
# Phasors arguments. They will work as normal with non-Complex arguments.
# fexp flog
#


#
# Numeric is a built-in class on which Fixnum, Bignum, etc., are based. Here
# some methods are added so that all number types can be treated to some extent
# as phase vectors.
#
class Numeric
#
# Returns a Phasor <tt>(0,<i>self</i>)</tt>.
#
def phase
Phasor(0, self)
end
#
# Convert scalar to phasor
#
def to_phasor
Phasor(abs, angle)
end
end


#
# Redefine some Complex methods to co-operate with Phasors
#
class Complex < Numeric
undef step if defined?(self.step)
#
# Attempts to coerce +other+ to a Complex number.
# Redefined to return Phasors if other is a Phasor.
#
def coerce(other)
if Complex.generic?(other)
return Complex.new!(other), self
elsif other.class == Phasor
return other, Phasor.new(self.abs, self.angle)
else
super
end
end
#
# Redefined to accurately compare phasors to complex numbers
# (using abs and angle instead of real and imag)
#
def == (other)
if other.kind_of?(Complex) or Complex.generic?(other)
self.abs == other.abs and self.angle == other.angle
else
other == self
end
end
end


#
# Creates a Phasor. +a+ and +b+ should be Numeric. The result will be
# <tt>a∠b</tt>.
#
def Phasor(a, b = 0)
b = b % (Phasor::C) #if defined?(Unify) # Note! Spiral drawing does not work with Unify
Phasor.new(a, b)
end


# === Description
#
# Complex number class using polar coordinates.
#
# The representation of a complex number by its polar coordinates is
# called the polar form of the complex number. These can also
# be viewed as phase vectors (phasors) and represented with angle notation
# or with exponential notation using Euler's formula.
#
# === Definition
#
# Alternatively to the cartesian representation z = x + iy,
# the complex number z can be specified by polar coordinates.
# The polar coordinates are:
# r = |z| ≥ 0, called the absolute value or modulus, and
# φ = arg(z), called the argument or the angle of z.
#
# For more information, see:
# http://en.wikipedia.org/wiki/Complex_number#Polar_form and
# http://en.wikipedia.org/wiki/Polar_coordinates#Complex_numbers
#
class Phasor < Numeric

# step could be used to step between angles (& abs too!) by either
# a) requiring argument to be a scalar (easier) and trace a straight line between points.
# b) giving a phasor/complex number as argument and tracing a circle trough 3 points.
undef step if defined?(self.step)

def scalar?
false
end

def Phasor.scalar?(other) # :nodoc:
other.kind_of?(Integer) or
other.kind_of?(Float) or
(defined?(Rational) and other.kind_of?(Rational))
end
#
# Creates a +Phasor+ <tt>a</tt>∠<tt>b</tt>.
#
def Phasor.new!(a, b = 0)
new(a,b)
end
def initialize(a, b = 0)
raise TypeError, "non numeric 1st arg `#{a.inspect}'" if !a.kind_of? Numeric
raise TypeError, "`#{a.inspect}' for 1st arg" if a.kind_of? Complex
raise TypeError, "non numeric 2nd arg `#{b.inspect}'" if !b.kind_of? Numeric
raise TypeError, "`#{b.inspect}' for 2nd arg" if b.kind_of? Complex
@abs = a
@angle = b
#self.real; self.image # to make Polar#abs & other methods which read @real & @image to work
self
end

def real
@real ||= @abs * Math.cos(@angle)
end

def image
@image ||= @abs * Math.sin(@angle)
end
alias imag image

def abs
@abs ||= Math.hypot(@real, @image)
end
alias amp abs
def angle
@angle ||= Math.atan2!(@image, @real)
end
alias arg angle
alias phase angle

#
# Addition with scalar or phase vector.
#
# Calculated with trigonometric functions, because this
# is faster than round-trip conversion to complex numbers.
#
def + (other)
if other.kind_of?(Complex)
gamma = other.angle - @angle
abs = Math.sqrt!( self.abs2 + other.abs2 + (2 * @abs * other.abs * Math.cos!(gamma)) )
angle = Math.asin!( other.abs / (abs / Math.sin!(Math::PI - gamma)) ) + @angle
Phasor(abs, angle)
elsif Phasor.scalar?(other)
gamma = other.angle - @angle
abs = Math.sqrt!( self.abs2 + other ** 2 + (2 * @abs * other.abs * Math.cos!(gamma)) )
angle = Math.asin!( other.abs / (abs / Math.sin!(Math::PI - gamma)) ) + @angle
Phasor(abs, angle)
else
x , y = other.coerce(self)
x + y
end
end

#
# Subtraction with scalar or phase vector.
#
# Calculated with trigonometric functions, because this
# is faster than round-trip conversion to complex numbers.
#
# FIXME:
# Substracting phasors, where image is negative, does not give correct result.
#
def - (other)
if other.kind_of?(Complex)
gamma = other.angle - @angle
abs = Math.sqrt!( self.abs2 + other.abs2 - (2 * @abs * other.abs * Math.cos!(gamma)) )
angle = @angle - Math.asin!( other.abs / (abs / Math.sin!(gamma)) ).abs
if gamma < 0 then angle = Math::PI + angle end
Phasor(abs, angle)
elsif Phasor.scalar?(other)
gamma = other.angle - @angle
abs = Math.sqrt!( self.abs2 + other ** 2 - (2 * @abs * other.abs * Math.cos!(gamma)) )
angle = @angle - Math.asin!( other.abs / (abs / Math.sin!(gamma)) ).abs
if gamma < 0 then angle = Math::PI + angle end # Not sure if this is always correct!
Phasor(abs, angle)
else
x , y = other.coerce(self)
x - y
end
end

#
# Multiplication with scalar or phase vector.
#
def * (other)
if other.kind_of?(Complex)
abs = @abs * other.abs
angle = @angle + other.angle
Phasor(abs, angle)
elsif Phasor.scalar?(other)
Phasor(@abs * other, @angle)
else
x , y = other.coerce(self)
x * y
end
end

#
# Division by scalar or phase vector.
#
def / (other)
if other.kind_of?(Complex)
abs = @abs / other.abs
angle = @angle - other.angle
Phasor(abs, angle)
elsif Phasor.scalar?(other)
Phasor(@abs / other, @angle)
else
x , y = other.coerce(self)
x / y
end
end

#
# Exponentiation with scalar or phase vector.
#
def ** (other)
if other == 0
return Phasor(1)
end
if other.kind_of?(Complex)
Math.fexp( (Math.log!(@abs) + @angle * Phasor::I) * other )
elsif Phasor.scalar?(other)
Phasor(@abs ** other, @angle * other) # Does not work if either angle == 0!
else
x , y = other.coerce(self)
x ** y
end
end
#
# Remainder after division by a scalar or phase vector.
#
def % (other)
if other.kind_of?(Complex)
Complex(@real % other.real, @image % other.image).to_phasor
elsif Phasor.scalar?(other)
Complex(@real % other, @image % other).to_phasor
else
x , y = other.coerce(self)
x % y
end
end
#
# Principal nth root of scalar or phase vector.
#
def ^ (other)
self ** (1.0/other)
end
#
# Nth roots of unity
#
def Phasor.nth_roots (n, amp = 1, rotations = 1)
acc = Phasor::C/n
(0...n * rotations).map { |k| Phasor(amp, k * acc) }
end

#
# Square of the absolute value.
#
def abs2
@abs*@abs
end

#
# Convert a phasor to complex number.
#
def to_complex
Complex(real, imag)
end
#
# Complex conjugate (<tt>p * p.conjugate == p.abs**2</tt>).
#
def conjugate
Phasor(@abs, - @angle)
end
alias conj conjugate
#
# Compares the absolute values of the two numbers.
#
def <=> (other)
self.abs <=> other.abs
end
#
# Phasor is a kind of <tt>Complex Number</tt>.
#
def kind_of?(klass)
klass == Complex ? true : super;
end
#
# Test for numerical equality (<tt>a == a + 0.phase</tt>).
#
def == (other)
if other.kind_of?(Complex) or Phasor.scalar?(other)
@abs == other.abs and @angle == other.angle
else
other == self
end
end
#
# Attempts to coerce +other+ to a phase vector.
#
def coerce(other)
if other.kind_of?(Complex) or Phasor.scalar?(other)
return Phasor.new!(other.abs, other.angle), self
else
super
end
end

#
# Angle notation of the phase vector.
#
def to_s
if @abs != 0
@abs.to_s + "∠" + @angle.to_s
else
"∠" + @angle.to_s
end
end
#
# Returns a hash code for the phase vector.
#
def hash
@abs.hash ^ @angle.hash
end
#
# Returns "<tt>Phasor(<i>abs</i>, <i>angle</i>)</tt>".
#
def inspect
sprintf("Phasor(%s, %s)", abs, angle)
end
#
# Angle measure of full circle in radians.
#
C = Math::PI*2 # radians
# TODO: Enable use of other units than radians.
#C = 360 # degrees - Note! changing units does not work reliably!
#C = 1.0
#
# <b>I</b> is the imaginary number. It exists at point (0,1) on the complex plane.
#
I = Phasor(1, C/4.0)
#
# Absolute value (aka modulus):
# distance from the zero point on the complex plane.
#
attr :abs
#
# Argument (angle from (1,0) on the complex plane).
#
attr :angle
end

module Math
# Redefined to handle a Phasor argument.
def flog(z)
if Phasor.scalar?(z) and z >= 0
log!(z)
else
Complex(log!(z.abs), z.angle).to_phasor
# FIXME: Find formula to calculate with phasors
#Phasor( log!(z.abs * Math.cos!(z.angle)), z.abs * Math.sin!(z.angle) )
end
end

# Redefined to handle a Phasor argument.
def fexp(p)
if Phasor.scalar?(p) and p >= 0
exp!(p)
else
# http://en.wikipedia.org/wiki/Complex_exponential_function
Phasor( exp!(p.abs * Math.cos!(p.angle)), p.abs * Math.sin!(p.angle) )
end
end
module_function :flog
module_function :fexp
end
(1-1/2)