plot.surface <- function(bgrid,u,b0,max=TRUE,textb0=TRUE,xlab='b1',ylab='b2',...) {
# Description: draws countour plot of the expected utility at a fixed value of b0
# Input:
#   - bgrid: grid with b0,b1,b2 values as the element 'bgrid' returned by 'read.forwsim'
#   - u: expected utility for each point in the grid
#   - b0: only grid points with the indicated b0 will be used for plotting
#   - max: if set to TRUE the maximum will be indicated with a square. 
#   - textb0: if TRUE, the value of b0 will be added to the top right of the plot
#   - ...: optional arguments to 'contour'
# Return: nothing
# Side effects: creates a contour plot with u as a function of the first two columns in bgrid

if (missing(b0)) stop('b0 must be specified')
if (length(b0)>1) stop('b0 must be a vector of length 1')
if (nrow(bgrid)!=length(u)) stop('Dimensions of bgrid and u must match')

b1 <- as.numeric(names(table(bgrid[,2])))
b2 <- as.numeric(names(table(bgrid[,3])))

b <- bgrid[bgrid[,1]==b0,]
u <- u[bgrid[,1]==b0]
contour(b1,b2,matrix(u,nrow=length(b1),ncol=length(b2),byrow=T),xlab=xlab,ylab=ylab,...)
if (textb0) mtext(paste('b0=',round(b0,3)),3,adj=1)
if (max) points(b[u==max(u),2],b[u==max(u),3],pch=15,cex=2)

}
