initial code

crust.py

@@ -0,0 +1,31 @@
+from polytropes import monotrope
+from polytropes import polytrope
+import units as cgs
+
+
+##################################################
+#SLy (Skyrme) crust
+KSLy = [6.80110e-9, 1.06186e-6, 5.32697e1, 3.99874e-8] #Scaling constants
+GSLy = [1.58425, 1.28733, 0.62223, 1.35692] #polytropic indices
+RSLy = [1.e4, 2.44034e7, 3.78358e11, 2.62780e12 ] #transition depths
+
+tropes = []
+trans = []
+
+pm = None
+for (K, G, r) in zip(KSLy, GSLy, RSLy):
+    m = monotrope(K*cgs.c**2, G)
+    tropes.append( m )
+
+    #correct transition depths to avoid jumps
+    if not(pm == None):
+        rho_tr = (m.K / pm.K )**( 1.0/( pm.G - m.G ) )
+        #print rho_tr, np.log10(rho_tr), r, rho_tr/r
+    else:
+        rho_tr = r
+    pm = m
+
+    trans.append(rho_tr)
+
+#Create crust using polytrope class
+SLyCrust = polytrope(tropes, trans)

eoslib.py

@@ -0,0 +1,154 @@
+import units as cgs
+from math import pi
+
+from polytropes import monotrope
+from polytropes import polytrope
+
+#import numpy as np
+
+#Dictionary from Read et al 2009 
+# all M_max < 2Msun commented out
+eosLib = {
+#    'PAL6'  :[ 34.380,  2.227,  2.189,  2.159, 'npem' ],
+    'SLy'   :[ 34.384,  3.005,  2.988,  2.851, 'npem' ],
+#    'APR1'  :[ 33.943,  2.442,  3.256,  2.908, 'npem' ],
+#    'APR2'  :[ 34.126,  2.643,  3.014,  2.945, 'npem' ],
+    'APR3'  :[ 34.392,  3.166,  3.573,  3.281, 'npem' ],
+    'APR4'  :[ 34.269,  2.830,  3.445,  3.348, 'npem' ],
+#    'FPS'   :[ 34.283,  2.985,  2.863,  2.600, 'npem' ],
+    'WFF1'  :[ 34.031,  2.519,  3.791,  3.660, 'npem' ],
+    'WFF2'  :[ 34.233,  2.888,  3.475,  3.517, 'npem' ],
+#    'WFF3'  :[ 34.283,  3.329,  2.952,  2.589, 'npem' ],
+#    'BBB2'  :[ 34.331,  3.418,  2.835,  2.832, 'npem' ],
+#    'BPAL12':[ 34.358,  2.209,  2.201,  2.176, 'npem' ],
+    'ENG'   :[ 34.437,  3.514,  3.130,  3.168, 'npem' ],
+    'MPA1'  :[ 34.495,  3.446,  3.572,  2.887, 'npem' ],
+    'MS1'   :[ 34.858,  3.224,  3.033,  1.325, 'npem' ],
+#    'MS2'   :[ 34.605,  2.447,  2.184,  1.855, 'npem' ],
+    'MS1b'  :[ 34.855,  3.456,  3.011,  1.425, 'npem' ],
+#    'PS'    :[ 34.671,  2.216,  1.640,  2.365, 'meson' ],
+#    'GS1a'  :[ 34.504,  2.350,  1.267,  2.421, 'meson' ],
+#    'GS2a'  :[ 34.642,  2.519,  1.571,  2.314, 'meson' ],
+#    'BGN1H1':[ 34.623,  3.258,  1.472,  2.464, 'hyperon' ],
+#    'GNH3'  :[ 34.648,  2.664,  2.194,  2.304, 'hyperon' ],
+#    'H1'    :[ 34.564,  2.595,  1.845,  1.897, 'hyperon' ],
+#    'H2'    :[ 34.617,  2.775,  1.855,  1.858, 'hyperon' ],
+#    'H3'    :[ 34.646,  2.787,  1.951,  1.901, 'hyperon' ],
+    'H4'    :[ 34.669,  2.909,  2.246,  2.144, 'hyperon' ],
+#    'H5'    :[ 34.609,  2.793,  1.974,  1.915, 'hyperon' ],
+#    'H6a'   :[ 34.593,  2.637,  2.121,  2.064, 'hyperon' ],
+#    'H7'    :[ 34.559,  2.621,  2.048,  2.006, 'hyperon' ],
+#    'PCL2'  :[ 34.507,  2.554,  1.880,  1.977, 'hyperon' ],
+#    'ALF1'  :[ 34.055,  2.013,  3.389,  2.033, 'quark' ],
+    'ALF2'  :[ 34.616,  4.070,  2.411,  1.890, 'quark' ],
+#    'ALF3'  :[ 34.283,  2.883,  2.653,  1.952, 'quark' ],
+#    'ALF4'  :[ 34.314,  3.009,  3.438,  1.803, 'quark' ]
+}
+
+
+#EoS class using dense eos from Read et al (2009) 
+def get_eos(key):
+
+    #read eos table for parameters
+    ll = eosLib[ key ]
+
+    #dense eos starting pressure
+    p1 = 10**ll[0]
+
+    #polytrope indices
+    g1 = ll[1] 
+    g2 = ll[2]
+    g3 = ll[3]
+
+    #transition densities
+    r1 = 2.8e14
+    r2 = 10**14.7
+    r3 = 10**15.0
+
+    #scaling constants
+    K1 = p1/(r2**g1)
+    K2 = K1 * r2**(g1-g2)
+    K3 = K2 * r3**(g2-g3)
+
+    tropes = [monotrope(K1, g1),
+              monotrope(K2, g2),
+              monotrope(K3, g3) ]
+    trans = [r1, r2, r3]
+
+    dense_eos = polytrope(tropes, trans)
+
+    return dense_eos
+
+# Smoothly glue core to SLy crust
+# for polytropic eos we can just unpack 
+# and repack the piecewise presentation
+def glue_crust_and_core(crust, core):
+
+    #unpack crust and core
+    tropes_crust = crust.tropes
+    trans_crust  = crust.transitions
+
+    tropes_core = core.tropes
+    trans_core  = core.transitions
+
+    #find transition depth
+    rho_tr = (tropes_core[0].K / tropes_crust[-1].K )**( 1.0/( tropes_crust[-1].G - tropes_core[0].G ) )
+    #print "Transition from core to crust at", rho_tr, np.log10(rho_tr), crust.edens_inv( crust.pressure( rho_tr ) )/cgs.GeVfm_per_dynecm
+    trans_core[0] = rho_tr
+    #trans_crust[-1] = rho_tr
+
+    #repack
+    tropes = tropes_crust + tropes_core
+    trans  = trans_crust  + trans_core
+
+    for trope in tropes:
+        trope.a = 0.0
+
+    eos = polytrope( tropes, trans )
+    return eos
+
+
+
+
+
+
+##################################################
+# Some simple phenomenological mono/bitrope EoSs
+def simple_eos():
+
+    K1 = 3.99873692e-8 
+    G1 = 1.35692395 
+    r1 = 0
+
+    K2 = 2.23872092e-10
+    G2 = 3 
+    #a = 0.010350691 * c *c
+    r2 = 1.4172900e14
+
+    m1 = monotrope(K1*(cgs.c**2), G1)
+    m2 = monotrope(K2, G2)
+    eos = polytrope([m1, m2], [r1, r2])
+
+    return eos
+
+def simple_eos2():
+    Gamma0 = 5.0/3.0 
+    K0 = (3.0*pi**2)**(2.0/3.0)*cgs.hbar**2/(5.0*cgs.mn**(8.0/3.0))
+    Gamma1 = 3.0 # high densities: stiffer equation of state
+    #Gamma1 = 2.5 # high densities: softer equation of state
+    rho1 = 5e14
+    P1 = K0*rho1**Gamma0
+    K1 = P1/rho1**Gamma1
+
+    m1 = monotrope(K0, Gamma0)
+    m2 = monotrope(K1, Gamma1)
+    eos = polytrope([m1, m2], [0.0, rho1])
+
+    return eos
+
+def simpler_eos():
+    K = 1.982e-6
+    G = 2.75
+    m = monotrope(K, G)
+    eos = polytrope([m], [0.0])
+    return eos

