Here is a version that has my rudimentary attempt to add a conformal map.
Code:
// SetHSV() and GetColor() based on a
// c++ program from :
// http://commons.wikimedia.org/wiki/File:Color_complex_plot.jpg
// by Claudio Rocchini
// http://en.wikipedia.org/wiki/Domain_coloring
PlotRes(r);
SetHSV(h,s,v);
GetColor(v);
EXPORT Plot()
BEGIN
PRINT();
RECT();
PlotRes(16);
PlotRes(8);
PlotRes(4);
PlotRes(2);
PlotRes(1);
WHILE 1 DO
FREEZE();
END;
END;
EXPORT PlotRes(r)
BEGIN
LOCAL xi,yi,co,va;
xi:=(Xmax-Xmin)/320;
yi:=(Ymax-Ymin)/240;
FOR X FROM 0 TO 320-r STEP r DO
FOR Y FROM 0 TO 240-r STEP r DO
IF ((Xmin+xi*X)==0) THEN
//avoid typical undefined points at origin
va:=F1(Xmin+xi*(X+0.001) + i*(Ymin+yi*Y));
ELSE
va:=F1(Xmin+xi*X + i*(Ymin+yi*Y));
END;
//LOG coloring
//co:=RGB(MIN(255,LN(1+ABS(IM(va)))*128),MIN(255,LN(1+ABS(RE(va)))*128),MIN(255,MAX(0,(LN(1+ABS(va))-2)*128)));
//Domain coloring
co:= GetColor(va);
RECT_P(G0,X,Y,X+r-1,Y+r-1,co);
END;
END;
END;
SetHSV(h, s, v)
BEGIN
LOCAL r, g, b;
LOCAL z, f, p, q, t, i;
IF(s==0) THEN
r:=v;
g:=v;
b:=v;
ELSE
IF(h==1) THEN h := 0; END;
z := FLOOR(h*6);
i := IP(z);
f := h*6 - z;
p := v*(1-s);
q := v*(1-s*f);
t := v*(1-s*(1-f));
CASE
IF i==0 THEN r:=v; g:=t; b:=p; END;
IF i==1 THEN r:=q; g:=v; b:=p; END;
IF i==2 THEN r:=p; g:=v; b:=t; END;
IF i==3 THEN r:=p; g:=q; b:=v; END;
IF i==4 THEN r:=t; g:=p; b:=v; END;
IF i==5 THEN r:=v; g:=p; b:=q; END;
END;
END;
r :=MIN(255,IP(256*r));
g :=MIN(255,IP(256*g));
b :=MIN(255,IP(256*b));
RETURN RGB(r,g,b);
END;
GetColor(v)
BEGIN
LOCAL a:=0;
LOCAL m,ranges,rangee,k,sat,val;
IF v≠0 THEN a:=ARG(v); END;
WHILE (a<0) DO a := a+ (2*π); END;
a := a/(2*π);
// RE Conformal mapping
m := ABS(RE(v));
ranges := 0;
rangee := 1;
WHILE(m>rangee) DO
ranges := rangee;
rangee := rangee * e;
END;
k:=(m-ranges)/(rangee-ranges);
IF (k<0.5) THEN
sat:=k*2;
ELSE
sat:=1 -(k -0.5) *2;
END;
val := sat;
sat := 1 - (1-sat)^3; sat := 0.4 + sat*0.6;
val := 1 - val;
val := 1 - (1-val)^3;
val := 0.6 + val*0.4;
IF (val > 0.9999) OR (sat >0.9999) THEN return SetHSV(a,sat,val); END;
//IM Conformal mapping
m := ABS(IM(v));
ranges := 0;
rangee := 1;
WHILE(m>rangee) DO
ranges := rangee;
rangee := rangee * e;
END;
k:=(m-ranges)/(rangee-ranges);
IF (k<0.5) THEN
sat:=k*2;
ELSE
sat:=1 -(k -0.5) *2;
END;
val := sat;
sat := 1 - (1-sat)^3; sat := 0.4 + sat*0.6;
val := 1 - val;
val := 1 - (1-val)^3;
val := 0.6 + val*0.4;
IF (val > 0.9999) OR (sat >0.9999) THEN return SetHSV(a,sat,val); END;
//Domain Coloring
m := ABS(v);
ranges := 0;
rangee := 1;
WHILE(m>rangee) DO
ranges := rangee;
rangee := rangee * e;
END;
k:=(m-ranges)/(rangee-ranges);
IF (k<0.5) THEN
sat:=k*2;
ELSE
sat:=1 -(k -0.5) *2;
END;
val := sat;
sat := 1 - (1-sat)^3; sat := 0.4 + sat*0.6;
val := 1 - val;
val := 1 - (1-val)^3;
val := 0.6 + val*0.4;
return SetHSV(a,sat,val);
END;
Here is the example plot of SIN(X) (reminder: X treated like Z, which is x + iy) :
One more example for (X^2-1)*(X-2-i)^2/(X^2+2+2*i):
As a bonus, here is attempt at a plot of Julia set for fc=Z^2+c, c=-0.8+0.156i. I plotted it by setting F2=X^2-0.8+0.156*i, and F1=F2(F2(F2(F2(.....(F2(X))...)))). This is not optimized for fractals at all so plotting this is quite slow, but given enough time is doable: