# 1. データをリセットして点の数を設定
rm(list = ls())
N <- 200000

# 2. 箱の中にランダムな点を打つ
# Xは 0 から π、 Yは 0 から 1 
x_points <- runif(N, min = 0, max = pi)
y_points <- runif(N, min = 0, max = 1)

# 3. 点が「うねうねの山の下側」にあるか判定
# Rでは「eのX乗」は exp(X) と書きます
inside <- y_points <= (sin(x_points) * exp(-x_points))

# 4. 箱の中に入った割合を計算
ratio <- sum(inside) / N

# 5. 【どう見ても面積πの箱】を掛け算
box_area <- pi * 1
estimated_area <- ratio * box_area

# 6. 正しい理論値を計算
true_area <- (1 + exp(-pi)) / 2

# 結果の表示
cat("======================================\n")
cat(" モンテカルロ法:", estimated_area, "\n") # 0.432付近になります
cat(" 正しい理論値  :", true_area, "\n")       # 0.4321691...
cat("======================================\n")

plot_n <- 5000
colors <- ifelse(inside[1:plot_n], "blue", "red")

# グラフの描画
plot(x_points[1:plot_n], y_points[1:plot_n], col = colors, pch = 20,
     main = "減衰振動の面積計算", xlab = "X", ylab = "Y")

# 実際の数式通りの線を黒で引く
curve(sin(x) * exp(-x), add = TRUE, col = "black", lwd = 3, from = 0, to = pi)