# Westermann Morphospace functions 


# for Ritterbush and Bottjer, 2012. Westermann Morphospace displays ammonoid shell shape and hypothetical paleoecology. Paleobiology, 38(3), 424-446.

# Last update: 9/21/13

# for questions or updates, email Kathleen at ritterbu@uchicago.edu

Westyspace<-function(data,identifier,obj.name,title,color="black",symbol=1,size=1){
	
# We will build a complex object that contains several matricies.  The name of the complex object is given as an argument in the function.

# First, we generate a table of values for the shape parameters.

parameters<- matrix(nrow=nrow(data), ncol=3)

rownames(parameters)<-data[,identifier]
colnames(parameters)<-c("U","w","Th")

parameters[,1]<-(data[,"Umbilical.Diameter"]/data[,"Diameter"])

parameters[,2]<-(data[,"a"]/(data[,"Diameter"]-data[,"a"]-data[,"Umbilical.Diameter"]))

parameters[,3]<-(data[,"b"]/data[,"Diameter"])



# Next, we scale each of these parameter values.

scaled<-matrix(nrow=nrow(parameters), ncol=3)

rownames(scaled)<-rownames(parameters)
colnames(scaled)<-c("U.scaled","w.scaled","Th.scaled")

scaled[,1]<-(parameters[,1]/.52)
scaled[,2]<-((parameters[,2]-1)/.77)
scaled[,3]<-((parameters[,3]-.14)/.54)



# Finally, we normalize these scaled values into percentages of the entire shape characterization.

percents<-matrix(nrow=nrow(scaled),ncol=3)

rownames(percents)<-rownames(scaled)
colnames(percents)<-c("Serpenticone","Oxycone","Sphaerocone")

percents[,1]<-(scaled[,1]/(scaled[,1]+scaled[,2]+scaled[,3]))
percents[,2]<-(scaled[,2]/(scaled[,1]+scaled[,2]+scaled[,3]))
percents[,3]<-(scaled[,3]/(scaled[,1]+scaled[,2]+scaled[,3]))


# Create x,y coordinates to plot the ternary diagram
coordinates<-matrix(nrow=nrow(percents), ncol=2)

for(i in 1:nrow(percents)){
	coordinates[i,1]<-percents[i,3]+(percents[i,2]/2)
	}
	
for(i in 1:nrow(percents)){
	coordinates[i,2]<-(1-(percents[i,2]*sqrt(3)/2))
	
	}
	
row.names(coordinates)<-row.names(percents)


# Indicate corners of the triangle to draft the frame
#create frame:
corner.guys<-matrix(nrow=3, ncol=3)

corner.guys[1,]<-c(1,0,0)
corner.guys[2,]<-c(0,1,0)
corner.guys[3,]<-c(0,0,1)

row.names(corner.guys)<-c("serpenticone", "oxycone", "sphearocone")

corner.xy<-matrix(nrow=3,ncol=2)

for(i in 1:nrow(corner.guys)){
	corner.xy[i,1]<-corner.guys[i,3]+(corner.guys[i,2]/2)
	}
	
for(i in 1:nrow(corner.guys)){
	corner.xy[i,2]<-1-(corner.guys[i,2]*sqrt(3)/2)
	}


# create the object
assign(obj.name,list(Parameters=parameters,Scaled=scaled,Percents=percents,Coordinates=coordinates), envir=.GlobalEnv)

plot( coordinates,frame.plot=FALSE, ann=FALSE,axes=FALSE, xlim=c(-.1,1.1), ylim=c(-.1,1.1), col=color,pch=symbol,cex=size)

lines(corner.xy)
lines(corner.xy[1,],corner.xy[3,])
text(x=c(0,.5,1),y=c(1.05,0.1,1.05),labels=c(row.names(corner.guys)))
title(main=title)
	
	
identify(coordinates,labels=row.names(coordinates), cex=.5)

	}
	
	
	##### Export format: This will create only the data points and the triangle around them.  They can be copied/edited into publication figures with Adobe Illustrator or other vector graphics program.

	
#plot(coordinates,frame.plot=FALSE, ann=FALSE,axes=FALSE, xlim=c(0,1), ylim=c(0,1))
#lines(corner.xy)
#lines(corner.xy[1,],corner.xy[3,])





