• WikiのようなProjective Varietyのような絵(こちら)が書きたかったが…

my.pr <- function(xy,x0,y0,y1){
	k <- 1/(xy[,2]-y1+1)
	new.y <- y0 * k
	new.x <- x0 + (xy[,1]-x0) * k
	return(cbind(new.x,new.y))
}


x <- y <- seq(from=-10,to=10,length=1000)
xyz<- expand.grid(x,y)


a <- xyz[,2]^2-(xyz[,1]^3-xyz[,1]+1)
#a <- xyz[,2]^2-(xyz[,1]^3+1)

a.s <- which(abs(a) < 0.1)

plot(xyz[a.s,1],xyz[a.s,2])
x0 <- 0
y0 <- 10
y1 <- min(xyz)-1

new.xy <- my.pr(xyz[a.s,1:2],x0,y0,y1)
xlim <- range(new.xy[,1])
ylim <- range(new.xy[,2])
plot(new.xy,xlim=xlim,ylim=ylim,col=2)


x1 <- y1 <- min(xyz)-1
x2 <- y2 <- seq(from=-10,to=10,length=50)

xy2<- expand.grid(x2,y2)

new.xy2 <- my.pr(xy2,x0,y0,y1)
par(new=TRUE)
plot(new.xy2,cex=0.1,xlim=xlim,ylim=ylim)


yy <- seq(from=min(xy2),to=max(xy2),length=100)
xx <- rep(1,length(yy))

new.xy3 <- my.pr(cbind(xx,yy),x0,y0,y1)
par(new=TRUE)
plot(new.xy3,xlim=xlim,ylim=ylim,type="l",col=3)




library(rgl)
x <- y <- z <- seq(from=-2,to=2,length=100)
xyz<- expand.grid(x,y,z)


b <- xyz[,2]^2*xyz[,3] -(xyz[,1]^3-xyz[,1]*xyz[,3]^2+xyz[,3]^3)
b.s <- which(abs(b) < 0.01)

plot3d(xyz[b.s,])

plot(xyz[b.s,1]/xyz[b.s,3],xyz[b.s,2]/xyz[b.s,3])