<?php

// 潮汐力の計算
function calculateTidalForce($G, $M, $d, $a, $theta) {
    // 水平成分
    $fTh = 3/2 * ($G * $M / pow($d, 3)) * $a * pow(sin($theta), 2);
    
    // 鉛直成分
    $fTv = 3 * ($G * $M / pow($d, 3)) * $a * (pow(cos($theta), 2) - 1/3);
    
    return ['horizontal' => $fTh, 'vertical' => $fTv];
}

// 潮汐ポテンシャルの計算
function calculateTidalPotential($G, $M, $d, $a, $theta) {
    return 3/2 * ($G * $M / pow($d, 3)) * pow($a, 2) * (pow(cos($theta), 2) - 1/3);
}

// 平衡潮汐の計算
function calculateEquilibriumTide($g, $M, $E, $a, $c, $d, $theta) {
    $D = 3/4 * $g * ($M / $E) * pow($a / $c, 3) * $a;
    return 2 * $D * pow($c / $d, 3) * (pow(cos($theta), 2) - 1/3);
}

// 観測点での潮汐ポテンシャルの計算
function calculateObservationPointPotential($D, $c, $d, $phi, $delta, $T) {
    $cosTheta = cos($phi) * cos($delta) * cos($T) + sin($phi) * sin($delta);
    return $D * pow($c / $d, 3) * (
        pow(cos($phi), 2) * pow(cos($delta), 2) * pow(cos($T), 2) +
        pow(sin($phi), 2) * pow(sin($delta), 2) * cos($T) +
        3 * (pow(sin($phi), 2) - 1/3) * (pow(sin($delta), 2) - 1/3)
    );
}

// 使用例
$G = 6.674e-11; // 万有引力定数
$M = 7.34767309e22; // 月の質量
$E = 5.97219e24; // 地球の質量
$d = 384400000; // 地球と月の平均距離
$a = 6371000; // 地球の半径
$g = 9.81; // 重力加速度
$c = $d; // 地球と月の平均距離（簡略化のため同じ値を使用）
$theta = deg2rad(45); // 例として45度を使用
$phi = deg2rad(35); // 観測点の緯度（例：35度）
$delta = deg2rad(5); // 天体の赤緯（例：5度）
$T = deg2rad(30); // 時角（例：30度）

$tidalForce = calculateTidalForce($G, $M, $d, $a, $theta);
$tidalPotential = calculateTidalPotential($G, $M, $d, $a, $theta);
$equilibriumTide = calculateEquilibriumTide($g, $M, $E, $a, $c, $d, $theta);
$observationPointPotential = calculateObservationPointPotential($D, $c, $d, $phi, $delta, $T);

echo "Tidal Force (Horizontal): " . $tidalForce['horizontal'] . " N\n";
echo "Tidal Force (Vertical): " . $tidalForce['vertical'] . " N\n";
echo "Tidal Potential: " . $tidalPotential . " J/kg\n";
echo "Equilibrium Tide: " . $equilibriumTide . " m\n";
echo "Observation Point Potential: " . $observationPointPotential . " J/kg\n";

<?php

/**
 * 時角を計算する
 * 
 * @param int $year 年
 * @param int $month 月
 * @param int $day 日
 * @param int $hour 時
 * @param int $minute 分
 * @param int $second 秒
 * @param float $longitude 経度（東経は正、西経は負）
 * @return float 時角（度数法）
 */
function calculateHourAngle($year, $month, $day, $hour, $minute, $second, $longitude) {
    // 修正ユリウス日を計算
    $mjd = calculateModifiedJulianDate($year, $month, $day, $hour, $minute, $second);
    
    // グリニッジ恒星時を計算
    $gst = calculateGreenwichSiderealTime($mjd);
    
    // 地方恒星時を計算
    $lst = calculateLocalSiderealTime($gst, $longitude);
    
    // 時角を計算（0から360度の範囲に正規化）
    $hourAngle = fmod($lst * 15, 360);
    if ($hourAngle < 0) {
        $hourAngle += 360;
    }
    
    return $hourAngle;
}

/**
 * 修正ユリウス日を計算する
 */
function calculateModifiedJulianDate($y, $m, $d, $h, $mi, $s) {
    if ($m <= 2) {
        $y -= 1;
        $m += 12;
    }
    $ret = (int)(365.25 * $y) + (int)($y / 400) - (int)($y / 100);
    $ret += (int)(30.59 * ($m - 2)) + $d - 678912;
    $ret += $h / 24 + $mi / 1440 + $s / 86400 - 0.375;
    return $ret;
}

/**
 * グリニッジ恒星時を計算する
 */
function calculateGreenwichSiderealTime($mjd) {
    $t = ($mjd - 51544.5) / 36525;
    $gst = 280.46061837 + 360.98564736629 * ($mjd - 51544.5) + 0.000387933 * $t * $t - $t * $t * $t / 38710000;
    return fmod($gst, 360) / 15; // 時間単位に変換
}

/**
 * 地方恒星時を計算する
 */
function calculateLocalSiderealTime($gst, $longitude) {
    $lst = $gst + $longitude / 15;
    if ($lst < 0) {
        $lst += 24;
    } elseif ($lst >= 24) {
        $lst -= 24;
    }
    return $lst;
}

// 使用例
$year = 2024;
$month = 9;
$day = 11;
$hour = 21;
$minute = 7;
$second = 0;
$longitude = 139.75; // 東京の経度

$hourAngle = calculateHourAngle($year, $month, $day, $hour, $minute, $second, $longitude);
echo "時角: " . $hourAngle . "度";