fitting not reporting rprs correctly?
Opened this issue · 0 comments
I am trying to write a simple code to inject transits in fake data and then recover the planet parameters. However, when I use FitTransit I think it is reporting the rprs incorrectly. Specifically the input rprs and the fitted rprs don't match. I have based my code mainly on the example, and have tried to make as simple as possible, only fitting rprs while holding the other parameters fixed. I run through 5 different depth transits, and only the first in the loop appears correct. Output plots show data in black, input model in orange, and fitted model in blue. For me the list of input rprs are (0.1, 0.11, 0.12, 0.13, 0.14) and the fitted rprs are (0.1, 0.14, 0.19, 0.23, 0.26). What am I missing here? Thanks,
import ktransit
import matplotlib.pyplot as plt
import numpy as np
import random
from ktransit import FitTransit
#make the model planet transit
exptime = 0.0188
time = np.arange(0.25,1.75,exptime)
M = ktransit.LCModel()
#input variables
star_rho = 1.5
star_zpt = 0.0
planet_T0 = 1.0
planet_period = 1.0
planet_rprs = np.arange(0.1,.15,0.01)
#do this for a bunch of rprs
for counter in range(0,5): #n_aor
M.add_star(rho=star_rho, # mean stellar density in cgs units
zpt=star_zpt) # a photometric zeropoint, incase the normalisation was wonky
M.add_planet( T0=planet_T0, # a transit mid-time
period=planet_period, # an orbital period in days
rprs=planet_rprs[counter] ) # planet stellar radius ratio
M.add_data( time=time, # timestamps to evaluate the model on
itime=np.zeros_like(time)+exptime ) # integration time of each timestamp
tmod = M.transitmodel
plt.plot(time,tmod, color = 'orange', label = str(planet_rprs[counter]))
#make some noisey data
pmap_data = np.random.uniform(-0.005, 0.005, len(time))
#add the transit to the data
fake = tmod + pmap_data
plt.scatter(time, fake)
#and now: can I recover the rprs parameter?
fitT = FitTransit()
fitT.add_guess_star(rho=star_rho, zpt = star_zpt)
fitT.add_guess_planet(
period=planet_period,
T0=planet_T0,
rprs=planet_rprs[counter])
fitT.add_data(time=time, flux=fake)#, ferr=ferr)
vary_star = [] # free stellar parameters 'rho','zpt'
vary_planet = ([ 'rprs' # free planetary parameters#'T0', 'rprs'
]) # free planet parameters are the same for every planet you model
fitT.free_parameters(vary_star, vary_planet)
fitT.do_fit() # run the fitting
fitT.print_results() # print some results
#overplot fitted transit model
plt.plot(time, fitT.transitmodel, label = str(fitT.fitresultplanets['pnum0']['rprs']) , color = 'blue')
plt.legend()
plt.show()