import numpy as np
import matplotlib.pyplot as plt


def Rho_dtic(TC):
    """ valid between -45C and +25C"""

    return 0.7987 - 0.0036*(TC - 25.)

def C2K(C):
    return C + 273.15
def K2C(K):
    return K - 273.15

def C2F(C):
    return (9./5.)*C + 32.
def F2C(F):
    return (5./9.)*(F-32.)

def C2R(C):
    return (9./5.)*(C+273.15)
def R2C(R):
    return (5./9.)*R - 273.15

def F2K(F):
    return (5./9.)*(F-32.) + 273.15
def K2F(K):
    return  (9./5.)*(K-273.15) + 32.



C = 5. * np.arange(15) -45.

C_lowextrap = -85.
C_ext = np.linspace(C_lowextrap, C[0], 5)

print "rho at -30F: ", Rho_dtic((-30.-32.)*(5./9.))

rho = Rho_dtic(C)
rho_ext  = Rho_dtic(C_ext)

K, K_ext = C2K(C), C2K(C_ext)
F, F_ext = C2F(C), C2F(C_ext)
R, R_ext = C2R(C), C2R(C_ext)

K_limits = np.array([50, 300])
C_limits = K2C(K_limits)
F_limits = K2F(K_limits)
R_limits = C2R(K2C(K_limits))

put_rat = 0.03
K_put = K_limits.max() - put_rat*(K_limits.max()-K_limits.min())
C_put = C_limits.max() - put_rat*(C_limits.max()-C_limits.min())
F_put = F_limits.max() - put_rat*(F_limits.max()-F_limits.min())
R_put = R_limits.max() - put_rat*(R_limits.max()-R_limits.min())

R_info = np.array([302, 335, 350, 355, 380, 420])
C_info = R2C(R_info)
K_info = C2K(R2C(R_info))
F_info = C2F(R2C(R_info))
info = ['hard wax', 'solid wax', 'soft wax', 'gel',
        'gel, thick liquid', 'liquid w/ some solids']

# NOW LOX!!!!!!

"""
# A10.....
R_A10 = np.array([98.0, 110, 120, 130, 140, 150, 160, 165])  # R is absolute Farenheit
C_A10 = R2C(R_A18)
K_A10 = C2K(R2C(R_A18))
F_A10 = C2F(R2C(R_A18))
lb_A10 = np.array([81.5, 79.5, 78.2, 76.7, 74.9, 73.2, 71.5, 70.8])
f = 0.01602  # g/cm^3 per lb/ft^3
rho_A10 = f*lb_A18
"""

# A02.....
K_A02 = np.array([68.0, 70, 74, 78, 80, 82, 86, 89])
C_A02 = K2C(K_A02)
F_A02 = K2F(K_A02)
R_A02 = C2R(K2C(K_A02))
rho_A02 = np.array([1.2489, 1.2393, 1.2200, 1.2008,
                   1.1911,1.1815, 1.1623, 1.1479 ])

# A18.....
A, B, n = 0.43533, 0.28772, 0.29240
Tc, Tmin, Tmax = 154.58, 54.35, 154.58

K_A18 = 55. + 5.*np.arange(8)
C_A18 = K2C(K_A18)
F_A18 = K2F(K_A18)
R_A18 = C2R(K2C(K_A18))

def Rho_L(T):
    return A * (B**(-(1 - T/Tc)**n))

rho_A18 = Rho_L(K_A18)

F_NK33_LOX = -310.
F_NK33_RP1 =  -30.

C_NK33_LOX = F2C(F_NK33_LOX)
K_NK33_LOX = C2K(F2C(F_NK33_LOX))
R_NK33_LOX = C2R(F2C(F_NK33_LOX))

C_NK33_RP1 = F2C(F_NK33_RP1)
K_NK33_RP1 = C2K(F2C(F_NK33_RP1))
R_NK33_RP1 = C2R(F2C(F_NK33_RP1))

rho_NK33_LOX = Rho_L(C2K(F2C(F_NK33_LOX)))
rho_NK33_RP1 = Rho_dtic(F2C(F_NK33_RP1))

NK33_LOX_label = "NK-33  LOX  -310F"
NK33_RP1_label = "NK-33  RP1   -30F"
C_NK33_LOXlabel   = -180.
C_NK33_RP1label   =  -70.
F_NK33_LOXlabel   = C2F(C_NK33_LOXlabel)
F_NK33_RP1label   = C2F(C_NK33_RP1label)
K_NK33_LOXlabel   = C2K(C_NK33_LOXlabel)
K_NK33_RP1label   = C2K(C_NK33_RP1label)
R_NK33_LOXlabel   = C2R(C_NK33_LOXlabel)
R_NK33_RP1label   = C2R(C_NK33_RP1label)


fig = plt.figure(figsize=[9,10])

rho_limits = np.array([0.75, 1.35])
h = rho_limits.min() + 0.01
h_put = rho_limits.max() - 0.02
#bbox = {'fc': '0.8', 'pad': 0}
thing     = {'ha': 'center', 'va': 'bottom', 'bbox': None,
             'fontsize':15}
