import gpkit
import gpkit.interactive
gpkit.interactive.init_printing()

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

var = gpkit.Variable
vec = gpkit.VectorVariable

dCpdy = var("\\frac{dC_P}{dy}")
a = var("a")
b = var("b")
lam = var("\\lambda")
s = var("s")
xi = var("\\xi")
F = var("F")
y = var("y")
eps = var("\\epsilon")

gp = gpkit.GP(1/dCpdy,
              [1 >= a + b,
               2*a*b >= s*lam*y + s*xi,
               xi**2 >= (lam*y)**2 + 4*a*b,
               s*b >= dCpdy/(8*lam*y**2*F) + eps*a*b],
              {F: 1, y: 0.75, eps: 1.0/30.0})
gp

print gp.solve().table()

B = 3
N = 10

f = var("f")
pi = var("\\pi", np.pi)
DCp = var("(\Delta C_P)")
xi = vec(N, "\\xi")
zeta = vec(N, "\\zeta")
y = vec(N, "y")
dy = vec(N, "(\Delta y)")
lam = var("\\lambda", ('sweep', np.linspace(0.01, 8, 20)))

constraints = [1 >= a + b,
               2*a*b >= s*lam*y + s*xi,
               xi**2 >= (lam*y)**2 + 4*a*b,
               s*b >= DCp/(8*lam*dy**2*F) + eps*a*b]

constraints += [ y[-1] + dy[-1]/2 <= 1,
                [y[i] + dy[i]/2 + dy[i+1]/2 <= y[i+1] for i in range(N-1)],
                 #B*zeta[0]/(2*pi) + dy[0]/2 <= y[0],#switch0
                 dy[0]/2 <= y[0],#switch0
                 1/y >= 1 + 2*f*s/(B*a),
                 1 >= F**1.56 + 1.655*F**5.51*f**-2.76   ]

sol = gpkit.GP(1/(N*DCp), constraints, {eps: ('sweep', [0.2, 0.1, 0.05, 0.033])}).solve('mosek')
print sol.table()

B = 3
N = 20
Nlams = 20

f = var("f")
pi = var("\\pi", np.pi)
DCp = var("(\Delta C_P)")
xi = vec(N, "\\xi")
zeta = vec(N, "\\zeta")
y = vec(N, "y")
dy = vec(N, "(\Delta y)")
lam = var("\\lambda", ('sweep', np.linspace(0.01, 8, Nlams)))

constraints = [1 >= a + b,
               2*a*b >= s*lam*y + s*xi,
               xi**2 >= (lam*y)**2 + 4*a*b,
               s*b >= DCp/(8*lam*dy**2*F*N/1.7) + eps*a*b]

constraints += [ y[-1] + dy[-1]/2 <= 1,
                [y[i] + dy[i]/2 + dy[i+1]/2 <= y[i+1] for i in range(N-1)],
                 #B*zeta[0]/(2*pi) + dy[0]/2 <= y[0],#switch0
                 dy[0]/2 <= y[0],#switch0
                 1/y >= 1 + 2*f*s/(B*a),
                 1 >= F**1.56 + 1.655*F**5.51*f**-2.76   ]

sol = gpkit.GP(1/(N*DCp), constraints, {eps: ('sweep', [0.2, 0.1, 0.05, 0.033])}).solve('mosek')

plt.plot(sol(lam).reshape(4, Nlams).T, sol(N*DCp).reshape(4, Nlams).T, '--', linewidth=1.5)
plt.gca().set_color_cycle(None)
Fsub = sol(F)
lamsub = sol(lam)
epssub = sol(eps)

##

lams = []
NDCps = []
for j in range(Fsub.size):
    lam = var("\\lambda", lamsub[j])

    constraints = [1 >= a + b,
                   2*a*b >= s*lam*y + s*xi,
                   xi**2 >= (lam*y)**2 + 4*a*b,
                   s*b >= DCp/(8*lam*dy**2*F*N/1.7) + eps*a*b]

    constraints += [ y[-1] + dy[-1]/2 <= 1,
                    [y[i] + dy[i]/2 + dy[i+1]/2 <= y[i+1] for i in range(N-1)],
                     #B*zeta[0]/(2*pi) + dy[0]/2 <= y[0],#switch0
                     dy[0]/2 <= y[0],#switch0
                    ]

    sol = gpkit.GP(1/(N*DCp), constraints, {F: Fsub[j], eps: epssub[j]}).solve('mosek', printing=False)
    
    lams.append(sol(lam))
    NDCps.append(sol(N*DCp))
    
lams = np.array(lams)
NDCps = np.array(NDCps)
plt.plot(lams.reshape(4, Nlams).T, NDCps.reshape(4, Nlams).T, linewidth=1.5)
plt.plot([0, 8], [16.0/27.0, 16.0/27.0], 'k')
plt.ylim([0, 0.7])
plt.xlabel('Tip Speed Ratio (lambda)')
plt.ylabel('C_P')

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()
    styles = "<style>\n%s\n</style>" % (open(os.path.join(base,'profile_default/static/custom/custom.css'),'r').read())
    return HTML(styles)
css_styling()