ich stehe vor folgendem Problem:
Ich habe den folgenden Python File erstellt, womit ich die Parameter eines Planeten herausfinden kann, nun möchte aus dieser Abfolge von Befehlen eine Funktion definieren, die ich in anderen Files als Modul importieren kann. Und zwar möchte ich diese Funktion für eine Vielzahl von Planeten durchführen lassen.
Den ausgewählten Planeten rufe ich mit nobj=x auf (Zeile 10):
Vielen Dank für jede Art von Hilfestellung, viele Grüße
vdv84
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=4
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
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)
mpl.errorbar(xnorm, ynorm, yerr=yerr, fmt='bp', ms=2)
mpl.plot(time, y, 'rp', ms=2)
#mpl.hist(ynorm-y)
mpl.show()
