/*********************************************/
/* SAS program to accompany
   "Animate snowfall in SAS" by Rick Wicklin
   Published 14DEC2016 on The DO Loop blog
   http://blogs.sas.com/content/iml/2016/12/14/animate-snow-sas.html
*/
/* Define routines to construct a Koch snowflake.
   See  http://blogs.sas.com/content/iml/2016/12/12/koch-snowflake.html
*/
title;footnote;
proc iml;
/* divide a line segment (2 pts) into 4 segments (5 pts).  
   Create middle pt by rotating through -60 degrees */ 
start KochDivide(A, E);
   segs = j(5, 2, .);                      /* allocate matrix to hold 4 shorter segs */
   v = (E-A) / 3;                          /* vector 1/3 as long as orig */
   segs[{1 2 4 5}, ] = A + v @ T(0:3);     /* endpts of new segs */
   /* Now compute middle point. Use ATAN2 to find direction angle. */
   theta = -constant("pi")/3 + atan2(v[2], v[1]); /* change direction angle by pi/3 */
   w = cos(theta) || sin(theta);           /* vector to middle point */
   segs[3,] = segs[2,] + norm(v)*w;
   return segs;
finish;

start KochPoly(P0, iters=5);
P = P0;
do j=1 to iters;
   N = nrow(P) - 1;        /* old number of segments */
   newP = j(4*N+1, 2);     /* new number of segments + 1 */
   do i=1 to N;            /* for each segment */
      idx = (4*i-3):(4*i+1);    /* rows for new segments */
      newP[idx, ] = KochDivide(P[i,], P[i+1,]); /* generate new segments */
   end;
   P = newP;               /* update polygon and iterate */
end;
return P;
finish;
store module=(KochPoly KochDivide);
QUIT;

/* Main Program:
   Create frames for animated GIF that has falling (Koch) snowflakes. 
*/
%let numFrames = 30;                /* number of images in GIF */

proc iml;
pi = constant('pi');
call randseed(123);
load module=(KochPoly KochDivide);
/* create equilateral triangle as base for snowflake */
angles = -pi/6 // pi/2 // 7*pi/6;   /* angles for equilateral triangle */
Tri = cos(angles) || sin(angles);   /* vertices of equilateral triangle */
Tri = Tri // Tri[1,];               /* close the polygon */
P = KochPoly(Tri, 3);               /* create the Koch snowflake */

/* define points to use as centers of the snowflakes */
Nx = 3; Ny = 3;
center = expandgrid( ((1:Nx)-0.5)/ Nx, ((1:Ny)-0.5)/ Ny) 
         + 0.25*(randfun(Nx*Ny // 2, "Uniform")-0.5);
/* duplicate centers and translate vertically so they are initially off screen */
center = center // (center + {0 1});

/* constants and matrices used in loop */
D = diag({0.1 0.1});       /* scaling matrix for snowflakes */
R = j(2,2);                /* allocate rotation matrix */
maxT = &numFrames;         /* number of time points in the sequence */
dt = 2*pi/maxT;            /* step size is multiple of 2*pi for periodic functions */
results = j(nrow(P), 4, .);
create Anim from results[c={"Frame" ID X Y}];
do i = 1 to maxT;
   results[,1] = i;        /* BY variable */
   t = i*dt;               /* current time */
   Xoffset = 0.02*sin(t);  /* small periodic mation in X (gentle wind) */
   Yoffset = -i/maxT;      /* steady falling in Y */
   c = center + (Xoffset || Yoffset);
   a = 3*t/6;              /* rotation angle multiple of pi/6 */
   R[1,1] = cos(a);  R[1,2] = -sin(a);  R[2,1] = sin(a);  R[2,2] = cos(a);
   glyph = P * D * R;      /* scale and rotate snowflake  */
   do j = 1 to nrow(center);          /* for each center */
      results[,2] = j;                /* ID variable */
      results[,3:4] = c[j,] + glyph;  /* translate snowflake */
      append from results;
   end;
end;
close;
QUIT;

/* add "snowy ground" polygon */
data BG;
bgID = "Snow"; 
do Frame  = 1 to &numFrames;
   bgX = 0; bgY = 0; output;
   bgX = 1; bgY = 0; output;
   bgX = 1; bgY = 0.33; output;
   bgX = 0; bgY = 0.33; output;
end;
run;

/* add message of greeting */
data Msg;
do Frame  = 1 to &numFrames;
   xc = 0.5; yc=0.5;
   txt = "Happy Holidays|Peace on Earth"; /* use '|' as split character */
   output;
end;
run;

/* Merge. Each BY-group level has to have a copy of the ground and text */
data All;
merge Anim BG Msg;
by Frame;
if first.Frame then cnt=0;
cnt + 1;
if cnt>4 then do; bgX=.; bgY=.; end;
if cnt>1 then do; xc=.; yc=.; txt=""; end;
drop cnt;
run;

/* set up ODS to write animated GIF:
   http://blogs.sas.com/content/iml/2016/08/22/animation-by-statement-proc-sgplot.html
*/
ods graphics / imagefmt=GIF width=4in height=4in antialiasmax=60000;  /* size of each image GIF */
options papersize=('4 in', '4 in')                    /* set size for images */
        nodate nonumber nobyline                      /* do not show date, time, or frame number */
        animduration=0.15 animloop=yes noanimoverlay  /* animation details */
        printerpath=gif animation=start;              /* start recording images to GIF */
ods printer file='C:\AnimGif\Snow\Anim.gif';          /* images saved into animated GIF */
ods html select none;                                 /* suppress screen output */

/* write each frame */
proc sgplot data=All noautolegend;
by Frame;
styleattrs wallcolor=CXD6EAF8;                                               * light blue "sky";
polygon x=bgX y=bgY ID=bgID / outline fill fillattrs=(color=whitesmoke);     * white ground;
text x=xc y=yc text=txt / splitchar='|' splitpolicy=splitalways          
          position=center textattrs=(size=32 weight=bold color=ForestGreen); * green text;
polygon x=x y=y ID=ID / outline fill fillattrs=(color=snow) transparency=0.30;
xaxis min=0 max=1 offsetmin=0 offsetmax=0 display=none;
yaxis min=0 max=1 offsetmin=0 offsetmax=0 display=none;
run;

ods html select all;                                  /* restore screen output */ 
options printerpath=gif animation=stop;               /* stop recording images */
ods printer close;                                    /* close the animated GIF file */