「インスタンスの開始時にそれらを計算する」または「それらの構造を事前計算してハードコードされたPython構造に出力する」ことを示唆する回答は、1年分の日の出を保存することによって必要となる365倍の乗数、または2000倍の乗数を無視しているようです。計算は、インスタンスの開始時に行われます。pyEphemを使用すると、2000の日の出と日の入りの計算に2秒以上かかります。2000の場所の1年間の日の出と日の入りをソースコードに保存すると、20メガバイト以上が使用される可能性があります。数値を効率的にピクルス化するには、2 * 365 * 2000 * 8=11,680,000バイトが必要です。
より速く、より良く機能するアプローチは、ある場所の時間の最小二乗モデルを他の場所の時間の観点から設定することです。これにより、以下で説明するように、使用される合計スペースを約70分の1に削減できます。
まず、ポイントAとBが同じ緯度にあり、高度と地平線のパラメータが類似している場合、Aでの日の出はBでの日の出に対して一定の時間オフセットで発生します。たとえば、AがBの西15度の場合、日の出はAではBよりも1時間遅れています。次に、ポイントA、B、Cが同じ経度で低緯度にある場合、あるポイントでの日の出時刻は、他の2つの線形の組み合わせとしてかなり正確に計算できます。高緯度または精度を高めるために、いくつかの時間曲線の線形結合を使用できます。第3に、春分日である3月20日のA地点の日の出時刻を正規化ポイントとして使用できるため、すべての計算を同じ緯度に正規化できます。
次の表は、4つの時間曲線の線形結合を使用した場合の精度の結果を示しています。赤道から最大46°離れた経度の場合、結果は約0.5秒以内にとどまります。48°から60°の場合、結果は5秒以内に留まります。64°では、結果に最大2分誤差が生じる可能性があり、65°では、最大約6分になる可能性があります。しかし、これらの時間はおそらくほとんどの実用的な目的には十分です。66°では、pyEphemがスローする例外を処理しないため、以下に示すプログラムが機能しなくなることに注意してください。「AlwaysUpError：'Sun'は2013/6/1407:20:15でまだ地平線上にあります」が発生しますが、66°は北極圏の66.5622°Nを下回っています。
プログラムを変更して、必要な数の時間曲線を使用するようにするのは簡単です（lata = ...プログラムのさまざまなステートメントを参照）。必要な精度が得られますが、より多くの曲線と係数を格納する必要があります。もちろん、モデルは時間曲線のサブセットを使用するように変更できます。たとえば、10個の曲線を保存し、任意のターゲット緯度に最も近い緯度の4個に基づいて計算を行うことができます。ただし、このデモプログラムでは、このような改良は行われていません。
Lat. 0.0: Error range: -0.000000 to 0.000000 seconds
Lat. 5.0: Error range: -0.370571 to 0.424092 seconds
Lat. 10.0: Error range: -0.486193 to 0.557997 seconds
Lat. 15.0: Error range: -0.414288 to 0.477041 seconds
Lat. 20.0: Error range: -0.213614 to 0.247057 seconds
Lat. 25.0: Error range: -0.065826 to 0.056358 seconds
Lat. 30.0: Error range: -0.382425 to 0.323623 seconds
Lat. 35.0: Error range: -0.585914 to 0.488351 seconds
Lat. 40.0: Error range: -0.490303 to 0.400563 seconds
Lat. 45.0: Error range: -0.164706 to 0.207415 seconds
Lat. 47.0: Error range: -0.590103 to 0.756647 seconds
Lat. 48.0: Error range: -0.852844 to 1.102608 seconds
Lat. 50.0: Error range: -1.478688 to 1.940351 seconds
Lat. 55.0: Error range: -3.342506 to 4.696076 seconds
Lat. 60.0: Error range: -0.000002 to 0.000003 seconds
Lat. 61.0: Error range: -7.012057 to 4.273954 seconds
Lat. 62.0: Error range: -21.374033 to 12.347188 seconds
Lat. 63.0: Error range: -51.872753 to 27.853411 seconds
Lat. 64.0: Error range: -124.000365 to 59.661029 seconds
Lat. 65.0: Error range: -351.425224 to 139.656187 seconds
上記のアプローチを使用して、2000の場所ごとに、5つの浮動小数点数を格納する必要があります。3月20日の日の出時刻と、4つの時間曲線の4つの乗数係数です。（前述の70分の1の削減は、365の数値ではなく、場所ごとに5つの数値を格納することによるものです。）各時間曲線について、365の数値が格納され、エントリiは3月20日の日の出の時間差です。4つの時間曲線を格納すると、それらの2000を格納するのと同じ量の1/500のスペースが使用されるため、曲線の格納スペースは乗数係数のスペースによって支配されます。
scipy.optimize.leastsqを使用して係数を解決するプログラムを提供する前に、ipythonインタープリターで精度テーブルを作成し、エラーを視覚化するためのプロットを描画するために使用できる2つのコードスニペットを次に示します。
import sunrise as sr
for lat in range(0, 65, 5):
sr.lsr(lat, -110, 2013, 4)
上記は、前に示したエラーテーブルのほとんどを生成します。の3番目のパラメーターlsrが呼び出されdaySkip、値4はlsr、テストを高速化するために4日ごと（つまり、1年のうち約90日のみ）で機能します。使用sr.lsr(lat, -110, 2013, 1)すると同様の結果が得られますが、4倍の時間がかかります。
sr.plotData(15,1./(24*3600))
上記は、sunrise.plotDataにすべてをプロットするように指示しています（日の出データの近似値、モデルの結果の近似値、秒単位でスケーリングされた残差、および基本曲線）。
プログラムを以下に示します。主に北半球の経度でテストされていることに注意してください。時間曲線が十分に対称である場合、プログラムはそのまま南半球の経度を処理します。エラーが大きすぎる場合は、南半球の経度を基本曲線に追加するか、赤道の南にある別の曲線セットを使用するようにモデルを変更できます。このプログラムでは日没は計算されないことに注意してください。日没の場合は、next_setting(ephem.Sun())呼び出しに類似したprevious_rising(ephem.Sun())呼び出しを追加し、さらに4つの時間曲線を保存します。
#!/usr/bin/python
import ephem, numpy, scipy, scipy.optimize
# Make a set of observers (observation points)
def observers(lata, lona):
def makeIter(x):
if hasattr(x, '__iter__'):
return x
return [x]
lata, lona = makeIter(lata), makeIter(lona)
arr = []
for lat in lata:
for lon in lona:
o = ephem.Observer()
o.lat, o.lon, o.elevation, o.temp = str(lat), str(lon), 1400, 0
arr.append(o)
return tuple(arr)
# Make a year of data for an observer, equinox-relative
def riseData(observr, year, skip):
yy = ephem.Date('{} 00:00'.format(year))
rr = numpy.arange(0.0, 366.0, skip)
springEquinox = 78
observr.date = ephem.Date(yy + springEquinox)
seDelta = observr.previous_rising(ephem.Sun()) - yy - springEquinox + 1
for i, day in enumerate(range(0, 366, skip)):
observr.date = ephem.Date(yy + day)
t = observr.previous_rising(ephem.Sun()) - yy - day + 1 - seDelta
rr[i] = t
return numpy.array(rr)
# Make a set of time curves
def makeRarSet(lata, lona, year, daySkip):
rar=[]
for o in observers(lata, lona):
r = riseData(o, year, daySkip)
rar.append(r)
x = numpy.arange(0., 366., daySkip)
return (x, rar)
# data() is an object that stores curves + results
def data(s):
return data.ss[s]
# Initialize data.ss
def setData(lata, lona, year, daySkip):
x, rar = makeRarSet(lata, lona, year, daySkip)
data.ss = rar
# Compute y values from model, given vector x and given params in p
def yModel(x, p):
t = numpy.zeros(len(x))
for i in range(len(p)):
t += p[i] * data(i)
return t
# Compute residuals, given params in p and data in x, y vectors.
# x = independent var, y = dependent = observations
def residuals(p, y, x):
err = y - yModel(x, p)
return err
# Compute least squares result
def lsr(lat, lon, year, daySkip):
latStep = 13.
lata = numpy.arange(0., 66.4, latStep)
lata = [ 88 * (1 - 1.2**-i) for i in range(8)]
l, lata, lstep, ldown = 0, [], 20, 3
l, lata, lstep, ldown = 0, [], 24, 4
while l < 65:
lata.append(l); l += lstep; lstep -= ldown
#print 'lata =', lata
setData(lata, lon, year, daySkip)
x, ya = makeRarSet(lat, lon, year, daySkip)
x, za = makeRarSet(lat, 10+lon, year, daySkip)
data.ss.append(za[0])
y = ya[0]
pini = [(0 if abs(lat-l)>latStep else 0.5) for l in lata]
pars = scipy.optimize.leastsq(residuals, pini, args=(y, x))
data.x, data.y, data.pv = x, y, yModel(x, pars[0])
data.par, data.err = pars, residuals(pars[0], y, x)
#print 'pars[0] = ', pars[0]
variance = numpy.inner(data.err, data.err)/len(y)
#print 'variance:', variance
sec = 1.0/(24*60*60)
emin, emax = min(data.err), max(data.err)
print ('Lat. {:4.1f}: Error range: {:.6f} to {:.6f} seconds'.format(lat, emin/sec, emax/sec))
def plotData(iopt, emul):
import matplotlib.pyplot as plt
plt.clf()
x = data.x
if iopt == 0:
iopt = 15
emul = 1
if iopt & 1:
plt.plot(x, data.y)
plt.plot(x, data.y + 0.001)
plt.plot(x, data.y - 0.001)
if iopt & 2:
plt.plot(x, data.pv)
if iopt & 4:
plt.plot(x, emul*data.err)
if iopt & 8:
for ya in data.ss:
plt.plot(x, ya)
plt.show()