mrtommyb / ktransit

A simple exoplanet transit modeling tool in python
GNU General Public License v3.0
52 stars 20 forks source link

fitting not reporting rprs correctly? #11

Open jkrick opened 6 years ago

jkrick commented 6 years ago

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