How to use the phoebe.u function in phoebe

def test_sun(plot=False):
    b = phoebe.default_star(starA='sun')

    b.set_value('teff', 1.0*u.solTeff)
    b.set_value('rpole', 1.0*u.solRad)
    b.set_value('mass', 1.0*u.solMass)
    b.set_value('period', 24.47*u.d)

    b.set_value('incl', 23.5*u.deg)
    b.set_value('distance', 1.0*u.AU)

    assert(b.get_value('teff', u.K)==5772)
    assert(b.get_value('rpole', u.solRad)==1.0)
    assert(b.get_value('mass', u.solMass)==1.0)
    assert(b.get_value('period', u.d)==24.47)
    assert(b.get_value('incl', u.deg)==23.5)
    assert(b.get_value('distance', u.m)==1.0*u.AU.to(u.m))

    b.add_dataset('lc', pblum=1*u.solLum)

set to True.

    """

    b.run_delayed_constraints()

    kwargs.setdefault('check_visible', False)
    kwargs.setdefault('check_default', False)

    computeps = b.get_compute(compute, force_ps=True, **_skip_filter_checks)
    ltte = computeps.get_value('ltte', **kwargs)

    # make sure times is an array and not a list
    times = np.array(times)

    vgamma = b.get_value('vgamma', context='system', unit=u.solRad/u.d, **_skip_filter_checks)
    t0 = b.get_value('t0', context='system', unit=u.d, **_skip_filter_checks)

    hier = b.hierarchy
    starrefs = hier.get_stars()
    orbitrefs = hier.get_orbits()
    s = b.filter(context='component', **_skip_filter_checks)

    periods, eccs, smas, t0_perpasses, per0s, long_ans, incls, dpdts, \
    deccdts, dperdts, components = [],[],[],[],[],[],[],[],[],[],[]


    for component in starrefs:

        # we need to build a list of all orbitlabels underwhich this component
        # belongs.  For a simple binary this is just the parent, but for hierarchical
        # systems we need to get the labels of the outer-orbits as well

def _true_anom_to_phase(true_anom, period, ecc, per0):
    """
    TODO: add documentation
    """
    phshift = 0

    mean_anom = true_anom - (ecc*sin(true_anom))*u.deg

    Phi = (mean_anom + per0) / (360*u.deg) - 1./4

    # phase = Phi - (phshift - 0.25 + per0/(360*u.deg)) * period
    phase = (Phi*u.d - (phshift - 0.25 + per0/(360*u.deg)) * period)*(u.cycle/u.d)

    return phase

hier = b.hierarchy

    for i,comp in enumerate(components):
        c_times = np.array([])
        c_rvs = np.array([])
        c_sigmas = np.array([])

        for dataset in datasets:
            rvc_ps = b.get_dataset(dataset=dataset, component=comp, **_skip_filter_checks)
            rvc_rvs = rvc_ps.get_value(qualifier='rvs', unit=u.km/u.s, **_skip_filter_checks)
            if not len(rvc_rvs):
                # then no observations here
                continue

            rvc_times = rvc_ps.get_value(qualifier='times', unit=u.d, **_skip_filter_checks)
            if len(rvc_rvs) != len(rvc_times):
                raise ValueError("rv@{}@{} does not match length of times@{}@{}".format(comp, dataset, comp, dataset))

            rvc_sigmas = rvc_ps.get_value(qualifier='sigmas', unit=u.km/u.s, **_skip_filter_checks)
            if not len(rvc_sigmas):
                rvc_sigmas = np.full_like(rvc_rvs, fill_value=np.nan)

            mask_enabled = b.get_value(qualifier='mask_enabled', dataset=dataset, default=False, **_skip_filter_checks)
            if mask and mask_enabled:
                mask_phases = b.get_value(qualifier='mask_phases', dataset=dataset, **_skip_filter_checks)
                if len(mask_phases):
                    logger.warning("applying mask_phases - set mask_enabled=False to disable")
                    mask_t0 = b.get_value(qualifier='phases_t0', dataset=dataset, **_skip_filter_checks)
                    phases_for_mask = b.to_phase(rvc_times, component=None, t0=mask_t0)

                    inds = phase_mask_inds(phases_for_mask, mask_phases)

# now pull general compute options
        if compute is not None:
            if isinstance(compute, str):
                compute_ps = b.get_compute(compute, check_visible=False)
            else:
                # then hopefully compute is the parameterset
                compute_ps = compute
            eclipse_method = compute_ps.get_value(qualifier='eclipse_method', **kwargs)
            horizon_method = compute_ps.get_value(qualifier='horizon_method', check_visible=False, **kwargs)
            dynamics_method = compute_ps.get_value(qualifier='dynamics_method', **kwargs)
            irrad_method = compute_ps.get_value(qualifier='irrad_method', **kwargs)
            boosting_method = compute_ps.get_value(qualifier='boosting_method', **kwargs)
            if conf.devel:
                mesh_init_phi = compute_ps.get_value(qualifier='mesh_init_phi', unit=u.rad, **kwargs)
            else:
                mesh_init_phi = 0.0
        else:
            eclipse_method = 'native'
            horizon_method = 'boolean'
            dynamics_method = 'keplerian'
            irrad_method = 'none'
            boosting_method = 'none'
            mesh_init_phi = 0.0

        # NOTE: here we use globals()[Classname] because getattr doesn't work in
        # the current module - now this doesn't really make sense since we only
        # support stars, but eventually the classname could be Disk, Spot, etc
        if 'dynamics_method' in kwargs.keys():
            # already set as default above
            _dump = kwargs.pop('dynamics_method')

primary component.
    * `mean_anom` (float/quantity, optional): mean anomaly.
    * `incl` (float/quantity, optional): orbital inclination.
    * `q` (float, optional): mass ratio.
    * `sma` (float/quantity, optional): semi-major axis of the orbit.
    * `long_an` (float/quantity, optional): longitude of the ascending node.

    Returns
    --------
    * (, list): ParameterSet of all newly created
         objects and a list of all necessary
        constraints.
    """
    params = []

    params += [FloatParameter(qualifier='period', timederiv='dpdt', value=kwargs.get('period', 1.0), default_unit=u.d, limits=(1e-6,None), description='Orbital period')]
    params += [FloatParameter(qualifier='freq', value=kwargs.get('freq', 2*np.pi/3.0), default_unit=u.rad/u.d, advanced=True, description='Orbital frequency')]
    params += [FloatParameter(qualifier='dpdt', value=kwargs.get('dpdt', 0.0), default_unit=u.s/u.yr, advanced=True, description='Time derivative of orbital period')]
    params += [FloatParameter(qualifier='per0', timederiv='dperdt', value=kwargs.get('per0', 0.0), default_unit=u.deg, description='Argument of periastron')]
    params += [FloatParameter(qualifier='dperdt', value=kwargs.get('dperdt', 0.0), default_unit=u.deg/u.yr, advanced=True, description='Periastron change')]
    params += [FloatParameter(qualifier='ecc', timederiv='deccdt', value=kwargs.get('ecc', 0.0), default_unit=u.dimensionless_unscaled, limits=(0.0,0.999999), description='Eccentricity')]
    # if conf.devel:
        # NOTE: if adding this back in, will need to update the t0_* constraints in builtin.py and re-enable in parameters.HierarchyParameter.is_time_dependent
        # params += [FloatParameter(qualifier='deccdt', value=kwargs.get('deccdt', 0.0), default_unit=u.dimensionless_unscaled/u.d, advanced=True, description='Eccentricity change')]
    params += [FloatParameter(qualifier='t0_perpass', value=kwargs.get('t0_perpass', 0.0), default_unit=u.d, description='Zeropoint date at periastron passage of the primary component')]  # TODO: d vs JD
    params += [FloatParameter(qualifier='t0_supconj', value=kwargs.get('t0_supconj', 0.0), default_unit=u.d, description='Zeropoint date at superior conjunction of the primary component')]  # TODO: d vs JD
    params += [FloatParameter(qualifier='t0_ref', value=kwargs.get('t0_ref', 0.0), default_unit=u.d, description='Zeropoint date at reference point for the primary component')]
    params += [FloatParameter(qualifier='mean_anom', value=kwargs.get('mean_anom', 0.0), default_unit=u.deg, advanced=True, description='Mean anomaly at t0@system')]
    #params += [FloatParameter(qualifier='ph_perpass', value=kwargs.get('ph_perpass', 0.0), default_unit=u.cycle, description='Phase at periastron passage')]
    #params += [FloatParameter(qualifier='ph_supconj', value=kwargs.get('ph_supconj', 0.0), default_unit=u.cycle, description='Phase at superior conjunction')]
    #params += [FloatParameter(qualifier='ph_infconj', value=kwargs.get('ph_infconj', 0.0), default_unit=u.cycle, description='Phase at inferior conjunction')]
    #params += [FloatParameter(qualifier='t0_ph0', value=kwargs.get('t0_ph0', 0.0), default_unit=u.d, description='Zeropoint to anchor at phase=0.0')] # TODO: d vs JD

