Skip to content

Instantly share code, notes, and snippets.

@dutta-alankar

dutta-alankar/field-length-sim.py

Last active January 28, 2026 16:34
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save dutta-alankar/a2d96ba972adad5ed754a03939a67352 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save dutta-alankar/a2d96ba972adad5ed754a03939a67352 to your computer and use it in GitHub Desktop.
Download ZIP
Field length calculation in sims
Raw
field-length-sim.py
#!/bin/python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 27 11:54:48 2026
@author: alankar.
Usage: time python lF-res.py
"""
import numpy as np
from scipy.interpolate import interp1d
import matplotlib.colors
import matplotlib.gridspec
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.lines import Line2D
from labellines import labelLine, labelLines
from decimal import Decimal
import os
import pathlib
import re
import numpy as np
from scipy.interpolate import interp1d
from scipy.stats import binned_statistic
import h5py
dark = False
extension = "pdf"
transparent = False
## Plot Styling
## Plot Styling
matplotlib.rcParams["xtick.direction"] = "in"
matplotlib.rcParams["ytick.direction"] = "in"
matplotlib.rcParams["xtick.top"] = False
matplotlib.rcParams["ytick.right"] = True
matplotlib.rcParams["xtick.minor.visible"] = True
matplotlib.rcParams["ytick.minor.visible"] = True
matplotlib.rcParams["axes.grid"] = True
matplotlib.rcParams["grid.linestyle"] = ":"
matplotlib.rcParams["grid.linewidth"] = 0.8
matplotlib.rcParams["grid.color"] = "gray"
matplotlib.rcParams["grid.alpha"] = 0.3
matplotlib.rcParams["lines.dash_capstyle"] = "round"
matplotlib.rcParams["lines.solid_capstyle"] = "round"
matplotlib.rcParams["legend.handletextpad"] = 0.4
matplotlib.rcParams["axes.linewidth"] = 1.0
matplotlib.rcParams["lines.linewidth"] = 3.5
matplotlib.rcParams["ytick.major.width"] = 1.2
matplotlib.rcParams["xtick.major.width"] = 1.2
matplotlib.rcParams["ytick.minor.width"] = 1.0
matplotlib.rcParams["xtick.minor.width"] = 1.0
matplotlib.rcParams["ytick.major.size"] = 11.0
matplotlib.rcParams["xtick.major.size"] = 11.0
matplotlib.rcParams["ytick.minor.size"] = 5.0
matplotlib.rcParams["xtick.minor.size"] = 5.0
matplotlib.rcParams["xtick.major.pad"] = 10.0
matplotlib.rcParams["xtick.minor.pad"] = 10.0
matplotlib.rcParams["ytick.major.pad"] = 6.0
matplotlib.rcParams["ytick.minor.pad"] = 6.0
matplotlib.rcParams["xtick.labelsize"] = 26.0
matplotlib.rcParams["ytick.labelsize"] = 26.0
matplotlib.rcParams["axes.titlesize"] = 24.0
matplotlib.rcParams["axes.labelsize"] = 28.0
matplotlib.rcParams["axes.labelpad"] = 8.0
plt.rcParams["font.size"] = 28
matplotlib.rcParams["legend.handlelength"] = 2
matplotlib.rcParams['legend.title_fontsize'] = '22'
# matplotlib.rcParams["figure.dpi"] = 200
matplotlib.rcParams["axes.axisbelow"] = True
matplotlib.rcParams["figure.figsize"] = (13,10)
if dark:
plt.style.use('dark_background')
colors = ["black", "yellowgreen", "steelblue", "darkorchid", "plum", "goldenrod", "crimson"]
cmap_name = "gist_rainbow" #earth_r"
# Source - https://stackoverflow.com/a
# Posted by jsalonen, modified by community. See post 'Timeline' for change history
# Retrieved 2026-01-28, License - CC BY-SA 3.0
def fexp(number):
(sign, digits, exponent) = Decimal(number).as_tuple()
return len(digits) + exponent - 1
def fman(number):
return Decimal(number).scaleb(-fexp(number)).normalize()
pc = 3.086e+18
mu = 0.609
kB = 1.380649e-16
mp = 1.6726e-24
km = 1.000e+05
s = 1.0
XH = 0.75
n_hot = 1.0e-03 # cgs cm^-3
UNIT_LENGTH = 100*pc
UNIT_DENSITY = n_hot*mu*mp
UNIT_VELOCITY = 1*km/s
gamma = 5/3.
chi = 100
Tcl = 1.0e+04
kappa_sp = 5.7e-07*(8.0e+04)**2.5 # cgs (erg/cm/s/K) at 8e4 K
kappa_sp = kappa_sp * (1.0e+06/8.0e+04)**2.5 # at 10^6 K
eta_sp = 0.11 # *(8.0e+04/1.0e+06)**2.5 # cgs (g/cm/s) at 10^6 K, same as visc_mu
Lbox = 1.0 # code
mach_shear = 0.5 # relative velocity mach number
n_cl = chi*n_hot
v_sh = np.sqrt(gamma*kB*chi*Tcl/(mu*mp))*mach_shear/UNIT_VELOCITY # code
bump_cool = 5.0e+00
Lambd_solar = interp1d(*np.loadtxt("cooltable-orig.dat",
unpack=True),
fill_value="extrapolate")
temperature = np.logspace(np.log10(Tcl), np.log10(Tcl*chi), 100)
cool_factor_min = np.min(temperature/(bump_cool*Lambd_solar(temperature)))
cool_factor_max = np.max(temperature/(bump_cool*Lambd_solar(temperature)))
# print("debug: ", cool_factor_max, cool_factor_min)
kappa_sp = 5.7e-07*(8.0e+04)**2.5 # cgs (erg/cm/s/K) at 8e4 K
kappa_sp = kappa_sp * (1.0e+06/8.0e+04)**2.5 # at 10^6 K
eta_sp = 0.11 # *(8.0e+04/1.0e+06)**2.5 # cgs (g/cm/s) at 10^6 K
custom_res = True
def field_length(fk, Nres, Nwav, color):
fnu = 0.5
kappa = fk*(mu*XH*n_hot*Lbox*UNIT_LENGTH/Nwav)**2/cool_factor_max
eta = fnu*np.sqrt(chi)*(n_hot*mu*mp/200)*(v_sh*Lbox/Nwav)*UNIT_VELOCITY*UNIT_LENGTH # cgs; only v_sh & L_box in code units
deltax_kappa = np.sqrt(fk*cool_factor_min/cool_factor_max)*Lbox/(Nres*Nwav)/chi # code
Npoints = int(np.power(2, np.ceil(np.log(max(Lbox/deltax_kappa, Nres*Nwav/fnu))/np.log(2))))
kappa_code = kappa*((mu*mp/kB)/(UNIT_DENSITY*UNIT_VELOCITY*UNIT_LENGTH))
eta_code = eta/(UNIT_DENSITY*UNIT_VELOCITY*UNIT_LENGTH)
print(f"kappa = {kappa_code:.2e} code = {kappa:.2e} erg/cm/s/K")
print(f"kappa/kappa_sp({(chi*Tcl):.1e} K) = {(kappa/kappa_sp):.2e}")
v_turb = v_sh/np.sqrt(chi)
t_cool_mix = (gamma/(gamma-1))*(n_hot*kB*chi*Tcl)/((np.sqrt(chi)*n_hot*mu*XH)**2*(bump_cool*Lambd_solar(np.sqrt(chi)*Tcl))) # cgs
# Nwav should not be in t_turb
t_turb = Lbox*UNIT_LENGTH/(v_turb*UNIT_VELOCITY) # cgs
Da = t_turb/t_cool_mix
nu = eta/(0.5*n_hot*(1+chi)*mu*mp) # cgs
Re = v_sh*Lbox*UNIT_VELOCITY*UNIT_LENGTH/(Nwav*nu)
till = 5.0
tKH = np.sqrt(chi)*Lbox/v_sh # code
if custom_res:
Np_used = 256
deltax_used = Lbox/Np_used # code
else:
Np_used = Npoints
deltax_used = deltax_kappa
temperature = np.logspace(4, 6, 100)
lF_all_phase = ((Lbox/Nwav) *
(temperature/(chi*Tcl)) *
np.sqrt(fk*temperature/(bump_cool*Lambd_solar(temperature))/cool_factor_max)) # code
print(f"L_box/l_F_max = {(Lbox/np.max(lF_all_phase)):.2e}")
# print(f"N_lambda/sqrt(fk) = {(Nwav/np.sqrt(fk)):.2e}")
print(f"L_box/l_F_min = {(Lbox/np.min(lF_all_phase)):.2e}")
plt.plot(temperature, lF_all_phase/deltax_used, color=color, label=rf"${fman(kappa/kappa_sp):.2f}\times10^{{ {int(fexp(kappa/kappa_sp))} }}$")
plt.xlim(xmin=np.min(temperature), xmax=np.max(temperature))
lines = plt.hlines([Nres, Lbox/deltax_used], xmin=np.min(temperature), xmax=np.max(temperature),
color="k", linestyle=":")
# Extract segments (returns list of [[x0, y0], [x1, y1]])
segs = lines.get_segments()
first_seg = segs[0]
# Create an invisible Line2D object with the same coordinates
fake_line = Line2D(first_seg[:, 0], first_seg[:, 1], visible=False)
plt.gca().add_line(fake_line)
x_data = fake_line.get_xdata()
x_min, x_max = x_data.min(), x_data.max()
target_x = x_min + 0.85 * (x_max - x_min)
labelLine(
fake_line,
target_x,
label=rf"$\ell_F/\Delta={Nres}$",
align=True,
yoffset=0.01,
ha="right",
fontsize=16,
backgroundcolor="none",
color="k"
)
last_seg = segs[1]
# Create an invisible Line2D object with the same coordinates
fake_line = Line2D(last_seg[:, 0], last_seg[:, 1], visible=False)
plt.gca().add_line(fake_line)
x_data = fake_line.get_xdata()
x_min, x_max = x_data.min(), x_data.max()
target_x = x_min + 0.10 * (x_max - x_min)
labelLine(
fake_line,
target_x,
label=rf"$\ell_F/\Delta={Np_used}$",
align=True,
yoffset=0.01,
ha="right",
fontsize=16,
backgroundcolor="none",
color="k"
)
return kappa, Np_used
fk = 8.0e-06
Nres = 4
Nwav = 0.02
kappa, Np_used = field_length(fk, Nres, Nwav, colors[6])
fk = 8.0e-05
Nres = 4
Nwav = 0.02
kappa, Np_used = field_length(fk, Nres, Nwav, colors[1])
fk = 8.0e-04
Nres = 4
Nwav = 0.02
kappa, Np_used = field_length(fk, Nres, Nwav, colors[2])
fk = 8.0e-03
Nres = 4
Nwav = 0.02
kappa, Np_used = field_length(fk, Nres, Nwav, colors[3])
fk = 8.0e-02
Nres = 4
Nwav = 0.02
kappa, Np_used = field_length(fk, Nres, Nwav, colors[5])
lines = plt.gca().get_lines()
for line in lines:
pattern = r"(\d+\.?\d*)\s*\\times\s*10\^\{\s*(-?\d+)\s*\}"
match = re.search(pattern, line.get_label())
if match:
coeff = match.group(1)
exp = match.group(2)
# Combine into Python's scientific notation format (e.g., "3.29e-2")
value = float(f"{coeff}e{exp}")
# print(f"Original: {line.get_label()} => Result: {value}")
else:
continue
kappa_line = value*kappa_sp #cgs
# print(rf"$\kappa={fman(kappa_line):.2f} \times 10^{{ {int(fexp(kappa_line))} }} {{\rm erg/cm/s/K}}$")
x_data = line.get_xdata()
x_min, x_max = x_data.min(), x_data.max()
target_x = x_min + 0.65 * (x_max - x_min)
labelLine(
line,
target_x,
label=rf"$\kappa={fman(kappa_line):.1f} \times 10^{{ {int(fexp(kappa_line))} }} {{\rm erg/cm/s/K}}$",
align=True,
yoffset=0.01,
ha="right",
fontsize=14,
backgroundcolor="none",
)
ymin, ymax = plt.gca().get_ylim()
plt.fill_between(temperature,
np.full(temperature.shape, 1.0),
np.full(temperature.shape, ymin),
color="gray", alpha=0.5)
# plt.ylim(ymin=ymin, ymax=ymax)
plt.legend(loc="upper left", prop={"size": 16}, ncol=1,
title=rf"$\frac{{\kappa}}{{\kappa_{{\rm sp}}({fman(chi*Tcl):.1f} \times 10^{{ {int(fexp(chi*Tcl))} }} K)}}$",
framealpha=0.5, shadow=False, fancybox=True)
plt.xscale("log")
plt.yscale("log")
plt.xlabel(r"Temperature [K]")
plt.ylabel(r"$\ell_{\rm F}/\Delta$")
plt.title(rf"$L_{{\rm box}} = {(UNIT_LENGTH/pc):.1f}\ \rm pc$, $n_{{\rm hot}} = {fman(n_hot):.1f} \times 10^{{ {int(fexp(n_hot))} }}\ {{\rm cm^{{-3}}}}$, $\Lambda_{{\rm bump}}={bump_cool:.1f}$, $L_{{\rm box}}/\Delta = {Np_used}$")
print(f"./lF-T-many-kappa-N_{Np_used}.{extension}")
plt.savefig(f"./lF-T-many-kappa-N_{Np_used}.{extension}",
transparent=transparent,
bbox_inches="tight")
Display the source blob
Display the rendered blob
Raw
lF-T-many-kappa-N_256.svg
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment