Skip to content

Instantly share code, notes, and snippets.

@oliverphilcox
Last active March 16, 2021 16:17
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save oliverphilcox/559f086f1bf63b23d55c508b2f47bad3 to your computer and use it in GitHub Desktop.
Compute the explicit solution of the spherical collapse equations, using the method of Slepian & Philcox (2021).
import numpy as np
def spherical_collapse(time_array, N_fft=32, eps=1e-4):
"""Compute the spherical collapse radius r as a function of time.
This is done with the method of Slepian & Philcox (2021), via contour integration and FFTs.
Note that the radius and time are given in dimensionless units.
Args:
time_array (np.ndarray): Array of times in units of the free-fall time.
Keyword Args:
N_fft (int): Number of FFT points, default: 32.
eps (float): Small parameter used to define contours, default 1e-4.
Returns:
np.ndarray: Array of r(t)/r_0 where r_0 is the initial radius.
"""
# Input checks
assert N_fft>0, "Need at least one FFT grid-point!"
assert (eps>0)&(eps<0.1), "Epsilon parameter must be small, else we might miss root!"
# Define cycloid function
def f(z, t):
return z + 0.5*np.sin(2.*z)-np.pi/2.*t
# Define contour integrand
def g(x, t):
return 1./f(np.pi/4.+(np.pi/4.-eps)*np.exp(2.0j*np.pi*x), t)
# Create FFT array and evaluate g(x;t)
x_array = np.arange(0.,1.,1./N_fft)[:,np.newaxis]
gx_array = g(x_array, time_array[np.newaxis,:])
# Perform FFT
ft_gx_array = np.fft.fft(gx_array,axis=0)
# Compute z0
z0_array = np.pi/4.+(np.pi/4.-eps)*np.real_if_close(ft_gx_array[-2]/ft_gx_array[-1])
# Compute output
radius_output = np.cos(z0_array)**2.
return radius_output
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment