Homework 5 solutions

We want to create a table the stores the values of our pressure at different densities and temperatures. Let’s start by defining the constants and the analytic function for the pressure:

import numpy as np

# CGS constants

k = 1.38e-16
a = 7.56e-15
mu_I = 1.26
mu_e = 1.15
K = 1.e13
m_u = 1.67e-24

Now we’ll write our pressure function. Note: this is an approximation of what the real pressure is inside a star. In particular, we’re assuming that the electrons behave as a non-relativistic degenerate gas, but at high temperatures and low densities, they are probably better described as an ideal gas.

def pressure(rho, T):
    return k * rho * T / (mu_I * m_u) + (1./3.) * a * T**4 + K * rho**(5./3.) * mu_e**(-5./3.)

Now lets write a function to tabulate this on a grid of density and temperature

def make_table(n_rho, n_T, rho_min=1.0, rho_max=1.e4, T_min=1.e6, T_max=1.e8):
    
    # create the density and temperature points where the data is stored
    
    rho_v = np.logspace(np.log10(rho_min), np.log10(rho_max), n_rho)
    T_v = np.logspace(np.log10(T_min), np.log10(T_max), n_T)
    
    # now create the table and fill it
    
    data = np.zeros((n_rho, n_T))
    
    for i, rho in enumerate(rho_v):
        for j, T in enumerate(T_v):
            data[i, j] = pressure(rho, T)
            
    return rho_v, T_v, data

Let’s create this table

n_rho = 10
n_T = 10

rho_v, T_v, data = make_table(n_rho, n_T)

Let’s look at a point just to see if it seems reasonable

i = 2
j = 5

print(f"rho = {rho_v[i]}, T = {T_v[j]}, P = {data[i, j]:20.10g}")

rho = 7.742636826811269, T = 12915496.650148828, P =      6.868491677e+15

Next we need to figure out what our interpolating function look like. We want to write the pressure as:

\[P(\rho, T) = a (\rho - \rho_i)(T - T_j) + b (\rho - \rho_i) + c (T - T_j) + d\]

where we get the coefficients from the 4 table values that surround the evaluation point.

This gives us the constraints:

\[\begin{align*} P_{i,j} &= d \\ P_{i+1,j} &= b \Delta \rho + d \\ P_{i,j+1} &= c \Delta T + d \\ P_{i+1,j+1} &= a \Delta \rho \Delta T + b \Delta \rho + c \Delta T + d \end{align*}\]

which is easily solved to give:

\[\begin{align*} a &= \frac{P_{i+1,j+1} - P_{i+1,j} - P_{i,j+1} + P_{i,j}}{\Delta \rho \Delta T} \\ b &= \frac{P_{i+1,j} - P_{i,j}}{\Delta \rho} \\ c &= \frac{P_{i,j+1} - P_{i,j}}{\Delta T} \\ d &= P_{i,j} \end{align*}\]

Now we can write our interpolation function

def interpolate_pressure(rho0, T0, rho_v, T_v, data):
    """given a point (rho0, T0), and the table data: rho_v, T_v, data,
    return an estimate of the pressure via interpolation"""
    
    # first we need to find the index i and j such that rho_v[i] < rho0 <= rho_v[i+1]
    
    i = np.argwhere(rho_v > rho0)[0][0] - 1
    assert rho_v[i] < rho0 <= rho_v[i+1]
    
    j = np.argwhere(T_v > T0)[0][0] - 1
    assert T_v[j] < T0 <= T_v[j+1]
        
    # now compute the coefficients
    drho = rho_v[i+1] - rho_v[i]
    dT = T_v[j+1] - T_v[j]
    
    a = (data[i+1, j+1] - data[i+1, j] - data[i, j+1] + data[i, j]) / (drho * dT)
    b = (data[i+1, j] - data[i, j]) / drho
    c = (data[i, j+1] - data[i, j]) / dT
    d = data[i, j]
    
    return a * (rho0 - rho_v[i]) * (T0 - T_v[j]) + b * (rho0 - rho_v[i]) + c * (T0 - T_v[j]) + d

interpolate_pressure(150, 1.5e7, rho_v, T_v, data)

1.8263573046485e+17