Complex Plotting App

[danielmewes]

I was interested in exploring the complex roots of polynomials a while back.
Unfortunately the Advanced Graphing app doesn't like imaginary numbers in its equations ( http://h30434.www3.hp.com/t5/Other-HP-Co...-p/4171298 ), so I was left with writing my own program for visualizing such functions.

Specifically, given a function f(x) : C -> C, I wanted to get a two dimensional plot over some domain [a, b] + i[c, d] where different colors correspond to different output values.

The forum didn't allow me to attach the hpapp file, but you can download it here:
http://dmewes.com/~daniel/Complex%20Plot.hpapp

To use it:
1. Go to Symb, enter F1(X)= (for example X^3-1)
2. Optionally: Go to Plot->Setup to set the boundaries for the real part (X Rng) and imaginary part (Y Rng) of the parameter.
3. Press Plot to start plotting
4. Press On to leave the plot

A few additional notes:
* The X axis corresponds to the real part of the input parameter, the Y axis to the imaginary part.
* The current shader uses a logarithmic scale, and ignores the sign of the real and imaginary parts. You can change the line `co:=RGB...` in the code to use a different shader.
* With the default shader, zero appears as black, (mostly) imaginary numbers appear red, and (mostly) real ones appear green. Results with both an imaginary and real part get shades of yellow.
* Because plotting every pixel is relatively slow, the app starts by plotting a low-resolution plot, and then refines it down to pixel resolution if you keep it running.


I've also attached a photo and a screenshot showing the plots of
1. X^5 - 1
2. sin(0.5*X^3 + X) + cos(X^2)
   


Program code (replace i with the proper imaginary i symbol in `va:=F1...`):

Code:


PlotRes(r);
EXPORT Plot()
BEGIN
 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
   va:=F1(Xmin+xi*X + i*(Ymin+yi*Y));
   co:=RGB(MIN(255,LN(1+ABS(IM(va)))*128),MIN(255,LN(1+ABS(RE(va)))*128),MIN(2​55,MAX(0,(LN(1+ABS(va))-2)*128)));
   RECT_P(G0,X,Y,X+r-1,Y+r-1,co);
  END;
 END;
END;

[dwgg]

Thank you for this program. I was also looking into plotting complex functions, and one other option I found is to use the 3D contour plot feature of the Graph3D program, using RE() and IM(), and setting the rotation angle to x=90, y=0, and z=0.. however your program looks better in plotting complex functions.

I also modified and added Domain Coloring to your program, by porting the code in http://commons.wikimedia.org/wiki/File:C...x_plot.jpg. I think Domain Coloring could be more useful and it also looks prettier.

Here is a screenshot of domain coloring plot showing sin(0.5*X^3 + X) + cos(X^2):
   

Below is the modified code:

Code:


// SetHSV() and GetColor() based on
// 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,MA​X(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*π);
    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;
 
    sat := 1 - (1-sat)^3; sat := 0.4 + sat*0.6;
 
    IF (k<0.5) THEN
        val:=k*2;
    ELSE
        val:=1 -(k -0.5) *2;
    END;

    val := 1 - val;
    val := 1 - (1-val)^3; 
    val := 0.6 + val*0.4;

    return SetHSV(a,sat,val);
END;

[Han]

Both versions are very cool! As for using Graph3D to do complex plots...

Graph3D can be modified very easily to do 3D complex plots. The main question is what schemes are people most interested in. For example, one might be to use the x-axis and y-axis for the input z=x+yi and z-axis for Re(f(z)) while the color of the point is defined (in some way) by Im(f(z)).

Another option may be to toggle between showing Re(f(z)) and Im(f(z)) on the z-axis. As already suggested, defining V1 as f(z), V2 as Re(f(z)) and V3 as Im(f(z)) and plotting only V2 and V3 simultaneously would be a quick way to do this.

Are there other schemes that are more preferrable/common? Are there any advantages to using 3D complex graphs?

[dwgg]

Hi Han,

With Graph3D, I am already using something similar to what you're suggesting as second option ("Re(f(z)) and Im(f(z)) on the z-axis"). I did minor modification on the Graph3D source code ( namely, disabling function select for 3D contour, then modifying g3D_drawcontour to allow drawing contour plots for multiple functions). What I get is something similar to Conformal map. Here is an example for SIN(X+i*Y):
   

