%let name=mandelbrot;

/* 
Set your current-working-directory (to read/write files), if you need to ...
%let rc=%sysfunc(dlgcdir('c:\someplace\public_html')); 
*/
filename odsout '.';

/*
SAS/Graph version of
http://rosettacode.org/wiki/Mandelbrot_set#Fortran
*/

%let x_min=-2.5;
%let x_max=1.5;
%let y_min=-1.5;
%let y_max=1.5;

%let i_points=800;
%let j_points=600;

data mandelbrot (keep = j i colorvar);

width=&x_max-&x_min;
height=&y_max-&y_min;
x_centre=&x_min+(width/2);
y_centre=&y_min+(height/2);

/*
n_max is the number of levels of color to use.
You'll need to assign a color for each one in your annotated points.
*/
n_max=41;

/*
Although the plot looks like smooth/continuous color, it is really
composed of lots of individual 'points' - if you use enough points,
and pack them together dense enough, then they will look like solid
colors. (i = x direction, and j = y direction)
*/
i_points=&i_points;
j_points=&j_points;

dx_di=width/i_points;
dy_dj=height/j_points;
x_offset=x_centre - 0.5*(i_points+1)*dx_di;
y_offset=y_centre - 0.5*(j_points+1)*dy_dj;

do j = 1 to j_points;
 y_0 = y_offset + dy_dj * j;

 do i = 1 to i_points;
  x_0 = x_offset + dx_di * i;
  x=0.0;
  y=0.0;
  n=0;

  do while (1=1);
   x_sqr=x**2;
   y_sqr=y**2;
   /* outside */
   if (x_sqr+y_sqr > 4.0) then do;
    colorvar=(n_max-n)+1;
    output;
    leave;
   end;
   /* inside */
   if (n eq n_max) then do;
    colorvar=(n_max-n)+1;
    output;
    leave;
   end;
   y=y_0 + 2.0*x*y;
   x=x_0 + x_sqr - y_sqr;
   n+1;
  end;

 end;
end;
run;



ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" 
 (title="Mandelbrot Plot") 
 style=htmlblue;

ods graphics / 
 noscale /* if you don't use this option, the text will be resized */
 imagefmt=png imagename="&name"
 width=800px height=600px noborder; 

title1 c=gray33 h=20pt "Mandelbrot Plot w/ SAS Proc SGplot";

proc sgplot data=mandelbrot;
scatter y=j x=i / markerattrs=(symbol=circlefilled size=1px)
 colorresponse=colorvar
 colormodel=(black white);
yaxis offsetmin=0 offsetmax=0 labelpos=top;
xaxis offsetmin=0 offsetmax=0;
run;

proc sgplot data=mandelbrot;
scatter y=j x=i / markerattrs=(symbol=circlefilled size=1px)
 colorresponse=colorvar 
 colormodel=(
   black
   CXCC99FF CX9966CC CX663399 CX003366
   CXCC99FF CX9966CC CX663399 CX003366
   CXCC99FF CX9966CC CX663399 CX003366
   CXCC99FF CX9966CC CX663399 CX003366
   CXCC99FF CX9966CC CX663399 CX003366
   CXFF99CC CXCC6699 CX993366 CX660033
   CXFFCC99 CXCC9966 CX996633 CX663300
   CXCCFF99 CX99CC66 CX669933 CX336600
   CX99FFCC CX66CC99 CX339966 CX006633
   CX99CCFF CX6699CC CX336699 CX003366
   );
yaxis offsetmin=0 offsetmax=0 labelpos=top;
xaxis offsetmin=0 offsetmax=0;
run;

ods graphics / width=550px height=550px; 

proc sgplot data=mandelbrot noautolegend;
scatter y=j x=i / markerattrs=(symbol=circlefilled size=1px)
 colorresponse=colorvar
 colormodel=(black white);
yaxis values=(50 to 550 by 50) offsetmin=0 offsetmax=0 display=none;
xaxis values=(100 to 650 by 50) offsetmin=0 offsetmax=0 display=none;
run;

title;
proc print data=mandelbrot (obs=50) noobs
 style(data)={font_size=11pt}
 style(header)={font_size=11pt};
run;

quit;
ODS HTML CLOSE;
ODS LISTING;