doch das eingerückte gehört zum ende der schleife. ich wollte jetzt nur nicht die komplette schleife hier rein posten weil die über 200 zeilen geht, ich dachte das wird vielleicht zu unübersichtlich, aber ich kann auch nochmal die komplette schleife hier rein posten.
hier sind meine exceptions in den zeilen 248 - 251 zu finden in externen modulen und darauffolgend ein continue.
Code: Alles auswählen
kepid=[11446443, 8554498, 10666592, 3861595, 10748390, 8191762]
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", "a")
KeplerID = str(nobj)
#print "KeplerID", KeplerID
#print "ma.parameterSummary(per)", ma["per"]
Period = str(ma["per"])
#print "Period", Period
Inklination = str(ma["i"])
Halbachse = str(ma["a"])
TransitT0 = str(ma["T0"])
RpRs = str(ma["p"])
lines_of_text=("KeplerID = ", KeplerID, ", Period = ", Period, ", Inklination = ", Inklination, ", Große Halbachse in Sternradien = ", Halbachse, ", TransitT0 = ", TransitT0, ", Rp/Rs = ", RpRs, "\n")
#fo.writelines(lines_of_text)
##exit()
f = open("foo.py", "a")
f.writelines(lines_of_text)
#mpl.errorbar(xnorm, ynorm, yerr=yerr, fmt='bp', ms=2)
#mpl.plot(time, y, 'rp', ms=2)
##mpl.hist(ynorm-y)
#mpl.show()
import pyaErrTemplate
import pyaOtherErrors
import warn
import pyaValErrs
continue
f.close()