How to use the phoebe.units.conversions function in phoebe

twigs_atm = mybundle.search('atm')
    for atm in twigs_atm:
        if atm.split('@')[0] == 'value':
            continue
        if mybundle[atm] in ['blackbody','kurucz']:
            mybundle[atm] = 'blackbody_uniform_none_teff.fits'
            passband_twig = 'passband@{}'.format("@".join(atm.split('@')[1:]))
            if passband_twig in mybundle and mybundle[passband_twig] == 'JOHNSON.V':
                mybundle[passband_twig] = 'johnson_v.ptf'
    
    mybundle['pblum@secondary'] = phb.getpar('phoebe_plum2')
    mybundle.run_compute(label='from_legacy', irradiation_alg='point_source')
    #mybundle.get_system().compute(animate=True)
    lc_ph2 = mybundle['flux@{}@lcsyn'.format(dataref)]
    
    U = phb2.units.conversions.Unit
    R1 = U(mybundle.get_system()[0].params['component'].request_value('r_pole'), 'm')
    T1 = U(mybundle['teff@primary'], 'K')
    R2 = U(mybundle.get_system()[1].params['component'].request_value('r_pole'), 'm')
    T2 = U(mybundle['teff@secondary'], 'K')
    sigma = U('sigma')
    L1 = (4*np.pi*R1**2*sigma*T1**4).convert('Lsol')
    L2 = (4*np.pi*R2**2*sigma*T2**4).convert('Lsol')
    print("Numerical bolometric luminosity (primary) = {} Lsol".format(phb2.convert('erg/s', 'Lsol',mybundle['primary'].luminosity())))
    print("Numerical bolometric luminosity (secondary) = {} Lsol".format(phb2.convert('erg/s', 'Lsol',mybundle['secondary'].luminosity())))
    print("Eq. sphere bolometric luminosity (primary) = {}".format(L1))
    print("Eq. sphere bolometric luminosity (secondary) = {}".format(L2))
    print("Numerical passband luminosity (primary) = {} Lsol".format(phb2.convert('erg/s', 'Lsol',mybundle['primary'].luminosity(ref='LC'))))
    print("Numerical passband luminosity (secondary) = {} Lsol".format(phb2.convert('erg/s', 'Lsol',mybundle['secondary'].luminosity(ref='LC'))))
    
    
    # Passband luminosities:

# First load data
            filename = self.get_value('filename')
            columns = self.get_value('columns')
            #-- check if the data are already in here, and only reload when
            #   they are not, or when force is True
            if not force and (self['columns'][0] in self and len(self[self['columns'][0]])>0):
                return False
            data_columns = np.loadtxt(filename).T
            for i,col in enumerate(data_columns):
                if not columns[i] in self:
                    self.add(dict(qualifier=columns[i],value=col,description=''.format(filename)))
                else:
                    if 'user_units' in self and columns[i] in self['user_units']:
                        to_unit = self.get_parameter(columns[i]).get_unit()
                        from_unit = self['user_units'][columns[i]]
                        self[columns[i]] = conversions.convert(from_unit, to_unit, col)
                    else:
                        self[columns[i]] = col
           
            # Then try to load header:
            with open(filename, 'r') as ff:
                while True:
                    line = ff.readline()
                    
                    # End of file
                    if not line:
                        break
                   
                    # End of header
                    elif not line[0] == '#':
                        break

