おっぱい関数

「おっぱい関数」を材料として3次元表示についてまとめる。

「おっぱい関数」は元々誰が作ったか不詳だが，検索すると
Julia + Excel　http://www.hirax.net/diaryweb/2012/03/01.html
Maxima　https://plus.google.com/u/0/117731699005945352808/posts/Yswj5tswtNd
などが引っかかる。変形したものもある。
Python　http://kemeconoajito.blog88.fc2.com/blog-entry-175.html
R によるものは，http://siguniang.wordpress.com/2011/11/09/3d-surface-plot-with-latticergl/
にある。

library(lattice)
 
g <- function(x, y) {
 1/8 * (6 * exp(-((2/3 * abs(x) - 1)^2 + (2/3 * y)^2) - 1/3 * (2/3 * y + 1/2)^3) +
        2/3 * exp(-2.818^11 * ((abs(2/3 * x) - 1)^2 + (2/3 * y)^2)^2) +
        2/3 * y - (2/3 * x)^4)
}
m <- 50
x <- seq(-3, 3, length.out=m)
y <- seq(-3, 3, length.out=m)
grid   <- expand.grid(x=x, y=y)
grid$z <- g(grid$x, grid$y)
wireframe(z ~ x * y, grid, shade=TRUE)
library(rgl)
grid$col <- ifelse(sqrt((abs(grid$x) - 3/2) ^ 2 + grid$y ^2 ) < 0.45, 'hotpink', 'lightpink')# nipple
open3d()
rgl.surface(x, y, grid$z, color=grid$col)

この関数部分を小生が書き直した

n <- 100
x0 <- y0 <- seq(-3, 3, length=n+1)
xy <- expand.grid(x=x0, y=y0)
x <-  xy[,1]
y <-  xy[,2]
t0 <- (2 * abs(x) / 3 - 1) ^ 2
t1 <- 2 * y / 3
t2 <- t1 ^ 2
t4 <- -t0 - t2 - (t1 + 0.5) ^ 3 / 3
z <- 0.125 * (6 * exp(t4) +
     2 * exp(-exp(11) * (t0 + t2) ^ 2) / 3 +
     2.5 * t1 -16 * x ^ 4 / 81)
library(rgl)
plot3d(x, y, z, col = "burlywood1", aspect = c(1, 1, 0.6), axes=FALSE, xlab="", ylab="", zlab="")

# 以下の，ふたつの方法もやってみる
open3d()
rgl.surface(x0, y0, z, color="burlywood1")
library(lattice)
wireframe(z ~ x * y, xy, shade=TRUE, aspect=c(1, 0.5), xlab="", ylab="", zlab="")

# 以下で出力される csv ファイルを Excel で 3d 表示
write(matrix(z, n+1), file="oppai2.csv", ncolumns=n+1, sep=",")

変形したものも書いておく。
http://kemeconoajito.blog88.fc2.com/blog-entry-175.html

n <- 100
xy <- expand.grid(seq(-3, 3, length=n+1), seq(-3, 3, length=n+1))
x <-  xy[,1]
y <-  xy[,2]
t0 <- (2 * abs(x) / 3 - 1) ^ 2
t1 <- 2 * y / 3
t2 <- t1 ^ 2
t3 <- (t0 + t2) ^ 2
t4 <- -t0 - t2 -  (t1 + 0.5) ^ 3 / 3
z <- 0.125 * (7 * exp(t4) + 1.5 * exp(-exp(10) * t3) +
     exp(-exp(5) * t3) + 2.5 * t1 - 16 * x ^ 4 / 81)
library(rgl)
plot3d(x, y, z, col = "burlywood1", aspect = c(1, 1, 0.6), axes=FALSE, xlab="", ylab="", zlab="")

ついでに見つかった，2次元の「おっぱい関数」も R で書いたものを載せておく。

# https://twitter.com/math_neko/status/227672102487068672/photo/1
y <- seq(.Machine$double.eps, 1.2, length=300)
x <- 3*y*log(y)-1/36*exp(-(36*y-36/exp(1))^4)
plot(1, 1, type="n", xlim=c(-1.2, 1), ylim=c(-0.1, 1))
lines(x, y)

# http://goo.gl/rtZS7
n <- 5000
curve(sqrt(1-(x-1)^2), n=n, xlim=c(-2.5, 2.5), ylim=c(-0.1,  1.25), asp=1)
curve(sqrt(1-(x+1)^2), n=n, add=TRUE)
curve(sqrt(0.01-(x-1)^2)+1, n=n, add=TRUE)
curve(sqrt(0.01-(x+1)^2)+1, n=n, add=TRUE)

# http://i.imgur.com/IHAhR.jpg
t <- seq(-1.4, 1.5, length=200)
y <- 2*t^3-0.5*t-1
x <- 0.9*(exp(-t^4+t^2+t)+0.2*(exp(-100*((y-0.2)^2)))-0.5*exp(-(y-3.5)^2)-0.15*exp(-(y-5)^2))
plot(x, y, type="l", asp=1)