How to use gpkit - 10 common examples

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)

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')