[FIXED] Get different results when loop using more than one value

Issue

I'm trying to write a light curve simulation for transit, however, when use more than one value for the inclination ('ibound' in the code below), there will be some weird horizontal line in the final figure, but it will not happen if only one value is used each time, e.g.ibound=[80,82] is not ok, but ibound=[80] or ibound=[82] gives the correct result.

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()

Solution

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()

if you change your block on this, will be what you wanted. why it combines all the lines already described

Figure

Answered By - Dmitrij