spherical_triangle_area <- function(R, a, b, c) {
# 直交座標への変換
to_cartesian <- function(lat, lon) {
c(cos(lat) * cos(lon), cos(lat) * sin(lon), sin(lat))
}

p <- to_cartesian(a[1], a[2])
q <- to_cartesian(b[1], b[2])
r <- to_cartesian(c[1], c[2])

# 辺のベクトルを計算
u_ab <- q - p
u_bc <- r - q
u_ca <- p - r

# 内角を計算
angle_between <- function(u, v) {
acos(sum(u * v) / (sqrt(sum(u^2)) * sqrt(sum(v^2))))
}

alpha <- angle_between(u_ab, -u_ca)
beta <- angle_between(u_bc, -u_ab)
gamma <- angle_between(u_ca, -u_bc)

# 超過角を計算
E <- alpha + beta + gamma - pi

# 球面三角形の面積を計算
R^2 * E
}