data triangle ;
  input  id $ x1  y1  x2  y2  x3  y3 ;             
  cards ;
  a -1 1 3 2 1.5 3
 
  ;
  run ;


proc iml worksize=512;
 reset fuzz;
base1=1; 
base2=2; 
use triangle ;
read all into R[colname=coords rowname=id]; 

N=nrow(R);
m=ncol(R);
p=3;   /* actual number of coordinate pairs */
nland=j(N,1)*p;
notu=j(N,1)*n;
 X=j(n,p); /* saves space for necessary arrays */
 Y=j(n,p);
 W=j(n,p);
 Z=j(n,p);
 do i=1 to p;  /* separates out coordinates in x and y arrays */
     x[,i]=r[,(i-1)*2+1];
     y[,i]=r[,(i-1)*2+2];
 end;
 
 call gstart ;
 call gopen("triangle");
 call gwindow({-2 -2, 6 6}); 
 xbox = { -2 2 2 -2 } ;
 ybox = { -2 -2  6 6 } ;
 call gxaxis({-2,0},6,6);
   call gyaxis({0,-2},6,6);

 call gpoly(x,y,,,'X','red');
 call gshow ;
 call gclose ;

xmn=X-x[,:]*j(1,p);  /* three steps to compute centroid size=c_size */
ymn=Y-y[,:]*j(1,p);
c_size=sqrt(xmn[,##]+ymn[,##]);
free xmn ymn;
 /*print r[format=4.0 colname=var rowname=species]; */
 X=X-x[,base1]*j(1,p);    /* translate all specimens to base1 as */
 Y=Y-y[,base1]*j(1,p);           /* the 0,0 point */

call gstart ;
 call gopen("tritran");
 call gwindow({-2 -2, 6 6}); 
 xbox = { -2 2 2 -2 } ;
 ybox = { -2 -2  6 6 } ;
 call gxaxis({-2,0},6,6);
   call gyaxis({0,-2},6,6);

 call gpoly(x,y,,,'X','red');
 call gshow ;
 call gclose ;



 den=sqrt(Y[,base2]##2+X[,base2]##2);
do j=1 to N;
     sintheta=Y[j,base2]/den[j,1];
     costheta=X[j,base2]/den[j,1];
     w[j,]=costheta*X[j,]+sintheta*Y[j,];
     z[j,]=-sintheta*X[j,]+costheta*Y[j,];
end;

call gstart ;
 call gopen("trirot");
 call gwindow({-2 -2, 6 6}); 
 xbox = { -2 2 2 -2 } ;
 ybox = { -2 -2  6 6 } ;
 call gxaxis({-2,0},6,6);
   call gyaxis({0,-2},6,6);

 call gpoly(w,z,,,'X','red');
 call gshow ;
 call gclose ;



size=w[,base2];   /* set aside baseline length */
do J=1 to n;
   dbs=1/size[j,1];
   w[j,]=w[j,]*dbs;
   z[j,]=z[j,]*dbs;
end;
do j=1 to p;     /* put var back in r size adjusted */
 r[,(j-1)*2+1]=w[,j];
 r[,(j-1)*2+2]=z[,j];
end;
Uni=j(n,2);       /* next 5 lines compute uniform components */
wmn=w[:,];
zmn=z[:,];
denom=zmn[,##];
uni[,1]=((W-j(n,1)*wmn)*zmn`)/denom;
uni[,2]=((Z-j(n,1)*zmn)*zmn`)/denom;
coords={ 'nland' 'notu' 'size' 'c_size' 'uni_x' 'uni_y'}||coords;
R=nland||notu||size||c_size||uni||R;  /* create array with all variables */
print R[];
create new3 from R[colname=coords rowname = id];
append from R[rowname =id];

do i=1 to p;  
     w[,i]=w[,i] - 0.5;
   
 end;
print w[]; print z[];
call gstart ;
call gopen("scaled");
 call gwindow({-2 -2, 6 6}); 

 call gxaxis({-2,0},6,6);
   call gyaxis({0,-2},6,6);

 call gpoly(w,z,,,'X','red');
 call gshow ;
 call gclose ;
 

 call gstart ;
 call gopen("newwin") ;
 call gwindow({-2 -2, 2 2}); 
 call gxaxis({-2,0},4,4,1);
   call gyaxis({0,-2},4,4,1);

 call gpoly(w,z,,,'X','red');
 call gshow ;
 call gclose ;

run;
quit ;