@acer 

very nice code I wanted exactly this kind of coloring 

I used this article from wikipedia because there is a nice formula for coloring but I have done this only in plot3d like below:

plot3d(abs(x^2 - y^2 + (2*x*y)*I), color = [argument(x^2 - y^2 + (2*x*y)*I), 2*arctan(abs(x^2 - y^2 + (2*x*y)*I))/Pi, 1, colortype = HSV]);

thank you for solving this problem in 2d, do you know maybe what is this colortype function equivalent when we write implicitplot3d?

implicitplot3d(abs(cos(x) + cos(y) + cos(z) + cos(x)*cos(z)*cos(y)*I) = 1, x = -3 .. 3, y = -3 .. 3, z = -3 .. 3, grid = [41, 41, 41], style = surface, color = [sin(x*y), cos(x*y), 0.8, colortype = HSV])

@Kitonum 

I thought about this

https://en.wikipedia.org/wiki/Domain_coloring 

maybe there is some kind of color function in Maple where we can implement code of RGB coloring 

@Carl Love yes, I forgot about it but now it's okay

@acer 

I added two parametrs and I added display to show contours but it doesn't work do you know what is wrong with this modification of your code?

restart;
(m,n) := 200, 200:
(aa,bb,cc,dd) := -3, 3, -3, 3:
M := Array(1..m, 1..n,
           (i,j) -> evalhf( I*(bb-(i-1)*(bb-aa)/(m-1))
                            + (cc+(j-1)*(dd-cc)/(n-1)) ),
           datatype=complex[8], order=Fortran_order):
 

f := (z) -> ((z^2 + a)/(z^3 + b))^3;

P:=proc(MM,param1,param2) local Z,A,temp,img;
  uses ImageTools, plots;
  Z := map[evalhf, inplace](subs(a=param1,b=param2,eval(f)), copy(M)):
  A := Array(1..m, 1..n, 1..3, 1.0, datatype=float[8], order=Fortran_order):
  temp := Array(map[evalhf](z->argument(z), Z), datatype=float[8]):
  A[..,..,1] := map[evalhf,inplace](`*`, FitIntensity(temp,inplace), 360):
  img := HSVtoRGB(A):
  display(plot([[0,0]],x=cc..dd,y=aa..bb,axes=boxed,background=img,
       gridlines=false,size=[400,400]),contourplot(subs(a = param1,b=param2, z = x + y*I, eval(f(z))), x = -3 .. 3, y = -3 .. 3));
end proc:
 

P(M,0.5,2);
 

@Carl Love it's look almost the same (mirror image) but what if function is more complicated you can treat data as variable?

@acer 

restart;
with(plots);
f := z -> z^5 + z^4 - z + 1;
A := plot3d(1, theta = 0 .. 2*Pi, phi = Pi/2 .. Pi, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-f(x + y*I))], axes = none, style = surface, grid = [201, 201]);
B := plot3d(1, phi = 0 .. 2*Pi, theta = 0 .. Pi/2, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-f(I/(y + x*I)))], axes = none, style = surface, grid = [201, 201]);
display(A, B);

g := z -> GAMMA(z - I)^2/GAMMA(z + I)^2;

A1 := plot3d(1, theta = 0 .. 2*Pi, phi = Pi/2 .. Pi, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-g(x + y*I))], axes = none, style = surface, grid = [201, 201]);
B1 := plot3d(1, phi = 0 .. 2*Pi, theta = 0 .. Pi/2, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-g(I/(y + x*I)))], axes = none, style = surface, grid = [201, 201]);
display(A1, B1);
 

@acer 

Thank you, I have done something like this:

S1 := plot3d(1, theta = 0 .. 2*Pi, phi = 0 .. Pi/2, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-g(-I*(y + x*I)))]);
S2 := plot3d(1, theta = 0 .. 2*Pi, phi = Pi/2 .. Pi, coords = spherical, colorscheme = ["xyzcoloring", (x, y, z) -> argument(-f(x + y*I))]);
display(S1, S2, orientation = [-10, 57, 52]);
 