compute = kwargs.get('compute')
        starrefs = kwargs.get('starrefs')
        orbitrefs = kwargs.get('orbitrefs')
        step_size = kwargs.get('step_size')
        orbit_error = kwargs.get('orbit_error')
        time0 = kwargs.get('time0')

        # write the input file
        # TODO: need to use TemporaryFiles to be MPI safe
        fi = open('_tmp_pd_inp', 'w')
        fi.write('{} {}\n'.format(len(starrefs), time0))
        fi.write('{} {}\n'.format(step_size, orbit_error))
        fi.write('\n')
        fi.write(' '.join([str(b.get_value(qualifier='mass', component=star,
                context='component', unit=u.solMass) * c.G.to('AU3 / (Msun d2)').value)
                for star in starrefs])+'\n') # GM

        fi.write(' '.join([str(b.get_value(qualifier='requiv', component=star,
                context='component', unit=u.AU))
                for star in starrefs])+'\n')

        if info['kind'] == 'lc':
            # TODO: this will make two meshing calls, let's create and extract from the dictionary instead, or use set_value=True
            pblums = [b.get_value(qualifier='pblum', dataset=info['dataset'], component=starref, unit=u.W, check_visible=False) for starref in starrefs]

            u1s, u2s = [], []
            for star in starrefs:
                if b.get_value(qualifier='ld_func', component=star, dataset=info['dataset'], context='dataset') == 'quadratic':
                    ld_coeffs = b.get_value(qualifier='ld_coeffs', component=star, dataset=info['dataset'], context='dataset', check_visible=False)
                else:
                    # TODO: can we still interpolate for quadratic manually using b.compute_ld_coeffs?

