%let name=mandelbrot;
filename odsout '.';

/*
From example at the end of Robert Ray's TSSQL/TSPL example.
/u/sasrsr/expo05/frac.sas

Also used in his TSSQL/TSPL training video:
http://video.na.sas.com/Video/video_flash.aspx?index=sas9_2_training_session_sept08_1
*/


/* Setup table that contains one row for each coordinate to compute */
proc tspl;
  table inp(overwrite=yes);
    dcl double p q r;
    dcl integer maxiterate;

    method init();
      dcl int n m;
      dcl int i j k;
      dcl double pmin pmax qmin qmax;

      n = 1024;
      m = 1024;
      pmin = -1.5; pmax = -0.5;
      qmin = -0.5; qmax =  0.5;

      /* A different location on the fractal */
      /***
      pmin = -1.336; pmax = -1.0625;
      qmin =  0.110; qmax =  0.5;
      ***/

      r = 100.0;
      maxiterate = 50;

      do k = 1 to n*m;
        i = k/m;
        j = mod(k,m);
        p = i*(pmax-pmin)/(n-1)+pmin;
        q = j*(qmax-qmin)/(m-1)+qmin;
        output;
      end;
    end;
  run;
quit;

/* Create Thread Program */
proc tspl;
  thread fracthread;
    dcl int mesh;
    dcl double x y rr nx ny;

    keep p q mesh;

    method run();
      set inp;
   
      x = p;
      y = q;
      rr = r**2;

      mesh = 0;
      do while (mesh < maxiterate and (x**2+y**2 < rr));
        nx = x**2 - y**2 + p;
        ny = 2.0*x*y + q;
        x = nx;
        y = ny;
        mesh = mesh+1;
      end;
    end;
  run;
quit;

/* Execute Thread Program */
proc tspl;
  table mandelbrot(overwrite=yes);
    dcl thread fracthread frac;
    method run();
      set from frac threads=8;
    end;
  run;
quit;

/* Colors for the contour plot */
pattern1  v=solid c=CX003366;
pattern2  v=solid c=CX336699;
pattern3  v=solid c=CX6699CC;
pattern4  v=solid c=CX99CCFF;
pattern5  v=solid c=CX006633;
pattern6  v=solid c=CX339966;
pattern7  v=solid c=CX66CC99;
pattern8  v=solid c=CX99FFCC;
pattern9  v=solid c=CX336600;
pattern10 v=solid c=CX669933;
pattern11 v=solid c=CX99CC66;
pattern12 v=solid c=CXCCFF99;
pattern13 v=solid c=CX663300;
pattern14 v=solid c=CX996633;
pattern15 v=solid c=CXCC9966;
pattern16 v=solid c=CXFFCC99;
pattern17 v=solid c=CX660033;
pattern18 v=solid c=CX993366;
pattern19 v=solid c=CXCC6699;
pattern20 v=solid c=CXFF99CC;
pattern21 v=solid c=CX003366;
pattern22 v=solid c=CX663399;
pattern23 v=solid c=CX9966CC;
pattern24 v=solid c=CXCC99FF;
pattern25 v=solid c=CX003366;
pattern26 v=solid c=CX663399;
pattern27 v=solid c=CX9966CC;
pattern28 v=solid c=CXCC99FF;
pattern29 v=solid c=CX003366;
pattern30 v=solid c=CX663399;
pattern31 v=solid c=CX9966CC;
pattern32 v=solid c=CXCC99FF;
pattern33 v=solid c=CX003366;
pattern34 v=solid c=CX663399;
pattern35 v=solid c=CX9966CC;
pattern36 v=solid c=CXCC99FF;
pattern37 v=solid c=CX003366;
pattern38 v=solid c=CX663399;
pattern39 v=solid c=CX9966CC;
pattern40 v=solid c=CXCC99FF;
pattern41 v=solid c=black;


goptions device=png;
goptions noborder;
 
ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" (title="SAS/Graph Mandelbrot Plot (proc gcontour)") style=minimal;

goptions gunit=pct htitle=4 htext=2.5 ftitle="arial/bold" ftext="arial";

title "Mandelbrot Plot";
proc gcontour data=mandelbrot;
  plot q*p=mesh /
    nolegend
    pattern
    join
    levels = 5 to 45
 des="" name="&name";  
run; 

quit;
ODS HTML CLOSE;
ODS LISTING;