Mau*_*itz 1 math geometry r astronomy swift3
我正在尝试实现这个解决方案来计算 Swift3 中太阳的位置。然后我将其包装在另一个函数中，该函数简单地循环一天，从午夜每 10 分钟步进一次，直到 23:50。
我并不真正理解 R，并且我没有完全理解答案的一些细节，特别是似乎是带有方括号的某种 if/clamp 函数。当我感到困惑时，我尽力与Python版本进行比较。否则，唯一的区别是由于使用了NSDate，它简化了顶部的一些代码。
我得到的一些值似乎是正确的，当我绘制结果时，我可以看到曲线的基础。然而，一次调用（例如上午 7 点）和下一次调用（7 点 10 分）的结果截然不同。
我强烈怀疑我在钳位方面做错了什么，并且输入中的微小变化会以不同的方式进行修改/截断并影响输出。但我看不出来。任何了解这个算法的人都可以帮忙吗？
这是我得到的输出示例：
2017-06-21 00:10:00 +0000 -16.0713262209521 31.7135341633943
2017-06-21 00:20:00 +0000 61.9971433936385 129.193513530349
2017-06-21 00:30:00 +0000 22.5263575559266 78.5445189561018
2017-06-21 00:40:00 +0000 29.5973897349096 275.081637736092
2017-06-21 00:50:00 +0000 41.9552795956374 262.989819486864
正如您所看到的，它在迭代之间剧烈波动。地球不是那样转的！我的代码如下，这个版本只是将结果发送到日志：
class func julianDayFromDate(_ date: Date) -> Double {
let ti = date.timeIntervalSince1970
return ((ti / 86400.0) + 2440587)
}
class func sunPath(lat: Double, lon: Double, size: CGSize) -> UIImage {
var utzCal = Calendar(identifier: .gregorian)
utzCal.timeZone = TimeZone(secondsFromGMT: 0)!
let year = utzCal.component(.year, from: Date())
let june = DateComponents(calendar: utzCal, year: year, month: 6, day: 21).date!
// now we loop for every 10 minutes (2 degrees) and plot those points
for time in stride(from:0, to:(24 * 60), by: 10) {
let calcdate = june.addingTimeInterval(Double(time) * 60.0)
let (alt, az) = sun(date: calcdate, lat: lat, lon: lon)
print(calcdate, alt, az)
}
class func sun(date: Date, lat: Double, lon: Double) -> (altitude: Double, azimuth: Double) {
// these come in handy
let twopi = Double.pi * 2
let deg2rad = Double.pi / 180.0
// latitude to radians
let lat_radians = lat * deg2rad
// the Astronomer's Almanac method used here is based on Epoch 2000, so we need to
// convert the date into that format. We start by calculating "n", the number of
// days since 1 January 2000
let n = julianDayFromDate(date) - 2451545.0
// it continues by calculating the position in ecliptic coordinates,
// starting with the mean longitude of the sun in degrees, corrected for aberation
var meanlong_degrees = 280.460 + (0.9856474 * n)
meanlong_degrees = meanlong_degrees.truncatingRemainder(dividingBy: 360.0)
// and the mean anomaly in degrees
var meananomaly_degrees = 357.528 + (0.9856003 * n)
meananomaly_degrees = meananomaly_degrees.truncatingRemainder(dividingBy: 360.0)
let meananomaly_radians = meananomaly_degrees * deg2rad
// and finally, the eliptic longitude in degrees
var elipticlong_degrees = meanlong_degrees + (1.915 * sin(meananomaly_radians)) + (0.020 * sin(2 * meananomaly_radians))
elipticlong_degrees = elipticlong_degrees.truncatingRemainder(dividingBy: 360.0)
let elipticlong_radians = elipticlong_degrees * deg2rad
// now we want to convert that to equatorial coordinates
let obliquity_degrees = 23.439 - (0.0000004 * n)
let obliquity_radians = obliquity_degrees * deg2rad
// right ascention in radians
let num = cos(obliquity_radians) * sin(elipticlong_radians)
let den = cos(elipticlong_radians)
var ra_radians = atan(num / den)
ra_radians = ra_radians.truncatingRemainder(dividingBy: Double.pi)
if den < 0 {
ra_radians = ra_radians + Double.pi
} else if num < 0 {
ra_radians = ra_radians + twopi
}
// declination is simpler...
let dec_radians = asin(sin(obliquity_radians) * sin(elipticlong_radians))
// and from there, to local coordinates
// start with the UTZ sidereal time
let cal = Calendar.current
let h = Double(cal.component(.hour, from: date))
let m = Double(cal.component(.minute, from: date))
let f: Double
if h == 0 && m == 0 {
f = 0.0
} else if h == 0 {
f = 60.0 / m
} else if h == 0 {
f = 24.0 / h
} else {
f = (24.0 / h) + (60.0 / m)
}
var utz_sidereal_time = 6.697375 + 0.0657098242 * n + f
utz_sidereal_time = utz_sidereal_time.truncatingRemainder(dividingBy: 24.0)
// then convert that to local sidereal time
var localtime = utz_sidereal_time + lon / 15.0
localtime = localtime.truncatingRemainder(dividingBy: 24.0)
var localtime_radians = localtime * 15.0 * deg2rad
localtime_radians = localtime.truncatingRemainder(dividingBy: Double.pi)
// hour angle in radians
var hourangle_radians = localtime_radians - ra_radians
hourangle_radians = hourangle_radians.truncatingRemainder(dividingBy: twopi)
// get elevation in degrees
let elevation_radians = (asin(sin(dec_radians) * sin(lat_radians) + cos(dec_radians) * cos(lat_radians) * cos(hourangle_radians)))
let elevation_degrees = elevation_radians / deg2rad
// and azimuth
let azimuth_radians = asin( -cos(dec_radians) * sin(hourangle_radians) / cos(elevation_radians))
// now clamp the output
let azimuth_degrees: Double
if (sin(dec_radians) - sin(elevation_radians) * sin(lat_radians) < 0) {
azimuth_degrees = (Double.pi - azimuth_radians) / deg2rad
} else if (sin(azimuth_radians) < 0) {
azimuth_degrees = (azimuth_radians + twopi) / deg2rad
} else {
azimuth_degrees = azimuth_radians / deg2rad
}
return (elevation_degrees, azimuth_degrees)
}
好吧，在下载了 OSX 的 R 解释器后，发现它没有调试器，发现有多种方法可以进行打印，但都有各自的注意事项，等等，我发现了我正在寻找的问题。它确实错误地限制了其中一个值。这是一个有效的 Swift3 版本，它应该很容易转换为任何类似 C 的语言，并且比原始版本更容易阅读。您必须提供您自己的前两个函数版本，这些函数适用于目标平台的日期格式。这truncatingRemainer是某人的大胆想法，即 Double 上不应该有 % 运算符，这是正常的MOD。
// convinience method to return a unit-epoch data from a julian date
class func dateFromJulianDay(_ julianDay: Double) -> Date {
let unixTime = (julianDay - 2440587) * 86400.0
return Date(timeIntervalSince1970: unixTime)
}
class func julianDayFromDate(_ date: Date) -> Double {
//==let JD = Integer(365.25 * (Y + 4716)) + Integer(30.6001 * (M +1)) +
let ti = date.timeIntervalSince1970
return ((ti / 86400.0) + 2440587.5)
}
// calculate the elevation and azimuth of the sun for a given date and location
class func sun(date: Date, lat: Double, lon: Double) -> (altitude: Double, azimuth: Double) {
// these come in handy
let twopi = Double.pi * 2
let deg2rad = Double.pi / 180.0
// latitude to radians
let lat_radians = lat * deg2rad
// the Astronomer's Almanac method used here is based on Epoch 2000, so we need to
// convert the date into that format. We start by calculating "n", the number of
// days since 1 January 2000. So if your date format is 1970-based, convert that
// a pure julian date and pass that in. If your date is 2000-based, then
// just let n = date
let n = julianDayFromDate(date) - 2451545.0
// it continues by calculating the position in ecliptic coordinates,
// starting with the mean longitude of the sun in degrees, corrected for aberation
var meanlong_degrees = 280.460 + (0.9856474 * n)
meanlong_degrees = meanlong_degrees.truncatingRemainder(dividingBy: 360.0)
// and the mean anomaly in degrees
var meananomaly_degrees = 357.528 + (0.9856003 * n)
meananomaly_degrees = meananomaly_degrees.truncatingRemainder(dividingBy: 360.0)
let meananomaly_radians = meananomaly_degrees * deg2rad
// and finally, the eliptic longitude in degrees
var elipticlong_degrees = meanlong_degrees + (1.915 * sin(meananomaly_radians)) + (0.020 * sin(2 * meananomaly_radians))
elipticlong_degrees = elipticlong_degrees.truncatingRemainder(dividingBy: 360.0)
let elipticlong_radians = elipticlong_degrees * deg2rad
// now we want to convert that to equatorial coordinates
let obliquity_degrees = 23.439 - (0.0000004 * n)
let obliquity_radians = obliquity_degrees * deg2rad
// right ascention in radians
let num = cos(obliquity_radians) * sin(elipticlong_radians)
let den = cos(elipticlong_radians)
var ra_radians = atan(num / den)
ra_radians = ra_radians.truncatingRemainder(dividingBy: Double.pi)
if den < 0 {
ra_radians = ra_radians + Double.pi
} else if num < 0 {
ra_radians = ra_radians + twopi
}
// declination is simpler...
let dec_radians = asin(sin(obliquity_radians) * sin(elipticlong_radians))
// and from there, to local coordinates
// start with the UTZ sidereal time, which is probably a lot easier in non-Swift languages
var utzCal = Calendar(identifier: .gregorian)
utzCal.timeZone = TimeZone(secondsFromGMT: 0)!
let h = Double(utzCal.component(.hour, from: date))
let m = Double(utzCal.component(.minute, from: date))
let f: Double // universal time in hours and decimals (not days!)
if h == 0 && m == 0 {
f = 0.0
} else if h == 0 {
f = m / 60.0
} else if m == 0 {
f = h
} else {
f = h + (m / 60.0)
}
var utz_sidereal_time = 6.697375 + 0.0657098242 * n + f
utz_sidereal_time = utz_sidereal_time.truncatingRemainder(dividingBy: 24.0)
// then convert that to local sidereal time
var localtime = utz_sidereal_time + lon / 15.0
localtime = localtime.truncatingRemainder(dividingBy: 24.0)
let localtime_radians = localtime * 15.0 * deg2rad
// hour angle in radians
var hourangle_radians = localtime_radians - ra_radians
hourangle_radians = hourangle_radians.truncatingRemainder(dividingBy: twopi)
// get elevation in degrees
let elevation_radians = (asin(sin(dec_radians) * sin(lat_radians) + cos(dec_radians) * cos(lat_radians) * cos(hourangle_radians)))
let elevation_degrees = elevation_radians / deg2rad
// and azimuth
let azimuth_radians = asin( -cos(dec_radians) * sin(hourangle_radians) / cos(elevation_radians))
// now clamp the output
let azimuth_degrees: Double
if (sin(dec_radians) - sin(elevation_radians) * sin(lat_radians) < 0) {
azimuth_degrees = (Double.pi - azimuth_radians) / deg2rad
} else if (sin(azimuth_radians) < 0) {
azimuth_degrees = (azimuth_radians + twopi) / deg2rad
} else {
azimuth_degrees = azimuth_radians / deg2rad
}
// all done!
return (elevation_degrees, azimuth_degrees)
}