gmapプロシジャでコッホ雪片を描く

Koch_curve.sas

data koch0;
    id = 1; x = cos(constant("pi") * 1 / 6); y = sin(constant("pi") * 1 / 6); density = 0; output;
    id = 1; x = cos(constant("pi") * 5 / 6); y = sin(constant("pi") * 5 / 6); density = 0; output;
    id = 1; x = cos(constant("pi") * 9 / 6); y = sin(constant("pi") * 9 / 6); density = 0; output;
run;


%macro koch(in1=, out1=, dens=);

    proc datasets library = work nolist;
        delete result;
    run;
    quit;
    
    data __m1;
        set &in1. nobs = nobs;
        if _n_ eq 1 then call symputx("nobs", nobs);
        n = _n_;
    run;
    
    %do i = 1 %to &nobs.;
    
        data _null_;
            set __m1(firstobs=&i. obs=&i.);
            call symputx("x1", x);
            call symputx("y1", y);
        run;
        
        %if &i. ne &nobs. %then %do;
            data _null_;
                set __m1(firstobs=%eval(&i.+1) obs=%eval(&i.+1));
                call symputx("x2", x);
                call symputx("y2", y);
            run;
        %end;
        %else %do;
            data _null_;
                set __m1(obs=1);
                call symputx("x2", x);
                call symputx("y2", y);
            run;
        %end;

        data __m2;
            retain id 1 density &dens.;
            
            x = (&x2. + &x1. * 2) / 3;
            y = (&y2. + &y1. * 2) / 3;
            n = &i. + 0.3;
            output;
            
            x = (&x2. * 2 + &x1.) / 3;
            y = (&y2. * 2 + &y1.) / 3;
            n = &i. + 0.7;
            output;
            
            x = cos(-constant("pi") / 3) * ((&x2. * 2 + &x1.) / 3 - (&x2. + &x1. * 2) / 3)
              - sin(-constant("pi") / 3) * ((&y2. * 2 + &y1.) / 3 - (&y2. + &y1. * 2) / 3) + (&x2. + &x1. * 2) / 3;
            y = sin(-constant("pi") / 3) * ((&x2. * 2 + &x1.) / 3 - (&x2. + &x1. * 2) / 3)
              + cos(-constant("pi") / 3) * ((&y2. * 2 + &y1.) / 3 - (&y2. + &y1. * 2) / 3) + (&y2. + &y1. * 2) / 3;
            n = &i. + 0.5;
            output;
        run;
        
        proc append base = result data = __m2;
        run;

    %end;

    proc append base = result data = __m1;
    run;
    
    proc sort data = result out = &out1.(drop=n);
        by n;
    run;

%mend koch;

%koch(in1=koch0, out1=koch1, dens=1);
%koch(in1=koch1, out1=koch2, dens=2);
%koch(in1=koch2, out1=koch3, dens=3);
%koch(in1=koch3, out1=koch4, dens=4);
%koch(in1=koch4, out1=koch5, dens=5);
%koch(in1=koch5, out1=koch6, dens=6);


%macro draw_koch_map(dens=);

    proc gmap map = koch6 data = koch0 density = &dens.;
        id id;
        choro id / statistic = first nolegend;
    run;
    quit;
    
%mend draw_koch_map;

%draw_koch_map(dens=0);
%draw_koch_map(dens=1);
%draw_koch_map(dens=2);
%draw_koch_map(dens=3);
%draw_koch_map(dens=4);
%draw_koch_map(dens=5);
%draw_koch_map(dens=6);