label_line.py

@@ -0,0 +1,51 @@
+import numpy as np
+import math
+import matplotlib.pyplot as plt
+
+def label_line(line, label_text, near_i=None, near_x=None, near_y=None, rotation_offset=0, offset=(0,0)):
+    """call 
+        l, = plt.loglog(x, y)
+        label_line(l, "text", near_x=0.32)
+    """
+    def put_label(i):
+        """put label at given index"""
+        i = min(i, len(x)-2)
+        dx = sx[i+1] - sx[i]
+        dy = sy[i+1] - sy[i]
+        rotation = np.rad2deg(math.atan2(dy, dx)) + rotation_offset
+        pos = [(x[i] + x[i+1])/2. + offset[0], (y[i] + y[i+1])/2 + offset[1]]
+        txt = plt.text(pos[0], pos[1], label_text, size=5, rotation=rotation, color = line.get_color(),
+        ha="center", va="center", bbox = dict(ec='1',fc='1',pad=0))
+
+        return txt
+
+
+    x = line.get_xdata()
+    y = line.get_ydata()
+    ax = line.get_axes()
+    if ax.get_xscale() == 'log':
+        sx = np.log10(x)    # screen space
+    else:
+        sx = x
+    if ax.get_yscale() == 'log':
+        sy = np.log10(y)
+    else:
+        sy = y
+
+    # find index
+    if near_i is not None:
+        i = near_i
+        if i < 0: # sanitize negative i
+            i = len(x) + i
+        put_label(i)
+    elif near_x is not None:
+        for i in range(len(x)-2):
+            if (x[i] < near_x and x[i+1] >= near_x) or (x[i+1] < near_x and x[i] >= near_x):
+                put_label(i)
+    elif near_y is not None:
+        for i in range(len(y)-2):
+            if (y[i] < near_y and y[i+1] >= near_y) or (y[i+1] < near_y and y[i] >= near_y):
+                put_label(i)
+    else:
+        raise ValueError("Need one of near_i, near_x, near_y")
+

