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()
#print "ac", ac


kepid=[11446443, 8554498]

for nobj in kepid:
  #print "ac[nobj]", ac[nobj] #anders zu readKepler_3
  nplanets = len(ac)
  print "nplanets", nplanets

  klc=ka.KeplerMultiLC(nobj)
  lcd=klc.getLCData(kid=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=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
  a=c[0]["a\\rs"]
  print "a in RS", a
  b=c[0]["b"]
  print "b", b
  RS=c[0]["rs"]
  print "RS in R*", RS
  ibo=np.arccos(b/a)
  print "inc in Bogenmaß", ibo
  inc=ibo/(2*np.pi)*360.
  print "inc in Grad", inc
  p=c[0]["rp\\rs"]
  print "p=PR/SR", p
  dur=c[0]["dur"]
  print "dur in h", dur



  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)
    #regbest=np.polyfit(xwotwoa, ywotwoa, 2)
    #fitbest=np.polyval(regbest, x)
    #mpl.plot(xwotwoa, fitbest, "r-")
    #mpl.plot(x, fit, "r-")
    #mpl.show()
      
    #fit für aktive Sterne
    #xwota=x
    #ywota=y
    indtransits=[]
    for xtact1 in xtransit:
      indtransitandcon=np.where(np.logical_and(x >= xtact1-1.5*xdelta, x <= xtact1+1.5*xdelta))[0]
      indtransits.extend(indtransitandcon)
      #xwota=np.delete(xwota, indtransitandcon)
      #ywota=np.delete(ywota, indtransitandcon)
      #xtransitandcon=np.append(x, indtransitandcon)
      #ytransitandcon=np.append(y, indtransitandcon)
      #
      #mpl.show()
    xtransitandcon=x[indtransits]
    ytransitandcon=y[indtransits]
  #mpl.plot(xtransitandcon, ytransitandcon, "b")
  #mpl.show()
  #mpl.plot(x, y)
  #mpl.show()
  #exit()
    
  xonlycon=xtransitandcon
  yonlycon=ytransitandcon
  for xtact2 in xtransit:
    indonlytran=np.where(np.logical_and(xonlycon >= xtact2-xdelta, xonlycon <= xtact2+xdelta))
    xonlycon=np.delete(xonlycon, indonlytran)
    yonlycon=np.delete(yonlycon, indonlytran)
    #regact=np.polyfit(xonlycon, yonlycon, 2)
    #fitact=np.polyval(regact, xonlycon)
    #yonlycon /= fitact
    #mpl.plot(xonlycon, yonlycon, "gp")
    #mpl.show()
  #exit()

  xnorm=xtransitandcon
  ynorm=ytransitandcon
  yerr = np.zeros(ynorm.size)

  for xtact3 in xtransit:
    contind=np.where(np.logical_and(xonlycon > xtact3-3*xdelta, xonlycon < xtact3+3*xdelta))
    if len(contind[0]) == 0: continue
    regact=np.polyfit(xonlycon[contind], yonlycon[contind], 2)

    ind=np.where(np.logical_and(xtransitandcon > xtact3-3*xdelta, xtransitandcon < xtact3+3*xdelta))  
    fitact=np.polyval(regact, xtransitandcon[ind])  
    ynorm[ind] /= fitact
    #mpl.plot(xtransitandcon, ynorm)
    #mpl.show()
    
    res=ynorm[contind]-1.0
    std=np.std(res)
    yerr[ind] += std
    
    #indau=np.where(np.logical_or(res >= 1.5*std, res <= -1.5*std))[0]
    #indnorm=np.delete(ind, indau)
    #xonlyconwoa=np.delete(xonlycon, indau)
    #yonlyconwoa=np.delete(yonlycon, indau)
    #regactwoa=np.polyfit(xonlyconwoa, yonlyconwoa, 2)
    #fitactwoa=np.polyval(regactwoa, xnorm[indnorm])
    #ynorm[indnorm] /= fitactwoa
    
  #mpl.errorbar(xnorm, ynorm, yerr=yerr, fmt='bp', ms=2)
  #mpl.show()
  #exit()


  #mpl.plot(x, y)
  #mpl.show()
  #exit()

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








  # ... and now the forTrans module
  from PyAstronomy.modelSuite import forTrans as ft
  # Load a simple timer for demonstration purposes
  import time as timerMod

  # Create MandelAgolLC object
  ma = ft.MandelAgolLC()

  # Set parameters
  ma["per"] = P
  #ma["i"] = 90.0
  ma["i"] = inc
  ma["a"] = a
  ma["T0"] = float(t0)
  ma["p"] = p #0.3
  ma["linLimb"] = 0.5
  ma["quadLimb"] = 0.5
  ma["b"] = 0.

  # Choose some (large) time axis
  time = xnorm
  #time = np.arange(xnorm.min(),xnorm.max(),.01)

  # ... and calculate model
  y = ma.evaluate(time)
  #y += np.random.normal(0.,.01,time.size)
  #yerr=np.ones(time.size)*.01

  # Let's see what happened ...
  #mpl.errorbar(time, y, yerr=yerr, fmt='kp', markersize=2)
  #mpl.show()

  ma.thaw(["per","T0","linLimb","quadLimb","a","i","p"])
  ma.fit(time,ynorm,yerr=yerr, maxiter=1e4, maxfun=1e4, xtol=1e-8, ftol=1e-8)
  ma.parameterSummary()
  fit=ma.evaluate(time)

  fo = open("foo.txt", "wb")
  fo.write(ma["per"])
  exit()
  
  #mpl.errorbar(xnorm, ynorm, yerr=yerr, fmt='bp', ms=2)
  #mpl.plot(time, y, 'rp', ms=2)
  ##mpl.hist(ynorm-y)
  #mpl.show()