How to use the gpkit.Variable function in gpkit

def __init__(self):

        CL     = Variable('C_L', '-', 'Lift coefficient')
        CLmax  = Variable('C_{L_{max}}', '-', 'Max lift coefficient')
        CD     = Variable('C_D', '-', 'Drag coefficient')
        D      = Variable('D', 'N', 'Total aircraft drag (cruise)')
        Dfuse  = Variable('D_{fuse}', 'N', 'Fuselage drag')
        Dht    = Variable('D_{ht}', 'N', 'Horizontal tail drag')
        Dvt    = Variable('D_{vt}', 'N', 'Vertical tail drag')
        Dwing  = Variable('D_{wing}', 'N', 'Wing drag')
        LD     = Variable('\\frac{L}{D}', '-', 'Lift/drag ratio')
        Lh     = Variable('L_h', 'N', 'Horizontal tail downforce')
        Lw     = Variable('L_w', 'N', 'Wing lift')
        M      = Variable('M', '-', 'Cruise Mach number')
        R      = Variable('Range', 'nautical_miles', 'Range')
        Sw     = Variable('S_w', 'm**2', 'Wing reference area')
        Te     = Variable('T_e', 'N', 'Engine thrust at takeoff')
        TSFC   = Variable('TSFC', 'lb/lbf/hr',
                          'Thrust specific fuel consumption')
        V      = Variable('V_{\\infty}', 'm/s', 'Cruise velocity')
        VTO    = Variable('V_{TO}', 'm/s', 'Takeoff speed')
        W      = Variable('W', 'N', 'Total aircraft weight')
        Weng   = Variable('W_{eng}', 'N', 'Engine weight')
        Wfuel  = Variable('W_{fuel}', 'N', 'Fuel weight')
        Wfuse  = Variable('W_{fuse}', 'N', 'Fuselage weight')
        Wht    = Variable('W_{ht}', 'N', 'Horizontal tail weight')

def SimPleAC():
    "Creates SimpleAC model"
    # Env. constants
    g = Variable("g", 9.81, "m/s^2", "gravitational acceleration")
    mu = Variable("\\mu", 1.775e-5, "kg/m/s", "viscosity of air")
    rho = Variable("\\rho", 1.23, "kg/m^3", "density of air")
    rho_f = Variable("\\rho_f", 817, "kg/m^3", "density of fuel")

    # Non-dimensional constants
    C_Lmax = Variable("C_{L,max}", 1.6, "-", "max CL with flaps down")
    e = Variable("e", 0.92, "-", "Oswald efficiency factor")
    k = Variable("k", 1.17, "-", "form factor")
    N_ult = Variable("N_{ult}", 3.3, "-", "ultimate load factor")
    S_wetratio = Variable("(\\frac{S}{S_{wet}})", 2.075, "-",
                          "wetted area ratio")
    tau = Variable("\\tau", 0.12, "-", "airfoil thickness to chord ratio")
    W_W_coeff1 = Variable("W_{W_{coeff1}}", 2e-5, "1/m",
                          "wing weight coefficent 1")  # 12e-5 originally
    W_W_coeff2 = Variable("W_{W_{coeff2}}", 60, "Pa",
                          "wing weight coefficent 2")

    # Dimensional constants
    Range = Variable("Range", 3000, "km", "aircraft range")
    TSFC = Variable("TSFC", 0.6, "1/hr", "thrust specific fuel consumption")
    V_min = Variable("V_{min}", 25, "m/s", "takeoff speed")
    W_0 = Variable("W_0", 6250, "N", "aircraft weight excluding wing")

    # Free Variables
    LoD = Variable("L/D", "-", "lift-to-drag ratio")
    D = Variable("D", "N", "total drag force")
    V = Variable("V", "m/s", "cruising speed")
    W = Variable("W", "N", "total aircraft weight")

aircraft = Aircraft(Nclimb, Ncruise, enginestate)

        # vectorize
        with Vectorize(Nclimb):
            climb = ClimbSegment(aircraft, Nclimb)

        with Vectorize(Ncruise):
            cruise = CruiseSegment(aircraft, Nclimb)

        statelinking = StateLinking(climb.state, cruise.state, enginestate, Nclimb, Ncruise)

        # declare new variables
        W_ftotal = Variable('W_{f_{total}}', 'lbf', 'Total Fuel Weight')
        W_fclimb = Variable('W_{f_{climb}}', 'lbf',
                            'Fuel Weight Burned in Climb')
        W_fcruise = Variable('W_{f_{cruise}}', 'lbf',
                             'Fuel Weight Burned in Cruise')
        W_total = Variable('W_{total}', 'lbf', 'Total Aircraft Weight')