polytropes.py

@@ -0,0 +1,106 @@
+import units as cgs
+
+#Monotropic eos
+class monotrope:
+
+    #transition continuity constant 
+    a = 0.0
+    cgsunits = cgs.c**2
+
+    def __init__(self, K, G):
+        self.K = K / self.cgsunits
+        self.G = G
+        self.n = 1.0/(G - 1)
+
+    #pressure P(rho)
+    def pressure(self, rho):
+        return self.cgsunits * self.K * rho**self.G
+
+    #energy density mu(rho)
+    def edens(self, rho):
+        return (1.0 + self.a)*rho + (self.K/(self.G - 1)) * rho**self.G
+
+    #for inverse functions lets define rho(P)
+    def rho(self, press):
+        if press < 0.0:
+            return 0.0
+        return ( press/self.cgsunits/self.K )**(1 / self.G)
+
+
+# Piecewise polytropes
+class polytrope:
+
+    def __init__(self, tropes, trans, prev_trope = None ):
+        self.tropes      = tropes
+        self.transitions = trans
+
+        self.prs  = []
+        self.eds  = []
+
+        for (trope, transition) in zip(self.tropes, self.transitions):
+
+            if not( prev_trope == None ):
+                trope.a = self._ai( prev_trope, trope, transition )
+            else:
+                transition = 0.0
+
+
+            ed = trope.edens(transition) 
+            pr = trope.pressure(transition)
+
+            self.prs.append( pr )
+            self.eds.append( ed )
+
+            prev_ed = ed
+            prev_tr = transition
+            prev_trope = trope
+
+
+    def _ai(self, pm, m, tr):
+        return pm.a + (pm.K/(pm.G - 1))*tr**(pm.G-1) - (m.K/(m.G - 1))*tr**(m.G-1)
+
+    def _find_interval_given_density(self, rho):
+        if rho <= self.transitions[0]:
+            return self.tropes[0]
+
+        for q in range( len(self.transitions) - 1 ):
+            if self.transitions[q] <= rho < self.transitions[q+1]:
+                return self.tropes[q]
+
+        return self.tropes[-1]
+
+    #inverted equations as a function of pressure
+    def _find_interval_given_pressure(self, press):
+        if press <= self.prs[0]:
+            return self.tropes[0]
+
+        for q in range( len(self.prs) - 1):
+            if self.prs[q] <= press < self.prs[q+1]:
+                return self.tropes[q]
+
+        return self.tropes[-1]
+
+
+    ################################################## 
+    def pressure(self, rho):
+        trope = self._find_interval_given_density(rho)
+        return trope.pressure(rho)
+
+    #vectorized version
+    def pressures(self, rhos):
+        press = []
+        for rho in rhos:
+            pr = self.pressure(rho)
+            press.append( pr )
+        return press
+
+    def edens_inv(self, press):
+        trope = self._find_interval_given_pressure(press)
+        rho = trope.rho(press)
+        return trope.edens(rho)
+
+    def rho(self, press):
+        trope = self._find_interval_given_pressure(press)
+        return trope.rho(press)
+
+