and now I see two problems: first is quality of rendering this color function, second is scaling and translate the S2 because rotation is ok and easy to implement we just multiply by some complex number. Do you have any idea how to improve quality and scaling? 

I have also found https://www.wolfram.com/language/12/complex-visualization/riemann-sphere.html?product=language nice post which is about this idea.

@acer 

I used code from https://www.mapleprimes.com/questions/226790-Is-There-Any-Maple-Code-For-The-domain

restart;
with(ImageTools);
m, n := 600, 600;
                        m, n := 600, 600

a, b, c, d := -1.2, 1.2, -1.2, 1.2;
               a, b, c, d := -1.2, 1.2, -1.2, 1.2

M := Array(1 .. m, 1 .. n, (i, j) -> evalhf((b - (i - 1)*(b - a)/(m - 1))*I + c + (j - 1)*(d - c)/(n - 1)), datatype = complex[8], order = Fortran_order);
f := z -> z^5 - z^4 - z + 1;
  f := proc (z) options operator, arrow; z^5-z^4-z+1 end proc

Z := map[evalc, inplace](f, copy(M));
A := Array(1 .. m, 1 .. n, 1 .. 3, 1.0, datatype = float[8], order = Fortran_order);
temp := Array(map[evalhf](z -> argument(-z), Z), datatype = float[8]);
A[() .. (), () .. (), 1] := map[evalhf, inplace](`*`, FitIntensity(temp, inplace), 360);
img := HSVtoRGB(A);
G := plot([[0, 0]], x = c .. d, y = a .. b, axes = boxed, background = img, gridlines = false);
 

mapping 1/z turn inside out this phase portrait so function 1/z^5-1/z^4+1/z+1 give you (after rotation) second phase portrait. 

these plots represent stratographic projection so when you have sphere with upper hemisphere with this phase plot and you project this to plane you get this image but in this situatuion we want to inverse procedure.

@acer 

is it possible to stretch this part of image inside the circle on the hemisphere? Because that mapping transform whole square but I guess it should be only this image inside the circle.

@dharr 

I added lacking part of this equation but after "evaluating" nothing happens

@Carl Love with(plots);
a := -Pi;
b := Pi;
implicitplot3d(1/2*(sin(2*x)*cos(y)*sin(z) + sin(2*y)*cos(z)*sin(x) + sin(2*z)*cos(x)*sin(y)) - 1/2*(cos(2*x)*cos(2*y) + cos(2*y)*cos(2*z) + cos(2*z)*cos(2*x)) + 0.15 = 0, x = a .. b, y = a .. b, z = 1/2*a .. 3/2*b, grid = [70, 70, 70], style = patchnogrid, shading = xy, lightmodel = light2);

only change style to patch

@vv Thank you a lot, this is really intresting solution. I like that it's not only for this but for much larger familiy of piecewise periodic functions.

@acer 

I make contours by contourplot and take absolute value and argument to make these lines next I put all outputs in display.

contourplot(abs((x + y*I - 1)/((x + y*I)^2 + x + y*I + 1)), x = -2 .. 2, y = -2 .. 2, contours = [1/64, 1/32, 1/16, 1/8, 1/4, 1/2, 1, 3/2, 2, 3, 4, 5, 6, 7], numpoints = 10000, color = black, thickness = 1/2)

contourplot(argument((x + y*I - 1)/((x + y*I)^2 + x + y*I + 1)), x = -2 .. 2, y = -2 .. 2, grid = [500, 500], color = black, thickness = 1/2)

if you have easier way to do this contours you may tell me, the worst part is about counting argument because lines are not smooth so I have to make large grid.

@acer 

Your code is exactly what I was looking for. Below I have your output. So you have just put minus sign before function. I didn't think that could help. Thank you once again.