%pylab inline
%config InlineBackend.figure_format = 'retina'

import gpkit
import gpkit.interactive
gpkit.interactive.init_printing()
mon = gpkit.mon
vec = gpkit.vecmon

# Constants #
rho = mon("\\rho", 1.23, "kg/m^3", "Density of Air")
W   = mon("W", "N", "Weight of Vehicle")
xi =  vec(2,"\\xi", "N", "Thrust per bin")

# Scalars #
A     = mon("A", "m^2", "Disk Area")
Omega = mon("\\Omega", "rpm", "Rotor RPM")
R     = mon("R", "m", "Rotor Radius")
P     = mon("P", "W", "Induced Power")

# Vectors #
r   = vec(2, "r", "-", "Non-dimensional Radius")
dr  = vec(2, "\\Delta r", "-", "Non-dimensional Radius Step")
Vi  = vec(2, "V_i", "m/s", "Induced Velocities")
dCT = vec(2, "dC_T", "-", "Incremental Thrust Coefficient of Each Bin")
dCP = vec(2, "dC_P", "-", "Incremental Power Coefficient of Each Bin")
dP  = vec(2, "dP", "W", "Incremental Power of Each Bin")

def BEMT_GP(N=5, Weight=1e4):
    Weight = float(Weight)
    global xi, r, dr, Vi, dCT, dCP, dP
    
    xi =  vec(N, "\\xi", "N", "Thrust per bin")
    
    constants = {W: Weight,
                 xi: Weight/float(N) * ones(N),
                 #xi: 2*Weight/float(N)**2 * array(range(1,N+1))
                 }
    
    # Vector reasignment #
    r   = vec(N, "r", "-", "Non-dimensional Radius")
    dr  = vec(N, "\\Delta r", "-", "Non-dimensional Radius Step")
    Vi  = vec(N, "V_i", "m/s", "Induced Velocities")
    dCT = vec(N, "dC_T", "-", "Incremental Thrust Coefficient of Each Bin")
    dCP = vec(N, "dC_P", "-", "Incremental Power Coefficient of Each Bin")
    dP  = vec(N, "dP", "W", "Incremental Power of Each Bin")

    phys_constraints = (A == pi*R**2,
                        R <= 8 * gpkit.units.m,
                        Omega <= 280 * gpkit.units.rpm )

    r_constraints=([r[j] >= dr[:j].sum() for j in range(1,N)],
                   1 >= dr.sum(),
                   r[0] == dr[0]/2,
                   [r[j] >= r[j-1] + .5*dr[j-1] + .5*dr[j] for j in range(1,N)],
                   r[-1] <= 1)

    main_constraints = (xi == rho*A*(Omega*R)**2*dCT,
                        0.25 == Vi**2*r*dr/(dCT*(Omega*R)**2),
                        0.25 == Vi**3*r*dr/(dCP*(Omega*R)**3),                  
                        dP == rho*A*(Omega*R)**3*dCP,
                        P >= dP.sum())
    objective = P
    eqns = phys_constraints + r_constraints + main_constraints
    return gpkit.GP(objective, eqns, constants)

gp = BEMT_GP()

gp

gp.solve()["local_model"][0]

from IPython.display import Math

Math('0.5404\\frac{\\left({\\xi}_{1}^{0.0144} {\\xi}_{5}^{-0.0318} \\prod_{1}^5 {\\xi}_{i}^{0.3035} \\right)}{\\rho^{0.5}}')

gp.solution["variables"]

gp.solution["sensitivities"]["variables"]["\\xi"]

p = xi/rho
p

def solve_and_plot_Vi(bins, weight):
    gp = BEMT_GP(bins, weight)
    sol = gp.solve(solver='mosek', printing=False)

    figure(figsize=(12,6))
    plot(sol(R*r), sol(xi/(R*dr)), 'r.')
    ylim(0, weight/4)
    ylabel('Thrust density [N/m]')
    xlabel('Radial position [m]')
    show()

from IPython.html.widgets import interactive

interactive(solve_and_plot_Vi, bins=(1,20), weight=(1e3, 1e5, 1e3))