#transformations worksheet
#Translated from transformations.mws worksheet by DS Wethey 22 Aug 2011
#

Ltransform<-function(A,seg,ADD=F)
{
if(ADD){par(new=T)}

plot(seg,type="l",xlim=c(-10,10),ylim=c(-10,10),lwd=3,asp=1)
	new<-matrix(ncol=2,nrow=nrow(seg))
	for (ii in 1:nrow(seg))
		{new[ii,]<-A%*%seg[ii,]}
lines(new,col="red",lwd=1)
abline(v=0);abline(h=0)
par(new=F)
}

Atransform<-function(A,seg,ADD=F)
{
if(ADD){par(new=T)}

plot(seg,type="l",xlim=c(-10,10),ylim=c(-10,10),lwd=3,asp=1)
	new<-matrix(ncol=2,nrow=nrow(seg))
	for(ii in 1:nrow(seg))
		{new[ii,]<-(A[,1:2]%*%seg[ii,])+A[,3]}
	lines(new,col="red",lwd=1)
	abline(v=0);abline(h=0)

par(new=F)
}

A<-matrix(nrow=2,ncol=2,c(1,3,3,1),byrow=T);A #symmetric matrix
B<-matrix(nrow=2,ncol=2,c(1,2,0,1),byrow=T);B #shear
C<-matrix(nrow=2,ncol=2,c(cos(pi/3),-sin(pi/3),sin(pi/3),cos(pi/3)),byrow=T);C #rotation

SQ<-matrix(ncol=2,c(0, 0, 1, 0, 1, 1, 0, 1 , 0, 0),byrow=T) #a box
plot(SQ,xlim=c(-2,2),ylim=c(-2,2),asp=1)
polygon(SQ,col="green")

plot(quasicircle,xlim=c(-2,2),ylim=c(-2,2),asp=1)
quasicircle<-matrix(ncol=2,c(cos((0:72)*pi/36), sin((0:72)*pi/36)))
polygon(quasicircle,col="yellow")

seg1<-  matrix(ncol=2,byrow=T,c(0, 0, 3, -2)); plot( seg1,type="l");  #  a segment
seg2<-  matrix(ncol=2,byrow=T,c(0, 0 , 1, 1));  plot( seg2,type="l" ); # another segment
seg3<-  matrix(ncol=2,byrow=T,c(0, 0, -1, 1)); plot( seg3,type="l" ); #  yet another segment
seg4<-  matrix(ncol=2,byrow=T,c(0, 0, 0.5, 0.2)); plot( seg4,type="l" ); # and just one more

plot.new()
Ltransform(A,seg1)

plot.new()
Ltransform(B,seg1)

plot.new()
Ltransform(C,seg1)

plot.new()
Ltransform(A,seg2)
# where did seg2 go?  

plot.new()
Ltransform(A,seg3)

plot.new()
Ltransform(A,SQ)
Ltransform(A,seg2,T)

plot.new()
Ltransform(A,quasicircle)
Ltransform(A,seg2,T)
Ltransform(A,matrix(ncol=2,byrow=2,c(0,0,cos(3*pi/4),sin(3*pi/4))),T)

#Can you see which parts of the square went where?  Once again, what happened to the original segment 2?
plot.new()
Ltransform(B,SQ)
Ltransform(B,seg2,T)

#Can you see why this is called a shear?
plot.new()
Ltransform(C,SQ)
Ltransform(C,seg2,T)

#Now it is your turn! Define new matrices, new vectors (segmens), and see what the transformations 
#do to both the vectors and the unit square


#Now lets look at affine transformations.
#In this case instead of taking [x,y] to [ax+by, cx+dy], which takes [0,0] to itself,
#  we take [x,y] to [ax+by+r,cx+dy+s] which takes [0,0] to [r,s]
#Since there are 6 parameters involved, you need to specify a 2x3 matrix

F<-matrix(nrow=2,ncol=3,byrow=T,c(1,-0.5,1,-0.5,1,0))
plot.new()
Atransform(F,seg1)

plot.new()
Atransform(F,seg3)

plot.new()
Atransform(F,seg4)

Atransform(F,SQ)
Atransform(F,quasicircle)

#Eigenvalues and eigenvectors
#A matrix is entered by giving the number of rows, number of columns, the order of entries,
#       and a list of the entries from left to right, top row to bottom row. 
E<-matrix(nrow=2,ncol=2,byrow=T,c(-4,5/2,-5/2,1))
eigen(E)

#The eigen command calculates eigenvalues and eigenvectors for the matrix

V<-eigen(E)
V

#The eigenvectors can be seen by
V$vectors

#The eigenvalues can be seen by
V$values

#extract individual eigenvectors using the column select command:
v1<-V$vectors[,1] ;v1
v2<-V$vectors[,2] ;v2

#Check for yourself that this is correct!
#Is it true that A v1 = (-3/2) v1?
#Let's try using our transform procedure to do a graphical check
seg<-matrix(nrow=2,ncol=2,byrow=T,c(0,0,1,1))
plot.new()
Ltransform(E,seg)

#Try this for multiples of (1,1) -- what do you notice?

#Here's another example, this time with two distinct complex eigenvalues.  
#Oftentimes we don't care about specific complex eigenvalues; 
#   the mere absence of real eigenvalues is all we need to know.  
#   Other times we only care about features of the complex eigenvalues, 
#   such as whether they are pure imaginary, or if they have positive or 
#   negative real parts.  

F<- matrix(nrow=2, ncol=2, byrow=T,c(4, 1, -2, 3));   eigen(F);

#The next one has two distinct pure imaginary eigenvalues.
 G<- matrix(nrow=2, ncol=2, byrow=T,c(-0.4, -4, 2, 0.4));   eigen(G);

#The next has two distict real eigenvalues.
 H<-  matrix(nrow=2, ncol=2, byrow=T,c(2, 3, 4, 1));   eigen(H);
 
 EV<-eigen(H);
 
 v1<-EV$vectors[,1];v1
 v2<-EV$vectors[,2];v2
 
#we can rescale v1 and v2 by dividing them by their first elements:

v1a<-v1/v1[1]; v1a
v2a<-v2/v2[1]; v2a
  
seg1<- matrix(nrow=2,ncol=2,byrow=T,c(0, 0, v1a));  seg2<-matrix(nrow=2,ncol=2,byrow=T,c(0,0,v2a));
 
plot.new()
Ltransform(H,seg1)
plot.new()
Ltransform(H,seg2)