And here is Domain Coloring plot using the same function (entering SIN(X) with my above code):
   

I think what will be good to have is putting these two together, so we get something like below (example from https://gandhiviswanathan.wordpress.com/...functions/):
   

Then we can examine zeros, poles, and conformality all together on the same plot (as well as making the graph look nice).

[Mark Hardman]

(02-01-2015 06:42 AM)dwgg Wrote:  

I already hate myself for saying this but...

These are perfect for Super Bowl Sunday!

[dwgg]

Only if this plot could have persuaded the Seahawks to do running play instead of doing passing on that last one! (all discussions regarding Superbowl XLIX will ultimately lead to this...)

[rprosperi]

(02-02-2015 03:43 PM)dwgg Wrote:  Only if this plot could have persuaded the Seahawks to do running play instead of doing passing on that last one! (all discussions regarding Superbowl XLIX will ultimately lead to this...)

A perfect play IMHO. And yes, I am from the Boston area...

How on earth did they make that call?! At least the head coach admitted it was his mistake.

[dwgg]

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,MA​X(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:
   

[danielmewes]

Pretty cool extensions dwgg

[dwgg]

A correction regarding "conformal mapping" I mentioned above.. it seems conformal map refers to something else entirely (mapping one space to another while preserving angle) and what I added on is simply drawing contours of function.. sorry for the confusion. I don't know a lot about complex analysis and got confused since many pages that talks about conformal map had domain coloring graph with contours.

BTW, here is a minor bug fix to put the Y axis on correct direction (positive on the top).

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*(240-r-Y)));
   ELSE
    va:=F1(Xmin+xi*X + i*(Ymin+yi*(240-r-Y)));
   END;

//LOG coloring
   //co:=RGB(MIN(255,LN(1+ABS(IM(va)))*128),MIN(255,LN(1+ABS(RE(va)))*128),MIN(255,MA​X(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 Contour
    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 Contour
    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;

Also, here is an animation I made of ((Z^2+c)*(Z-2-i)^2)/(X^2+2+2i) from c=-16 to c=9

https://www.youtube.com/watch?v=_HP7MA7UP3E

[Han]

Very, very cool!

I was just thinking about how you could speed up the program. There is a lot of redundancy in the program. In fact, calling PlotRes(n) for any value n>1 is unnecessary as far as the final plot is concerned. I understand that they are there to give users a sense the program is refining the plot over time.

At each new resolution, we only need to update 3/4 of each grid section from the previous draw. The current version redraws all 4 divisions.

I wonder if we can improve the speed by making use of this observation.

[Han]

Usage: Plot(r) where the initial resolution is a grid of squares of size \( 2^r \times 2^r \). On the emulator, this modified version takes half the time vs. the original. On the calculator, the modified version takes 75% of the time of the original.

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

SetHSV(h,s,v);
GetColor(v);
EvalF();

EXPORT Plot(r)
BEGIN
  local x1,x2,y1,y2,co;
  local dx:=(Xmax-Xmin)/320;
  local dy:=(Ymax-Ymin)/240;
  local z1;
  local a,b,d,k,x,y;

  d:=2^r;
  FOR x FROM 0 TO 320-d STEP d DO
    FOR y FROM 0 TO 240-d STEP d DO
      z1:=Xmin+x*dx+i*(Ymin+y*dy);
      co:=EvalF(z1);
      RECT_P(G0,x,y,x+d-1,y+d-1,co); 
    END;
  END;

IF r THEN
  FOR k FROM 1 TO r DO
    d:=2^(r-k);
    FOR x FROM 0 TO 160/d-1 DO
      FOR y FROM 0 TO 120/d-1 DO
        a:=x*2*d; b:=y*2*d+d;
        z1:=Xmin+a*dx+i*(Ymin+b*dy);
        co:=EvalF(z1);
        RECT_P(G0,a,b,a+d-1,b+d-1,co);
        z1:=z1+d*dx;
        co:=EvalF(z1);
        RECT_P(G0,a+d,b,a+2*d-1,b+d-1,co);
        z1:=z1-i*d*dy;
        co:=EvalF(z1);
        RECT_P(G0,a+d,b-d,a+2*d-1,b-1,co);
      END;
    END;
  END;
END;

  FREEZE;
//  WAIT(-1);
END;

EvalF(z)
BEGIN
  IF RE(z) THEN
    RETURN(GetColor(F1(z)));
  ELSE
    RETURN(GetColor(F1(z+.001)));
  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;

