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

 
1
# coding: utf-8
2
#
3
#   phasor.rb - 
4
#   	$Release Version: 0.1 $
5
#   	$Revision: 0.4 $
6
#   	$Date: 2009/03/26 07:34:05 $
7
#   	by Peter Hillerström
8
#
9
# ----
10
#
11
# === Author
12
#
13
# Written by Peter Hillerström.
14
# 
15
# === Attributions and license
16
#
17
# This class is adapted from the Ruby Complex number class 
18
# written by Keiju ISHITSUKA (SHL Japan Inc.)
19
#
20
# === Usage
21
#
22
# phasor.rb implements the Phasor class for complex numbers in polar form.
23
# Some methods in other Numeric classes (including Complex) are redefined or 
24
# added to allow greater interoperability with Complex numbers and Phasors.
25
#
26
# Phase vectors can be created in the following manner:
27
# - <tt>Phasor(radius, theta)</tt>
28
#   
29
# Additionally, note the following:
30
# - <tt>Phasor::I</tt> (the mathematical constant <i>i</i>)
31
# - <tt>Numeric#phase</tt> (e.g. <tt>5.phase -> Phasor(0, 5.0)</tt>)
32
#
33
# The following +Math+ module methods are defined to do calculations using 
34
# Phasors arguments. They will work as normal with non-Complex arguments.
35
#    fexp flog
36
# 
37

    
38

    
39
#
40
# Numeric is a built-in class on which Fixnum, Bignum, etc., are based.  Here
41
# some methods are added so that all number types can be treated to some extent
42
# as phase vectors.
43
#
44
class Numeric
45
  #
46
  # Returns a Phasor <tt>(0,<i>self</i>)</tt>.
47
  #
48
  def phase
49
    Phasor(0, self)
50
  end
51
  
52
  #
53
  # Convert scalar to phasor
54
  #
55
  def to_phasor
56
    Phasor(abs, angle)
57
  end
58
end
59

    
60

    
61
#
62
# Redefine some Complex methods to co-operate with Phasors
63
#
64
class Complex < Numeric
65
  
66
  undef step if defined?(self.step)
67
  
68
  #
69
  # Attempts to coerce +other+ to a Complex number.
70
  # Redefined to return Phasors if other is a Phasor.
71
  #
72
  def coerce(other)
73
    if Complex.generic?(other)
74
      return Complex.new!(other), self
75
    elsif other.class == Phasor
76
      return other, Phasor.new(self.abs, self.angle)
77
    else
78
      super
79
    end
80
  end
81
  
82
  #
83
  # Redefined to accurately compare phasors to complex numbers
84
  # (using abs and angle instead of real and imag)
85
  #
86
  def == (other)
87
    if other.kind_of?(Complex) or Complex.generic?(other)
88
      self.abs == other.abs and self.angle == other.angle
89
    else
90
      other == self
91
    end
92
  end
93
    
94
end
95

    
96

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

    
106

    
107
# === Description
108
# 
109
# Complex number class using polar coordinates.
110
# 
111
# The representation of a complex number by its polar coordinates is 
112
# called the polar form of the complex number. These can also 
113
# be viewed as phase vectors (phasors) and represented with angle notation 
114
# or with exponential notation using Euler's formula.
115
# 
116
# === Definition
117
#
118
# Alternatively to the cartesian representation z = x + iy,
119
# the complex number z can be specified by polar coordinates.
120
# The polar coordinates are:
121
# r = |z| ≥ 0, called the absolute value or modulus, and
122
# φ = arg(z), called the argument or the angle of z. 
123
# 
124
# For more information, see:
125
# http://en.wikipedia.org/wiki/Complex_number#Polar_form and
126
# http://en.wikipedia.org/wiki/Polar_coordinates#Complex_numbers
127
#
128
class Phasor < Numeric
129

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

    
135
  def scalar?
136
    false
137
  end
138

    
139
  def Phasor.scalar?(other) # :nodoc:
140
    other.kind_of?(Integer) or
141
    other.kind_of?(Float) or
142
    (defined?(Rational) and other.kind_of?(Rational))
143
  end
144
  
145
  #
146
  # Creates a +Phasor+ <tt>a</tt>∠<tt>b</tt>.
147
  #
148
  def Phasor.new!(a, b = 0)
149
    new(a,b)
150
  end
151
  
152
  def initialize(a, b = 0)
153
    raise TypeError, "non numeric 1st arg `#{a.inspect}'" if !a.kind_of? Numeric
154
    raise TypeError, "`#{a.inspect}' for 1st arg" if a.kind_of? Complex
155
    raise TypeError, "non numeric 2nd arg `#{b.inspect}'" if !b.kind_of? Numeric
156
    raise TypeError, "`#{b.inspect}' for 2nd arg" if b.kind_of? Complex
157
    @abs = a
158
    @angle = b
159
    #self.real; self.image # to make Polar#abs & other methods which read @real & @image to work
160
    self
161
  end
162

    
163
  def real
164
    @real ||= @abs * Math.cos(@angle)
165
  end
166

    
167
  def image
168
    @image ||= @abs * Math.sin(@angle)
169
  end
170
  alias imag image
171

    
172
  def abs
173
    @abs ||= Math.hypot(@real, @image)
174
  end
175
  alias amp abs
176
  
177
  def angle
178
    @angle ||= Math.atan2!(@image, @real)
179
  end
180
  alias arg angle
181
  alias phase angle
182

    
183
  #
184
  # Addition with scalar or phase vector.
185
  # 
186
  # Calculated with trigonometric functions, because this
187
  # is faster than round-trip conversion to complex numbers.
188
  #
189
  def + (other)
190
    if other.kind_of?(Complex)
191
      gamma = other.angle - @angle
192
      abs = Math.sqrt!( self.abs2 + other.abs2 + (2 * @abs * other.abs * Math.cos!(gamma)) )
193
      angle = Math.asin!( other.abs / (abs / Math.sin!(Math::PI - gamma)) ) + @angle
194
      Phasor(abs, angle)
195
    elsif Phasor.scalar?(other)
196
      gamma = other.angle - @angle
197
      abs = Math.sqrt!( self.abs2 + other ** 2 + (2 * @abs * other.abs * Math.cos!(gamma)) )
198
      angle = Math.asin!( other.abs / (abs / Math.sin!(Math::PI - gamma)) ) + @angle
199
      Phasor(abs, angle)
200
    else
201
      x , y = other.coerce(self)
202
      x + y
203
    end
204
  end
205

    
206
  #
207
  # Subtraction with scalar or phase vector.
208
  # 
209
  # Calculated with trigonometric functions, because this
210
  # is faster than round-trip conversion to complex numbers.
211
  #
212
  # FIXME:
213
  # Substracting phasors, where image is negative, does not give correct result.
214
  #
215
  def - (other)
216
    if other.kind_of?(Complex)
217
      gamma = other.angle - @angle
218
      abs = Math.sqrt!( self.abs2 + other.abs2 - (2 * @abs * other.abs * Math.cos!(gamma)) )
219
      angle = @angle - Math.asin!( other.abs / (abs / Math.sin!(gamma)) ).abs
220
      if gamma < 0 then angle = Math::PI + angle end
221
      Phasor(abs, angle)
222
    elsif Phasor.scalar?(other)
223
      gamma = other.angle - @angle
224
      abs = Math.sqrt!( self.abs2 + other ** 2 - (2 * @abs * other.abs * Math.cos!(gamma)) )
225
      angle = @angle - Math.asin!( other.abs / (abs / Math.sin!(gamma)) ).abs
226
      if gamma < 0 then angle = Math::PI + angle end # Not sure if this is always correct!
227
      Phasor(abs, angle)
228
    else
229
      x , y = other.coerce(self)
230
      x - y
231
    end
232
  end
233

    
234
  #
235
  # Multiplication with scalar or phase vector.
236
  #
237
  def * (other)
238
    if other.kind_of?(Complex)
239
      abs = @abs * other.abs
240
      angle = @angle + other.angle
241
      Phasor(abs, angle)
242
    elsif Phasor.scalar?(other)
243
      Phasor(@abs * other, @angle)
244
    else
245
      x , y = other.coerce(self)
246
      x * y
247
    end
248
  end
249

    
250
  #
251
  # Division by scalar or phase vector.
252
  #
253
  def / (other)
254
    if other.kind_of?(Complex)
255
      abs = @abs / other.abs
256
      angle = @angle - other.angle
257
      Phasor(abs, angle)
258
    elsif Phasor.scalar?(other)
259
      Phasor(@abs / other, @angle)
260
    else
261
      x , y = other.coerce(self)
262
      x / y
263
    end
264
  end
265

    
266
  #
267
  # Exponentiation with scalar or phase vector.
268
  #
269
  def ** (other)
270
    if other == 0
271
      return Phasor(1)
272
    end
273
    if other.kind_of?(Complex)
274
      Math.fexp( (Math.log!(@abs) + @angle * Phasor::I) * other )
275
    elsif Phasor.scalar?(other)
276
      Phasor(@abs ** other, @angle * other) # Does not work if either angle == 0!
277
    else
278
      x , y = other.coerce(self)
279
      x ** y
280
    end
281
  end
282
  
283
  #
