/*********************************************/ /* 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 */