thing_put = {'ha': 'right',  'va': 'top', 'bbox': None,
             'fontsize':18}
thing_LOX = {'ha': 'left',  'va': 'bottom', 'bbox': None,
             'fontsize':14}
thing_RP1 = {'ha': 'right',  'va': 'bottom', 'bbox': None,
             'fontsize':14}
ticklab_fontsize = 15
lws, lwd, lwn = 2, 1.5, 1.0

# plt.subplot(4,1,1)
ax = plt.subplot(4, 1, 1)
for item in (ax.get_xticklabels() + ax.get_yticklabels()):
    item.set_fontsize(ticklab_fontsize)
plt.plot(C, rho, lw=lws)
plt.plot(C_ext, rho_ext, '--r', lw=lwd)
plt.xlim(C_limits)
plt.ylim(rho_limits)
plt.minorticks_on()
plt.text(C_put, h_put, "Centigrade", thing_put)
for a, b in zip(C_info, info):
    plt.text(a, h, b, thing, rotation=90)
plt.plot(C_A02, rho_A02, 'or')
plt.plot(C_A18, rho_A18, lw=lws)
plt.plot(C_limits, [rho_NK33_LOX]*2, '-k', lw=lwn)
plt.plot(C_limits, [rho_NK33_RP1]*2, '--k', lw=lwn)
plt.text(C_NK33_LOXlabel, rho_NK33_LOX+0.02, NK33_LOX_label, thing_LOX)
plt.text(C_NK33_RP1label, rho_NK33_RP1+0.02, NK33_RP1_label, thing_RP1)






ax = plt.subplot(4, 1, 2)
for item in (ax.get_xticklabels() + ax.get_yticklabels()):
    item.set_fontsize(ticklab_fontsize)
plt.plot(K, rho, lw=lws)
plt.plot(K_ext, rho_ext, '--r', lw=lwd)
plt.xlim(K_limits)
plt.ylim(rho_limits)
plt.minorticks_on()
plt.text(K_put, h_put, "Kelvin", thing_put)
for a, b in zip(K_info, info):
    plt.text(a, h, b, thing, rotation=90)
plt.plot(K_A02, rho_A02, 'or')
plt.plot(K_A18, rho_A18, lw=lws)
plt.plot(K_limits, [rho_NK33_LOX]*2, '-k', lw=lwn)
plt.plot(K_limits, [rho_NK33_RP1]*2, '--k', lw=lwn)
plt.text(K_NK33_LOXlabel, rho_NK33_LOX+0.02, NK33_LOX_label, thing_LOX)
plt.text(K_NK33_RP1label, rho_NK33_RP1+0.02, NK33_RP1_label, thing_RP1)



ax = plt.subplot(4, 1, 3)
for item in (ax.get_xticklabels() + ax.get_yticklabels()):
    item.set_fontsize(ticklab_fontsize)
plt.plot(F, rho, lw=lws)
plt.plot(F_ext, rho_ext, '--r', lw=lwd)
plt.xlim(F_limits)
plt.ylim(rho_limits)
plt.minorticks_on()
plt.text(F_put, h_put, "Farenheit", thing_put)
for a, b in zip(F_info, info):
    plt.text(a, h, b, thing, rotation=90)
plt.plot(F_A02, rho_A02, 'or')
plt.plot(F_A18, rho_A18, lw=lws)
plt.plot(F_limits, [rho_NK33_LOX]*2, '-k', lw=lwn)
plt.plot(F_limits, [rho_NK33_RP1]*2, '--k', lw=lwn)
plt.text(F_NK33_LOXlabel, rho_NK33_LOX+0.02, NK33_LOX_label, thing_LOX)
plt.text(F_NK33_RP1label, rho_NK33_RP1+0.02, NK33_RP1_label, thing_RP1)



ax = plt.subplot(4, 1, 4)
for item in (ax.get_xticklabels() + ax.get_yticklabels()):
    item.set_fontsize(ticklab_fontsize)
plt.plot(R, rho, lw=lws)
plt.plot(R_ext, rho_ext, '--r', lw=lwd)
plt.xlim(R_limits)
plt.ylim(rho_limits)
plt.minorticks_on()
plt.text(R_put, h_put, "Rankine", thing_put)
for a, b in zip(R_info, info):
    plt.text(a, h, b, thing, rotation=90)
plt.plot(R_A02, rho_A02, 'or')
plt.plot(R_A18, rho_A18, lw=lws)
plt.plot(R_limits, [rho_NK33_LOX]*2, '-k', lw=lwn)
plt.plot(R_limits, [rho_NK33_RP1]*2, '--k', lw=lwn)
plt.text(R_NK33_LOXlabel, rho_NK33_LOX+0.02, NK33_LOX_label, thing_LOX)
plt.text(R_NK33_RP1label, rho_NK33_RP1+0.02, NK33_RP1_label, thing_RP1)



fig.subplots_adjust(left=0.06, right=0.96,
                    bottom=0.04, top=0.98,
                    hspace=0.20 )
plt.savefig('Rocket Wax')
plt.show()