题目描述

题目计算日出后指定时刻太阳圆盘可见部分的面积比例。这是一个涉及一个行星和一个恒星的二体问题。

假设行星为完美球体，半径为 $3950$ 英里。恒星（太阳）位于距行星中心 $92,900,000$ 英里处，被视为一个圆盘，其平面垂直于行星中心与太阳圆盘中心的连线。太阳圆盘半径为 $432,000$ 英里。行星均匀自转，每 $24$ 小时转一圈。太阳始终位于行星赤道的正上方。

时间从日出时刻（太阳第一缕光线到达赤道上给定参考点）开始测量，以秒为单位。对于每个输入的时间值，计算此时太阳圆盘照亮该参考点的面积比例。

输入格式

输入包含多个浮点数，每行一个，表示以秒为单位的时间值。时间值在 $0$ 到 $600$ 之间（含）。输入以文件结束符（$\text{EOF}$）终止。

输出格式

对于每个输入的时间值，输出一行一个浮点数，表示太阳圆盘可见部分的面积比例。如果比例为零，输出值应在 $[-0.001,0.001]$ 区间内。非零答案的误差需在 $0.1\%$ 以内。

样例

输入

0.0
600.0

输出

0.000000
1.000000

题目分析

本题的核心是计算太阳圆盘被行星遮挡部分的面积比例。由于太阳的角直径远大于行星的视直径，日出时太阳并非完全被遮挡，而是逐渐露出。

几何模型

设：

• 行星半径 $R_e=3950$ 英里
• 太阳半径 $R_s=432000$ 英里
• 行星中心到太阳中心距离 $D=92900000$ 英里
• 地球自转周期 $T=24$ 小时 $=86400$ 秒

日出时刻，太阳圆盘的上边缘恰好与行星表面在参考点处相切。此时，从行星中心到参考点与到太阳圆盘上边缘的视线之间的夹角为 $\alpha$。

角度 $\alpha$ 的计算

在日出时刻，由几何关系可得：
$$L=\sqrt{D^2+R_s^2-R_e^2}$$
该距离表示从行星中心到太阳圆盘上边缘切点的距离。

太阳圆盘对行星中心的半张角为 ${\theta }_s=\arctan (R_s/D)$，行星对太阳圆盘上边缘的切线半角为 ${\theta }_t=\arctan (L/R_e)$。因此：
$$\alpha ={\theta }_s+{\theta }_t=\arctan \left(\frac{R_s}{D}\right)+\arctan \left(\frac{L}{R_e}\right)$$

时间与角度的关系

行星自转角速度为：
$$\omega =\frac{2\pi }{86400}\text{ 弧度/秒}$$
设从日出时刻起经过 $t$ 秒，行星自转过的角度为 $\omega t$。此时，参考点相对于太阳中心的方向角为：
$$\beta =\alpha -\omega t$$

太阳圆盘可见部分

设 $h$ 为行星遮挡边缘与太阳圆盘中心的垂直距离。若行星表面在参考点处的切线与太阳中心方向的夹角为 $\beta$，则有：
$$h=R_e-D\cos \beta /\sin \beta$$

但更简洁的推导基于圆被直线切割的面积公式。设太阳圆盘被一条距离圆心为 $|h|$ 的直线切割，则切割后较小部分的面积比例为：
$$f(|h|)=\frac{1}{\pi }\arccos \left(\frac{|h|}{R_s}\right)-\frac{|h|}{\pi R_s^2}\sqrt{R_s^2-h^2}$$

• 若 $h\ge R_s$，太阳完全被遮挡，比例为 $0$。
• 若 $h\le -R_s$，太阳完全可见，比例为 $1$。
• 若 $-R_s<h<0$，可见部分为圆盘的较大部分，比例 $=1-f(|h|)$。
• 若 $0\le h<R_s$，可见部分为圆盘的较小部分，比例 $=f(h)$。

公式推导说明

由于 $D\gg R_e,R_s$，太阳的视张角很小，但本题要求精确计算，因此不能使用小角度近似。上述公式基于严格的几何关系。

复杂度分析

每组输入仅需常数次浮点运算，时间复杂度 $O(1)$。

代码实现

// Sunrise
// UVa ID: 415
// Verdict: Accepted
// Submission Date: 2016-08-05
// UVa Run Time: 0.000s
//
// 版权所有（C）2016，邱秋。metaphysis # yeah dot net

#include <bits/stdc++.h>

using namespace std;

const double pi = 2 * acos(0);

int main(int argc, char *argv[])
{
    cin.tie(0); cout.tie(0); ios::sync_with_stdio(false);

    double Sr = 432000, Er = 3950, D = 92900000;
    
    double L = sqrt(D * D + Sr * Sr - Er * Er);
    double alpha = atan(Sr / D) + atan(L / Er);
    
    double seconds;
    while (cin >> seconds)
    {
        double h, beta = alpha - 2.0 * pi * seconds / 3600.0 / 24.0;

        h = (Er - D * cos(beta)) / sin(beta);
        
        if (h >= Sr)
        {
            cout << "0.000000\n";
            continue;
        }
        
        if (h <= -Sr)
        {
            cout << "1.000000\n";
            continue;
        }
        
        double percentage = acos(fabs(h) / Sr) / pi - fabs(h) * sqrt(Sr * Sr - h * h) / (pi * Sr * Sr);
        if (h < 0)
            percentage = 1.0 - percentage;
        cout << fixed << setprecision(6) << percentage << '\n';
    }
    
    return 0;
}