Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
"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')))
# star1['radius'] = orbit['c2pars'].get_value_with_unit('radius')
B{Warning!} Changing C{sma} is now impossible:
>>> ps['sma'] = 10.
>>> print ps['sma'],ps['incl'],ps['asini']
17.6776695297 45.0 12.5
"""
#clear_memoization(self)
#-- clean up the contraint given by the user
splitted = constraint.split('=')
qualifier,expression = splitted[0],'='.join(splitted[1:])
qualifier = qualifier.split('{')[1].split('}')[0].strip()
self.constraints[qualifier] = expression
if include_as_parameter:
self.add(Parameter(qualifier=qualifier,value=0.))
elif isinstance(include_as_parameter,str):
self.add(Parameter(qualifier=qualifier,value=0.,unit=include_as_parameter))
if do_run_constraints and '{' in expression: # don't run constraints if it's only a number!
self.run_constraints()
"""
# If is in None, allow building the feedback with generic parameters
if init is None:
with open(self._emcee_file, 'r') as ff:
while True:
line = ff.readline()
if not line:
raise IOError("Invalid file")
if line[0] == '#':
continue
n_pars = len(line.strip().split())-2
break
adjustables = []
for i in range(n_pars):
adjustables.append(parameters.Parameter(qualifier=str(i), cast_type=float, value=0))
adjustables[-1]._unique_label = str(i)
return adjustables
if isinstance(init, list):
return init
# is this a bundle?
if hasattr(init, 'twigs'):
system = init.get_system()
# else it is a system
else:
system = init
try:
key, val = l.split('=')
except:
logger.error("line " + l[:-1] + " could not be parsed ('{}' probably not Phoebe legacy file)".format(inputfile))
raise IOError("Cannot parse phoebe file '{}': line '{}'".format(inputfile, l[:-1]))
key = key.strip()
# Start with parameters that determine container sizes:
if key == 'phoebe_lcno':
all_lcobs = [datasets.LCDataSet(user_columns=['time','flux','sigma'], user_units=dict(time='JD',flux='W/m3',sigma='W/m3')) for i in range(int(val))]
all_lcdeps[0] = [parameters.ParameterSet('lcdep', ld_coeffs=[0.5,0.5]) for i in range(int(val))]
all_lcdeps[1] = [parameters.ParameterSet('lcdep', ld_coeffs=[0.5,0.5]) for i in range(int(val))]
# We need to add extinction parameters to all lcdeps
ext_par = parameters.Parameter('extinction', context='lcdep', cast_type=float)
for i in range(int(val)):
all_lcdeps[0][i].add(ext_par.copy())
all_lcdeps[1][i].add(ext_par.copy())
continue
if key == 'phoebe_rvno':
all_rvobs = [datasets.RVDataSet(user_columns=['time','rv','sigma'], user_units=dict(time='JD',rv='km/s',sigma='km/s')) for i in range(int(val))]
rv_component_spec = [None for i in range(int(val))]
all_rvdeps[0] = [parameters.ParameterSet('rvdep', ld_coeffs=[0.5,0.5]) for i in range(int(val))]
all_rvdeps[1] = [parameters.ParameterSet('rvdep', ld_coeffs=[0.5,0.5]) for i in range(int(val))]
continue
#-- if this is an rvdep or lcdep, check which index it has
# and remove it from the string:
pattern = re.search('\[(\d)\]', key)
separate = re.search('\:', val)
df = 0
output = freqanalyse.spectrum_2D(time, wavelength[0],
flux, f0=f0, fn=fn, df=df, threads=6,
method='scargle', subs_av=True,
full_output=True, harmonics=0)
const, e_const = output['pars']['const'], output['pars']['e_const']
ampl, e_ampl = output['pars']['ampl'], output['pars']['e_ampl']
freq, e_freq = output['pars']['freq'], output['pars']['e_freq']
phase, e_phase = output['pars']['phase'], output['pars']['e_phase']
self.add(parameters.Parameter(qualifier='const', value=const, cast_type=np.array))
self.add(parameters.Parameter(qualifier='sigma_const', value=e_const, cast_type=np.array))
self.add(parameters.Parameter(qualifier='ampl', value=ampl, cast_type=np.array))
self.add(parameters.Parameter(qualifier='sigma_ampl', value=e_ampl, cast_type=np.array))
self.add(parameters.Parameter(qualifier='freq', value=freq, cast_type=np.array))
self.add(parameters.Parameter(qualifier='sigma_freq', value=e_freq, cast_type=np.array))
self.add(parameters.Parameter(qualifier='phase', value=phase, cast_type=np.array))
self.add(parameters.Parameter(qualifier='sigma_phase', value=e_phase, cast_type=np.array))
self.add(parameters.Parameter(qualifier='avprof', value=output['avprof'], cast_type=np.array))
self.add(parameters.Parameter(qualifier='model', value=output['model'], cast_type=np.array))
K1 = q*sma/(1+q)*2*np.pi/P*np.sqrt(1-e**2)
if K2 is None:
period = orbit.request_value('period','SI')
ecc = orbit['ecc']
q = orbit['q']
sma = orbit.request_value('sma','SI')
K2 = sma/(1+q)*2*np.pi/P*np.sqrt(1-e**2)
if not 'K1' in orbit:
orbit.add(parameters.Parameter(qualifier='K1',value=K1,
**kwargs))
else:
orbit['K1'] = K1
if not 'K2' in orbit:
orbit.add(parameters.Parameter(qualifier='K2',value=K2,
**kwargs))
else:
orbit['K2'] = K2
orbit.add_constraint('{q} = {K1}/{K2}')
orbit.add_constraint('{sma} = {K2}*{period}/(2*np.pi)*np.sqrt(1-{ecc}**2)*(1+{q})')
logger.info("orbit '{}': q and sma constrained by K1 and K2".format(orbit['label']))
return orbit.get_parameter('K1'), orbit.get_parameter('K2')
kwargs.setdefault('llim',0.0)
kwargs.setdefault('ulim',1.0)
kwargs.setdefault('frame','phoebe')
kwargs.setdefault('cast_type',float)
kwargs.setdefault('repr','%f')
star.pop_constraint('vrotcrit', None)
if vrotcrit is None:
radius = star.get_value('radius','m')
mass = star.get_value('mass','kg')
rotperiod = star.get_value('rotperiod','s')
if not 'vrotcrit' in star:
star.add(parameters.Parameter(qualifier='vrotcrit',
value=vrotcrit, **kwargs))
else:
star['vrotcrit'] = vrotcrit
if derive == 'rotperiod':
star.pop_constraint('rotperiod', None)
# Rotperiod = omega/omega_crit * omega_crit, derived from v/v_crit
star.add_constraint('{rotperiod} = 2*np.pi / (np.cos(3*np.arccos(0.5*{vrotcrit}) - np.pi) * np.sqrt( (8*constants.GG*{mass})/(27*{radius}**3)))')
logger.info("star '{}': 'rotperiod' constrained by 'vrotcrit'".format(star['label']))
star.run_constraints()
return star.get_parameter('vrotcrit')
kwargs.setdefault('description','Polar effective temperature')
kwargs.setdefault('llim',0)
kwargs.setdefault('ulim',1e20)
kwargs.setdefault('unit','K')
kwargs.setdefault('frame','phoebe')
kwargs.setdefault('cast_type',float)
kwargs.setdefault('repr','%f')
star.pop_constraint('teffpolar',None)
#-- set default value
if teffpolar is None:
teffpolar = star['teff']
if not 'teffpolar' in star:
star.add(parameters.Parameter(qualifier='teffpolar',
value=teffpolar,**kwargs))
else:
star['teffpolar'] = teffpolar
star.add_constraint('{teff} = {teffpolar}')
logger.info("star '{}': 'teff' redefined to be equal to 'teffpolar'".format(star['label']))
return star.get_parameter('teffpolar')
See, for example, [Esteves2013]_ or [Lopez-Morales2007]_ for a more detailed
explanation.
"""
if kwargs and 'albgeom' in ps:
raise ValueError("You cannot give extra kwargs to add_albgeom if albgeom already exist")
kwargs.setdefault('description',"Geometric albedo")
kwargs.setdefault('context',ps.context)
kwargs.setdefault('frame','phoebe')
kwargs.setdefault('cast_type',float)
kwargs.setdefault('repr','%f')
#-- remove any constraints on albgeom and add the parameter
ps.pop_constraint('albgeom',None)
if not 'albgeom' in ps:
ps.add(parameters.Parameter(qualifier='albgeom', value=albgeom,**kwargs))
else:
ps['albgeom'] = albgeom
ps.add_constraint('{alb} = 1.5*{albgeom}')
def add(self, parameter, with_qualifier=None, force=False):
"""
Add a parameter to the class instance.
If the parameter qualifier is already in the container, it will be
silently overwritten only if :envvar:`force=False`.
C{parameter} can either be a Parameter instance, or a dictionary with
stuff like C{qualifier}, C{value} etc.
"""
# Maybe we gave a dictionary with the properties instead of a Parameter.
# in that case, convert the dict to a Parameter object first
if not isinstance(parameter, Parameter):
parameter = Parameter(**parameter)
# Now add it to the container
if parameter.qualifier in self.container and not force:
raise KeyError('{0} already exists'.format(parameter.qualifier))
if with_qualifier is None:
with_qualifier = parameter.qualifier
self.container[with_qualifier] = parameter
# force context
self.container[with_qualifier].set_context(self.context)
#self.__dict__[parameter.qualifier] = parameter.get_value()