使用多个值循环时得到不同的结果

Get different results when loop using more than one value

我正在尝试为运输编写一个光变曲线模拟,但是,当使用多个倾角值时(下面代码中的'ibound'),会出现一些奇怪的水平线最后的数字,但是如果每次只使用一个值就不会发生,e.g.ibound=[80,82] 是不行的,但是 ibound=[80] 或 ibound=[82] 给出了正确的结果。

import numpy as np
import math
import matplotlib.pyplot as plt

r_p=1.314*(7.149e7)                 
r_s=0.826*(6.957e8)                 
Area_p=np.pi*r_p**2.
P=1.3*86400.
omega=2.*np.pi/P
a=0.0226*(1.496e11)


def r_cen(t):
  return a*np.sqrt( (np.sin(omega*t)**2) + (
      (np.cos(math.radians(i))**2) *
      (np.cos(omega*t)**2) ) )

def A_case2(rcen):
  Aobs = 0.
  nstep = 200
  for step in range(0,nstep):
    r=rcen-r_p+(float(step)*2.*r_p/float(nstep))
    if r>r_s:
      break
    else:
      theta = ((r**2.)+(rcen**2.)-(r_p**2.)) / (2.*r*rcen)
      if theta > 1.:
          theta = 1.
      Aobs=Aobs+(2.*r*np.arccos(theta)*(r_p*2./float(nstep)))
  return Aobs

LC=[]
phase = []
time = []
A=[]
ibound=[80,82]
for i in ibound:
  for t in range(-int(P/5),int(P/5)):
    rcen=r_cen(float(t))
    phase.append(float(t)/P)
    time.append(float(t))
    if rcen>r_s+r_p:
        A = 0.
    elif r_s>rcen+r_p:
        A = Area_p
    else:
        A = A_case2(rcen)
    LC.append(1.-(A/(np.pi*r_s**2)))
    

yaxis = LC
xaxis = phase
plt.plot(xaxis,yaxis,'r',linewidth=1)
plt.xlim([-0.1,0.1])
plt.xlabel('phase')
plt.ylabel('normalised counts')
plt.show()

它是 plt.plot() 函数的神器。当您将所有结果放入同一个列表时,绘图仪将尝试使用一条连续线连接所有数据点。所以第一遍的结束数据点将连接到第二遍的第一个点。

如果更改绘图参数,您可以清楚地看到这一点,

plt.plot(xaxis, yaxis, 'or', linewidth=1)

这给出了,

如您所见,绘图仪没有尝试用连续线连接数据点,因此没有额外的水平线。

您的 xaxis 值在中途有一个不连续点,因此 xaxis 看起来像(假设步长为 0.1)[-0.2, -0.1...0.1, 0.2, -0.2, -0.1...0.1, 0.2]。正如 DrBwts 所说,这就是 plt.plot 连接不连续点的位置,从而形成一条水平线。

如果您希望绘图仪仍然连接数据点(如果您希望图表在减小线宽后看起来不太起伏,这很有帮助),您可以替换代码中的最后一部分(来自 yaxis = LC起)如下:

yaxis = LC
xaxis = phase

num_elements_per_list = len(phase)//2

x1 = xaxis[:num_elements_per_list]
x2 = xaxis[num_elements_per_list:]

y1 = yaxis[:num_elements_per_list]
y2 = yaxis[num_elements_per_list:]

plt.plot(x1, y1, 'r', linewidth=1)
plt.plot(x2, y2, 'r', linewidth=1)
plt.xlim([-0.1,0.1])

# This line is just to make the x axis a bit prettier
plt.xticks(np.arange(-0.1, 0.101, 0.05))

plt.xlabel('Phase')
plt.ylabel('Normalised Counts')

plt.show()

这会生成以下图:

python

ibound=[80,82,85]  # the size of the list is not important
for i in ibound:
    LC=[]
    phase = []
    time = []
    A=[]
    for t in range(-int(P/5),int(P/5)):
        rcen=r_cen(float(t))
        phase.append(float(t)/P)
        time.append(float(t))
        if rcen>r_s+r_p:
            A = 0.
        elif r_s>rcen+r_p:
            A = Area_p
        else:
            A = A_case2(rcen)
        LC.append(1.-(A/(np.pi*r_s**2)))
        
    yaxis = LC
    xaxis = phase
    plt.plot(xaxis, yaxis, 'r', linewidth=0.5)
    plt.xlim([-0.1, 0.1])
    plt.xlabel('phase')
    plt.ylabel('normalised counts')
plt.show()

如果你改变你的块,这将是你想要的。 为什么它结合了已经描述的所有行

Figure