##        W_dry = Variable('W_{dry}', 'lbf', 'Dry Aircraft Weight')
        CruiseAlt = Variable('CruiseAlt', 'ft', 'Cruise Altitude [feet]')
        ReqRng = Variable('ReqRng', 'nautical_miles', 'Required Cruise Range')

        h = climb.state['h']
        hftClimb = climb.state['hft']
        dhft = climb.climbP['dhft']
        hftCruise = cruise.state['hft']

        # make overall constraints
        constraints = []

        constraints.extend([
            # weight constraints

#pressure ratios
        pihpt = Variable('\pi_{HPT}', '-', 'HPT Pressure Ratio')
        pilpt = Variable('\pi_{LPT}', '-', 'LPT Pressure Ratio')
        
        #turbine efficiences
        etaht = Variable('\eta_{ht}', 0.9, '-', 'Polytropic Efficiency of HPT')
        etalt = Variable('\eta_{lt}', 0.9, '-', 'Polytropic Efficiency of LPT')

        #flow faction f
        f = Variable('f', '-', 'Fuel Air Mass Flow Fraction')

        #BPR
        alpha = Variable('alpha', '-', 'By Pass Ratio')

        #relavent pressure ratio
        pitn = Variable('\pi_{tn}', '-', 'Turbine Nozzle Pressure Ratio')

        with SignomialsEnabled():
            constraints = [
                #HPT shafter power balance
                #SIGNOMIAL   
                SignomialEquality((1+f)*(ht41-ht45),ht3 - ht25),    #B.161

                #LPT shaft power balance
                #SIGNOMIAL  
                SignomialEquality((1+f)*(ht49 - ht45),-((ht25-ht18)+alpha*(ht21 - ht2))),    #B.165

                #HPT Exit states (station 4.5)
                Pt45 == pihpt * Pt41,
                pihpt == (Tt45/Tt41)**(4.60517),      #turbine efficiency is 0.9
                ht45 == Cpt1 * Tt45,

def __init__(self, ht, fuse, wing, fuseP, wingP, state, **kwargs):
        self.ht = ht
        self.fuse = fuse
        self.wing = wing
        self.fuseP = fuseP
        self.wingP = wingP
        
        #variables
        alpha   = Variable('\\alpha', '-', 'Horizontal tail angle of attack')
        D       = Variable('D_{ht}', 'N', 'Horizontal tail drag')
        Lh      = Variable('L_h', 'N', 'Horizontal tail downforce')
        SM      = VectorVariable(2, 'S.M.', '-', 'Stability margin')
        Rec     = Variable('Re_{c_h}', '-',
                           'Cruise Reynolds number (Horizontal tail)')
        CLah    = Variable('C_{L_{ah}}', '-', 'Lift curve slope (htail)')
        CLah0   = Variable('C_{L_{ah_0}}', '-',
                           'Isolated lift curve slope (htail)')
        CLaw    = Variable('C_{L_{aw}}', '-', 'Lift curve slope (wing)')

        
        CLh     = Variable('C_{L_h}', '-', 'Lift coefficient (htail)')
        dxlead  = Variable('\\Delta x_{{lead}_h}', 'm',
                           'Distance from CG to horizontal tail leading edge')
        dxtrail = Variable('\\Delta x_{{trail}_h}', 'm',
                           'Distance from CG to horizontal tail trailing edge')
        dxw     = VectorVariable(2, '\\Delta x_w', 'm',
                           'Distance from aerodynamic centre to CG')

def setup(self, aircraft, state, **kwargs):
        #submodels
        self.aircraft = aircraft
        self.aircraftP = AircraftP(aircraft, state)
        self.wingP = self.aircraftP.wingP
        self.fuseP = self.aircraftP.fuseP
                                  
        #variable definitions
        theta = Variable('\\theta', '-', 'Aircraft Climb Angle')
        excessP = Variable('excessP', 'W', 'Excess Power During Climb')
        RC = Variable('RC', 'feet/min', 'Rate of Climb/Decent')
        dhft = Variable('dhft', 'feet', 'Change in Altitude Per Climb Segment [feet]')
        RngClimb = Variable('RngClimb', 'nautical_miles', 'Down Range Covered in Each Climb Segment')

        #constraints
        constraints = []
        
        constraints.extend([ 
            RC == excessP/self.aircraftP['W_{avg}'],
            RC >= 500*units('ft/min'),
            
            #make the small angle approximation and compute theta
            theta * state['V']  == RC,
           
            dhft == self.aircraftP['tmin'] * RC,

k = Variable("k", 1.2, "-", "form factor")
    e = Variable("e", 0.95, "-", "Oswald efficiency factor")
    mu = Variable("\\mu", 1.78e-5, "kg/m/s", "viscosity of air")
    pi = Variable("\\pi", np.pi, "-", "half of the circle constant")
    rho = Variable("\\rho", 1.23, "kg/m^3", "density of air")
    tau = Variable("\\tau", 0.12, "-", "airfoil thickness to chord ratio")
    N_ult = Variable("N_{ult}", 3.8, "-", "ultimate load factor")
    V_min = Variable("V_{min}", 22, "m/s", "takeoff speed")
    C_Lmax = Variable("C_{L,max}", 1.5, "-", "max CL with flaps down")
    S_wetratio = Variable("(\\frac{S}{S_{wet}})", 2.05, "-",
                          "wetted area ratio")
    W_W_coeff1 = Variable("W_{W_{coeff1}}", 8.71e-5, "1/m",
                          "Wing Weight Coefficent 1")
    W_W_coeff2 = Variable("W_{W_{coeff2}}", 45.24, "Pa",
                          "Wing Weight Coefficent 2")
    CDA0 = Variable("(CDA0)", 0.031, "m^2", "fuselage drag area")
    W_0 = Variable("W_0", 4940.0, "N", "aircraft weight excluding wing")

    # Free Variables
    D = Variable("D", "N", "total drag force")
    A = Variable("A", "-", "aspect ratio")
    S = Variable("S", "m^2", "total wing area")
    V = Variable("V", "m/s", "cruising speed")
    W = Variable("W", "N", "total aircraft weight")
    Re = Variable("Re", "-", "Reynold's number")
    C_D = Variable("C_D", "-", "Drag coefficient of wing")
    C_L = Variable("C_L", "-", "Lift coefficent of wing")
    C_f = Variable("C_f", "-", "skin friction coefficient")
    W_w = Variable("W_w", "N", "wing weight")

    constraints = []

taucone  = Variable('\\tau_{cone}', 'N/m^2', 'Shear stress in cone')
        taufloor = Variable('\\tau_{floor}', 'N/m^2',
                            'Max allowable shear web stress')
        tcone    = Variable('t_{cone}', 'm', 'Cone thickness')
        tshell   = Variable('t_{shell}', 'm', 'Shell thickness')
        tskin    = Variable('t_{skin}', 'm', 'Skin thickness')
        waisle   = Variable('w_{aisle}', 'm', 'Aisle width')
        wfloor   = Variable('w_{floor}', 'm', 'Floor width')
        wfuse    = Variable('w_{fuse}', 'm', 'Fuselage width')
        wseat    = Variable('w_{seat}', 'm', 'Seat width')
        wsys     = Variable('w_{sys}', 'm',
                            'Width between cabin and skin for systems')
        xCGfu    = Variable('x_{CG_{fu}}', 'm', 'x-location of fuselage CG')
        xVbulk   = Variable('xVbulk', 'm^4', 'Volume moment of bulkhead')
        xVcyl    = Variable('xVcyl', 'm^4', 'Volume moment of cylinder')
        xVnose   = Variable('xVnose', 'm^4', 'Volume moment of nose')
        xWapu    = Variable('xWapu', 'N*m', 'Moment of APU')
        xWcone   = Variable('xWcone', 'N*m', 'Moment of cone')
        xWfix    = Variable('xWfix', 'N*m', 'Moment of fixed weights')
        xWfloor  = Variable('xWfloor', 'N*m', 'Moment of floor weight')
        xWfuse   = Variable('xWfuse', 'N*m', 'Fuselage moment')
        xWinsul  = Variable('xWinsul', 'N*m', 'Moment of insulation material')
        xWpadd   = Variable('xWpadd', 'N*m', 'Moment of misc weights')
        xWseat   = Variable('xWseat', 'N*m', 'Moment of seats')
        xWshell  = Variable('xWshell', 'N*m', 'Mass moment of shell')
        xWskin   = Variable('xWskin', 'N*m', 'Mass moment of skin')
        xWwindow = Variable('xWwindow', 'N*m', 'Mass moment of windows')
        x_upswp  = Variable('x_{up}', 'm', 'Fuselage upsweep point')
        xapu     = Variable('xapu', 'ft', 'x-location of APU')
        xconend  = Variable('xconend', 'm', 'x-location of cone end')
        xfix     = Variable('xfix', 'm', 'x-location of fixed weight')
        xshell1  = Variable('x_{shell1}', 'm', 'Start of cylinder section')

self.aircraft = aircraft
        self.wingP = aircraft.wing.dynamic(state)
        self.fuseP = aircraft.fuse.dynamic(state)
        self.engineP = aircraft.engine.dynamic(state)
        self.Pmodels = [self.wingP, self.fuseP, self.engineP]

        #variable definitions
        Vstall   = Variable('V_{stall}', 'knots', 'Aircraft Stall Speed')
        D        = Variable('D', 'N', 'Total Aircraft Drag')
        W_avg    = Variable('W_{avg}', 'N', 'Geometric Average of Segment Start and End Weight')
        W_start  = Variable('W_{start}', 'N', 'Segment Start Weight')
        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}']]),

Nw                  = Variable("N_{w}",                   "ft",       "Nacelle Width")
Nz                  = Variable("N_{z}",                   "-",        "Ultimate Load Factor; = 1.5 x Limit Load Factor")
Rkva                = Variable("R_{kva}",                 "kV*A",     "System Electrical Rating (Typically 40-60) ")
Scs                 = Variable("S_{cs}",                  "ft^2",     "Total Area of Control Surfaces")
Scsw                = Variable("S_{csw}",                 "ft^2",     "Control Surface Area (Wing Mounted")
Se                  = Variable("S_{e}",                   "ft^2",     "Elevator Area")
Sf                  = Variable("S_{f}",                   "ft^2",     "Fuselage Wetted Area")
S_floor             = Variable("S_{floor}",               "ft^2",     "Cargo Floor Area")
Sht                 = Variable("S_{ht}",                  "ft^2",     "Horizontal Tail Area")
Sn                  = Variable("S_{n}",                   "ft^2",     "Nacelle Wetted Area")
Svt                 = Variable("S_{vt}",                  "ft^2",     "Vertical Tail Area")
Sw                  = Variable("S_{w}",                   "ft^2",     "Trapezoidal Wing Area")
t_c_root_w          = Variable("\\frac{t}{c}_{root_w}",   "-",        "Thickness to Chord Ratio of Main Wing at the Root")
t_c_root_v          = Variable("\\frac{t}{c}_{root_v}",   "-",        "Thickness to Chord Ratio of Vertical Tail at the Root")
Vi                  = Variable("V_{i}",                   "gal",      "Volume of Integral Tanks")
Vp                  = Variable("V_{p}",                   "gal",      "Volume of Self Sealing (Protected) Tanks")
Vpr                 = Variable("V_{pr}",                  "ft^3",     "Volume of Pressurized Section")
V_stall             = Variable("V_{stall}",               "m/s",      "Stall Speed")
Vt                  = Variable("V_{t}",                   "gal",      "Total Fuel Volume")
W_APU_uninstalled   = Variable("W_{APU_{uninstalled}}",   "lb",       "Uninstalled APU Weight")
Wc                  = Variable("W_{c}",                   "lb",       "Maximum Cargo Weight")
Wdg                 = Variable("W_{dg}",                  "lb",       "Flight Design Gross Weight")
Wec                 = Variable("W_{ec}",                  "lb",       "Weight of Engine and Contents (per nacelle); \\approx 2.331 W_en^0.901 K_p K_tr")
Wl                  = Variable("W_{l}",                   "lb",       "Landing Design Gross Weight")
Wuav                = Variable("W_{uav}",                 "-",        "Uninstalled Avionics Weight (Typically 800-1400")
Lambda_w            = Variable("\\Lambda_{w}",            "-",        "Main Wing Sweep at 0.25 MAC")
Lambda_ht           = Variable("\\Lambda_{ht}",           "-",        "Horizontal Tail Sweep at 0.25 MAC")
Lambda_vt           = Variable("\\Lambda_{vt}",           "-",        "Vertical Tail Sweep at 0.25 MAC")
lambda_w            = Variable("\\lambda_{w}",            "-",        "Main Wing Taper Ratio")

W_wing                            = Variable("W_{wing}",                   "lbf", "Weight of Main Wing")
W_horizontal_tail                 = Variable("W_{horizontal-tail}",        "lbf", "Weight of Horizontal Tail")