Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
b.add_solver('estimator.rv_periodogram')
b.add_solver('estimator.lc_geometry')
b.add_solver('estimator.rv_geometry')
b.add_solver('estimator.ebai')
b.add_solver('optimizer.nelder_mead')
b.add_solver('optimizer.differential_evolution')
b.add_solver('sampler.emcee')
b.add_solver('sampler.dynesty')
# TODO: include constraint_func? Shouldn't matter since they're not in twigs
should_be_forbidden = b.qualifiers + b.contexts + b.kinds + [c.split('@')[0] for c in b.get_parameter('columns').choices]
if verbose:
for l in should_be_forbidden:
if l not in phoebe.parameters.parameters._forbidden_labels:
print(l)
for l in should_be_forbidden:
assert(l in phoebe.parameters.parameters._forbidden_labels)
phoebe.reset_settings()
return b
ref='mylc')
# create bodies
starA = universe.BinaryRocheStar(Apars, mesh=meshpars, orbit=orbitpars,
pbdep=[pbdep_lc,pbdep_rv_a], obs=[obs_rv_a])
starB = universe.BinaryRocheStar(Bpars, mesh=meshpars, orbit=orbitpars,
pbdep=[pbdep_lc,pbdep_rv_b], obs=[obs_rv_b])
system = universe.BinaryBag([starA,starB], orbit=orbitpars, label='V380Cyg',
position=globs,obs=[obs_lc])
# compute
#obs_lc.set_adjust('offset',True)
#obs_lc.set_adjust('scale',True)
tools.group([obs_rv_a, obs_rv_b], 'rv', scale=False, offset=True)
compute_options = parameters.ParameterSet(context='compute')
system.compute(params=compute_options, mpi=None)
lcobs = system.get_obs(category='lc', ref='mylc').asarray()
lcsyn = system.get_synthetic(category='lc', ref='mylc').asarray()
rvobs1 = system[0].get_obs(category='rv', ref='myrv_comp1').asarray()
rvsyn1 = system[0].get_synthetic(category='rv', ref='myrv_comp1').asarray()
rvobs2 = system[1].get_obs(category='rv', ref='myrv_comp2').asarray()
rvsyn2 = system[1].get_synthetic(category='rv', ref='myrv_comp2').asarray()
if not debug:
assert(np.mean((lcobs['flux']-(lcsyn['flux']*lcobs['scale'] - lcobs['offset']))**2/lcobs['sigma']**2)<0.00017)
assert(np.mean((rvobs1['rv']-(rvsyn1['rv'] - rvobs1['offset']))**2/rvobs1['sigma']**2)<2.14)
assert(np.mean((rvobs2['rv']-(rvsyn2['rv'] - rvobs2['offset']))**2/rvobs2['sigma']**2)<2.41)
else:
def test_pblums_l3(debug=False):
"""
Pblum and l3 scaling: synthetic V380
"""
return None
parfile = os.path.join(exdir, 'V380_Cyg_circular.par')
lcfile = os.path.join(exdir, 'mylc.lc')
rv1file = os.path.join(exdir, 'myrv_comp1.rv')
rv2file = os.path.join(exdir, 'myrv_comp2.rv')
Apars,Bpars,orbitpars,globs = create.from_library(parfile, create_body=False)
meshpars = parameters.ParameterSet(context='mesh:marching', delta=0.2)
# load datasets
obs_rv_a, pbdep_rv_a = phoebe.parse_rv(rv1file,
passband='JOHNSON.V',
columns=['time','rv','sigma'],
components=[None,'myprimary','myprimary'],
ref='myrv_comp1')
obs_rv_b, pbdep_rv_b = phoebe.parse_rv(rv2file,
passband='JOHNSON.V',
columns=['time','rv','sigma'],
components=[None,'mysecondary','mysecondary'],
ref='myrv_comp2')
obs_lc, pbdep_lc = phoebe.parse_lc(lcfile,
passband='JOHNSON.V',
"Possible solution: decrease total mass or period.\n"
"If you want to keep the mass ratio, you might set "
"M1={} Msol, M2={} Msol").format(asini, sma, new_mass1+1e-10, new_mass2+1e-10))
logger.info('Calculated system asini from period, ecc and semi-amplitudes: asini={} Rsol'.format(asini))
tools.add_asini(orbit,asini,derive='incl',unit='Rsol')
logger.info("Calculated incl: {} deg".format(orbit.request_value('incl','deg')))
#-- when the vsini is fixed, we can only adjust the radius or the synchronicity!
for i,(comp,star,vsini) in enumerate(zip([comp1,comp2],[star1,star2],[vsini1,vsini2])):
if vsini is None:
continue
radius = star.get_value_with_unit('radius')
#-- adjust the radius if possible
if star.get_adjust('radius'):
comp.add(parameters.Parameter(qualifier='radius',value=radius[0],unit=radius[1],description='Approximate stellar radius',adjust=star.get_adjust('radius')))
comp.add(parameters.Parameter(qualifier='vsini',value=vsini,unit='km/s',description='Approximate projected equatorial velocity',adjust=(vsini is None)))
orbit.add_constraint('{{c{:d}pars.radius}} = {{period}}*{{c{:d}pars[vsini]}}/(2*np.pi*np.sin({{incl}}))'.format(i+1,i+1))
orbit.run_constraints()
logger.info('Star {:d}: radius will be derived from vsini={:.3g} {:s} (+ period and incl)'.format(i+1,*comp.get_value_with_unit('vsini')))
#-- or else adjust the synchronicity
else:
radius_km = conversions.convert(radius[1],'km',radius[0])
syncpar = period_sec*vsini/(2*np.pi*radius_km*np.sin(incl[0]))
comp['syncpar'] = syncpar
logger.info('Star {:d}: Synchronicity parameters is set to synpar={:.4f} match vsini={:.3g} km/s'.format(i+1,syncpar,vsini))
#if 'radius' in orbit['c1pars']:
# logger.info('Radius 1 = {} {}'.format(*orbit['c1pars'].get_value_with_unit('radius')))
# star1['radius'] = orbit['c1pars'].get_value_with_unit('radius')
#if 'radius' in orbit['c2pars']:
# logger.info('Radius 2 = {} {}'.format(*orbit['c2pars'].get_value_with_unit('radius')))
starA['rotperiod'] = np.inf
starA['ld_func'] = 'claret'
starA['atm'] = 'kurucz'
starA['ld_coeffs'] = 'kurucz'
starB = parameters.ParameterSet(frame='phoebe',context='star',add_constraints=True)
starB['teff'] = 3000,'K'
starB['mass'] = 0.2413,'Msol'
starB['shape'] = 'sphere'
starB['radius'] = 0.2543,'Rsol'
starB['rotperiod'] = np.inf
starB['ld_func'] = 'linear'
starB['atm'] = 'blackbody'
starB['ld_coeffs'] = [0.5]
starC = parameters.ParameterSet(frame='phoebe',context='star',add_constraints=True)
starC['teff'] = 2800,'K'
starC['mass'] = 0.2127,'Msol'
starC['shape'] = 'sphere'
starC['radius'] = 0.2318,'Rsol'
starC['rotperiod'] = np.inf
starB['ld_func'] = 'linear'
starB['atm'] = 'blackbody'
starB['ld_coeffs'] = [0.5]
#-- inner binary
orbitBC = parameters.ParameterSet(frame='phoebe',context='orbit',add_constraints=True)
orbitBC['period'] = 1.76713,'d'
orbitBC['sma'] = 0.021986,'au'
orbitBC['ecc'] = 0.02334
orbitBC['per0'] = 89.52,'deg'
orbitBC['incl'] = 96.907-92.100,'deg'
def walk(mybundle):
for val,path in utils.traverse_memory(mybundle,
list_types=(Bundle, universe.Body, list,tuple),
dict_types=(dict, ),
parset_types=(parameters.ParameterSet, ),
get_label=(universe.Body, ),
get_context=(parameters.ParameterSet, ),
skip=()):
# First one is always root
path[0] = str(mybundle.__class__.__name__)
# All is left is to return it
yield path, val
# We're on the right track and can advance to find the
# next label in the structure!
#~ if name_of_this_level == structure_info[index]:
if isinstance(name_of_this_level,str) and fnmatch(name_of_this_level, structure_info[index]):
index += 1
# Keep looking if we didn't exhaust all specifications yet.
# If we're at the end already, this will avoid finding a
# variable at all (which is what we want)
if index < len(structure_info):
continue
# Now did we find it? We also need to check if the found parameter
# hasn't already been found (e.g. when two identical parameterSets
# are added).
if isinstance(val, parameters.Parameter):
#~ if val.get_qualifier() == qualifier and not val.get_unique_label() in found_labels:
if fnmatch(val.get_qualifier(), qualifier) and not val.get_unique_label() in found_labels:
# Special handling of orbits: you can't request orbital
# information of a component, only of the BodyBag.
if 'orbit' in val.get_context() and structure_info:
if not 'orbit' in structure_info[-1] and not (structure_info[-1] == path[-2]['label']):
continue
# Ignore parameters that are synthetics
if val.get_context()[-3:] == 'syn':
continue
found[build_twig_from_path(val,path)] = val
found_labels.append(val.get_unique_label())
def build_twig_from_path(thing, path):
"""
Build the twig for a parameter in a system.
@param thing: some object from a Bundle (Parameter, ParameterSet, Body...)
@type thing: some type
@param path: all levels of the hieararchy of the thing (as a list of objects)
@type path: list of objects
"""
twig = None
# We have different behaviour for twigs representing Parameters and
# ParameterSets
if isinstance(thing, parameters.Parameter):
# A twig for a parameterSet is built as:
# @[<label>@][@]
# what's the component:
labels = [ipath.get_label() for ipath in path if hasattr(ipath, 'get_label')]
component = labels[-1] if len(labels) else None
qualifier = thing.get_qualifier()
# Get the data label or parameter context
if isinstance(path[-2],str):
context = path[-2] # perhaps add path[-3] as well to get lcdep if there are clashes
else:
context = ''
if 'label' in path[-2]:
context += path[-2]['label'] + '@'</label>
kwargs.setdefault('description','Surface gravity')
kwargs.setdefault('unit',unit)
kwargs.setdefault('context',star.context)
kwargs.setdefault('adjust',False)
kwargs.setdefault('frame','phoebe')
kwargs.setdefault('cast_type',float)
kwargs.setdefault('repr','%f')
#-- default value
if surfgrav is None:
surfgrav = np.log10(constants.GG*star.get_value('mass','kg')/star.get_value('radius','m')**2) + 2.0
#-- remove any constraints on surfgrav and add the parameter
star.pop_constraint('surfgrav',None)
if not 'surfgrav' in star:
star.add(parameters.Parameter(qualifier='surfgrav',value=surfgrav,
**kwargs))
else:
star['surfgrav'] = surfgrav
#-- specify the dependent parameter
if derive is None:
star.pop_constraint('radius',None)
star.add_constraint('{surfgrav} = constants.GG*{mass}/{radius}**2')
logger.info("star '{}': 'surfgrav' constrained by 'mass' and 'radius'".format(star['label']))
elif derive=='mass':
star.pop_constraint('mass',None)
star.add_constraint('{mass} = {surfgrav}/constants.GG*{radius}**2')
logger.info("star '{}': '{}' constrained by 'surfgrav' and 'radius'".format(star['label'],derive))
elif derive=='radius':
star.pop_constraint('radius',None)
star.add_constraint('{radius} = np.sqrt((constants.GG*{mass})/{surfgrav})')
self.sections['compute'] = []
if not self.load_cfg('compute', basedir):
self.add_compute(label='preview',refl=False,heating=False,eclipse_alg='binary',subdiv_num=1)
self.add_compute(label='detailed',ltt=True,eclipse_alg='binary',boosting_alg='local')
#~ if devel:
if True:
self.sections['fitting'] = []
if not self.load_cfg('fitting', basedir):
self.add_fitting(context='fitting:lmfit', label='lmfit', computelabel='preview')
self.add_fitting(context='fitting:emcee', label='emcee', computelabel='preview')
if devel:
self.sections['gui'] = [] # currently assumes only 1 entry
if not self.load_cfg('gui', basedir):
self._add_to_section('gui', parameters.ParameterSet(context='gui'))
self.sections['logger'] = [] # currently assumes only 1 entry
if not self.load_cfg('logger', basedir):
self._add_to_section('logger', parameters.ParameterSet(context='logger'))
# now apply the logger settings
# if logger settings are changed, either reload the usersettings
# or call restart_logger()
self.restart_logger()
self._build_trunk()