Das deutsche Python-Forum

Code: Alles auswählen

# -*- coding: utf-8 -*-
from PyAstronomyKepler import keplerAccess as ka
import matplotlib.pylab as mpl
import numpy as np

kdb=ka.KeplerDB()
kpc=ka.KeplerPlanCan()
ac=kpc.getAllCandidates()

nobj=2
print "ac[nobj]", ac[nobj] #anders zu readKepler_3
nplanets = len(ac)
print "nplanets", nplanets

klc=ka.KeplerMultiLC(ac[nobj])
lcd=klc.getLCData(kid=ac[nobj], quarter=3, slc=True, cols=["TIME", "PDCSAP_FLUX"], verbose=False)
#anders zu readKepler_3

print len(lcd)
print lcd[0][:,1]

c=kpc.getCandidateInfo(kid=ac[nobj]) #anders zu readKepler_3
print c
print "bjdT0", c[0]["bjdT0"]
print "period", c[0]["period"]

t0=c[0]["bjdT0"]
print "t0", t0
P=c[0]["period"]
print "P", P


for i in range(len(lcd)):
  x=lcd[i][:,0]
  y=lcd[i][:,1]
  
  print "xmax, xmin", x.max(),x.min()
  ntransit=np.ceil((x[-1]-x[0])/P)
  print "ntransit", ntransit
  transit_offset=np.ceil((x.min()-t0)/P)
  print "transit_off", transit_offset
  xtransit=(np.arange(ntransit)+transit_offset)*P+t0
  xdelta=c[0]["dur"]/48. #anders zu readKepler_3 (dur in h, x in days, dur kompletter transit)
  print "xdelta", xdelta
  
  #a) 1. Weg Transits entfernen
  xwot=x
  ywot=y
  for xt in xtransit:
    ind=np.where(np.logical_and(xwot >= xt-xdelta, xwot <= xt+xdelta))[0]
    #mpl.plot(x[ind], y[ind], 'ks', markersize=3)
    xwot=np.delete(xwot,ind)
    ywot=np.delete(ywot,ind)
  print "xtransit", xtransit
  #anders zu readKepler_3
  #mpl.plot(xwot, ywot, 'ms', markersize=7)
  
  #b) 2. Weg Transits entfernen
  #xwot=x
  #ywot=y
  #for xt in xtransit:
    #ind=np.where(np.logical_or(xwot <= xt-xdelta, xwot >= xt+xdelta))[0]
    #xwot=xwot[ind]
    #ywot=ywot[ind]
  #mpl.plot(xwot, ywot, 'ms', markersize=7)
  
  #mpl.plot(x, y)
  #mpl.show()
  #continue
  #exit()
  
  #fit für inaktive Sterne
  #reg=np.polyfit(xwot, ywot, 2)
  #fit=np.polyval(reg, x)
  #fitwot=np.polyval(reg,xwot)
  #res=ywot-fitwot
  #std=np.std(res)
  #ausr=np.where(np.logical_or(res >= 3*std, res <= -3*std))[0]
  ##mpl.plot(xwot, ywot)
  ##mpl.plot(xwot[ausr], res[ausr], 'mp')
  #xwotwoa=np.delete(xwot, ausr)
  #ywotwoa=np.delete(ywot, ausr)
  ##mpl.plot(xwotwoa, ywotwoa)
  #regbest=np.polyfit(xwotwoa, ywotwoa, 2)
  #fitbest=np.polyval(regbest, x)
  #mpl.plot(xwotwoa, fitbest, "r-")
  #mpl.plot(x, fit, "r-")
  #mpl.show()
  #exit()
  
  
  #fit für aktive Sterne
  #xwota=x
  #ywota=y
  for xtact1 in xtransit:
    indtransitandcon=np.where(np.logical_and(x >= xtact1-1.5*xdelta, x <= xtact1+1.5*xdelta))[0]
    #xwota=np.delete(xwota, indtransitandcon)
    #ywota=np.delete(ywota, indtransitandcon)
    xtransitandcon=np.append(x, indtransitandcon)
    ytransitandcon=np.append(y, indtransitandcon)
  mpl.plot(xtransitandcon, ytransitandcon)
mpl.show()
exit()
  
  xonlycon=xtransitandcon
  yonlycon=ytransitandcon
  for xtact2 in xtransit:
    indonlytran=np.where(np.logical_or(xonlycon <= xtact2-xdelta, xonlycon >= xtact2+xdelta))
    xonlycon=np.delete(xonlycon, indonlytran)
    yonlycon=np.delete(yonlycon, indonlytran)
  mpl.plot(xonlycon, yonlycon)
mpl.show()
   
  regact=np.polyfit(xonlycon, yonlycon, 2)
  fitact=np.polyval(regact, xonlycon)
  resact=yonlycon-fitact
  stdact=np.std(resact)
  ausract=np.where(np.logical_or(resact >= 3*stdact, resact <= -3*stdact))[0]
  #mpl.plot(xwot[ausr], res[ausr], 'mp')
  xonlyconwoa=np.delete(xonlycon, ausract)
  yonlyconwoa=np.delete(yonlycon, ausract)
  #mpl.plot(xwotwoa, ywotwoa)
  regbestact=np.polyfit(xonlyconwoa, yonlyconwoa, 2)
  fitbestact=np.polyval(regbestact, xtransitandcon) 
  
  ytransitandcon /= fitbestact
  print "xtransitandcon", xtransitandcon
  print "ytransitandcon", ytransitandcon
  mpl.plot(xtransitandcon, ytransitandcon)
  
  #y /= fitbest
  #print "x", x
  #print "y", y
  
  #ywot /= np.median(ywot)
  #for k in in range(len(x))
    #ysteig=(y[-1]-y[0])/(x[-1]-x[0])*x
    #y /= ysteig
  
  #mpl.plot(x, y)
  #for j in range(len(xtransit)): mpl.axvline(xtransit[j])
  #mpl.plot(xtransit, ytransit, color='k')
  #mpl.xlabel('time')
  #mpl.ylabel('flux')
  #mpl.title('Lichtkurve von kid=ac[0]')
  #mpl.xlim(0, 10)
  #mpl.ylim(-1.5, 2.0)
mpl.show()