光变曲线计算学习

系外行星光变曲线文章［Mandel & Agol 2002］在 http://iopscience.iop.org/article/10.1086/345520/pdf,

在作者的算法介绍的网站 http://www.lpl.arizona.edu/~ianc/python/transit.html ，计算几个例子

transit.occultnonlin(z, p0, cn)的例子

# Reproduce Figure 2 of Mandel & Agol (2002):
from pylab import *
import transit
z = linspace(0, 1.2, 50)
cns = vstack((zeros(4), eye(4)))
figure()
for coef in cns:
    f = transit.occultnonlin(z, 0.1, coef)
    plot(z, f)

有误，会提示除零。

transit.py:1076: RuntimeWarning: divide by zero encountered in divide
  k = 0.5 * np.sqrt(-am1 / (z * p))
transit.py:1076: RuntimeWarning: invalid value encountered in sqrt
  k = 0.5 * np.sqrt(-am1 / (z * p))
transit.py:1150: RuntimeWarning: invalid value encountered in multiply
  special.hyp2f1(0.5, 0.5, 2.5 + nd4, 0.25/p2)).sum()
Traceback (most recent call last):
  File "<stdin>", line 2, in <module>
  File "transit.py", line 1154, in occultnonlin
    F[i08] =  -(1. / (np.pi*twoOmega)) * (N[:, i02] * cc/(nn + 4.) ).sum(0)
ValueError: NumPy boolean array indexing assignment cannot assign 8 input values to the 0 output values where the mask is true

但下面的transit.occultnonlin_small(z, p, cn)例子

# Reproduce Figure 2 of Mandel & Agol (2002):
from pylab import *
import transit
z = linspace(0, 1.2, 100)
cns = vstack((zeros(4), eye(4)))
figure()
for coef in cns:
    f = transit.occultnonlin_small(z, 0.1, coef)
    plot(z, f, '--')

可以运行。 先要从[source] 链接的网站http://www.lpl.arizona.edu/~ianc/python/_modules/transit.html#occultnonlin
下载网页版的代码。注意要去掉网页的注释,比如[doc]什么的，不然python会编译不过的。存储为transit.py

把给的例子略为修改一下，以便图形输出

import matplotlib;
import matplotlib.pyplot as plt;
# Reproduce Figure 2 of Mandel & Agol (2002):
from pylab import *
import transit
z = linspace(0, 1.2, 100)
cns = vstack((zeros(4), eye(4)))
figure()
for coef in cns:
    f = transit.occultnonlin_small(z, 0.1, coef)
    plot(z, f, '--')

plt.show();

就能产生［Mandel & Agol 2002］中的图2

figure_1.png

这是一个好的开始。