def set_unit(self, unit, convert=True):
        """
        Change the unit of a parameter.
        
        The values, lower and upper limits, step sizes and priors are also
        changed accordingly.
        
        @parameter unit: a physical unit
        @type unit: str, interpretable by L{conversions.convert}
        """
        #-- are we lazy? Possibly, somebody just set 'SI' or some other convention
        #   as a unit... in this case we have to derive the convention's version
        #   of the current unit... bastard!:
        if unit in conversions._conventions:
            unit = conversions.change_convention(unit,self.unit)
        
        if self.unit != unit:
            logger.debug("Converting parameter {} from {} to {}".format(self.qualifier,self.unit,unit))
            #-- the prior
            if convert:
                if hasattr(self,'prior'):
                    self.prior.convert(self.unit,unit)
                #-- the value
                new_value = conversions.convert(self.unit,unit,self.get_value())
                if hasattr(self,'llim'):
                    new_llim = conversions.convert(self.unit,unit,self.cast_type(self.llim))
                    self.llim = new_llim
                if hasattr(self,'ulim'):
                    new_ulim = conversions.convert(self.unit,unit,self.cast_type(self.ulim))
                    self.ulim = new_ulim

#-- 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')
    
    R1 = star1.get_value('radius','au')
    R2 = star2.get_value('radius','au')
    sma = orbit.get_value('sma','au'),'au'
    #-- we can also derive something on the radii from the eclipse duration, but
    #   only in the case for cicular orbits

old_dx = np.diff(old_x)
        x_shifted = np.hstack([old_x[:-1] - np.diff(old_x)/2.,old_x[-2:] + np.diff(old_x)[-1]/2.])
        #-- we definitely cannot go over the limits!
        x_shifted[0] = old_x[0]
        x_shifted[-1] = old_x[-1]
        old_dx = np.diff(x_shifted)
        
        if from_=='bounded':
            #old_x[old_xargs[1]] = args[1]
            #x_shifted[x_shiftedargs[1]] = args[1]
            new_x = transform_to_unbounded(old_x,*args)
            new_x_shifted = transform_to_unbounded(x_shifted,*args)
        else:
            new_x = conversions.convert(from_unit,to_unit,old_x)
            new_x_shifted = conversions.convert(from_unit,to_unit,x_shifted)
        new_dx = np.diff(new_x_shifted)
        new_y = old_dx/new_dx*old_y
        norm = np.trapz(new_y,x=new_x)
        new_y = new_y/norm
        self.distribution = 'histogram'
        self.distr_pars = dict(discrete=False,bins=new_x_shifted,prob=new_y)

values[name] = par.get_value()
        
        #-- now evaluate all the constraints, but convert the final values back
        #   to the original unit. We cannot use the "ps['qualifier'] = bla"
        #   method because we call L{run_constraints} in L{_setitem_}, causing
        #   infinite recursion.
        #-- also, the qualifier from the left hand side of the constraint doesn't
        #   need to be defined in the ParameterSet, so if it doesn't exist,
        #   just skip it
        if '.' in qualifier:
            return None
        
        value = eval(self.constraints[qualifier].format(**values))
        
        if unit:
            value = conversions.convert('SI',unit,value)
        
        return value

out = np.loadtxt(found, unpack=True)
            for col_in, col_file in zip(all_lcobs[i]['columns'], out):
                all_lcobs[i][col_in] = col_file
            
            if not all_lcobs[i]['sigma'].shape:
                all_lcobs[i]['sigma'] = len(all_lcobs[i])*[float(all_lcobs[i]['sigma'])]
            
            # Convert from mag to flux (perhaps with uncertainties)
            if all_lcobs[i]['user_columns'][1] == 'mag':
                passband = all_lcdeps[0][i]['passband']
                flux = all_lcobs[i]['flux']
                sigma = all_lcobs[i]['sigma']
                if len(sigma):
                    flux, sigma = conversions.convert('mag', 'W/m3', flux, sigma, passband=passband)
                else:
                    flux = conversions.convert('mag', 'W/m3', flux, passband=passband)
                all_lcobs[i]['flux'] = flux
                all_lcobs[i]['sigma'] = sigma
                
        all_lcobs[i]['filename'] = ""#all_lcobs[i]['ref'] + '_phoebe2.lc'
        
        
        if computehla:
            all_lcdeps[0][i]['pblum'] = -1
            all_lcobs[i].set_adjust('scale', True)
        
        if not usecla:
            all_lcdeps[1][i]['pblum'] = -1
        star1.add_pbdeps(all_lcdeps[0][i])
        star2.add_pbdeps(all_lcdeps[1][i])
    
    bodybag = universe.BodyBag([star1, star2], position=position, reddening=reddening)

# Overwrite phase default by setting of x_unit
    if x_unit is not None:
        x_unit_type = conversions.get_type(x_unit)
        if x_unit_type == 'angle':
            default_phased = True
        elif x_unit_type == 'time':
            default_phased = False
        else:
            raise ValueError(("Unallowed x_unit for plotting lc: {} is of type "
                              "{}, while only phase or time are "
                              "allowed").format(x_unit, x_unit_type))
    
    # Check y_unit type:
    if y_unit is not None:
        y_unit_type = conversions.get_type(y_unit)
        allowed = ['velocity']
        if not y_unit_type in allowed:
            raise ValueError(("Unallowed y_unit for plotting lc: {} is of type "
                              "{}, while only {} are allowed").format(y_unit,
                                   y_unit_type, ", ".join(allowed)))
    
    
    phased = kwargs.pop('phased', default_phased)
    ax = kwargs.pop('ax',plt.gca())
    scale = kwargs.pop('scale', 'obs')
        
    # Load synthetics: they need to be here
    loaded = syn.load(force=False)
    
    # Try to get the observations. They don't need to be loaded, we just need
    # the pblum and l3 values.

# the gravitional redshift is zero
    M = 0
    R = 1.
    if hasattr(the_system, 'bodies'):
        bodies = the_system.get_bodies()
    else:
        bodies = [the_system]
    
    rv_gravs = []
    for body in bodies:
        M = body.get_mass() * constants.Msol
        R = np.sqrt(body.mesh['_o_center'][:,0]**2 + \
                body.mesh['_o_center'][:,1]**2 + body.mesh['_o_center'][:,2]**2)
        R = R*constants.Rsol
        rv_grav = constants.GG*M/R/constants.cc
        rv_grav = conversions.convert('m/s', 'km/s', rv_grav)
        rv_gravs.append(rv_grav*np.ones(len(body.mesh)))
    rv_gravs = np.hstack(rv_gravs)
        
    #-- yes, it's that easy
    return rv_gravs

# the method
    ld_model = pbdep['ld_func']
    method = pbdep['method']
    profile = pbdep['profile']
    extend = pbdep.get('extend', 500.)
    keep = the_system.mesh['mu'] >= 0
    
    # Set the central wavelength "wc".
    wc = (wavelengths[0]+wavelengths[-1])/2.
    
    # Broaden the range of the wavelengths a bit, so that we are sure that also
    # neighbouring lines are taken into account. We clip the spectrum
    # aftwards. This "a bit" is taken to be 500 km/s, just to be sure.
    if extend > 0:
        w0 = conversions.convert('km/s', 'AA', -extend, wave=(wavelengths[0], 'AA'))
        wn = conversions.convert('km/s', 'AA', +extend, wave=(wavelengths[-1], 'AA'))
        step_before = (wavelengths[1]  - wavelengths[0])
        step_after =  (wavelengths[-1] - wavelengths[-2])
        wave_before = np.arange(w0, wavelengths[0], step_before)[:-1]
        wave_after = np.arange(wavelengths[-1], wn, step_after)[1:]
        wavelengths = np.hstack([wave_before, wavelengths, wave_after])    
    
    # If we're not seeing the star, we can easily compute the spectrum: it's
    # zero! Hihihi (hysterical laughter)!
    if not np.sum(keep):
        logger.info('Still need to compute (projected) intensity')
        the_system.intensity(ref=ref)
        
    # Check if there is any flux    
    the_system.projected_intensity(ref=ref, method='numerical')
    keep = the_system.mesh['proj_'+ref] > 0