(`times` is defined at the mid-exposure).
        Only applicable if `syn` is False and `is_lc` is True.

    Returns
    --------
    * ( or list, list): ParameterSet (if `as_ps`)
        or list of all newly created
         objects and a list of all necessary
        constraints.
    """

    params, constraints = [], []


    if is_lc:
        params += [FloatArrayParameter(qualifier='times', value=kwargs.get('times', []), default_unit=u.d, description='Model (synthetic) times' if syn else 'Observed times')]
        params += [FloatArrayParameter(qualifier='fluxes', value=_empty_array(kwargs, 'fluxes'), default_unit=u.W/u.m**2, description='Model (synthetic) flux' if syn else 'Observed flux')]

    if not syn:
        # TODO: should we move all limb-darkening to compute options since
        # not all backends support interp (and lookup is atm-dependent)
        params += [ChoiceParameter(qualifier='ld_mode', copy_for={'kind': ['star'], 'component': '*'}, component='_default',
                                   value=kwargs.get('ld_mode', 'interp'), choices=['interp', 'lookup', 'manual'],
                                   description='Mode to use for limb-darkening')]
        params += [ChoiceParameter(visible_if='ld_mode:lookup|manual', qualifier='ld_func',
                                   copy_for={'kind': ['star'], 'component': '*'}, component='_default',
                                   value=kwargs.get('ld_func', 'logarithmic'), choices=_ld_func_choices,
                                   description='Limb darkening model')]
        params += [ChoiceParameter(visible_if='ld_mode:lookup', qualifier='ld_coeffs_source',
                                   copy_for={'kind': ['star'], 'component': '*'}, component='_default',
                                   value=kwargs.get('ld_coeffs_source', 'auto'), choices=_ld_coeffs_source_choices,
                                   advanced=True,

spots_1=None, spots_2=None,
                                  flux_weighted=flux_weighted,
                                  verbose=1)


            # fill packets
            packetlist = []

            packetlist.append(_make_packet('times',
                                           info['times']*u.d,
                                           None,
                                           info))

            rvs = rvs1 if b.hierarchy.get_primary_or_secondary(info['component'])=='primary' else rvs2
            packetlist.append(_make_packet('rvs',
                                           rvs*u.km/u.s,
                                           None,
                                           info))
        else:
            raise TypeError("ellc only supports 'lc' and 'rv' datasets")

        return packetlist

# if 'cosbetas' in info['mesh_columns']:
                    # this_syn.set_value(time=time, qualifier='cosbetas', value=mesh['csbt{}'.format(cind)])

                # if 'areas' in info['mesh_columns']:
                    # TODO: compute and store areas from glump

                if 'loggs' in info['mesh_columns']:
                    packetlist.append(_make_packet('loggs',
                                                   mqtf(mesh['glog{}'.format(cind)]),
                                                   time,
                                                   info))


                if 'teffs' in info['mesh_columns']:
                    packetlist.append(_make_packet('teffs',
                                                   mqtf(mesh['tloc{}'.format(cind)])*u.K,
                                                   time,
                                                   info))

            # now we'll loop over the passband-datasets with requested columns
            for mesh_kind, mesh_dataset in zip(info['mesh_kinds'], info['mesh_datasets']):
                # print "*** legacy mesh pb", info['dataset'], mesh_dataset, mesh_kind
                if mesh_kind == 'lc':
                    lcind = lcinds[mesh_dataset]
                    for time in info['times']:
                        flux, mesh = phb1.lc((time,), lcind, True)

                        if 'abs_normal_intensities@{}'.format(mesh_dataset) in info['mesh_columns']:
                            packetlist.append(_make_packet('abs_normal_intensities',
                                                           mqtf(mesh['Inorm{}'.format(cind)])*u.erg*u.s**-1*u.cm**-3,
                                                           time,
                                                           info,