* 출처 : http://blogs.sas.com/content/iml/?s=koch-snowflake&searchsubmit=Find
* 결과 파일 : http://blogs.sas.com/content/iml/files/2016/12/kochsnow.gif
/*********************************************/
/* 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 */
Join us for SAS Innovate 2025, our biggest and most exciting global event of the year, in Orlando, FL, from May 6-9. Sign up by March 14 for just $795.