This script is used for calculating and plotting the stellar opacity data using interpolation methods. It involves reading the OPAL (Opacity Project At Livermore) opacity tables, converting and interpolating these data onto a specific grid expected by the DSEP (Dartmouth Stellar Evolution Program), and then calculating the fractional differences between interpolated and original data sets. The script also includes visualization of these calculations using matplotlib, showcasing the comparison between the given and calculated opacity values. It's tailored for researchers or students in astrophysics or related fields, focusing on aspects of stellar structure and evolution.

from pysep.opac.opal.defaults import GS98hz
from pysep.dm.filetypes import opacf_file
from pysep.opac.tops.api.convert import parse_RMO_TOPS_table_file
from pysep.opac.tops.api.convert import convert_rho_2_LogR
from pysep.opac.utils import get_target_log_R
from pysep.opac.utils import get_target_log_T
 
import matplotlib.pyplot as plt
import matplotlib as mpl
from scipy.interpolate import interp2d
import numpy as np
 
plt.rc('text', usetex=True)
plt.rc('font', family='Serif')
 
mpl.rcParams['xtick.minor.visible'] = True
mpl.rcParams['ytick.minor.visible'] = True
mpl.rcParams['xtick.direction'] = 'in'
mpl.rcParams['ytick.direction'] = 'in'
 
mpl.rcParams['xtick.top'] = True
mpl.rcParams['ytick.right'] = True
 
TOPS_Path = "/mnt/p/d/Astronomy/GraduateSchool/SecondYearProject/DSEP/runs/base/GS98_TOPS/GS98_opac/OP:77_0.70000004768_0.20029685231999994_0.0997031.dat"
 
rhoInit, logTInit, OPALTableInit = parse_RMO_TOPS_table_file(TOPS_Path)
 
 
# Interpolate raw TOPS data onto the LogR LogT grid DSEP expects
LogR, LogT, LogRMO = convert_rho_2_LogR(rhoInit, logTInit, OPALTableInit)
 
 
PostInter1 = (LogT, LogRMO[:, 13])
 
kappaFunc = interp2d(LogT, LogR, LogRMO.T, kind='cubic')
 
Rf = lambda T, rho: rho/((T*1e-6)**3)
 
targetLogT = get_target_log_T()
targetLogR = get_target_log_R()
 
rawValimg = list()
for logt in targetLogT:
    newList = list()
    t = 10**logt
    tidx = (abs(logTInit-logt)).argmin()
    for  logr in targetLogR:
        r = 10**logr
 
        rho = r*((t*1e-6)**3)
        idxRho = (abs(rhoInit - rho)).argmin()
 
        rhoB = max(idxRho-2, 0)
        rhoT = min(idxRho+2, rhoInit.shape[0])
 
        tB = max(tidx-2, 0)
        tT = min(tidx + 2, logTInit.shape[0])
        rawVal = OPALTableInit[rhoB:rhoT, tB:tT]
 
        kappa = interp2d(logTInit[tB:tT], rhoInit[rhoB:rhoT], rawVal)
        newList.append(kappa(logt, rho)[0])
    rawValimg.append(newList)
fig, ax = plt.subplots(1, 1, figsize=(10, 7))
inverseInterp = np.array(rawValimg)
 
fracDiff = (LogRMO-inverseInterp)/inverseInterp
lFracDiff = np.log10(fracDiff)
clFracDiff = np.nan_to_num(lFracDiff)
img = ax.imshow(lFracDiff, extent=(targetLogR.min(), targetLogR.max(), targetLogT.min(), targetLogT.max()), aspect=1.75)
 
ax.tick_params(axis="x", labelsize=22)
ax.tick_params(axis="y", labelsize=22)
ax.tick_params('both', length=12, width=2, which='major')
ax.tick_params('both', length=8, width=2, which='minor')
 
ax.axvline(x=-0.79, color='red')
cb = plt.colorbar(img)
cb.ax.tick_params(labelsize=17)
cb.ax.tick_params('y', length=5, width=1, which='major')
 
ax.set_xlabel('Log R [g cm$^{-3}$ K$^{-3}$]', fontsize=25)
ax.set_ylabel('Log T [K]', fontsize=25)
 
# cb.set_xlabel("Fractional Difference", fontsize=20)
 
 
plt.savefig("FractionalDifference.pdf", bbox_inches="tight")
print(lFracDiff.shape)
print((np.mean(fracDiff[:, 13])))
print((np.std(fracDiff[:, 13])))
print((np.mean(fracDiff[:, :])))
print((np.std(fracDiff[:, :])))
 
 
 
 
 
 
#
# opacity = np.zeros(shape=(logTInit.shape[0], rhoInit.shape[0]))
#
# opacityF = lambda x, rhoTarg: kappaFunc(x, Rf(10**x, rhoTarg))
#
# for i, rho in enumerate(rhoInit):
#     f = np.vectorize(lambda x: opacityF(x, rho))
#     opacity[:, i] = f(logTInit)
#
# n = 10
# print(opacity.shape, OPALTableInit.shape)
# print(list(enumerate(LogR)))
# print(LogT)
#
# print(LogRMO.shape)
# fig, ax = plt.subplots(1, 2, figsize=(20, 7))
#
# ax[0].imshow(OPALTableInit, extent=(logTInit.min(), logTInit.max(), rhoInit.min(), rhoInit.max()), aspect=2e-7)
# ax[0].plot(logTInit, 10**(-1.5)*((10**(logTInit))*1e-6)**3, color='red')
# ax[1].imshow(LogRMO.T, extent=(LogT.min(), LogT.max(), LogR.min(), LogR.max()))
# plt.show()