ich möchte, dass mein Python-Programm die berrechneten Werte über die Funktion "write to file" in eine Extradatei schreibt.
Ich benutze dazu ein Modul dass Messwerte fitten kann (Zeilen 205-235). Diese gefitteten Werte möchte ich in diese Extradatei schreiben (Zeilen 237-239). Es wird mir zwar kein Fehler ausgeworfen, jedoch bekomme ich als fertige Datei dann nur Buchstabensalat, obwohl ich natürlich Zahlenwerte erhalten möchte.
Dazu mein Quellcode:
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()
Für jeden Ratschlag bin ich sehr dankbar.
Viele Grüße
vdv84