import numpy as np
from gpkit.shortcuts import *
import gpkit.interactive
%matplotlib inline

"""
Explanation: <img src="fuellogo.svg" style="float:left; padding-right:1em;" width=150 />
AIRPLANE FUEL
Minimize fuel burn for a plane that can sprint and land quickly.
Set up the modelling environment
First we'll to import GPkit and turn on $\LaTeX$ printing for GPkit variables and equations.
End of explanation
"""


N_lift         = Var("N_{lift}", 6.0, "-", "Wing loading multiplier")
pi             = Var("\\pi", np.pi, "-", "Half of the circle constant")
sigma_max      = Var("\\sigma_{max}", 250e6, "Pa", "Allowable stress, 6061-T6")
sigma_maxshear = Var("\\sigma_{max,shear}", 167e6, "Pa", "Allowable shear stress")
g              = Var("g", 9.8, "m/s^2", "Gravitational constant")
w              = Var("w", 0.5, "-", "Wing-box width/chord")
r_h            = Var("r_h", 0.75, "-", "Wing strut taper parameter")
f_wadd         = Var("f_{wadd}", 2, "-", "Wing added weight fraction")
W_fixed        = Var("W_{fixed}", 14.7e3, "N", "Fixed weight")
C_Lmax         = Var("C_{L,max}", 1.5, "-", "Maximum C_L, flaps down")
rho            = Var("\\rho", 0.91, "kg/m^3", "Air density, 3000m")
rho_sl         = Var("\\rho_{sl}", 1.23, "kg/m^3", "Air density, sea level")
rho_alum       = Var("\\rho_{alum}", 2700, "kg/m^3", "Density of aluminum")
mu             = Var("\\mu", 1.69e-5, "kg/m/s", "Dynamic viscosity, 3000m")
e              = Var("e", 0.95, "-", "Wing spanwise efficiency")
A_prop         = Var("A_{prop}", 0.785, "m^2", "Propeller disk area")
eta_eng        = Var("\\eta_{eng}", 0.35, "-", "Engine efficiency")
eta_v          = Var("\\eta_v", 0.85, "-", "Propeller viscous efficiency")
h_fuel         = Var("h_{fuel}", 42e6, "J/kg", "fuel heating value")
V_sprint_reqt  = Var("V_{sprintreqt}", 150, "m/s", "sprint speed requirement")
W_pay          = Var("W_{pay}", 500*9.81, "N")
R_min          = Var("R_{min}", 1e6, "m", "Minimum airplane range")
V_stallmax     = Var("V_{stall,max}", 40, "m/s", "Stall speed")
# sweep variables
R_min          = Var("R_{min}", 5e6, "m", "Minimum airplane range")
V_stallmax     = Var("V_{stall,max}", 40, "m/s", "Stall speed")

"""
Explanation: declare constants
End of explanation
"""


V        = Vec(3, "V", "m/s", "Flight speed")
C_L      = Vec(3, "C_L", "-", "Wing lift coefficent")
C_D      = Vec(3, "C_D", "-", "Wing drag coefficent")
C_Dfuse  = Vec(3, "C_{D_{fuse}}", "-", "Fuselage drag coefficent")
C_Dp     = Vec(3, "C_{D_p}", "-", "drag model parameter")
C_Di     = Vec(3, "C_{D_i}", "-", "drag model parameter")
T        = Vec(3, "T", "N", "Thrust force")
Re       = Vec(3, "Re", "-", "Reynold's number")
W        = Vec(3, "W", "N", "Aircraft weight")
eta_i    = Vec(3, "\\eta_i", "-", "Aircraft efficiency")
eta_prop = Vec(3, "\\eta_{prop}", "-")
eta_0    = Vec(3, "\\eta_0", "-")
W_fuel   = Vec(2, "W_{fuel}", "N", "Fuel weight")
z_bre    = Vec(2, "z_{bre}", "-")
S        = Var("S", "m^2", "Wing area")
R        = Var("R", "m", "Airplane range")
A        = Var("A", "-", "Aspect Ratio")
I_cap    = Var("I_{cap}", "m^4", "Spar cap area moment of inertia per unit chord")
M_rbar   = Var("\\bar{M}_r", "-")
P_max    = Var("P_{max}", "W")
V_stall  = Var("V_{stall}", "m/s")
nu       = Var("\\nu", "-")
p        = Var("p", "-")
q        = Var("q", "-")
tau      = Var("\\tau", "-")
t_cap    = Var("t_{cap}", "-")
t_web    = Var("t_{web}", "-")
W_cap    = Var("W_{cap}", "N")
W_zfw    = Var("W_{zfw}", "N", "Zero fuel weight")
W_eng    = Var("W_{eng}", "N")
W_mto    = Var("W_{mto}", "N", "Maximum takeoff weight")
W_pay    = Var("W_{pay}", "N")
W_tw     = Var("W_{tw}", "N")
W_web    = Var("W_{web}", "N")
W_wing   = Var("W_{wing}", "N")

