How to use the gpkit.units function in gpkit

def test():
    """
    Test each different engine
    """
    #dict of initial guesses
    x0 = {
        'W_{engine}': 1e4*units('N'),
        'P_{t_0}': 1e1*units('kPa'),
        'T_{t_0}': 1e3*units('K'),
        'h_{t_0}': 1e6*units('J/kg'),
        'P_{t_{1.8}}': 1e1*units('kPa'),
        'T_{t_{1.8}}': 1e3*units('K'),
        'h_{t_{1.8}}': 1e6*units('J/kg'),
        'P_{t_2}': 1e1*units('kPa'),
        'T_{t_2}': 1e3*units('K'),
        'h_{t_2}': 1e6*units('J/kg'),
        'P_{t_2.1}': 1e3*units('K'),
        'T_{t_2.1}': 1e3*units('K'),
        'h_{t_2.1}': 1e6*units('J/kg'),
        'P_{t_{2.5}}': 1e3*units('kPa'),
        'T_{t_{2.5}}': 1e3*units('K'),
        'h_{t_{2.5}}': 1e6*units('J/kg'),
        'P_{t_3}': 1e4*units('kPa'),

def test():
    m = Mission(SimPleAC(),4)
    m.substitutions.update({
        'h_{cruise_m}'   :5000*units('m'),
        'Range_m'        :3000*units('km'),
        'W_{p_m}'        :3000*units('N'),
        '\\rho_{p_m}'    :1500*units('kg/m^3'),
        'C_m'            :120*units('1/hr'),
        'V_{min_m}'      :35*units('m/s'),
        'T/O factor_m'   :2,
    })
    m.cost = m['W_{f_m}']*units('1/N') + m['C_m']*m['t_m']
    sol = m.localsolve(verbosity=0)

span_max = Variable('b_{max}', 'm', 'Max Wing Span')

        K = Variable('K', '-', 'K for Parametric Drag Model')
        e = Variable('e', '-', 'Oswald Span Efficiency Factor')

        xw     = Variable('x_w', 'm', 'Position of wing aerodynamic center')
        cwma    = Variable('\\bar{c}_w', 'm',
                          'Mean aerodynamic chord (wing)')
        Vne     = Variable('V_{ne}', 'm/s', 'Never exceed velocity')
        
        constraints = []

        constraints.extend([
            #wing weight constraint
            #based off of a raymer weight and 737 data from TASOPT output file
            (S/(124.58*units('m^2')))**.65 == W_wing/(105384.1524*units('N')),

            #compute wing span and aspect ratio, subject to a span constraint
            AR == (span**2)/S,
            #AR == 9,

            span <= span_max,

            #compute K for the aircraft
            K == (pi * e * AR)**-1,

            xw == xw,
            
            Vne == Vne,
            cwma == cwma,
            ])

z_CG_0  = Variable('z_{CG}', 'm',
                           'CG height relative to bottom of fuselage')
        zwing   = Variable('z_{wing}', 'm',
                           'Height of wing relative to base of fuselage')


        with SignomialsEnabled():

            objective = W_lg

            constraints = [
                           # Track and Base geometry definitions
                           TCS([l_n+zwing+y_m*tan_gam>=l_m], reltol=1E-3), #[SP]
                           T == 2*y_m,
                           TCS([x_n + B <= x_m]),
                           x_n >= 5*units.m, # nose gear after nose

                           # Geometric constraints relating gear placement with
                           # fore/aft CG locations
                           TCS([dxn + x_n >= xcg], reltol=1E-3), # [SP]
                           TCS([dxm + xcg >= x_m], reltol=1E-4), # [SP]
                           # TODO forward and aft CG

                           # Maximum static loads through main and nose gears
                           L_n == W*dxm/B,
                           L_m == W*dxn/B,

                           # Dynamic braking load through nose gear
                           # (assumes deceleration of 10 ft/s^2)
                           L_n_dyn >= 0.31*((z_CG_0+l_m)/B)*W,

                           # For steering don't want too much or too little

'T_{t_6': 1e3*units('K'),
        'h_{t_6}': 1e6*units('J/kg'),
        'h_6': 1e6*units('J/kg'),
        'F_8': 1e2 * units('kN'),
        'F_6': 1e2 * units('kN'),
        'F': 1e2 * units('kN'),
        'F_{sp}': 1e-1,
        'TSFC': 1e-1,
        'I_{sp}': 1e4*units('s'),
        'u_6': 1e3*units('m/s'),
        'u_8': 1e3*units('m/s'),
        'm_{core}': 1e2*units('kg/s'),
        'm_{fan}': 1e3*units('kg/s'),
        '\\alpha': 1e1,
        'alphap1': 1e1,
        'm_{total}': 1e3*units('kg/s'),
        'T_2': 1e3*units('K'),
        'P_2': 1e2*units('kPa'),
        'u_2': 1e3*units('m/s'),
        'h_{2}': 1e6*units('J/kg'),
        'T_{2.5}': 1e3*units('K'),
        'P_{2.5}': 1e2*units('kPa'),
        'u_{2.5}': 1e3*units('m/s'),
        'h_{2.5}': 1e6*units('J/kg'),
        'P_{7}': 1e2*units('kPa'),
        'T_{7}': 1e3*units('K'),
        'u_7': 1e3*units('m/s'),
        'P_{5}': 1e2*units('kPa'),
        'T_{5}': 1e3*units('K'),
        'u_5': 1e3*units('m/s'),
        'P_{atm}': 1e2*units('kPa'),
        'T_{atm}': 1e3*units('K'),

'\\tau_{floor}': 30000. / 0.000145, # [TAS] [Al]
                'W\'\'_{floor}': 60.,  # [TAS]
                'W\'\'_{insul}': 22.,  # [TAS]
                'W\'_{window}': 145.*3.*units('N/m'),  # [TAS]
                'V_{mn}': 133.76*units('m/s'), 'V_{ne}':143.92*units('m/s'),

                'D_{reduct}': 1,

                # Fuselage subs
                'f_{seat}': 0.1,
                'W\'_{seat}': 1.,  # Seat weight determined by weight fraction instead
                'W_{cargo}': 0.1*units('N'), # Cargo weight determined by W_{avg. pass_{total}}
                'W_{avg. pass_{total}}':215.*units('lbf'),
                'f_{string}': 0.35,

                'h_{floor}': 5.12*units('in'),
                'w_{db}': 0.93*units('m'),
                '\\Delta P_{over}': 8.382 * units('psi'),
                'SPR': 8.,

                # TASOPT Fuselage substitutions
                'l_{nose}': 29.*units('ft'),
                'f_{L_{total/wing}}': 1.195,

                # Power system and landing gear subs
                'f_{hpesys}': 0.01, # [TAS]

                # Fractional weights
                'f_{fadd}': 0.2,  # [TAS]
                'f_{frame}': 0.25,  # [Philippe]
                'f_{lugg,1}': 0.4,  # [Philippe]
                'f_{lugg,2}': 0.1,  # [Philippe]

'CG height relative to bottom of fuselage')
        zwing   = Variable('z_{wing}', 'm',
                           'Height of wing relative to base of fuselage')

        d_nac   = Variable('d_{nacelle}', 'm', 'Nacelle diameter')

        with SignomialsEnabled():

            objective = W_lg

            constraints = [
                           # Track and Base geometry definitions
                           TCS([l_n+zwing+y_m*tan_gam>=l_m], reltol=1E-3), #[SP]
                           T == 2*y_m,
                           TCS([x_n + B <= x_m]),
                           x_n >= 5*units.m, # nose gear after nose

                           # Longitudinal tip over (static)
                           tan_phi == tan_15,

                           # Lateral tip over in turn (dynamic)
                           # www.dept.aoe.vt.edu/~mason/Mason_f/M96SC03.pdf
                           # stricter constraint uses forward CG
                           # cos(arctan(y/x))) = x/sqrt(x^2 + y^2)
                           1 >= (z_CG_0 + l_m)**2 * (y_m**2 + B**2) /
                                (dxn * y_m * tan_psi)**2,
                           tan_psi <= tan_63,

                           # Tail strike: Longitudinal ground clearance in
                           # takeoff, landing (Raymer says 10-15 degrees)
                           # TODO?: 2 cases:(i) upsweep angle > rotation angle,
                           # (ii) upsweep angle < rotation ang

W_end = Variable('W_{end}', 'N', 'Segment End Weight')
        W_burn = Variable('W_{burn}', 'N', 'Segment Fuel Burn Weight')
        WLoadmax = Variable('W_{Load_{max}}', 'N/m^2', 'Max Wing Loading')
        WLoad = Variable('W_{Load}', 'N/m^2', 'Wing Loading')
        t = Variable('tmin', 'min', 'Segment Flight Time in Minutes')
        thours = Variable('thr', 'hour', 'Segment Flight Time in Hours')

        constraints = []

        constraints.extend([
            #speed must be greater than stall speed
            state['V'] >= Vstall,


            #Figure out how to delete
            Vstall == 120*units('kts'),
            WLoadmax == 6664 * units('N/m^2'),

            #compute the drag
            TCS([D >= self.wingP['D_{wing}'] + self.fuseP['D_{fuse}']]),

            #constraint CL and compute the wing loading
            W_avg == .5*self.wingP['C_{L}']*self.aircraft['S']*state.atm['\\rho']*state['V']**2,      
            WLoad == .5*self.wingP['C_{L}']*self.aircraft['S']*state.atm['\\rho']*state['V']**2/self.aircraft.wing['S'],

            #set average weight equal to the geometric avg of start and end weight
            W_avg == (W_start * W_end)**.5,

            #constrain the max wing loading
            WLoad <= WLoadmax,

            #compute fuel burn from TSFC

fan = 1.60474
        lpc  = 4.98
        hpc = 35./8.

        #Percent of velocity loss from BL smeared across entire fan
        BLIVloss = 0.088197
        #Min cruise mach number
        Mcruisemin = 0.72

        fBLI = 0.4
        wake_benefit = 0.02
        wake_percent = 0.33

        substitutions = {
                'N_{land}': 6.,
                'p_s': 81.*units('cm'),
                'n_{eng}': 2.,
                'n_{aisle}':2.,
                'W_{avg. pass}': 180.*units('lbf'),
                'W_{carry on}': 15.*units('lbf'),
                'W_{checked}':40.*units('lbf'),
                'W_{fix}': 3000.*units('lbf'),
                'w_{aisle}': 0.51*units('m'),
                'w_{seat}': 0.5*units('m'),
                'w_{sys}': 0.1*units('m'),
                'r_E': 1.,  # [TAS]
                'p_{\\lambda_{vt}}':1.6,
                '\\lambda_{cone}': 0.3,  # [TAS]
                '\\rho_{cone}': 2700.,#*units('kg/m^3'),  # [TAS]
                '\\rho_{bend}': 2700.,#*units('kg/m^3'),  # [TAS]
                '\\rho_{floor}': 2700.,#*units('kg/m^3'),  # [TAS]
                '\\rho_{skin}': 2700.,#*units('kg/m^3'),  # [TAS]

#delta_tip <= b/5])

        # Engine Model
        W_eng = Variable('W_{eng}', 'lbf', 'Engine weight')
        P_shaft_max_sl = Variable('P_{shaft max sl}', 5.4,'hp', 'Maximum power at sea level')
        P_shaft_max =  Variable('P_{shaft max}', 'hp', 'Maximum power at altitude')
        eta_t = Variable(r'\eta_t', 0.5, '-', 'percent throttle')
        eng_cnst = Variable('eng_{cnst}', 0.5538, '-',
                            'engine constant based off of power weight model')
        W_engtot = Variable('W_{eng-tot}', 'lbf', 'Installed engine weight')
        W_engref = Variable('W_{eng-ref}', 4.4107, 'lbf', 'Reference engine weight')
        P_shaftref = Variable('P_{shaft-ref}', 2.295, 'hp', 'reference shaft power')

        # Engine Weight Constraints
        constraints.extend([W_eng/W_engref >= 0.5538*(P_shaft/P_shaftref)**1.075,
                            W_engtot >= 2.572*W_eng**0.922*units('lbf')**0.078])

        # Weights and masses
        W_payload = Variable('W_{fix}',10,'lbf', 'Payload weight')
        W_fuel = Variable('W_{fuel}', 'lbf', 'Fuel weight')
        W_zfw = Variable('W_{zfw}', 'lbf', 'Zero fuel weight')
        W_avionics = Variable('W_{avionics}', 2, 'lbf', 'Avionics weight')
        W_tank = Variable('W_{tank}', 'lbf', 'Fuel tank weight')
        W_wing = Variable('W_{wing}', 'lbf', 'Total wing structural weight')
        W = Variable('W', 'lbf', 'Aircraft weight') 

        # Weight Constraints
        constraints.extend([W_wing >= m_skin*g + m_cap*g,
                            W_cent >= W_payload + W_engtot + W_avionics + m_fuse*g,
                            W_zfw >= W_cent + W_wing*1.2,
                            W >= W_fuel + W_zfw])
	# Fuel Volume