カッシーニ曲線

x<-y<-seq(from=-2,to=2,by=0.01)
xy<-expand.grid(x,y)
a<-4
b<-4
z<-(xy[,1]^2+xy[,2]^2)^2-2*a^2*(xy[,1]^2-xy[,2]^2)+a^4

image(x,y,matrix(z,length(x),length(y)))