Rで自作関数:get_sun_position

#' 太陽の高度と方位を計算する関数
#' 
#' @param lat 緯度（度、10進数）
#' @param lon 経度（度、10進数）
#' @param datetime 日時（POSIXct型）
#' @return 高度（elevation）と方位（azimuth）のリスト（単位は度）
get_sun_position <- function(lat, lon, datetime) {
  
  if (!inherits(datetime, "POSIXct")) {
    stop("エラー: datetimeはPOSIXct型である必要があります。")
  }
  
  # --- 内部ヘルパー関数 ---
  deg2rad <- function(deg) deg * pi / 180
  rad2deg <- function(rad) rad * 180 / pi
  
  # --- 1. 時間の計算 ---
  unix_time <- as.numeric(datetime)
  jd <- (unix_time / 86400.0) + 2440587.5
  jc <- (jd - 2451545.0) / 36525.0
  
  # --- 2. 地球と太陽の軌道要素 ---
  l0 <- (280.46646 + jc * (36000.76983 + jc * 0.0003032)) %% 360
  m <- 357.52911 + jc * (35999.05029 - 0.0001537 * jc)
  e <- 0.016708634 - jc * (0.000042037 + 0.0000001267 * jc)
  
  c_rad <- sin(deg2rad(m)) * (1.914602 - jc * (0.004817 + 0.000014 * jc)) + 
           sin(deg2rad(2 * m)) * (0.019993 - 0.000101 * jc) + 
           sin(deg2rad(3 * m)) * 0.000289
           
  lambda <- l0 + c_rad
  omega <- 125.04 - 1934.136 * jc
  lambda_app <- lambda - 0.00569 - 0.00478 * sin(deg2rad(omega))
  
  eps0 <- 23.0 + (26.0 + ((21.448 - jc * (46.815 + jc * (0.00059 - jc * 0.001813))) / 60.0)) / 60.0
  eps <- eps0 + 0.00256 * cos(deg2rad(omega))
  
  # --- 3. 太陽赤緯と均時差 ---
  declination <- rad2deg(asin(sin(deg2rad(eps)) * sin(deg2rad(lambda_app))))
  
  y <- tan(deg2rad(eps / 2.0))^2
  eot <- 4.0 * rad2deg(
    y * sin(deg2rad(2 * l0)) - 
    2.0 * e * sin(deg2rad(m)) + 
    4.0 * e * y * sin(deg2rad(m)) * cos(deg2rad(2 * l0)) - 
    0.5 * y^2 * sin(deg2rad(4 * l0)) - 
    1.25 * e^2 * sin(deg2rad(2 * m))
  )
  
  # --- 4. 真太陽時と時角 ---
  time_utc_min <- (unix_time %% 86400) / 60.0
  tst <- (time_utc_min + eot + 4.0 * lon) %% 1440
  ha <- (tst / 4.0) - 180.0
  
  # --- 5. 高度と方位の計算（球面三角法） ---
  lat_rad <- deg2rad(lat)
  dec_rad <- deg2rad(declination)
  ha_rad  <- deg2rad(ha)
  
  # 幾何学的な高度
  el_rad <- asin(sin(lat_rad) * sin(dec_rad) + cos(lat_rad) * cos(dec_rad) * cos(ha_rad))
  elevation_true <- rad2deg(el_rad)
  
  # 大気差の補正 (Atmospheric Refraction) を追加
  # 光が地球の大気を通過する際の屈折を補正します（地平線に近いほど影響大）
  refraction <- ifelse(
    elevation_true > -1,
    1.02 / tan(deg2rad(elevation_true + 10.3 / (elevation_true + 5.11))) / 60.0,
    0
  )
  elevation <- elevation_true + refraction
  
  # 太陽方位 (Azimuth) の計算
  # atan2() を使用することで、北(0度)から時計回りの角度を全ての象限で取得
  sin_az <- -sin(ha_rad) * cos(dec_rad) / cos(el_rad)
  cos_az <- (sin(dec_rad) - sin(lat_rad) * sin(el_rad)) / (cos(lat_rad) * cos(el_rad))
  
  az_rad <- atan2(sin_az, cos_az)
  azimuth <- (rad2deg(az_rad) + 360) %% 360
  
  return(list(
    elevation = elevation,
    azimuth   = azimuth
  ))
}

# 緯度・経度を設定（例: 東京）
my_lat <- 35.689
my_lon <- 139.691

# 日時を設定（必ずタイムゾーン "Asia/Tokyo" などを指定する）
# 例として、2026年4月3日の正午（12:00）を計算
my_datetime <- as.POSIXct("2026-04-03 12:00:00", tz = "Asia/Tokyo")

# 関数の実行
result <- get_sun_position(lat = my_lat, lon = my_lon, datetime = my_datetime)

# 結果の表示
cat(sprintf("日時: %s\n", my_datetime))
cat(sprintf("太陽の高度: %.2f 度\n", result$elevation))
cat(sprintf("太陽の方位: %.2f 度 (北=0, 東=90, 南=180, 西=270)\n", result$azimuth))