User:Andrew Nikitin/Phase portraits

Script

The book "Visual Complex Functions" by Elias Wegert describes phase portraits, a visualization technique for complex functions of one complex variables. Each function value is represented by a pixel with a hue (which is an angular quantity) equal to function value's phase angle. Author argues that the resulting picture provides a lot of useful information about function.

Additionally, two other colorspace components, saturation and light, may be used to further enhance the view with a grid-like appearance of lines of equal phase and equal magnitude. Author calls this improved palette "enhanced phase portraits".

NB. 2017-04-07
NB. Phase portrait for complex functions
NB. (after Wegert, Elias. Visual Complex Functions)

require 'viewmat'

NB.*rgbfh v represent hue as unscaled RGB triple
NB. y=hue, 0..1
NB. http://stackoverflow.com/questions/2353211/hsl-to-rgb-color-conversion
rgbfh=:(0 >. 1 <. _1 + [: | _3 + 6 | 0 4 2 + 6 * ])">

NB.*rgbfhsl v convert hsl color representation to rgb
NB. y= 3#0..1
NB. returns unscaled rgb 3#0..1
rgbfhsl=:(3 : 0)"1
  'h s l'=.y
  l+s*(_0.5+rgbfh h)*-.|<:+:l
)

NB. Phase portraits

NB. unit square y points out
sq=:[: |. [: j.~/~ (%~ i:)

NB.*ccEnhPh v enhanced phase (with lines of equal phase and magnitude)
NB. Phase is represented by a hue component.
NB. Magnitude blocks are logarithmic and have light gradient toward
NB. increase.
NB. Phase blocks are separated with saturation drops.
ccEnhPh=:3 : 0
  8 ccEnhPh y
:
  NB. :: is to deal with occasional NaN error
  'r p' =.(i.@<:@#@$ |: ]) *. :: (0 0"_)">y
  p=.2p1%~2p1|p
  r=.(* -.@(_&= +. __&=)) 2p1%~^.r

  NB. 5^~ defines how sharp the gradient is at the edge
  NB. use "l" to visualize magnitude gradient and isomagnitude lines
  l=.0.5 +0.25*5^~_1 2 p. 1|x*r
  NB. use "s" to show isophase lines
  s=.-.0.5* 5^~1| x*p
  rgb=.rgbfhsl 0|:p,s,:l
  256#.255<.<.256*rgb
  NB.256#.(2^8-p)*"1<.(2^p=.3 3 2)(<:@[ <. *)"1 rgb NB. MSX screen 8
)

Try it out

z, identity

Right halfplane has "warm" color phase, left has "cold". Rainbow color order follows positive (counterclockwise) direction.

  viewrgb 8 ccEnhPh sq 128

Exponent

  viewrgb 8 ccEnhPh ^ 3*sq 128

Simple rational function

Random monic polynomial divided by z. One pole (in the center), 6 zeroes. Note that the pole has opposite order of phases around it.

  viewrgb 8 ccEnhPh (%~ (1,~j./?0$~2 6)&p.)) 1.5*sq 128

Tangent

Tangent of simple rational function similar to above. :: 0: is needed to deal with NaN errors.

  viewrgb 8 ccEnhPh (3 o. (%~ (1,~j./?0$~2 6)&p.)) :: 0:">  1.5*sq 128

Cosine

Cosine of simple rational function above. :: 0: is needed to deal with NaN errors.

  viewrgb 8 ccEnhPh (2 o. (%~ (1,~j./?0$~2 6)&p.)) :: 0:">  1.5*sq 128

Butterworth filter

Evenly spaced zeroes on the lower half of the unit circle.

  viewrgb 8 ccEnhPh (p. 1;^j.1p1+1p1*(%~i.@>:)16)&p. 2*sq 128

Multiple palettes

It seems that enhanced phase palette provides the most information (and also is the prettiest), but other palettes are possible and may be used to enhance different aspects of the function.

Here is a script that implements some other palettes. I had to factor out common and repetitive parts, so it may be less clear what is going on. This is why I keep the original script (above) as it shows standalone and clean implementation of the most informative case, the enhanced phase portraits.

NB. 2017-04-12
NB. Phase portrait for complex functions
NB. (after Wegert, Elias. Visual Complex Functions)

require 'viewmat'

NB.*rgbfh v represent hue as unscaled RGB triple
NB. y=hue, 0..1
NB. http://stackoverflow.com/questions/2353211/hsl-to-rgb-color-conversion
rgbfh=:(0 >. 1 <. _1 + [: | _3 + 6 | 0 4 2 + 6 * ])">

NB.*rgbfhsl v convert hsl color representation to rgb
NB. y= 3#0..1
NB. returns unscaled rgb 3#0..1
rgbfhsl=:(3 : 0)"1
  'h s l'=.y
  l+s*(_0.5+rgbfh h)*-.|<:+:l
)

NB.*clr v scale vector of 0..1 values to 0..255 integers and combine
NB. then into an integer as required by viewrgb and writebmp
clr24=:256 #. 255 <. [: <. 256&*
clr8e=:256 #. (2^3 3 2) ((256%[) * <:@[ <. *)"1 ] NB. simulate color depth of MSX screen 8

clr=:clr24 f.

NB. phase portrait

NB.*sq v unit square y points out
NB. The matrix is reversed so that +i direction points "up" on plots
NB. (towards index decrease)
sq=:[: |. [: j.~/~ (%~ i:)

NB.*CcPalette a construct color paletter verb from a rule mapping phase
NB. and magnitude to RGB01.
NB. u=color palette, a verb that converts phase, magnitude to rgb
NB. u takes:
NB. x = grid density
NB. y = 2 component array, 0&{ = log magnitude, 1&{=phase
CcPalette=:1 : 0
  [:
:
  NB. :: is to deal with occasional NaN error
  'r p' =.(i.@<:@#@$ |: ]) *. :: (0 0"_)">y
  p=.2p1%~2p1|p
  r=.(* -.@(_&= +. __&=)) 2p1%~^.r
  clr x u r,:p
)

NB.*DefaultLeft a provide default left argument
DefaultLeft=:1 : ('8 u y';':';'x u y')


NB. palette construction (x=grid density)

NB. two sided gradient, close to 0.5 in the middle, sharply drops to 0
NB. or 1 at the ends
NB. 5^~ defines how sharp the gradient is at the edge
grdLight=:0.5 + 0.25 * 5 ^~ _1 + 2 * 1 | *

NB. one sided gradient
NB. stays 1 most of the time, sharply drops to 0.5 at one end
grdSat=:1 - 0.5 * 5 ^~ 1 | *

Cc=.CcPalette DefaultLeft

NB.*ccEnhPh v enhanced phase (with lines of equal phase and magnitude)
NB. Phase is represented by a hue component.
NB. Magnitude blocks are logarithmic and have light gradient toward
NB. increase.
NB. Phase blocks are separated with saturation drops.
ccEnhPh=:([: rgbfhsl 0 |:(1{]),(grdSat 1&{),:(grdLight 0&{)) Cc

NB.*ccPh v phase only (hue=phase)
ccPh=:([:rgbfh 1{]) CcPalette DefaultLeft

NB.*ccMagPh v phase and magnitude (no lines of equal phase)
ccMagPh=:([: rgbfhsl 0 |:(1{]),1,:(grdLight 0{])) Cc

gray=:3&#">

NB. high contrast phase and magnitude
ccEnhPhBW=:(3 #"> ~:&(<.@:+:@:(1&|))/)@:* Cc

NB.*ccEnhPhGray v grayscale phase and magnitude
ccEnhPhGray=:(3 #"> -:@:+&(1&|)/)@:* Cc

NB.*ccMagGray v grayscale magnitude
ccMagGray=:(3 #"> 1 | 0 { ])@:* Cc

NB.*ccPhGray v grayscale phase
ccPhGray=:(3 #"> 1 | 1 { ])@:* Cc

NB. shortcut; use as e.g.
NB.   (3&o.) pp sq 128
pp=:1 : '[: viewrgb [: ccEnhPh (u :: 0:)">'