Domain Coloring

Domain coloring , assigning a color in the 2D complex plane, is obtained using the surface method map_color_from_op. Using the functions from the Complex Number Representation, Radius and Angle, with the addition of a coloring function, the above 3D figure was constructed.

The 2D plot of domain coloring is obtained using a PlanarSurface object, viewed normal to the xy-plane. This is shown in the figure below.

When the 3D shaded surfaces are viewed normal to the xy-plane, shading enhances the region with higher slopes, as shown below:

import numpy as np
from matplotlib import pyplot as plt
import s3dlib.surface as s3d
import s3dlib.cmap_utilities as cmu

#.. Domain Coloring Example

# 1. Define function to examine .....................................
def complexFunc(z):
    # https://en.wikipedia.org/wiki/Domain_coloring
    i = 1j
    num = (z**2 - 1)*(z - 2 -i)**2
    dem = z**2 + 2 + 2*i
    return num/dem

def funcRT(xyz) :
    x,y,z = xyz
    c = np.array(x,dtype=complex)
    c.imag = np.array(y)
    f = complexFunc(c)
    r = np.abs(f)                      #.. [0,inf]
    Z = 2*np.arctan(r)/np.pi           #.. [0,1]
    theta = (1 + np.angle(f)/np.pi)/2  #.. [0,1]
    return Z,theta 

coor_radius  = lambda c: [c[0],c[1],funcRT(c)[0]]
coor_theta   = lambda c: [c[0],c[1],funcRT(c)[1]]
isContinuous = lambda c,N: np.abs(N.T[2]) > 0.02

def domainColor(xyz) :
    Z,theta = funcRT(xyz)
    h,s,v = hsl_to_hsv(theta,1,Z)
    H = (h+5/6)%1  # equivalent cmap starts at magenta(5/6), not red(0)
    return H,s,v

def hsl_to_hsv(h, s, l):
    # from: https://gist.github.com/mathebox/e0805f72e7db3269ec22
    # Modified for arrays, used Numpy instead of python math.
    v = (2*l + s*(1-np.abs(2*l-1)))/2
    s = 2*(v-l)/v
    return h, s, v

# 2. Setup and map surfaces .........................................
rez = 6

surfaces = [None]*2
surface = s3d.PlanarSurface(rez,color='grey').domain([-3,3],[-3,3])
surface.map_geom_from_op( coor_radius )
surface.map_color_from_op( domainColor, rgb=False )
surfaces[0]=surface

# note: higher rez due to surface discontinuity.
surface = s3d.PlanarSurface(rez+2,color='grey').domain([-3,3],[-3,3])
surface.map_geom_from_op( coor_theta )
normals = surface.facenormals(v3d=False)
surface.clip(lambda c : isContinuous(c,normals))
surface.map_color_from_op( domainColor, rgb=False )
surfaces[1]=surface

# 3. Construct figure, add surface, plot ............................
elev,azim = 20,150

r_labels = [ '0','1/8','1/4','1/2','1','2','4','8',r'$\infty$']
r_axvals = [0] + [ 2*np.arctan(2**i)/np.pi for i in range(-3,4)] + [1]
t_labels = [r'-$\pi$',r'-$\pi$/2','0',r'$\pi$/2',r'$\pi$']
t_axvals = [ v/4 for v in range(0,5)]
z_ticks = [r_axvals,t_axvals]
z_axlab = [r_labels,t_labels]
title = [ r'$\mathcal{R}$, magnitude', r'$\theta$, angle']

fig = plt.figure(figsize=plt.figaspect(0.5/1))
info = r" f(z) =  $\frac{(z^{2}-1)(z-2-i)^{2}}{z^{2}+2+2i}$  "
fig.text(0.5,1,info, ha='center', va='top', fontsize='x-large')
for i,surface in enumerate(surfaces) :
    ax = fig.add_subplot(121+i, projection='3d', aspect='equal')
    ax.view_init(elev,azim)
    #ax.view_init(89.9,-90.01)  # <-- 2D-view from top 
    ax.set_proj_type('ortho')
    ax.set_title('\n'+title[i],color='grey', fontsize='xx-large')
    ax.set(xlim=( -3,3 ), ylim=( -3,3 ), zlim=( 0,1), zticks=z_ticks[i],
                xlabel='Re(z)', ylabel='Im(z)' )
    ax.set_zticklabels(z_axlab[i])

    hlt  = 0.1 if i==0 else 0.7
    surface.shade(0.25,ax=ax,direction=[0,1,1]).hilite(hlt,direction=[0,1,1])
    ax.add_collection3d(surface)

fig.tight_layout(pad=0)

# =======================================================================
# Domain Color Map ======================================================

# 2. Setup and map surfaces .........................................
cmu.hue_cmap('m','+m',name='mycym')

surface = s3d.PlanarSurface(rez,cmap='mycym').domain([-3,3],[-3,3])
surface.map_color_from_op( domainColor, rgb=False )

# 3. Construct figure, add surface, plot ............................

cbval_a = [-1,-.5,0,0.5,1]

fig = plt.figure(figsize=plt.figaspect(1))
fig.text(0.5,.92,info, ha='center', va='top', fontsize='x-large')

ax = plt.axes(projection='3d')
ax.set(xlim=( -3,3 ), ylim=( -3,3 ), zlim=( 0,1 ), 
            zticks=[], xlabel='Re(z)', ylabel='Im(z)' )
ax.view_init(89.9,-90.01)
ax.set_proj_type('ortho')
cbar = plt.colorbar(surface.cBar_ScalarMappable, ax=ax, ticks=cbval_a, shrink=0.6, pad=-.1 )
cbar.ax.set_yticklabels(t_labels)
cbar.set_label(r'$\theta$, angle', rotation=90, labelpad = 0)

sm = plt.cm.ScalarMappable(cmap='binary_r')
sm.set_array([])
cbar2 = fig.colorbar(sm, ax=ax, ticks=r_axvals,shrink=0.7, pad=-.1, orientation='horizontal' )
cbar2.ax.set_xticklabels(r_labels)
cbar2.set_label(r'$\mathcal{R}$, magnitude')

ax.add_collection3d(surface)

# =======================================================================

fig.tight_layout(pad=0)

plt.show()