
#################################################
##
## This file extracts gravitational-wave strain amplitude
## from AMSS-NCKU numerical-relativity outputs and plots it.
## Author: Xiaoqu
## Dates: 2024/10/01 --- 2025/11/20
##
#################################################

import numpy                               ## numpy for array operations
import scipy                               ## scipy for interpolation and signal processing
import math
import matplotlib
matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
import matplotlib.pyplot    as     plt     ## matplotlib for plotting
import os                                  ## os for system/file operations

import AMSS_NCKU_Input as input_data

# plt.rcParams['text.usetex'] = True  ## enable LaTeX fonts in plots


####################################################################################



####################################################################################

## Compute Fourier transform for discrete time data in [t0, t1].
## Adapted from deepseek examples.

## Parameters:
## time : time series [0, t0]
## f : function value series
## apply_window : apply windowing (recommended True)
## zero_pad_factor : zero-padding factor (recommended 2-8)

## Returns:
## frequency : frequency axis (Hz)
## frequency_spectrum : complex frequency spectrum

def compute_frequency_spectrum(time, signal_f, apply_window=True, zero_pad_factor=4):

    ## Sampling parameters
    N  = len(time)
    dt = time[1] - time[0]        ## sample interval
    fs = 1/dt                     ## sampling frequency
    omega_s = fs * 2.0 * math.pi  ## sampling angular frequency

    ## Data preprocessing
    f_detrended = signal_f - numpy.mean(signal_f)  ## remove DC offset

    ## Windowing (reduce spectral leakage)
    if apply_window:
        # window = scipy.signal.windows.hann(N)  # Hann window
        window = scipy.signal.windows.tukey(N, alpha=0.1)  # or Tukey window
        f_windowed = f_detrended * window
    else:
        f_windowed = f_detrended
    
    # Zero-padding (improve frequency resolution)
    M = zero_pad_factor * N
    f_padded = numpy.zeros(M)
    f_padded[:N] = f_windowed
    
    # Execute FFT to obtain complex frequency spectrum

    # Use numpy.fft.fft to perform the fast Fourier transform
    # Basic signature: numpy.fft.fft(a, n=None, axis=-1, norm=None)
    # a: input 1-D array to be transformed
    # n: transform length; if n>len(a) the input is zero-padded, if n<len(a) it is truncated
    # axis: axis along which to compute the FFT (default: last axis)
    # norm: normalization mode (None, 'ortho', etc.)
    # numpy.fft.fft returns a complex array; the corresponding frequency axis spans -fs/2 to fs/2
    
    frequency_spectrum = numpy.fft.fft(f_padded, norm='ortho')    
    # frequency_spectrum = numpy.fft.fft(f_detrended, norm='ortho')  
    # frequency_spectrum = numpy.fft.fft(f_windowed, norm='ortho')     
    
    # Frequency axis generation

    # Use numpy.fft.fftfreq to build the discrete frequency axis:
    # numpy.fft.fftfreq(n, d=1.0)
    # n: signal length (number of samples)
    # d: sample spacing (time between samples)
    frequency = numpy.fft.fftfreq(M, dt)
    # frequency = numpy.fft.fftfreq(N, dt)

    # The above gives frequency in Hz; convert to angular frequency
    frequency_omega = 2.0 * math.pi * frequency
    
    # Amplitude correction (compensate for window loss)
    if apply_window:
        window_gain = numpy.mean(window)  # window gain
        frequency_spectrum /= window_gain
    
    ## return frequency_omega, frequency_spectrum
    return frequency, frequency_omega, frequency_spectrum

    ## Window types, properties and typical use cases:
    ## Hann    : wide main lobe, fast sidelobe decay       -> black-hole quasi-periodic oscillations (QPO)
    ## Tukey   : adjustable flat-top region               -> gravitational-wave chirp signals
    ## Blackman: optimal sidelobe suppression              -> high dynamic-range data

####################################################################################


####################################################################################

def frequency_filter_integration(omega, frequency_spectrum, omega0):

    '''
    ## modifiy the part |omega| < omega0
    omega_filter = numpy.where(omega < omega0, omega0, omega)
    '''

    ## Replace region |omega| < omega0 with signed omega0
    # Build replacements: positive omegas -> +omega0, negative -> -omega0
    replacements = numpy.where(omega >= 0, omega0, -omega0)

    # Apply replacement where abs(omega) < omega0
    omega_filter = numpy.where( numpy.abs(omega) < omega0, replacements, omega )

    ## Integrand for frequency-domain integration
    ## Note: convolution in time corresponds to multiplication in frequency
    frequency_integration = - frequency_spectrum / (omega_filter)**2

    return frequency_integration

####################################################################################


####################################################################################

## This function replaces |omega| < omega0 with signed omega0

def omega_filter(omega, omega0):

    ## Replace region |omega| < omega0 with signed omega0
    # Build replacements: positive omegas -> +omega0, negative -> -omega0
    replacements = numpy.where(omega >= 0, omega0, -omega0)
    
    # Apply replacements where |omega| < omega0
    # omega_filter = numpy.where(mask, replacements, omega)
    omega_filter = numpy.where( numpy.abs(omega) < omega0, replacements, omega )

    return omega_filter
    
####################################################################################


####################################################################################

## Inverse Fourier transform utility

## Inputs:
## omega : frequency axis (Hz)
## F_omega : complex frequency-domain data
## sampling_factor : sampling multiplier for reconstruction
## original_zero_pad_factor : original zero-padding factor used in forward FFT

## Returns:
## t : time axis
## reconstructed : reconstructed time-domain signal


def inverse_fourier_transform(omega, F_omega, sampling_factor=2, original_zero_pad_factor=4):
    
    # Compute sampling parameters
    N = len(F_omega)
    domega = omega[1] - omega[0]  # frequency resolution

    # Determine sampling rate
    # To avoid aliasing, ensure sampling rate >= 2 * max frequency (Nyquist)
    if sampling_factor > 2: 
        sampling_rate_omega = sampling_factor * omega.max()
    else:
        sampling_rate_omega = 2.0 * omega.max()
    # The input omega is angular frequency; convert to ordinary frequency
    ## dt = 1.0 / sampling_rate_omega
    frequency = omega / (2.0*math.pi)
    dt = 2 * math.pi / sampling_rate_omega  
    
    '''
    # DC component check
    if not numpy.isclose(omega[0], 0):
        warnings.warn("Frequency-domain data does not include zero-frequency component; DC offset may result")
    '''
    
    # Perform inverse FFT (use same normalization as forward transform)
    reconstructed_signal = numpy.fft.ifft(F_omega, norm='ortho') 
    # Note: numpy.fft.ifft already handles normalization for 'ortho'

    ## Build time axis
    ## If zero-padding was used originally, recover the unpadded length
    if (original_zero_pad_factor > 1):
        N0 = N // original_zero_pad_factor
        t = numpy.arange(0, N0*dt, dt)
        reconstructed_signal2 = reconstructed_signal[:N0]
    ## If no zero-padding
    else:
        t = numpy.arange(0, N*dt, dt)
        reconstructed_signal2 = reconstructed_signal
    
    # Handle real signals
    if numpy.allclose(numpy.imag(reconstructed_signal2), 0):
        reconstructed_signal3 = numpy.real(reconstructed_signal2)
    
    return t[:len(reconstructed_signal3)], reconstructed_signal3


####################################################################################


####################################################################################

# Instantaneous frequency estimation using analytic signal (Hilbert transform)

def instantaneous_frequency(signal, sampling_rate):
    """
    Compute instantaneous frequency of a signal.
    :param signal: input time-domain sampled signal
    :param sampling_rate: sampling rate
    :return: time array and instantaneous frequency array
    """
    analytic_signal = scipy.signal.hilbert(signal)
    phase = numpy.unwrap(numpy.angle(analytic_signal))
    time = numpy.arange(len(signal)) / sampling_rate
    frequency = numpy.gradient(phase, time) / (2 * numpy.pi)
    return time, frequency

def get_frequency_at_t1(signal, sampling_rate, t1):
    """
    Get instantaneous frequency at time t1
    :param signal: input time-domain sampled signal
    :param sampling_rate: sampling rate
    :param t1: target time
    :return: instantaneous frequency at t1
    """
    time, freq = instantaneous_frequency(signal, sampling_rate)
    index = numpy.argmin(numpy.abs(time - t1))
    return freq[index]
    
    
####################################################################################



####################################################################################

## Function to plot gravitational-wave waveform h

## Inputs:
## outdir             path to data directory
## figure_outdir      path to figure output directory
## detector_number_i  detector index
## total_mass         total system mass

def generate_gravitational_wave_amplitude_plot( outdir, figure_outdir, detector_number_i ):


    # build file path
    file0 = os.path.join(outdir, "bssn_psi4.dat")

    if ( detector_number_i == 0 ):
        print()
        print("Plotting the gravitational-wave strain amplitude h")
        print()
        print("The corresponding data file is", file0)
        print()

    print()
    print( "Plotting gravitational-wave data for detector no.", detector_number_i )

    
    # read entire data file, assume whitespace-delimited floats
    data = numpy.loadtxt(file0)
    
    # extract columns from psi4 file
    time                 = data[:,0]
    psi4_l2m2m_real      = data[:,1]
    psi4_l2m2m_imaginary = data[:,2]
    psi4_l2m1m_real      = data[:,3]
    psi4_l2m1m_imaginary = data[:,4]
    psi4_l2m0_real       = data[:,5]
    psi4_l2m0_imaginary  = data[:,6]
    psi4_l2m1_real       = data[:,7]
    psi4_l2m1_imaginary  = data[:,8]
    psi4_l2m2_real       = data[:,9]
    psi4_l2m2_imaginary  = data[:,10]
    
    # Note: file0 is just a filename; no file.open() was used, so nothing to close
    # file0.close()
    
    # Use integer division to compute length per detector
    length = len(time) // input_data.Detector_Number 
    
    time2                 = numpy.zeros( (input_data.Detector_Number, length) )
    psi4_l2m2m_real2      = numpy.zeros( (input_data.Detector_Number, length) )
    psi4_l2m2m_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) )
    psi4_l2m1m_real2      = numpy.zeros( (input_data.Detector_Number, length) )
    psi4_l2m1m_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) )
    psi4_l2m0_real2       = numpy.zeros( (input_data.Detector_Number, length) )
    psi4_l2m0_imaginary2  = numpy.zeros( (input_data.Detector_Number, length) )
    psi4_l2m1_real2       = numpy.zeros( (input_data.Detector_Number, length) )
    psi4_l2m1_imaginary2  = numpy.zeros( (input_data.Detector_Number, length) )
    psi4_l2m2_real2       = numpy.zeros( (input_data.Detector_Number, length) )
    psi4_l2m2_imaginary2  = numpy.zeros( (input_data.Detector_Number, length) )
    
    # Split data into per-detector-radius series
    for i in range(input_data.Detector_Number):
        for j in range(length):
            time2[i,j]                 = time[                 j*input_data.Detector_Number + i ]
            psi4_l2m2m_real2[i,j]      = psi4_l2m2m_real[      j*input_data.Detector_Number + i ]
            psi4_l2m2m_imaginary2[i,j] = psi4_l2m2m_imaginary[ j*input_data.Detector_Number + i ]
            psi4_l2m1m_real2[i,j]      = psi4_l2m1m_real[      j*input_data.Detector_Number + i ]
            psi4_l2m1m_imaginary2[i,j] = psi4_l2m1m_imaginary[ j*input_data.Detector_Number + i ]
            psi4_l2m0_real2[i,j]       = psi4_l2m0_real[       j*input_data.Detector_Number + i ]
            psi4_l2m0_imaginary2[i,j]  = psi4_l2m0_imaginary[  j*input_data.Detector_Number + i ]
            psi4_l2m1_real2[i,j]       = psi4_l2m1_real[       j*input_data.Detector_Number + i ]
            psi4_l2m1_imaginary2[i,j]  = psi4_l2m1_imaginary[  j*input_data.Detector_Number + i ]
            psi4_l2m2_real2[i,j]       = psi4_l2m2_real[       j*input_data.Detector_Number + i ]
            psi4_l2m2_imaginary2[i,j]  = psi4_l2m2_imaginary[  j*input_data.Detector_Number + i ]

    
    ## Compute discrete Fourier transforms of Psi4 data
    ## l=2 m=-2 spectrum
    psi4_l2m2m_real_frequency, psi4_l2m2m_real_omega, psi4_l2m2m_real_omega_spectrem                                                             \
        = compute_frequency_spectrum( time2[detector_number_i], psi4_l2m2m_real2[detector_number_i],      apply_window=True, zero_pad_factor=4 )
    psi4_l2m2m_imaginary_frequency, psi4_l2m2m_imaginary_omega, psi4_l2m2m_imaginary_omega_spectrem                                              \
        = compute_frequency_spectrum( time2[detector_number_i], psi4_l2m2m_imaginary2[detector_number_i], apply_window=True, zero_pad_factor=4 )
    ## l=2 m=-1 spectrum
    psi4_l2m1m_real_frequency, psi4_l2m1m_real_omega, psi4_l2m1m_real_omega_spectrem                                                             \
        = compute_frequency_spectrum( time2[detector_number_i], psi4_l2m1m_real2[detector_number_i],      apply_window=True, zero_pad_factor=4 )
    psi4_l2m1m_imaginary_frequency, psi4_l2m1m_imaginary_omega, psi4_l2m1m_imaginary_omega_spectrem                                              \
        = compute_frequency_spectrum( time2[detector_number_i], psi4_l2m1m_imaginary2[detector_number_i], apply_window=True, zero_pad_factor=4 )
    ## l=2 m=0 spectrum
    psi4_l2m0_real_frequency, psi4_l2m0_real_omega, psi4_l2m0_real_omega_spectrem                                                                \
        = compute_frequency_spectrum( time2[detector_number_i], psi4_l2m0_real2[detector_number_i],       apply_window=True, zero_pad_factor=4 )
    psi4_l2m0_imaginary_frequency, psi4_l2m0_imaginary_omega, psi4_l2m0_imaginary_omega_spectrem                                                 \
        = compute_frequency_spectrum( time2[detector_number_i], psi4_l2m0_imaginary2[detector_number_i],  apply_window=True, zero_pad_factor=4 )
    ## l=2 m=1 spectrum
    psi4_l2m1_real_frequency, psi4_l2m1_real_omega, psi4_l2m1_real_omega_spectrem                                                                \
        = compute_frequency_spectrum( time2[detector_number_i], psi4_l2m1_real2[detector_number_i],       apply_window=True, zero_pad_factor=4 )
    psi4_l2m1_imaginary_frequency, psi4_l2m1_imaginary_omega, psi4_l2m1_imaginary_omega_spectrem                                                 \
        = compute_frequency_spectrum( time2[detector_number_i], psi4_l2m1_imaginary2[detector_number_i],  apply_window=True, zero_pad_factor=4 )
    ## l=2 m=2 spectrum
    psi4_l2m2_real_frequency, psi4_l2m2_real_omega, psi4_l2m2_real_omega_spectrem                                                                \
        = compute_frequency_spectrum( time2[detector_number_i], psi4_l2m2_real2[detector_number_i],       apply_window=True, zero_pad_factor=4 )
    psi4_l2m2_imaginary_frequency, psi4_l2m2_imaginary_omega, psi4_l2m2_imaginary_omega_spectrem                                                 \
        = compute_frequency_spectrum( time2[detector_number_i], psi4_l2m2_imaginary2[detector_number_i],  apply_window=True, zero_pad_factor=4 )


    # Compute detector distance from input parameters
    Detector_Interval   = ( input_data.Detector_Rmax - input_data.Detector_Rmin ) / ( input_data.Detector_Number - 1 )
    Detector_Distance_R = input_data.Detector_Rmax - Detector_Interval * detector_number_i
    
    #################################################

    ## Set minimum cutoff frequency for frequency-domain integration

    ## Create output file to record frequency-domain cutoff values

    file_cut_path = os.path.join( figure_outdir, "frequency_cut.txt" )
    file_cut      = open( file_cut_path, "w" )
    
    ## Compute total mass of the system and output
    
    total_mass    = 0.0
    puncture_mass = numpy.zeros( input_data.puncture_number )
    
    ## For 'Ansorg-TwoPuncture' initial data: normalize masses of the first two black holes
    if ( input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
        mass_ratio_Q = input_data.parameter_BH[0,0] / input_data.parameter_BH[1,0]
        BBH_M1 = mass_ratio_Q / ( 1.0 + mass_ratio_Q )
        BBH_M2 = 1.0          / ( 1.0 + mass_ratio_Q )
        for k in range( input_data.puncture_number ):
            if ( k == 0 ):
                puncture_mass[k] = BBH_M1 
            elif( k == 1 ):
                puncture_mass[k] = BBH_M2 
            else: 
                puncture_mass[k] = input_data.parameter_BH[k,0]
            total_mass += puncture_mass[k]
     
    ## For other initial-data methods: read puncture masses from input
    else:
        for k in range( input_data.puncture_number ):
            puncture_mass[k] = input_data.parameter_BH[k,0]
            total_mass += puncture_mass[k]
            
    ## Output total mass
    print( file=file_cut )
    
    for k in range( input_data.puncture_number ):
        print( f" mass[{k}] = {puncture_mass[k]} ", file=file_cut )
 
    print( file=file_cut )
        
    print( f" total mass = {total_mass} ", file=file_cut )
    print(                                 file=file_cut )

    ## Compute and output pairwise puncture distances
    puncture_distance = numpy.zeros( (input_data.puncture_number, input_data.puncture_number) )
    puncture_position = input_data.position_BH

    ## Compute pairwise puncture separations
    for k1 in range(input_data.puncture_number):
        for k2 in range(input_data.puncture_number):
            if (k1 != k2):
                puncture_distance[k1,k2] = (   ( puncture_position[k1,0] - puncture_position[k2,0] )**2           \
                                             + ( puncture_position[k1,1] - puncture_position[k2,1] )**2           \
                                             + ( puncture_position[k1,2] - puncture_position[k2,2] )**2  )**(0.5)
                print( f" puncture distance r[{k1,k2}] = {puncture_distance[k1,k2]} ", file=file_cut )
                print(                                                                 file=file_cut )
            ## If k1 == k2, avoid zero-distance artifacts on the diagonal
            else:
                puncture_distance[k1,k2] = (   ( puncture_position[0,0] - puncture_position[0,0] )**2           \
                                             + ( puncture_position[0,1] - puncture_position[0,1] )**2           \
                                             + ( puncture_position[0,2] - puncture_position[0,2] )**2  )**(0.5)

    print( file=file_cut )

    ## Estimate orbital periods and frequencies using Newtonian approximation
    orbital_period    = numpy.zeros( (input_data.puncture_number, input_data.puncture_number) )
    orbital_frequency = numpy.zeros( (input_data.puncture_number, input_data.puncture_number) )

    ## Estimate maximum orbital frequency using Newtonian approximation
    frequency_max = ( numpy.max(puncture_distance) / numpy.min( puncture_mass ) )**(0.5)

    ## Estimate orbital period and frequency for each pair using Newtonian mechanics
    for k1 in range(input_data.puncture_number):
        for k2 in range(input_data.puncture_number):
            if (k1 != k2):
                orbital_period[k1,k2]    = 2.0 * math.pi * ( puncture_distance[k1,k2]**3 / ( puncture_mass[k1] + puncture_mass[k2] ) )**(0.5)
                orbital_frequency[k1,k2] = 1.0 / orbital_period[k1,k2]
                print( f" orbital period estimate:    T_orbital[{k1,k2}] = {orbital_period[k1,k2]} ",    file=file_cut )
                print( f" orbital frequency estimate: f_orbital[{k1,k2}] = {orbital_frequency[k1,k2]} ", file=file_cut )
                print(                                                                                   file=file_cut )
            else:
                orbital_frequency[k1,k2] = frequency_max
                orbital_period[k1,k2]    = 1.0 / orbital_frequency[k1,k2]

    print( file=file_cut )

    ## Set minimum frequency cutoff based on orbital estimate
    orbital_frequency_min       = numpy.min( orbital_frequency )
    gravitational_frequency_min = 2.0 * orbital_frequency_min     ## GW frequency ~ 2 * orbital frequency for quadrupole
    print( " Orbital frequency estimate:            f_orbital_min =", orbital_frequency_min, file=file_cut )
    print( " Gravitational Wave frequency estimate: f_GW_min      =", orbital_frequency_min, file=file_cut )
    print(                                                                                   file=file_cut )

    ## Set minimum frequency cutoff based on orbital estimate
    frequency_cut = gravitational_frequency_min
    omega_cut     = 2.0 * math.pi * frequency_cut
    print( " Frequency Cut estimate: frequency_cut =", frequency_cut, file=file_cut )
    print( " Omega Cut estimate:     omega_cut     =", omega_cut,     file=file_cut )
    print(                                                            file=file_cut )

    ## Manual cutoff setting (deprecated)
    ## omega_cut = 2.0 * math.pi / 100.0

    #################################################

    ## Set tortoise coordinate (r*) for waveform retarded-time correction
    tortoise_R = Detector_Distance_R + 2.0 * total_mass * math.log( Detector_Distance_R / (2.0*total_mass) - 1.0)
    
    ## For more than two punctures, tortoise coordinate is ambiguous; use detector radius
    if ( input_data.puncture_number > 2 ):
        tortoise_R = Detector_Distance_R

    ## Set cutoff based on instantaneous frequency of the Psi4 signal
    ## Abandoned due to large errors
    '''
    ## Set initial time
    t1 = tortoise_R

    ## Compute instantaneous frequency of Psi4 signals
    ## instantaneous_frequency_psi4_l2m2_real = instantaneous_frequency( psi4_l2m2_real2[detector_number_i], len(psi4_l2m2_real2[detector_number_i]) )
    instantaneous_frequency_psi4_l2m2m_real = get_frequency_at_t1( psi4_l2m2m_real2[detector_number_i], len(psi4_l2m2m_real2[detector_number_i]), t1 ) / (2.0*math.pi)
    instantaneous_frequency_psi4_l2m1m_real = get_frequency_at_t1( psi4_l2m1m_real2[detector_number_i], len(psi4_l2m1m_real2[detector_number_i]), t1 ) / (2.0*math.pi)
    instantaneous_frequency_psi4_l2m0_real  = get_frequency_at_t1( psi4_l2m0_real2[detector_number_i],  len(psi4_l2m0_real2[detector_number_i]),  t1 ) / (2.0*math.pi)
    instantaneous_frequency_psi4_l2m1_real  = get_frequency_at_t1( psi4_l2m1_real2[detector_number_i],  len(psi4_l2m1_real2[detector_number_i]),  t1 ) / (2.0*math.pi)
    instantaneous_frequency_psi4_l2m2_real  = get_frequency_at_t1( psi4_l2m2_real2[detector_number_i],  len(psi4_l2m2_real2[detector_number_i]),  t1 ) / (2.0*math.pi)
    print( f" Instantaneous frequency at t - r* = 0, l=2 m=-2 psi4_real = {instantaneous_frequency_psi4_l2m2m_real:.2f} 1/M" )
    print( f" Instantaneous frequency at t - r* = 0, l=2 m=-1 psi4_real = {instantaneous_frequency_psi4_l2m1m_real:.2f} 1/M" )
    print( f" Instantaneous frequency at t - r* = 0, l=2 m=0  psi4_real = {instantaneous_frequency_psi4_l2m0_real:.2f}  1/M" )
    print( f" Instantaneous frequency at t - r* = 0, l=2 m=1  psi4_real = {instantaneous_frequency_psi4_l2m1_real:.2f}  1/M" )
    print( f" Instantaneous frequency at t - r* = 0, l=2 m=2  psi4_real = {instantaneous_frequency_psi4_l2m2_real:.2f}  1/M" )

    ## Add frequency cutoffs based on instantaneous frequency
    frequency_cut_l2m2m = abs( instantaneous_frequency_psi4_l2m2m_real ) * 1.5
    frequency_cut_l2m1m = abs( instantaneous_frequency_psi4_l2m1m_real ) * 1.5
    frequency_cut_l2m0  = abs( instantaneous_frequency_psi4_l2m0_real  ) * 1.5
    frequency_cut_l2m1  = abs( instantaneous_frequency_psi4_l2m1_real  ) * 1.5
    frequency_cut_l2m2  = abs( instantaneous_frequency_psi4_l2m2_real  ) * 1.5
    
    ## Add frequency-domain filter conditions
    omega_cut_l2m2m = 2.0 * math.pi / frequency_cut_l2m2m
    omega_cut_l2m1m = 2.0 * math.pi / frequency_cut_l2m1m
    omega_cut_l2m0  = 2.0 * math.pi / frequency_cut_l2m0
    omega_cut_l2m1  = 2.0 * math.pi / frequency_cut_l2m1
    omega_cut_l2m2  = 2.0 * math.pi / frequency_cut_l2m2

    if (omega_cut_l2m0 < omega_cut_l2m2):
        omega_cut_l2m0 = omega_cut_l2m2
    if (omega_cut_l2m1 < omega_cut_l2m2):
        omega_cut_l2m1 = omega_cut_l2m2
    if (omega_cut_l2m1m < omega_cut_l2m2):
        omega_cut_l2m1m = omega_cut_l2m2
    if (omega_cut_l2m2m < omega_cut_l2m2):
        omega_cut_l2m2m = omega_cut_l2m2
    '''
    
    ## Obtain integrand for inverse Fourier transform
    psi4_l2m2m_real_omega_integration      = frequency_filter_integration( psi4_l2m2m_real_omega,      psi4_l2m2m_real_omega_spectrem,      omega_cut )
    psi4_l2m2m_imaginary_omega_integration = frequency_filter_integration( psi4_l2m2m_imaginary_omega, psi4_l2m2m_imaginary_omega_spectrem, omega_cut )
    psi4_l2m1m_real_omega_integration      = frequency_filter_integration( psi4_l2m1m_real_omega,      psi4_l2m1m_real_omega_spectrem,      omega_cut )
    psi4_l2m1m_imaginary_omega_integration = frequency_filter_integration( psi4_l2m1m_imaginary_omega, psi4_l2m1m_imaginary_omega_spectrem, omega_cut )
    psi4_l2m0_real_omega_integration       = frequency_filter_integration( psi4_l2m0_real_omega,       psi4_l2m0_real_omega_spectrem,       omega_cut )
    psi4_l2m0_imaginary_omega_integration  = frequency_filter_integration( psi4_l2m0_imaginary_omega,  psi4_l2m0_imaginary_omega_spectrem,  omega_cut )
    psi4_l2m1_real_omega_integration       = frequency_filter_integration( psi4_l2m1_real_omega,       psi4_l2m1_real_omega_spectrem,       omega_cut )
    psi4_l2m1_imaginary_omega_integration  = frequency_filter_integration( psi4_l2m1_imaginary_omega,  psi4_l2m1_imaginary_omega_spectrem,  omega_cut )
    psi4_l2m2_real_omega_integration       = frequency_filter_integration( psi4_l2m2_real_omega,       psi4_l2m2_real_omega_spectrem,       omega_cut )
    psi4_l2m2_imaginary_omega_integration  = frequency_filter_integration( psi4_l2m2_imaginary_omega,  psi4_l2m2_imaginary_omega_spectrem,  omega_cut )

    ## Perform inverse Fourier transform in frequency domain to obtain gravitational-wave strain amplitudes
    ## l=2 m=-2 amplitude
    time_grid_h_plus_l2m2m, GW_h_plus_l2m2m \
        = inverse_fourier_transform( psi4_l2m2m_real_omega, psi4_l2m2m_real_omega_integration, sampling_factor=2, original_zero_pad_factor=4 )  
    time_grid_h_cross_l2m2m, GW_h_cross_l2m2m \
        = inverse_fourier_transform( psi4_l2m2m_imaginary_omega, psi4_l2m2m_imaginary_omega_integration, sampling_factor=2, original_zero_pad_factor=4 )
    ## l=2 m=-1 amplitude
    time_grid_h_plus_l2m1m, GW_h_plus_l2m1m \
        = inverse_fourier_transform( psi4_l2m1m_real_omega, psi4_l2m1m_real_omega_integration, sampling_factor=2, original_zero_pad_factor=4 )  
    time_grid_h_cross_l2m1m, GW_h_cross_l2m1m \
        = inverse_fourier_transform( psi4_l2m1m_imaginary_omega, psi4_l2m1m_imaginary_omega_integration, sampling_factor=2, original_zero_pad_factor=4 )
    ## l=2 m=0 amplitude
    time_grid_h_plus_l2m0, GW_h_plus_l2m0 \
        = inverse_fourier_transform( psi4_l2m0_real_omega, psi4_l2m0_real_omega_integration, sampling_factor=2, original_zero_pad_factor=4 )  
    time_grid_h_cross_l2m0, GW_h_cross_l2m0 \
        = inverse_fourier_transform( psi4_l2m0_imaginary_omega, psi4_l2m0_imaginary_omega_integration, sampling_factor=2, original_zero_pad_factor=4 )
    ## l=2 m=1 amplitude
    time_grid_h_plus_l2m1, GW_h_plus_l2m1 \
        = inverse_fourier_transform( psi4_l2m1_real_omega, psi4_l2m1_real_omega_integration, sampling_factor=2, original_zero_pad_factor=4 )  
    time_grid_h_cross_l2m1, GW_h_cross_l2m1 \
        = inverse_fourier_transform( psi4_l2m1_imaginary_omega, psi4_l2m1_imaginary_omega_integration, sampling_factor=2, original_zero_pad_factor=4 )
    ## l=2 m=2 amplitude
    time_grid_h_plus_l2m2, GW_h_plus_l2m2 \
        = inverse_fourier_transform( psi4_l2m2_real_omega, psi4_l2m2_real_omega_integration, sampling_factor=2, original_zero_pad_factor=4 )  
    time_grid_h_cross_l2m2, GW_h_cross_l2m2 \
        = inverse_fourier_transform( psi4_l2m2_imaginary_omega, psi4_l2m2_imaginary_omega_integration, sampling_factor=2, original_zero_pad_factor=4 )  

    
    # Construct time grids for computing gravitational-wave strain h
    # time_max = max( time2[detector_number_i] )
    # time_grid = numpy.linspace( tortoise_R, time_max, 2000 )
    # time_grid_new = numpy.linspace( 0, time_max-tortoise_R, len(time_grid) )  ## subtract tortoise coordinate from times
    # l=2 m=-2
    time_grid_h_plus_l2m2m_new  = time_grid_h_plus_l2m2m - tortoise_R
    time_grid_h_cross_l2m2m_new = time_grid_h_cross_l2m2m - tortoise_R
    # l=2 m=-1
    time_grid_h_plus_l2m1m_new  = time_grid_h_plus_l2m1m - tortoise_R
    time_grid_h_cross_l2m1m_new = time_grid_h_cross_l2m1m - tortoise_R
    # l=2 m=0
    time_grid_h_plus_l2m0_new  = time_grid_h_plus_l2m0 - tortoise_R
    time_grid_h_cross_l2m0_new = time_grid_h_cross_l2m0 - tortoise_R
    # l=2 m=1
    time_grid_h_plus_l2m1_new  = time_grid_h_plus_l2m1 - tortoise_R
    time_grid_h_cross_l2m1_new = time_grid_h_cross_l2m1 - tortoise_R
    # l=2 m=2
    time_grid_h_plus_l2m2_new  = time_grid_h_plus_l2m2 - tortoise_R
    time_grid_h_cross_l2m2_new = time_grid_h_cross_l2m2 - tortoise_R

    plt.figure( figsize=(8,8) )                                   ## figsize controls figure size
    plt.title( f" Gravitational Wave h   Detector Distance = { Detector_Distance_R } ", fontsize=18 )   ## fontsize controls text size
    plt.plot( time_grid_h_plus_l2m0_new, GW_h_plus_l2m0,  \
              color='red',    label="l=2 m=0 h+",                  linewidth=2 )
    plt.plot( time_grid_h_cross_l2m0_new, GW_h_cross_l2m0, \
              color='orange', label="l=2 m=0 hx",  linestyle='--', linewidth=2 )
    plt.plot( time_grid_h_plus_l2m1_new, GW_h_plus_l2m1,  \
              color='green',  label="l=2 m=1 h+",                  linewidth=2 )
    plt.plot( time_grid_h_cross_l2m1_new, GW_h_cross_l2m1, \
              color='cyan',   label="l=2 m=1 hx",  linestyle='--', linewidth=2 )
    plt.plot( time_grid_h_plus_l2m2_new, GW_h_plus_l2m2,  \
              color='black',  label="l=2 m=2 h+",                  linewidth=2 )
    plt.plot( time_grid_h_cross_l2m2_new, GW_h_cross_l2m2, \
              color='gray',   label="l=2 m=2 hx",  linestyle='--', linewidth=2 )
    if ( input_data.puncture_number > 2 ):
        plt.xlabel( "T - R [M]",  fontsize=16     )
    else:
        plt.xlabel( "T - R* [M]", fontsize=16     )
    plt.ylabel( r"R*h",           fontsize=16     )
    plt.xlim( 0.0, max(time_grid_h_plus_l2m0_new) )
    plt.legend( loc='upper right'                 )
    plt.grid(   color='gray', linestyle='--', linewidth=0.5 )  # show grid lines
    plt.savefig( os.path.join(figure_outdir, "Gravitational_Wave_h_Detector_" + str(detector_number_i) + ".pdf") )
    

    print( "Gravitational-wave plot for detector no.", detector_number_i, "finished.")
    print()

    if ( detector_number_i == (input_data.Detector_Number-1) ):
        print( "All gravitational-wave strain amplitude plots finished." )
        print()
    
    '''
    # The following block performs direct time-domain integration of Psi4.
    # This method was deprecated because of insufficient accuracy.
    # h = int_{0}^{t} dt' int_{0}^{t"} Psi4(t") dt"

    # Interpolate per-detector data to obtain smooth functions
    # Use cubic spline interpolation
    psi4_l2m2m_real2_interpolation      = scipy.interpolate.interp1d( time2[detector_number_i], psi4_l2m2m_real2[detector_number_i],      kind='cubic' )
    psi4_l2m2m_imaginary2_interpolation = scipy.interpolate.interp1d( time2[detector_number_i], psi4_l2m2m_imaginary2[detector_number_i], kind='cubic' )
    psi4_l2m1m_real2_interpolation      = scipy.interpolate.interp1d( time2[detector_number_i], psi4_l2m1m_real2[detector_number_i],      kind='cubic' )
    psi4_l2m1m_imaginary2_interpolation = scipy.interpolate.interp1d( time2[detector_number_i], psi4_l2m1m_imaginary2[detector_number_i], kind='cubic' )
    psi4_l2m0_real2_interpolation       = scipy.interpolate.interp1d( time2[detector_number_i], psi4_l2m0_real2[detector_number_i],       kind='cubic' )
    psi4_l2m0_imaginary2_interpolation  = scipy.interpolate.interp1d( time2[detector_number_i], psi4_l2m0_imaginary2[detector_number_i],  kind='cubic' )
    psi4_l2m1_real2_interpolation       = scipy.interpolate.interp1d( time2[detector_number_i], psi4_l2m1_real2[detector_number_i],       kind='cubic' )
    psi4_l2m1_imaginary2_interpolation  = scipy.interpolate.interp1d( time2[detector_number_i], psi4_l2m1_imaginary2[detector_number_i],  kind='cubic' )
    psi4_l2m2_real2_interpolation       = scipy.interpolate.interp1d( time2[detector_number_i], psi4_l2m2_real2[detector_number_i],       kind='cubic' )
    psi4_l2m2_imaginary2_interpolation  = scipy.interpolate.interp1d( time2[detector_number_i], psi4_l2m2_imaginary2[detector_number_i],  kind='cubic' )

    # Compute detector distance from input parameters
    Detector_Interval   = ( input_data.Detector_Rmax - input_data.Detector_Rmin ) / ( input_data.Detector_Number - 1 )
    Detector_Distance_R = input_data.Detector_Rmax - Detector_Interval * detector_number_i

    # Set tortoise coordinate
    tortoise_R = Detector_Distance_R + 2.0 * total_mass * math.log( Detector_Distance_R / (2.0*total_mass) - 1.0)
    
    # Construct time grid for gravitational-wave amplitude h
    time_max = max( time2[detector_number_i] )
    time_grid = numpy.linspace( tortoise_R, time_max, 2000 )
    time_grid_new = numpy.linspace( 0, time_max-tortoise_R, len(time_grid) )  # subtract tortoise coordinate

    GW_h_plus_l2m2m  = numpy.zeros( len(time_grid) )
    GW_h_cross_l2m2m = numpy.zeros( len(time_grid) )
    GW_h_plus_l2m1m  = numpy.zeros( len(time_grid) )
    GW_h_cross_l2m1m = numpy.zeros( len(time_grid) )
    GW_h_plus_l2m0   = numpy.zeros( len(time_grid) )
    GW_h_cross_l2m0  = numpy.zeros( len(time_grid) )
    GW_h_plus_l2m1   = numpy.zeros( len(time_grid) )
    GW_h_cross_l2m1  = numpy.zeros( len(time_grid) )
    GW_h_plus_l2m2   = numpy.zeros( len(time_grid) )
    GW_h_cross_l2m2  = numpy.zeros( len(time_grid) )

    # Solve for h by double numerical integration: h = int_{0}^{t} dt' int_{0}^{t"} Psi4(t") dt" 
    # The double integral can be reordered and simplified to h = int_{0}^{t} (t-t") Psi4(t") dt"
    def GW_h_plus_l2m2m_integrand(t, tmax):
        return psi4_l2m2m_real2_interpolation(t) * (tmax-t)
    def GW_h_cross_l2m2m_integrand(t, tmax):
        return psi4_l2m2m_imaginary2_interpolation(t) * (tmax-t)
    def GW_h_plus_l2m1m_integrand(t, tmax):
        return psi4_l2m1m_real2_interpolation(t) * (tmax-t)
    def GW_h_cross_l2m1m_integrand(t, tmax):
        return psi4_l2m1m_imaginary2_interpolation(t) * (tmax-t)
    def GW_h_plus_l2m0_integrand(t, tmax):
        return psi4_l2m0_real2_interpolation(t) * (tmax-t)
    def GW_h_cross_l2m0_integrand(t, tmax):
        return psi4_l2m0_imaginary2_interpolation(t) * (tmax-t)
    def GW_h_plus_l2m1_integrand(t, tmax):
        return psi4_l2m1_real2_interpolation(t) * (tmax-t)
    def GW_h_cross_l2m1_integrand(t, tmax):
        return psi4_l2m1_imaginary2_interpolation(t) * (tmax-t)
    def GW_h_plus_l2m2_integrand(t, tmax):
        return psi4_l2m2_real2_interpolation(t) * (tmax-t)
    def GW_h_cross_l2m2_integrand(t, tmax):
        return psi4_l2m2_imaginary2_interpolation(t) * (tmax-t)
    
    # Compute gravitational-wave strains h+ and hx
    # Redefine integrand with a lambda so it becomes a single-variable function for integration

    for j in range( len(time_grid) ):

        print( " j = ", j )

        GW_h_plus_l2m2m_integrand2 = lambda t: GW_h_plus_l2m2m_integrand(t, time_grid[j])
        ## Note: scipy.integrate.quad returns a tuple (value, error)
        GW_h_plus_l2m2m[j], err0 = scipy.integrate.quad( GW_h_plus_l2m2m_integrand2, 0.0, time_grid[j], limit=600 )
                               # epsabs=1e-8,  # absolute tolerance
                               # limit=600 )    # increase number of subintervals

        GW_h_cross_l2m2m_integrand2 = lambda t: GW_h_cross_l2m2m_integrand(t, time_grid[j])
        GW_h_cross_l2m2m[j], err0 = scipy.integrate.quad( GW_h_cross_l2m2m_integrand2, 0.0, time_grid[j], limit=600 )

        GW_h_plus_l2m1m_integrand2 = lambda t: GW_h_plus_l2m1m_integrand(t, time_grid[j])
        GW_h_plus_l2m1m[j], err0 = scipy.integrate.quad( GW_h_plus_l2m1m_integrand2, 0.0, time_grid[j], limit=600 )

        GW_h_cross_l2m1m_integrand2 = lambda t: GW_h_cross_l2m1m_integrand(t, time_grid[j])
        GW_h_cross_l2m1m[j], err0 = scipy.integrate.quad( GW_h_cross_l2m1m_integrand2, 0.0, time_grid[j], limit=600 )

        GW_h_plus_l2m0_integrand2 = lambda t: GW_h_plus_l2m0_integrand(t, time_grid[j])
        GW_h_plus_l2m0[j], err0 = scipy.integrate.quad( GW_h_plus_l2m0_integrand2, 0.0, time_grid[j], limit=600 )

        GW_h_cross_l2m0_integrand2 = lambda t: GW_h_cross_l2m0_integrand(t, time_grid[j])
        GW_h_cross_l2m0[j], err0 = scipy.integrate.quad( GW_h_cross_l2m0_integrand2, 0.0, time_grid[j], limit=600 )

        GW_h_plus_l2m1_integrand2 = lambda t: GW_h_plus_l2m1_integrand(t, time_grid[j])
        GW_h_plus_l2m1[j], err0 = scipy.integrate.quad( GW_h_plus_l2m1_integrand2, 0.0, time_grid[j], limit=600 )

        GW_h_cross_l2m1_integrand2 = lambda t: GW_h_cross_l2m1_integrand(t, time_grid[j])
        GW_h_cross_l2m1[j], err0 = scipy.integrate.quad( GW_h_cross_l2m1_integrand2, 0.0, time_grid[j], limit=600 )

        GW_h_plus_l2m2_integrand2 = lambda t: GW_h_plus_l2m2_integrand(t, time_grid[j])
        GW_h_plus_l2m2[j], err0 = scipy.integrate.quad( GW_h_plus_l2m2_integrand2, 0.0, time_grid[j], limit=600 )

        GW_h_cross_l2m2_integrand2 = lambda t: GW_h_cross_l2m2_integrand(t, time_grid[j])
        GW_h_cross_l2m2[j], err0 = scipy.integrate.quad( GW_h_cross_l2m2_integrand2, 0.0, time_grid[j], limit=600 )
            
    # Computation of gravitational-wave amplitudes h+ and hx complete

    # Now perform plotting
    plt.figure( figsize=(8,8) )                                   ## figsize controls figure size
    plt.title( f" Gravitational Wave h   Detector Distance = { Detector_Distance_R } ", fontsize=18 )   ## fontsize controls text size
    plt.plot( time_grid_new, GW_h_plus_l2m0,  \
              color='red',    label="l=2 m=0 h+",                  linewidth=2 )
    plt.plot( time_grid_new, GW_h_cross_l2m0, \
              color='orange', label="l=2 m=0 hx",  linestyle='--', linewidth=2 )
    plt.plot( time_grid_new, GW_h_plus_l2m1,  \
              color='green',  label="l=2 m=1 h+",                  linewidth=2 )
    plt.plot( time_grid_new, GW_h_cross_l2m1, \
              color='cyan',   label="l=2 m=1 hx",  linestyle='--', linewidth=2 )
    plt.plot( time_grid_new, GW_h_plus_l2m2,  \
              color='black',  label="l=2 m=2 h+",                  linewidth=2 )
    plt.plot( time_grid_new, GW_h_cross_l2m2, \
              color='gray',   label="l=2 m=2 hx",  linestyle='--', linewidth=2 )
    plt.xlabel( "T [M]",          fontsize=16 )
    plt.ylabel( r"R*h",           fontsize=16 )
    plt.legend( loc='upper right'             )
    plt.savefig( os.path.join(figure_outdir, "Gravitational_Wave_h_Detector_" + str(detector_number_i) + ".pdf") )
    
    print()
    print( "Gravitational-wave plot for detector no.", detector_number_i, "finished.")
    print( "Plotting of gravitational-wave strain amplitude h completed.")
    print()
    '''

    return

####################################################################################



####################################################################################

## Standalone usage example
'''
## outdir = "./BBH_q=1"
outdir = "./3BH"
for i in range( input_data.Detector_Number ):
    generate_gravitational_wave_amplitude_plot(outdir, outdir, i)
'''
####################################################################################