[rprosperi]

(02-07-2015 06:02 AM)Han Wrote:  Usage: Plot(r) where the initial resolution is a grid of squares of size \( 2^r \times 2^r \). On the emulator, this modified version takes half the time vs. the original. On the calculator, the modified version takes 75% of the time of the original.

Han - Though I admit to understanding very little about this topic, visually it's very cool, so I thought I would try it, but can't get it to run.

I copy/paste it into the emulator and it compiles with no issues, but when I run it (from Home), I get "Error: No definition in Symbolic view".

Which App must be current to run this (I assumed Advanced Graphing) and what is the proper way to invoke it?

[Han]

Define a complex function but in terms of X in the function app.

[rprosperi]

(02-07-2015 06:24 PM)Han Wrote:  Define a complex function but in terms of X in the function app.

Thanks, that helped a lot; interesting to explore some very vivid and colorful graphics. The instructions are plain as day in the initial post; sorry for not searching better before posting.

Anyhow, thanks for teaching me some more Prime in an entertaining way.

[Eddie W. Shore]

All I got from this link: http://dmewes.com/~daniel/Complex%20Plot.hpapp is a blank app. Am I missing something? Thanks!

[rprosperi]

(02-12-2015 03:26 AM)Eddie W. Shore Wrote:  All I got from this link: http://dmewes.com/~daniel/Complex%20Plot.hpapp is a blank app. Am I missing something? Thanks!

Copy the code in Post #12 above and paste into Emulator (and then download to device per normal, etc.).

[dwgg]

(02-07-2015 06:02 AM)Han Wrote:  Usage: Plot(r) where the initial resolution is a grid of squares of size \( 2^r \times 2^r \). On the emulator, this modified version takes half the time vs. the original. On the calculator, the modified version takes 75% of the time of the original.

Thank you Han for the improvement. I tried the modification and it does get rid of redundancy and improve the performance of the program. Moreover, since we are only drawing uncalculated pixels there is almost no time penalty in doing progressive resolution compared to doing pixel level single sweep.

Also adding one comment regarding conformality (I am learning more about complex analysis as I am further working on this program), it looks like you can indeed check conformality by checking contour grids: Conformality breaks where the contour is breaking up. You can read further in below link:

http://www.mai.liu.se/~hanlu09/complex/d...icode.html

Comparing with complex graphs I see from above link, I think now that in terms of graph features this program is meeting many of what one would look for in complex graph.. I think what we could improve further on is as below:

