不知道有先進朋友有沒有寫過，日出日落類似的方程式，以下是我找到的類似文章，如果
有可以提供一下參考！謝謝...
這裡提供一個簡便的方法來計算日出與日沒時間。為了說明方便，犧牲了計算的精密度。
用此方法計算日出與日沒時間，大約會有4分鐘以內的誤差。
步驟一
定義日出時間的時角（Hour Angle）
Z＝90? 實際的日出，但天色已微亮
Z＝96? 民用的日出，一般為街道熄燈的參考時刻
Z＝102? 航海用日出時間
Z＝108? 天文用日出沒時間，此時天空全黑
步驟二
修正時角
由於太陽有0.5度的視角，因此它的中心點位於地平下1/4度時即可看到，因此實際日出的
時角必須加以修正，Z＝90.5?
此外因為大氣折射的影響，在地平線下34'的位置，會讓我們以為它正處於地平線的位置
。因此實際日出的時角必須再次加以修正，Z＝91.067?
至於其它的選項則不需修正日出時間的時角。
步驟三
觀測當日，太陽中午時刻的赤道經度與緯度 （α，δ） ( 請參考「太陽位置的計算」
網頁)
觀測者地理經緯度 （λ，ψ）
觀測當日的均時差 E ( 計算方法可參考「時間」的網頁 )
步驟四
τ＝cos-1{ [ cos(Z)－sin(δ) sin(ψ)] / [ cos(ψ) cos(δ)] } /15
日沒時間＝12＋τ＋（Υ－λ）/15－E/60 (小時)
日出時間＝日沒時間＋24－2τ
說明：Υ－λ之中的Υ是指地方時間的地理經度。例如台灣位於東經121?31'38''，Υ＝1
20陛Aλ＝121?31'38''。
範例
2000年9月1日的實際日出與日沒時間？
步驟一
選擇實際的日出 Z＝90?
步驟二
修正時角
由於太陽有0.5度的視角，因此它的中心點位於地平下1/4度時即可看到，因此實際日出的
時角必須加以修正，Z＝90.5?
此外因為大氣折射的影響，在地平線下34'的位置，會讓我們以為它正處於地平線的位置
。因此實際日出的時角必須再次加以修正，Z＝91.067?
步驟三
觀測當日，太陽中午時刻的赤道經度與緯度 （10h42m29s，8?10'59''）
觀測者地理經緯度 （121?31'38''，25? 4'40''） ( 以台北為例 )
觀測當日的均時差 E＝－0.21分
步驟四
τ＝cos-1{ [ cos(91.067?)－sin(8?10'59'') sin(25? 4'40'')] / [ cos(25? 4'40'')
cos(8?10'59'')] } /15
τ＝6.267 hour
日沒時間＝18h10m
日出時間＝ 5h38m
說明：由於我們計算太陽的出沒時刻是以當日中午的太陽位置為基礎，事實上太陽與均時
差值無時無刻不在天球上移動，因此我們可以根據初次計算出來的時間，再次修正太陽的
實際位置，然後在依此方法重複計算。約三次的修正後，可以得到不錯的精確度，時間誤
差應該在分鐘以內。結果如下：
日沒時間＝18h12m41s
日出時間＝ 5h34m42s
...................
Sunrise/Sunset Algorithm
Source:
Almanac for Computers, 1990
published by Nautical Almanac Office
United States Naval Observatory
Washington, DC 20392
Inputs:
day, month, year: date of sunrise/sunset
latitude, longitude: location for sunrise/sunset
zenith: Sun's zenith for sunrise/sunset
offical = 90 degrees 50'
civil = 96 degrees
nautical = 102 degrees
astronomical = 108 degrees
NOTE: longitude is positive for East and negative for West
NOTE: the algorithm assumes the use of a calculator with the
trig functions in "degree" (rather than "radian") mode. Most
programming languages assume radian arguments, requiring back
and forth convertions. The factor is 180/pi. So, for instance,
the equation RA = atan(0.91764 * tan(L)) would be coded as RA
= (180/pi)*atan(0.91764 * tan((pi/180)*L)) to give a degree
answer with a degree input for L.
1. first calculate the day of the year
N1 = floor(275 * month / 9)
N2 = floor((month + 9) / 12)
N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3))
N = N1 - (N2 * N3) + day - 30
2. convert the longitude to hour value and calculate an approximate time
lngHour = longitude / 15
if rising time is desired:
t = N + ((6 - lngHour) / 24)
if setting time is desired:
t = N + ((18 - lngHour) / 24)
3. calculate the Sun's mean anomaly
M = (0.9856 * t) - 3.289
4. calculate the Sun's true longitude
L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634
NOTE: L potentially needs to be adjusted into the range [0,360) by adding/sub
tracting 360
5a. calculate the Sun's right ascension
RA = atan(0.91764 * tan(L))
NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/su
btracting 360
5b. right ascension value needs to be in the same quadrant as L
Lquadrant = (floor( L/90)) * 90
RAquadrant = (floor(RA/90)) * 90
RA = RA + (Lquadrant - RAquadrant)
5c. right ascension value needs to be converted into hours
RA = RA / 15
6. calculate the Sun's declination
sinDec = 0.39782 * sin(L)
cosDec = cos(asin(sinDec))
7a. calculate the Sun's local hour angle
cosH = (cos(zenith) - (sinDec * sin(latitude))) / (cosDec * cos(latitude))
if (cosH > 1)
the sun never rises on this location (on the specified date)
if (cosH < -1)
the sun never sets on this location (on the specified date)
7b. finish calculating H and convert into hours
if if rising time is desired:
H = 360 - acos(cosH)
if setting time is desired:
H = acos(cosH)
H = H / 15
8. calculate local mean time of rising/setting
T = H + RA - (0.06571 * t) - 6.622
9. adjust back to UTC
UT = T - lngHour
NOTE: UT potentially needs to be adjusted into the range [0,24) by adding/sub
tracting 24
10. convert UT value to local time zone of latitude/longitude
localT = UT + localOffset........
請輸入專案類型(網站專案或者應用程式專案)：