"""
Explanation: declare free variables
End of explanation
"""


W == 0.5*rho*C_L*S*V**2

"""
Explanation: Check the vector constraints:
End of explanation
"""


steady_level_flight = (W == 0.5*rho*C_L*S*V**2,
                       T >= 0.5*rho*C_D*S*V**2,
                       Re == (rho/mu)*V*(S/A)**0.5)

landing_fc = (W_mto <= 0.5*rho_sl*V_stall**2*C_Lmax*S,
              V_stall <= V_stallmax)

sprint_fc = (P_max >= T[2]*V[2]/eta_0[2],
             V[2] >= V_sprint_reqt)

drag_model = (C_D >= (0.05/S)*gpkit.units.m**2 +C_Dp + C_L**2/(pi*e*A),
              1 >= (2.56*C_L**5.88/(Re**1.54*tau**3.32*C_Dp**2.62) +
                   3.8e-9*tau**6.23/(C_L**0.92*Re**1.38*C_Dp**9.57) +
                   2.2e-3*Re**0.14*tau**0.033/(C_L**0.01*C_Dp**0.73) +
                   1.19e4*C_L**9.78*tau**1.76/(Re*C_Dp**0.91) +
                   6.14e-6*C_L**6.53/(Re**0.99*tau**0.52*C_Dp**5.19)))

propulsive_efficiency = (eta_0 <= eta_eng*eta_prop,
                         eta_prop <= eta_i*eta_v,
                         4*eta_i + T*eta_i**2/(0.5*rho*V**2*A_prop) <= 4)

# 4th order taylor approximation for e^x
z_bre_sum = 0
for i in range(1,5):
    z_bre_sum += z_bre**i/np.math.factorial(i)

range_constraints = (R >= R_min,
                     z_bre >= g*R*T[:2]/(h_fuel*eta_0[:2]*W[:2]),
                     W_fuel/W[:2] >= z_bre_sum)

punits = gpkit.units.parse_expression('N/W^0.8083')
weight_relations = (W_pay >= 500*g*gpkit.units.kg,
                    W_tw >= W_fixed + W_pay + W_eng,
                    W_zfw >= W_tw + W_wing,
                    W_eng >= 0.0372*P_max**0.8083 * punits,
                    W_wing/f_wadd >= W_cap + W_web,
                    W[0] >= W_zfw + W_fuel[1],
                    W[1] >= W_zfw,
                    W_mto >= W[0] + W_fuel[0],
                    W[2] == W[0])

wunits = gpkit.units.m**-4
munits = gpkit.units.parse_expression('Pa*m**6')
wing_structural_model = (2*q >= 1 + p,
                         p >= 2.2,
                         tau <= 0.25,
                         M_rbar >= W_tw*A*p/(24*gpkit.units.N),
                         .92**2/2.*w*tau**2*t_cap >= I_cap * wunits + .92*w*tau*t_cap**2,
                         8 >= N_lift*M_rbar*A*q**2*tau/S/I_cap/sigma_max * munits,
                         12 >= A*W_tw*N_lift*q**2/tau/S/t_web/sigma_maxshear,
                         nu**3.94 >= .86*p**-2.38 + .14*p**0.56,
                         W_cap >= 8*rho_alum*g*w*t_cap*S**1.5*nu/3/A**.5,
                         W_web >= 8*rho_alum*g*r_h*tau*t_web*S**1.5*nu/3/A**.5
                         )