284
  # Remainder after division by a scalar or phase vector.
285
  #
286
  def % (other)
287
    if other.kind_of?(Complex)
288
      Complex(@real % other.real, @image % other.image).to_phasor
289
    elsif Phasor.scalar?(other)
290
      Complex(@real % other, @image % other).to_phasor
291
    else
292
      x , y = other.coerce(self)
293
      x % y
294
    end
295
  end
296
  
297
  #
298
  # Principal nth root of scalar or phase vector.
299
  #
300
  def ^ (other)
301
    self ** (1.0/other)
302
  end
303
  
304
  #
305
  # Nth roots of unity
306
  #
307
  def Phasor.nth_roots (n, amp = 1, rotations = 1)
308
    acc = Phasor::C/n
309
    (0...n * rotations).map { |k| Phasor(amp, k * acc) }
310
  end
311

    
312
  #
313
  # Square of the absolute value.
314
  #
315
  def abs2
316
    @abs*@abs
317
  end
318

    
319
  #
320
  # Convert a phasor to complex number.
321
  #
322
  def to_complex
323
    Complex(real, imag)
324
  end
325
  
326
  #
327
  # Complex conjugate (<tt>p * p.conjugate == p.abs**2</tt>).
328
  #
329
  def conjugate
330
    Phasor(@abs, - @angle)
331
  end
332
  alias conj conjugate
333
  
334
  #
335
  # Compares the absolute values of the two numbers.
336
  #
337
  def <=> (other)
338
    self.abs <=> other.abs
339
  end
340
  
341
  #
342
  # Phasor is a kind of <tt>Complex Number</tt>.
343
  #
344
  def kind_of?(klass)
345
    klass == Complex ? true : super;
346
  end
347
  
348
  #
349
  # Test for numerical equality (<tt>a == a + 0.phase</tt>).
350
  #
351
  def == (other)
352
    if other.kind_of?(Complex) or Phasor.scalar?(other)
353
      @abs == other.abs and @angle == other.angle
354
    else
355
      other == self
356
    end
357
  end
358
  
359
  #
360
  # Attempts to coerce +other+ to a phase vector.
361
  #
362
  def coerce(other)
363
    if other.kind_of?(Complex) or Phasor.scalar?(other)
364
      return Phasor.new!(other.abs, other.angle), self
365
    else
366
      super
367
    end
368
  end
369

    
370
  #
371
  # Angle notation of the phase vector.
372
  #
373
  def to_s
374
    if @abs != 0
375
      @abs.to_s + "∠" + @angle.to_s
376
    else
377
      "∠" + @angle.to_s
378
    end
379
  end
380
  
381
  #
382
  # Returns a hash code for the phase vector.
383
  #
384
  def hash
385
    @abs.hash ^ @angle.hash
386
  end
387
  
388
  #
389
  # Returns "<tt>Phasor(<i>abs</i>, <i>angle</i>)</tt>".
390
  #
391
  def inspect
392
    sprintf("Phasor(%s, %s)", abs, angle)
393
  end
394
  
395
  #
396
  # Angle measure of full circle in radians.
397
  #
398
  C = Math::PI*2 # radians
399
  # TODO: Enable use of other units than radians.
400
  #C = 360        # degrees - Note! changing units does not work reliably!
401
  #C = 1.0
402
  
403
  #
404
  # <b>I</b> is the imaginary number. It exists at point (0,1) on the complex plane.
405
  #
406
  I = Phasor(1, C/4.0)
407
  
408
  #
409
  # Absolute value (aka modulus):
410
  # distance from the zero point on the complex plane.
411
  #
412
  attr :abs
413
  
414
  #
415
  # Argument (angle from (1,0) on the complex plane).
416
  #
417
  attr :angle
418
end
419

    
420
module Math
421
  
422
  # Redefined to handle a Phasor argument.
423
  def flog(z)
424
    if Phasor.scalar?(z) and z >= 0
425
      log!(z)
426
    else
427
      Complex(log!(z.abs), z.angle).to_phasor
428
      # FIXME: Find formula to calculate with phasors
429
      #Phasor( log!(z.abs * Math.cos!(z.angle)), z.abs * Math.sin!(z.angle) )
430
    end
431
  end
432

    
433
  # Redefined to handle a Phasor argument.
434
  def fexp(p)
435
    if Phasor.scalar?(p) and p >= 0
436
      exp!(p)
437
    else
438
      # http://en.wikipedia.org/wiki/Complex_exponential_function
439
      Phasor( exp!(p.abs * Math.cos!(p.angle)), p.abs * Math.sin!(p.angle) )
440
    end
441
  end
442
  
443
  module_function :flog
444
  module_function :fexp
445
  
446
end