1. Any further speed improvement (if possible)
2. Capability for saving previous calculated plots and using it to redraw if function is not modifed (like Graph3D... so we don't have to keep redrawing when we come back)
3. Tracing feature (based on re and im coordinate print out z value)
4. implement option for edge detection filter for the contour lines so the contour lines are of uniform width (see further info in above link).

[Han]

(02-12-2015 06:29 AM)dwgg Wrote:  

(02-07-2015 06:02 AM)Han Wrote:  Usage: Plot(r) where the initial resolution is a grid of squares of size \( 2^r \times 2^r \). On the emulator, this modified version takes half the time vs. the original. On the calculator, the modified version takes 75% of the time of the original.

Thank you Han for the improvement. I tried the modification and it does get rid of redundancy and improve the performance of the program. Moreover, since we are only drawing uncalculated pixels there is almost no time penalty in doing progressive resolution compared to doing pixel level single sweep.

Also adding one comment regarding conformality (I am learning more about complex analysis as I am further working on this program), it looks like you can indeed check conformality by checking contour grids: Conformality breaks where the contour is breaking up. You can read further in below link:

http://www.mai.liu.se/~hanlu09/complex/d...icode.html

Comparing with complex graphs I see from above link, I think now that in terms of graph features this program is meeting many of what one would look for in complex graph.. I think what we could improve further on is as below:

1. Any further speed improvement (if possible)
2. Capability for saving previous calculated plots and using it to redraw if function is not modifed (like Graph3D... so we don't have to keep redrawing when we come back)
3. Tracing feature (based on re and im coordinate print out z value)
4. implement option for edge detection filter for the contour lines so the contour lines are of uniform width (see further info in above link).

Up late grading; thought I'd respond before bed…

1. I wonder if PIXON_P (as opposed to RECT_P) would be any faster for r=1 since we're only drawing a single pixel (fewer flops for sure since some arithmetic is done for computing the dimensions of the colored squares.

2. This can be done; not sure about data set sizes (in terms of storage) but you could perhaps just store the evaluation of the function at x+yi. The color can be calculated (slower but less storage needed -- up to you). You can even predefine an array of size 320x240 and simply save the output values in entry (a,b) -- see the variables in the loop -- as those are also the top left corners of each colored square.

3. Shouldn't be hard once #2 is added on.

4. This can be done using a dynamic list of lists. Each list within the list would correspond to a single contour level. You can expand each contour level list with concat; same goes for the container list. (Graph3D uses some similar for its cross-sections shown during trace mode). The idea I have in mind may be way to slow, though. That is, each new contour point is compared with existing points within its contour-level list. You would need to design an algorithm to decide whether or not a new contour point gets added based on the currently existing list of contour points (perhaps a simple distance calculation to see if the new point would create a new "vertex" for the contour line as opposed to simply making the contour line thicker). (I have not read the link so this may be a very poor idea to actually implement.)

Edit for #4: Or perhaps just compute the image of the lines x+0i and 0+yi under the function f and use LINE_P.

#5. Add a feature for different "color" schemes (I actually find the yellow -> red -> black to be quite pleasing to my eyes) and even allow users to use a 320x240 image for the mapping (as opposed to colors).

[dwgg]

Color scheme change (Red to Yellow for now by limiting hue space range to 1/6 the original.. for black to red to yellow I think we have to define some other color space other than HSV which I am using now) plus some minor modifications:

1. define r as local in Plot() since I prefer overriding the plot key as was in the original code by danielmewes.
2. Flip Y axis
3. tune rangee to improve contours and logarithmic abs value tracking mask (I think both r and initial rangeee values should be configurable.. later on maybe we will implement setup screen)
4. also fixed RE contour to grey and IM contour to white for this version.

Example plot of (X+2)^2*(X-1-2*i)*(X+i), in the range +/- 4, +/- 3i:
   

Code below:

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

SetHSV(h,s,v);
GetColor(v);
EvalF();

EXPORT Plot()
BEGIN
  local x1,x2,y1,y2,co;
  local dx:=(Xmax-Xmin)/320;
  local dy:=(Ymax-Ymin)/240;
  local z1;
  local a,b,b1,d,k,x,y;
  local r:=4;  //r is size of initial square which is 2^r by 2^r pixels

  d:=2^r;
  FOR x FROM 0 TO 320-d STEP d DO
    FOR y FROM 0 TO 240-d STEP d DO
      z1:=Xmin+x*dx+i*(Ymin+(240-d-y)*dy);
      co:=EvalF(z1);
      RECT_P(G0,x,y,x+d-1,y+d-1,co); 
    END;
  END;
IF r THEN
  FOR k FROM 1 TO r DO
    d:=2^(r-k);
    FOR x FROM 0 TO 160/d-1 DO
      FOR y FROM 0 TO 120/d-1 DO
        a:=x*2*d; b:=y*2*d; b1:=(120/d-1-y)*2*d+d;
        z1:=Xmin+a*dx+i*(Ymin+b1*dy);
        co:=EvalF(z1);
        RECT_P(G0,a,b,a+d-1,b+d-1,co);
        z1:=z1+d*dx;
        co:=EvalF(z1);
        RECT_P(G0,a+d,b,a+2*d-1,b+d-1,co);
        z1:=z1+i*d*dy;
        co:=EvalF(z1);
        RECT_P(G0,a+d,b-d,a+2*d-1,b-1,co);
      END;
    END;
  END;
END;

  FREEZE;
  WAIT(-1);
END;

EvalF(z)
BEGIN
  IF RE(z) THEN
    RETURN(GetColor(F1(z)));
  ELSE
    RETURN(GetColor(F1(z+.001)));
  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 Contour
    m := ABS(RE(v));
    ranges := 0;
    rangee := 0.5;

    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/6,sat,val);
      return RGB(200,200,200); 
    END;

//IM Contour
    m := ABS(IM(v));
    ranges := 0;
    rangee := 0.5;

    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/6,sat,val);
      return RGB(255,255,255); 
    END;

//Domain Coloring
    m := ABS(v);
    ranges := 0;
    rangee := 0.01;

    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/6,sat,val);
END;