eqns = (weight_relations + range_constraints + propulsive_efficiency
        + drag_model + steady_level_flight + landing_fc + sprint_fc + wing_structural_model)

m = gpkit.Model(W_fuel.sum(), eqns)

"""
Explanation: Form the optimization problem
In the 3-element vector variables, indices 0, 1, and 2 are the outbound, return and sprint flights.
End of explanation
"""


m.interact()

"""
Explanation: Design an airplane
End of explanation
"""


m.solution["localmodel"]

"""
Explanation: The "local model" is the power-law tangent to the Pareto frontier, gleaned from sensitivities.
End of explanation
"""


from gpkit.interactive.plotting import sensitivity_plot
_ = sensitivity_plot(m)

"""
Explanation: plot design frontiers
End of explanation
"""


from string import Template

fuelupdate_js = Template("""
var W_eng = $W_eng,
    lam = $lam

fuel.shearinner.scalex = 1-$tcap*10
fuel.shearinner.scaley = 1-$tweb*100
fuel.airfoil.scaley = $tau/0.13
fuel.fuse.scalex = $W_fus/24000
fuel.wing.scalex = $b/2/14
fuel.wing.scaley = $cr*1.21
""")

def fuelupdate_py(sol):
    varstrs = ("p", "S", "A", "t_{cap}", "t_{web}", "w",
               "\\tau", "W_{eng}", "W_{mto}", "W_{wing}")
    p, S, A, t_cap, t_web, w, tau, W_eng, W_mto, W_wing = sol.getvars(*varstrs)
    lam = 0.5*(p-1)
    return fuelupdate_js.substitute(lam = lam,
                                    b = (S*A)**0.5,
                                    cr = 2/(1+lam)*(S/A)**0.5,
                                    tcap = t_cap/tau,
                                    tweb = t_web/w,
                                    tau = tau,
                                    W_eng = W_eng,
                                    W_fus = W_mto - W_wing - W_eng)

fuelconstraint_js = """
fuel.engine1.scale = Math.pow(W_eng/3000, 2/3)
fuel.engine2.scale = Math.pow(W_eng/3000, 2/3)
fuel.engine1.y = 6*lam
fuel.engine2.y = 6*lam
fuel.wingrect.scaley = 1-lam
fuel.wingrect.y = -6 + 5*lam
fuel.wingtaper.scaley = lam
fuel.wingtaper.y = 5*lam
"""

gpkit.interactive.showcadtoon("fuel",
"position:absolute; height:0; right:0; top:24em;")

gpkit.interactive.ractorpy(m, fuelupdate_py,
         {"V_{stall,max}": (20, 50, 0.3),
          "R_{min}": (1e6, 1e7, 1e5),
          "V_{sprintreqt}": (100, 200, 1)},
         fuelconstraint_js)

gpkit.interactive.ractorjs("fuel", m, fuelupdate_py,
         {"V_{stall,max}": (20, 50, 3),
          "R_{min}": (1e6, 1e7, 2e6),
          "V_{sprintreqt}": (100, 200, 20)},
         fuelconstraint_js)

"""
Explanation: Interactive analysis
Let's investigate it with the cadtoons library. Running cadtoon.py flightconditions.svg in this folder creates an interactive SVG graphic for us.
First, import the functions to display HTML in iPython Notebook, and the ractivejs library.
End of explanation
"""


from IPython import utils
from IPython.core.display import HTML
import os
def css_styling():
    """Load default custom.css file from ipython profile"""
    base = utils.path.get_ipython_dir()
    csspath = os.path.join(base,'profile_default/static/custom/custom.css')
    styles = "<style>\n%s\n</style>" % (open(csspath,'r').read())
    return HTML(styles)
css_styling()

"""
Explanation: This concludes the aircraft example. Try playing around with the sliders up above until you're bored; then check out one of the other examples. Thanks for reading!
Import CSS for nbviewer
If you have a local iPython stylesheet installed, this will add it to the iPython Notebook:
End of explanation
"""


