Advertisement
Guest User

Untitled

a guest
Jul 18th, 2015
164
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 2.54 KB | None | 0 0
  1. #SimpleBinary: fitting
  2. import phoebe
  3. import numpy as np
  4. import matplotlib.pyplot as plt
  5.  
  6. logger = phoebe.get_basic_logger()
  7. star1 = phoebe.ParameterSet('star', mass=1.2, radius=1.1, teff=6321.,
  8.                             atm='kurucz', ld_coeffs='kurucz', ld_func='claret')
  9. star2 = phoebe.ParameterSet('star', mass=0.8, radius=0.9, teff=4123.,
  10.                             atm='kurucz', ld_coeffs='kurucz', ld_func='claret')
  11. comp1, comp2, orbit = phoebe.create.binary_from_stars(star1, star2, period=(10.,'d'))
  12. comp1['syncpar'] = 1.0
  13. comp2['syncpar'] = 1.0
  14. orbit['ecc'] = 0.1
  15. mesh = phoebe.ParameterSet('mesh:marching', delta=0.07)
  16. lcdep = phoebe.ParameterSet('lcdep', passband='KEPLER.V', atm='kurucz',
  17.                             ld_func='claret', ld_coeffs='kurucz', ref='my kepler lc')
  18. rvdep = phoebe.ParameterSet('rvdep', passband='JOHNSON.V', atm='kurucz',
  19.                             ld_func='claret', ld_coeffs='kurucz')
  20. #loading observations
  21. input_lc, lcdep_=phoebe.parse_lc('lc.obs')
  22. input_lc['ref']='my kepler lc'
  23. input_lc.get_parameter('pblum').set_adjust(True)
  24. input_lc.get_parameter('l3').set_adjust(True)
  25. phoebe.parameters.tools.add_esinw_ecosw(orbit)
  26. orbit['incl'] = 88.
  27. incl = orbit.get_parameter('incl')
  28. incl.set_adjust(True)
  29. incl.set_prior(distribution='uniform', lower=0.0, upper=180.)
  30. orbit['ecosw'] = 0.05
  31. ecosw = orbit.get_parameter('ecosw')
  32. ecosw.set_adjust(True)
  33. ecosw.set_prior(distribution='uniform', lower=-1, upper=1)
  34. orbit['esinw'] = 0.05
  35. esinw = orbit.get_parameter('esinw')
  36. esinw.set_adjust(True)
  37. esinw.set_prior(distribution='uniform', lower=-1, upper=1)
  38. fitparams = phoebe.ParameterSet('fitting:lmfit', method='nelder')
  39. star1 = phoebe.BinaryRocheStar(comp1, mesh=mesh, pbdep=[lcdep,rvdep])#, orbit=orbit)
  40. star2 = phoebe.BinaryRocheStar(comp2, mesh=mesh, pbdep=[lcdep,rvdep])#, orbit=orbit)
  41. system = phoebe.BinaryBag([star1, star2], orbit=orbit, obs=[input_lc])
  42. params = phoebe.ParameterSet('compute', eclipse_alg='binary')
  43. feedback = phoebe.fitting.run(system, params, fitparams=fitparams, mpi=None, accept=True)
  44. system.save('mymodel.phoebe')
  45. feedback.save('myfeedback.phoebe')
  46. system = phoebe.load_body('mymodel.phoebe')
  47. feedback = phoebe.load_ps('myfeedback.phoebe')
  48. syn = system.get_synthetic(category='lc', ref='my kepler lc').asarray()
  49. plt.figure()
  50. plt.plot(syn['time'], syn['flux'], 'ko-')
  51. syn1 = system[0].get_synthetic(category='rv').asarray()
  52. syn2 = system[1].get_synthetic(category='rv').asarray()
  53. plt.figure()
  54. plt.plot(syn1['time'], syn1['rv'], 'ko-')
  55. plt.plot(syn2['time'], syn2['rv'], 'ro-')
  56. plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement