Last active
March 16, 2021 16:17
Star
You must be signed in to star a gist
Compute the explicit solution of the spherical collapse equations, using the method of Slepian & Philcox (2021).
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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