Westypoints<-function(data,identifier,color="black",symbol=1,size=1){
	
# First, we generate a table of values for the shape parameters.

parameters<- matrix(nrow=nrow(data), ncol=3)

rownames(parameters)<-data[,identifier]
colnames(parameters)<-c("U","w","Th")

parameters[,1]<-(data[,"Umbilical.Diameter"]/data[,"Diameter"])

parameters[,2]<-(data[,"a"]/(data[,"Diameter"]-data[,"a"]-data[,"Umbilical.Diameter"]))

parameters[,3]<-(data[,"b"]/data[,"Diameter"])



# Next, we scale each of these parameter values.

scaled<-matrix(nrow=nrow(parameters), ncol=3)

rownames(scaled)<-rownames(parameters)
colnames(scaled)<-c("U.scaled","w.scaled","Th.scaled")

scaled[,1]<-(parameters[,1]/.52)
scaled[,2]<-(parameters[,2]-1)
scaled[,3]<-((parameters[,3]-.14)/.54)



# Finally, we normalize these scaled values into percentages of the entire shape characterization.

percents<-matrix(nrow=nrow(scaled),ncol=3)

rownames(percents)<-rownames(scaled)
colnames(percents)<-c("Serpenticone","Oxycone","Sphaerocone")

percents[,1]<-(scaled[,1]/(scaled[,1]+scaled[,2]+scaled[,3]))
percents[,2]<-(scaled[,2]/(scaled[,1]+scaled[,2]+scaled[,3]))
percents[,3]<-(scaled[,3]/(scaled[,1]+scaled[,2]+scaled[,3]))


# Create x,y coordinates to plot the ternary diagram
coordinates<-matrix(nrow=nrow(percents), ncol=2)

for(i in 1:nrow(percents)){
	coordinates[i,1]<-percents[i,3]+(percents[i,2]/2)
	}
	
for(i in 1:nrow(percents)){
	coordinates[i,2]<-(1-(percents[i,2]*sqrt(3)/2))
	
	}
	
row.names(coordinates)<-row.names(percents)

	
#save(parameters, file="parameters.RData", ascii=T)
#save(scaled, file="scaled.RData", ascii=T)
#save(percents, file="percents.RData", ascii=T)
#save(coordinates, file="coordinates.RData", ascii=T)	
	
points(coordinates, ann=FALSE, xlim=c(-.1,1.1), ylim=c(-.1,1.1), col=color,pch=symbol,cex=size)
#lines(corner.xy)
#lines(corner.xy[1,],corner.xy[3,])
#text(x=c(0,.5,1),y=c(1.05,0.1,1.05),labels=row.names(corner.xy))

	
	
identify(coordinates,labels=row.names(coordinates), cex=.5)

	}
	
	# These functions are used if the coordinates have already been generated:


Westyplot<-function(coordinates,title, color="black",symbol=1,size=1){
#create frame:
corner.guys<-matrix(nrow=3, ncol=3)

corner.guys[1,]<-c(1,0,0)
corner.guys[2,]<-c(0,1,0)
corner.guys[3,]<-c(0,0,1)

row.names(corner.guys)<-c("serpenticone", "oxycone", "sphearocone")

corner.xy<-matrix(nrow=3,ncol=2)

for(i in 1:nrow(corner.guys)){
	corner.xy[i,1]<-corner.guys[i,3]+(corner.guys[i,2]/2)
	}
	
for(i in 1:nrow(corner.guys)){
	corner.xy[i,2]<-1-(corner.guys[i,2]*sqrt(3)/2)
	}
	# plot data:
plot( coordinates,frame.plot=FALSE, ann=FALSE,axes=FALSE, xlim=c(-.1,1.1), ylim=c(-.1,1.1), col=color,pch=symbol,cex=size)

lines(corner.xy)
lines(corner.xy[1,],corner.xy[3,])
text(x=c(0,.5,1),y=c(1.05,0.1,1.05),labels=row.names(corner.xy))
title(main=title)
	
	
identify(coordinates,labels=row.names(coordinates), cex=.5)

	
}
	
WestypointsAdd<-function(coordinates, color="black",symbol=1,size=1){
	points(coordinates, ann=FALSE,xlim=c(-.1,1.1), ylim=c(-.1,1.1), col=color,pch=symbol,cex=size)
#lines(corner.xy)
#lines(corner.xy[1,],corner.xy[3,])
#text(x=c(0,.5,1),y=c(1.05,0.1,1.05),labels=row.names(corner.xy))
identify(coordinates,labels=row.names(coordinates), cex=.5)
	}