搞数学、IT的人往往都不能歌善舞，也许算法、逻辑与程序是他们最好的表达方式，因此在这分享个画图的程序。一切都源于笛卡尔的那个极坐标方程，他和公主的传说大家都知道吧！？（可惜是假的 哈哈）
还有很多可以美化的地方，请大家帮助一起完善。

/*注：程序为创建新型曲线的动画，在SAS中运行，输出结果为gif文件。
去C盘下点开gif文件用IE浏览器就可以看了，FireFox好像不行。
一切都源于极坐标方程“r=a*(1-sin(t))”，在此基础上对x、y缩放使曲线更像一个“心”。
在循环中，对“心”按半径进行展开，并随机平移、旋转和改变染色，用来获得动画效果。
具体参数的调整请参考程序中的注释。
*/

%let path=C:\HeartAnimation.gif;/*输出动画文件路径*/

%let N=20; /*心型展开的个数*/
%let RotateFreq=3; /*心形随机旋转的频率*/
%let AddFreq=5; /*动画重现频率*/

%let xPropotion=0.6;/*心形x轴缩放比例*/
%let yPropotion=1.5;/*心形y轴缩放比例*/
%let r0max=5;/*心形半径最大值*/
%let r0min=2;/*心形半径最小值*/
%let x0max=5;/*心形平移最大值*/
%let y0max=5;/*心形平移最小值*/
%let colorList=RED YELLOW BLUE PINK ORANGE;/*心形颜色列表。每种颜色之间必须以1个空格进行分割，否则即使颜色数会出错。*/

%let cback=SkyBlue;/*背景颜色*/
%let iteration=1;/*动画循环次数，0为不间断循环*/
%let delay=20;/*动画图片间隔时间*/

%macro Heart;
%do r=1 %to &r0.;
data Heart;
        pi=constant('pi');
        do t=-pi to pi by 0.1;
                r0=0.5*&r.*(1-sin(&t0.+t));
                x=&x0.+&xPropotion.*r0*cos(t);
                y=&y0.+&yPropotion.*r0*sin(t);
                s=&s.;
                add=0;
                r=&r.;
                color=&color.;
                output;
        end;
        keep x y s add r color;
    run;

    proc append base=PlotData data=Heart force;
    run;
%end;
%mend;

proc sql;/*初始化PlotData，准备append*/
   create table PlotData (x num, y num, s num, r num, add num, color num);
quit;

%macro PlotData;
%do s=1 %to &N.;
    data _null_;/*生成%Heart所需要的宏变量*/
                pi=constant('pi');
                colorNum=1+length("&colorList.")-length(compress("&colorList."));
        call symput('x0',2*&x0max.*(ranuni(0)-0.5));
        call symput('y0',2*&y0max.*(ranuni(0)-0.5));
        call symput('r0',&r0min.+floor((&r0max.-&r0min)*ranuni(0)));
        call symput('color',1+floor(colorNum*ranuni(0)));       
                if mod(&s.,&RotateFreq.)=0 then call symput('t0',2*pi*ranuni(0));               
                else call symput('t0',0);       
                if mod(&s.,&AddFreq.)=0 then call symput('add',1);               
                else call symput('add',0);
    run;

    %Heart;

        %if &add.=1 %then %do;
                data PlotData_Add;/*按周期重新之前的动画*/
                        set PlotData;
                        s=&s.;
                        add=1;
                run;

                proc append base=PlotData data=PlotData_Add force;
                run;
        %end;
%end;

data PlotData_Last;/*生成基本心形曲线图*/
        pi=constant('pi');
        do t=-pi to pi by 0.001;
                r0=1*(1-sin(&t0.+t));
                x=&r0max.*r0*cos(t);
                y=&y0max.+&r0max.*r0*sin(t);
                s=2*&N.+1;
                add=0;
                r=.;
                color=1;
                output;
        end;
        keep x y s r add color;
run;

proc append base=PlotData data=PlotData_Last force;
run;
%mend;
%PlotData;

filename out "&path.";
goptions
        reset=all
        device=gifanim
        gsfname=out
    gsfmode=replace
        ftext=swissb
        ctext=navy
        htext=1
    cback=&cback.
        iteration=&iteration.
        delay=&delay.
;

symbol1 value=# h=2;
goptions colors=(&colorList.);

%let xRange=%eval(1+&x0max.+&r0max.);
%let yRange=%eval(1+&y0max.+&r0max.);
axis1 order=(-&yRange. to &yRange. by 1) offset=(1,1) minor=(number=1) label=none ;
axis2 order=(-&xRange. to &xRange. by 1) offset=(1,1) minor=(number=1) label=none ;

proc sort data=PlotData;
        by s add r;
proc gplot data=PlotData uniform;
        plot y*x=color  /vaxis=axis1 haxis=axis2 CFRAME=SkyBlue;
        by s add r;
run;
quit;

data _null_;
        file out recfm=n mod;
        put '3B'x;
run;
filename out clear;